Aphids are one of the most harmful insect species for crops and are thus an important group to model in ALMaSS because they drive both management and natural enemy dynamics. Due to the enormous number of aphids, using an agent-based method to model aphids is not feasible. To deal with this, we propose creating an Aphid model using a spatial subpopulation model approach. Built on the dynamic landscape model in ALMaSS, subpopulations are assumed to occupy the landscape in a set of regular grid cells of equal size. The aphid population will be tracked in each grid using a stage-structured population model; population size and behaviour will depend on the habitat structure under the grid. They will develop, reproduce, and die based on the aphids’ physiological processes, the effects of predators, and the grids’ weather and vegetation host conditions. Aphid dispersal by flying and local movement to adjacent grids will interconnect the populations in the grid cells.

Acyrthosiphon pisum, ALMaSS, Aphis fabae, black bean aphid, grain aphid, Myzus persicae, pea aphid, peach aphid, Sitobion avenae, subpopulation

Aphids are a group of sap-sucking insects comprising 4000 species, among which 100 species colonise agriculturally important crops (Emden and Harrington 2007). Due to their rapid reproduction rate and high dispersal ability, aphids can cause considerable economic damage to agriculture (Hullé et al. 2010). Aphids mainly cause damage to plants in three ways.

1. Direct physical damage through feeding and diverting phloem-sap. Aphids penetrate the plant surface and feed on the nutrients in the phloem sap from the phloem vessels, while the injected saliva prevents the plant tissue from sealing (Forbes 1977); this diverts part of the plant’s nutrients to the aphid. If the number of aphids present per plant is large, it can cause a significant reduction in plant biomass and growth rate (Burd and Burton 1992).
2. Indirect damage by transmitting plant viruses. The predominant cause of crop damage by aphids is the transmission of plant viruses. Aphids can spread 275 of the >700 insect-vectored plant viruses (Francki et al. 1991; Nault 1997). Valenzuela and Hoffman (2014) reviewed the yield effects caused by aphids on grain crops in Australia. They found that feeding and virus injuries caused an overall economic cost of $241 million annually and $482 million yearly, respectively. Tatchel (1989) estimated that medium and severe aphid outbreaks in Britain potentially caused £70 million per year and £120 million per year, respectively, for wheat alone.
3. Indirect damage by honeydew production and honeydew-induced moulds. Honeydew produced by aphids can impact photosynthesis by impairing carbon dioxide exchange (Rabbinge et al. 1981). In addition, honeydew promotes the growth of fungi, further impairing photosynthesis.

Many factors have an impact on the size of the aphid population. The effect of climatic conditions, especially temperature and wind, and the presence of natural enemies have been shown to be key drivers affecting aphid populations and outbreaks (Wellings and Dixon 1987). For example, warmer winter temperatures can lead to aphid outbreaks in early-sown crops (Clement et al. 2010). Warmer spring/summer temperatures increase aphid fecundity and decrease development time, resulting in a higher intrinsic rate of increase. Wind and rainfall can also influence aphid outbreaks, as those factors can change the ability and timing of migration for aphids, thereby hindering their ability to spread to different crops (Clement 2006). Sowing timing has also been shown to impact aphid population dynamics. Early sown crops, possibly due to cold January/February temperatures, resulted in aphid populations crashing before flowering, reducing economic damage (McVean et al. 1999). On the other hand, late sowing resulted in aphid populations peaking around flowering. For host alternating aphid species, the primary host plant spatial distribution and density are key factors in determining, e.g., Aphis fabae, outbreak risks in Britain (Cammell et al. 1989). Another important factor affecting aphid populations is aphid control conducted by humans. Insecticide application is a common and fairly simple way of controlling aphid populations; however, this management strategy has drawbacks. Intensive application of insecticides over many years has led to widespread resistance among several species of aphids (Edwards et al. 2008; Bass et al. 2014; Bass et al. 2015). In some instances, insecticide application can increase aphid populations by damaging the aphid predator populations more than the aphids, thereby removing the biological control (Kerns and Gaylor 1993). Biological control is another more sustainable option for managing aphid outbreaks, achieved by introducing specialist predators or parasitic wasps to control aphid populations (Pervez and Omkar 2006; Boivin et al. 2012). As such, aphid populations do not exist in isolation but are important in supporting a range of generalist and specialist populations of predator and parasitoids.

In order to predict the aphid population, it is important to have a modelling tool to include all the listed factors influencing or potentially influencing aphid population dynamics. Duffy et al. (2017) produced a model to describe Sitobion avenae population dynamics on winter wheat. The model accurately described the seasonal population dynamics of S. avenae on winter wheat, making it a potentially valuable tool for aphid pest management in that specific context. However, their model describes the processes without dealing with space explicitly.

In this paper, we follow the Formal Model format proposed by (Topping et al. 2022). We aim to describe how to model four aphid species, Acyrthosiphon pisum, Aphis fabae, Myzus persicae and Sitobion avenae, using a subpopulation method (Duan and Topping 2024) to run in a realistic landscape environment (Topping and Duan 2024a) (usually 10 km X 10 km) in the Animal Landscape and Man Simulation System (ALMaSS) (Topping et al. 2003; Topping 2022).

The implementation of the four aphid models will be the first model using the subpopulation approach within the ALMaSS framework (Duan and Topping 2024); as such, it forms a test case for the new method. The model should represent the four aphid species’ population distribution and abundance in dynamic landscapes when implemented. The models should respond to changes in crop management and will be used to assess the effectiveness of control measures and to help understand the spatiotemporal dynamics of aphid populations. A primary aim of the models will be coupling the dynamic aphid models with aphid predators, e.g., ladybird, within the same system. We aim to use this coupled system to assess management impacts on natural enemies and evaluate aphid control using ladybirds. As such, the aphid model may form an important tool in evaluating integrated pest management at regional scales.

The aphid models in ALMaSS differ from the animal models developed in the system because it is not an agent-based model (Topping and Duan 2024b). Agent-based models require the understanding and knowledge of an individual’s characteristics and behaviours. Due to the enormous numbers of aphids and lack of individual-level knowledge, it is impossible to model aphids as individual agents, nor is it even possible to define super-individuals (Scheffer et al. 1995). Hence, to manage this practically, the aphid models will utilise a hybrid approach between stage-structured population and dynamic spatial models. We refer to this as the ALMaSS subpopulation model (Topping and Duan 2024b). For this aphid implementation, the model splits the landscape into equally sized (10 meters by 10 meters by default) relatively small units as a continuous grid. Each square is used to model a single subpopulation of aphids and describes the development, reproduction, and death of the aphids within that grid square per time step (one day). Conditions within each grid square are defined based on its context in the ALMaSS landscapes (e.g., host plant availability). All subpopulations can, in theory, interact based on dispersal, which is implemented by local movement and long-distance flight between grids.

The aphids are divided into a complex stage-structured population representation within each subpopulation. Various life stages are defined, and individuals of the same age form cohorts within these. These cohorts are reasonable since these individuals may be extremely similar, even genetically identical. Each day, the number of individuals passing from one day to the next is modelled, and when a cohort reaches the last day defined by degree days of that stage, it is transferred to the next life stage or death. Death can occur at any stage or age, but new individuals can only enter as eggs or newly produced nymphs. This modelling approach has the advantage that it still uses the power of the ALMaSS dynamic landscape modelling but also copes with many billions of aphids. The resolution of the population planned is at 10 × 10 m grid squares on landscapes of 10 × 10 km. Unlike the models to model aphids on a single host in a single field developed by previous researchers (e.g., Duffy, Fealy et al. 2017; Carter 1982), our model attempts to model aphids at a 100 km² scale, requiring the inclusion a variety of mechanisms previous models did not account for. Our approach offers a more holistic view of aphid population dynamics at a large scale, but it requires a higher level of knowledge about all mechanisms as they interact with each other at the landscape level. The approach allows for interactions between the aphids and other organisms (e.g., parasitoids) explicitly within the model; it can also be coupled with other models in the ALMaSS system to generate external interactions. Subpopulations may be limited not only by habitat conditions conferring population limits but also by interactions with natural enemies.

A subpopulation approach was taken to represent the billions of aphids in the simulated landscapes. However, this limits the detail in which individual interactions are specified; thus, the model will not represent feedback related to individual characteristics. For instance, (Llewellyn and Brown 1985) showed that adult weight positively correlates with fecundity, but this process cannot be easily represented in a subpopulation model. This, and similar limitations in detail, are a general drawback of the approach, although the subpopulation model created is still very detailed. Another limitation is that, as is typical when building a model of this detail, much of the data comes from laboratory and field studies carried out in different places, times, and methodologies.

Aphids exhibit considerable intraspecific variation in how they can utilise host plants, with distinct populations of the same species demonstrating varied performance across different host plants (Via 1991b; Nikolakakis et al. 2003; Gorur et al. 2005). A. pisum consists of several genetically different groups specialised for various host plants. However, the species shows a high degree of ecological specialisation among populations. Populations perform better on their respective “home” plants (e.g., higher fecundity, faster development and lower mortality) (Ferrari et al. 2008). A higher resistance to parasitoids and predators when on the “home” plant was demonstrated for A. pisum (Hufbauer and Via 1999). Since it is not possible to know the details of this specialisation, and often hosts in experiments and fields vary with location, this has some important consequences for the model:

• This ecological specialisation is not included in this version.
• A potential implication could be that populations spread faster and further than they naturally could since every aphid population in the model performs similarly on every possible host plant.
• The total population growth might be overestimated. Ferrari et al. (2008) showed that A. pisum clones had fecundity rates of up to 2.61 (+-0.13) times higher on the host plant that they were collected from than other host plants. Furthermore, since host plants in the model are not differentiated but treated as equal hosts, this might overestimate population growth each time a clone switches host but receives no growth penalty. To accurately reflect this, the origin of each aphid would need to be tracked, which would be too resource-intensive to implement.
• The model aphids described are, in large part, a reflection of the laboratory data used to simulate key performance parameters, such as fecundity, development time and longevity. This makes the assumption that, in the model, all populations perform as in these experiments on every type of host plant.

In addition to the ecological specialisation previously outlined, the model presently does not account for the intraspecific variation observed in aphid life cycles. Clones of A. pisum and S. avenae have been shown to have different responses to photoperiods, with some clones producing a combination of sexual morphs and virginoparae, while others exclusively produced virginoparae. Those differences reflect different overwintering strategies. However, our current model simplifies these dynamics and assumes a homogenous response across aphid populations and similar response between the four modelled species. The model’s current framework is designed to focus on macro-scale population dynamics and general trends in habitat utilisation. The limit in the detail of the ALMaSS landscape model also influences the aphid model via host plants. Host plants are assigned based on suitable general vegetation types, and each 1 m² is assigned a vegetation type. However, there is considerable variation within a habitat type (e.g., a field margin) in the natural world. This variation cannot be captured in the present model and may lead to over or underestimation of habitat suitability for the aphid. It will undoubtedly lead to a more homogeneous pattern of host plant distribution than would be the case in reality. A key issue for laboratory-based model parameters is that these were usually obtained under constant temperature conditions. The effect of fluctuating temperatures is thereby not accounted for in the simulation. For instance, there may be faster development at alternating temperatures in the lower thermal range compared to higher temperatures (Siddiqui et al. 1973). In fact, a key trait of insects is that ecothermic development is non-linear, and asymmetrical (Angilletta Jr 2009). This is a general problem unless there is enough information to implement a non-linear day-degree model, considering changing rates towards temperature extremes. If this information is available, a future version of the model could implement a thermal performance curve model to describe growth and temperature relationships (e.g. Angilletta 2006).

Additionally, aphids do not necessarily attain ambient air temperature, as they are in close contact with the host plant, bringing their temperature closer to that of the plant. The leaf temperature, in turn, depends on the growth stage of the plant and environmental conditions. This nuance adds another layer of complexity to accurately modelling aphid development and performance.

An overview of the aphid model, the components (movement, morph determination, development, reproduction, and mortality) driving the population dynamics and the variables influencing the components is shown in Fig. 1. Development describes the development rate for nymphs and longevity for adults (Influenced by temperature). Reproduction calculates the rate at which offspring is produced (Sexual/Asexual), determined by environmental temperature and age. Movement determines the direction, timing and distance that winged aphids can fly. Morph determination determines which type of nymphs are born (Male, oviparae, alate virginoparae or alate virginoparae), influenced by photoperiod, aphid density and host quality. Mortality calculates how many aphids die each day based on parasitoids, temperature and density. We present a detailed description of the processes in the following sections.

Figure 1. 

Process diagram illustrating the framework of the aphid model. Black arrows represent direct influence, grey indirect influence and dotted lines represent identified but not simulated influence.

The aphid life cycle exhibits remarkable phenotypic plasticity, characterised by the production of distinct morphological forms despite being genetically identical. Form diversity includes winged and wingless variants, which are pivotal for understanding aphid ecology and modelling their population dynamics (Ogawa and Miura 2014). To lay a foundation for the subsequent model description, we provide an overview of crucial life stages and their terminologies, adapted from Hille Ris Lambers (1966).

Fundatrix Viviparous parthenogenetic female developing from a fertilised egg developing in spring.

Virginopara Viviparous parthenogenetic female producing other parthenogenetic viviparae. Can be winged (alate virginopara) and unwinged (apterous virginopara)

Emigrant Alate virginopara that migrates from the primary host to the secondary host.

Immigrant Alate virginopara that migrates from the secondary host to the primary host.

Oviparous female/ovipara Female that lays eggs for overwintering.

Gynopara Vivipara that produces only oviparae.

Sexupara Vivipara that produces both sexual morphs (oviparae and males).

Sexuals Oviparous females and males.

Andropara Apterous virginopara that produces only males.

The general holocyclic aphid life cycle is shown in Fig. 2. Eggs hatch in early spring, giving rise to fundatrices. Fundatrices are the first generation of asexual females that give birth to nymphs, which will develop into winged (Alate) or unwinged (Apterous) asexual adults. Nymphs undergo four moults during development into adults, which will not be modelled for the sake of simplicity. Asexual reproduction continues throughout the summer, which mainly happens in crop hosts. Several environmental factors, such as crowding, host quality, and photoperiod, influence winged offspring formation. The most predominantly identified factors are crowding and host quality. As the aphid densities on a host plant increase, and host plant quality decreases, winged offspring will be formed to disperse to other hosts.

Figure 2. 

General aphid life cycle. Blue arrows represent development, red arrows stand for reproduction and green represent egg laying.

Aphids can be divided into classes when it comes to host plant choice, namely, non-host alternating (Monoecious), and host alternating (Heteroecious) species (Hille Ris Lambers 1966). A. pisum and S. avenae being monoecious and A. fabae and M. persicae being heteroecious species (Wöhrmann and Hales 1989; Emden and Harrington 2007). As autumn approaches, days become shorter and temperature declines. This triggers the formation of a single sexual generation, consisting of winged/wingless males and wingless sexual females (Ovipara) (Lamb and Pointing 1972; Tsitsipis and Mittler 1977; Hand and Wratten 1985; Margaritopoulos and Tsitsipis 2002). The sexual females produce fertilised eggs, which undergo diapause during the winter and hatch in early spring.

Winged asexual aphids (alate virginoparae) are produced due to several ecological conditions, e.g., most importantly, a decrease in host plant quality and crowding, as a means of dispersal and survival (Sutherland 1969). Winged aphids will fly from the current host in search of a new suitable host when the weather allows. Winged aphid morphs are less fecund than their wingless counterparts, highlighting the trade-off between dispersal and fecundity in polymorphism in aphids (Wratten 1977). Sexual morphs (Males and females) are produced in autumn in response to environmental cues, such as decreasing temperatures and short-day conditions (Tagu et al. 2005). This response to photoperiod is suppressed by a transgenerational timer, thereby avoiding sexual morph production in spring (Matsuda et al. 2017). Lees (1989) recorded male production, for A. pisum clones, at photoperiods 13–14 hours of daylight and ovipara production at 0–13 hours of daylight. Hand and Wratten (1985) recorded sexual morph production for brown S. avenae clones at photoperiods < 14 h. Both sources note a high variability between clones in the response to photoperiods and their progeny sequence.

The Development describes the rate at which nymphs develop into adults, how long eggs take to develop (see Table 2) and adult life expectancy. The development time is assumed to be linearly dependent on temperature; therefore degree days (DD) are used to describe these development times. In addition to that development of different morphs (winged morphs and sexual morphs) is determined by environmental factors. A thermal constant (k) is used to determine development times for nymphs and eggs. For adults, longevity (Age-specific survival) is determined by a degree-day model, assuming a linear relationship between temperature and longevity.

Eggs only develop after exposure to cold temperatures, preventing early hatching (Bonnemaison 1951). During the winter, eggs enter diapause, a period of slow, temperature-independent development, followed by a temperature-dependent development (Shingleton et al. 2003). Complete egg development takes approximately 100 days until it hatches (Miura et al. 2003).

Aphid nymph development has been shown to linearly correlate with temperature, within the lower and upper development threshold (LDT/UDT, in °C) (Bieri et al. 1983; Hutchison and Hogg 1985).

Adult aphid longevity (excluding mortalities caused by abiotic/biotic factors) depends on their experienced temperature (Campbell and Mackauer 1977). No data on the life expectancy of sexual morphs could be found.

Planned implementation

Development time will be determined by accumulated degree-days (ADD) for eggs, nymphs and adults. For each day the accumulated degree days for each morph and age group is calculated based on the hourly temperatures according to equation (2). For temperatures above the turning temperature (TT), the lower development threshold (LDT) is used to calculate the accumulated degree days. To account for the non-linear portion of the development close to the threshold a turning temperature is used, giving the aphids some development close and below the LDT. For adult morphs the TT is set significantly higher than the other morphs to avoid excessive longevities which would occur by using a degree day model to track longevity.

$DD=\left\{\begin{array}{c}0,T<0\\ T-LDT,T\ge TT\\ T-0,T<TT\end{array}\right.$ (2)

As information on the egg development of A. pisum, A. fabae, M. persica and S. avenae was not found in the literature, thermal requirements of Myzocallis coryli and Dysaphis plantaginea aphid species were used as a starting point to calibrate the thermal requirements for eggs in the aphid model (Messing and AliNiazee 1991; Graf et al. 2006). Messing and Aliniazee (1991) and Graf et al. (2006) found development times of 140- and 250-degree days respectively. From this data the equation in Table 3 was derived, assuming that before 0% hatch chance before 140 degree days and a 100% hatch chance at 250 degree days.

Table 3.

Download as

CSV

XLSX

Egg development functions and thresholds implemented for the four aphid species.

Aphid	Function for hatch chance, where x is the accumulated day degrees	Lower development threshold (°C)
A. pisum	0.0029x − 0.16	5
A. fabae	0.0029x − 0.16	5
M. persica	0.0029x − 0.16	5
S. avenae	0.0029x − 0.16	5

For nymphs, the four instars will not be modelled but will be grouped into one called ‘nymph’, representing the total development time of all instars. Data from Hutchinson and Hogg (1984) was used to calculate maximum degree days (MDD) for A. pisum nymphs by summarising the four instars maximum degree days. MDD for Apterae and Alate are based on average values from Hutchinson and Hogg (1984) (Table 4). Depending on the region where the model is applied different Lower development thresholds should be chosen. Accumulated degree days will be calculated according to equation 1 using the LDT.

Table 4.

Download as

CSV

XLSX

Development times and thresholds for the nymphal life stages of the four aphid species.

	Thermal constant (DD)	Lower development threshold (°C)	Source
A. pisum – winged	148.59+-4.48	2.73	Assumed
A. pisum – unwinged	140.25+-2.62	2.73(2.73)	(Hutchison and Hogg 1984)
A. fabae – winged	192	2.73	Assumed
A. fabae – unwinged	173+-15.72	2.73	(Akca et al. 2015)/Assumed
M. persica – winged	119.8+-1.3	4(4)	(Liu and Meng 1999)
M. persica – unwinged	133+-1.2	4(4)	(Liu and Meng 1999)
S. avenae – winged	142	2.73	Assumed
S. avenae – unwinged	128.99+-10.33	2.73 (4)	(Lykouressis 1985)

For adults, the effect of temperature on longevity within the favourable temperature range is assumed to be linear. A degree-day-based approach is used to determine the longevity of adult aphids. This physiological age is used in calculating both fecundity and longevity, as they are age dependent. Longevity curves from published studies were utilised to model species’ survival over time. These curves depict the proportion of an initial population that remains alive at specific intervals. We transformed the survival percentages into corresponding daily mortality probabilities to align these curves with our 25-degree day (DD) intervals. These probabilities reflect the daily likelihood of death for individuals within each age group. Then, at each 25-DD interval, we applied these mortality probabilities to the respective age group’s remaining population to calculate the fraction that would succumb. The survivorship curves forming the basis for the mortality rates are shown in Fig. 3.

Figure 3. 

Adult survival curves implemented for the four aphid species. Survival rate as a fraction (0–1) surviving of the initial population (Bieri et al. 1983; Sengonca et al. 1994; La Rossa et al. 2013; Akca et al. 2015).

Mortality describes how many individuals die daily due to biotic and abiotic factors. Mortality can be caused by extreme temperatures, extreme weather conditions, a decrease in plant quality, predation (e.g., Coccinellidae), diseases (e.g., fungal pathogens), and parasitism (e.g., parasitoid wasps). Harrison and Barlow (1973) reported LT50 at −9 °C and 34.5 °C to be 784.6 min and 693.6 min, respectively, for first-instar A. pisum nymphs. They also noted that temperature tolerance increased with increasing aphid age. Extreme weather conditions (apart from extreme temperatures) are not accounted for in the model, as the effect of weather on aphid population dynamics is more indirect. A decrease in host quality has been shown to increase mortality for the peach aphid Myzus persicae (Kift et al. 1998). Mortality by predation can be described in two ways, i.e., generalist and specialist predators. Generalist predators act as a constant rate of mortality since generalist predator populations do not depend on a particular prey species and are, therefore, not linked in their population dynamics. Specialist predators act as a delayed density-dependent mortality factor on the population. Snyder and Ives (2003) conducted a field study of pea aphids on alfalfa to determine a model for mortality rates caused by biotic factors (host quality, specialist predator, generalist predator, disease). Although they found the effects to be additive in their study, they suggested that if the study had been continued for a longer period, this would have changed. The overall impact of fungal pathogens on aphid population dynamics depends on the pathogen’s transmission rates and the host’s susceptibility. Transmission is influenced by a variety of factors including temperature, plant chemical production, light and most critically humidity (Steinkraus 2006). Mortality due to fungal pathogens will be indirectly modelled by a model adapted from Snyder and Ives (2003). Temperature, plant chemical production, light and humidity are not considered, as the model will only depend on aphid density.

Planned implementation

Egg mortality is assumed to be constant 70% mortality rate for all aphid species, based on the data presented in Leather (1993). This value was transformed into a daily mortality rate assuming an average egg development period of 100 days, which results in those eggs having a daily chance of 1.19% +-0.077 to die in the model. The downside of this approach is that egg survival rate in effect will be determined by average winter/spring temperatures, as low average winter/spring temperatures will cause the egg to develop for a longer period, thereby increasing the total amount of egg mortality. In contrast, higher winter temperatures would cause less egg mortality in addition to the earlier hatch date.

Nymph and adult morphs experience density-dependent and density-independent mortality. Density-independent mortality is applied as constant background mortality to simulate the effect of generalist predators. The effect of generalist predators is simulated by a constant 2.5% and 2% mortality rate for A. pisum /S. avenae and A. fabae/M. persicae respectively. Delayed density-dependent mortality is simulated by equation 3, adapted from Snyder (2003) as shown in Fig. 4. This simulates the delayed density-dependent mortality due to disease and decrease in host quality due to crowding and specialist predators. The parameters used for the four aphid species are listed in Table 5.

Figure 4. 

Delayed density dependent mortality model applied at 1200 gBiomass/m² (Snyder and Ives 2003).

Table 5.

Download as

CSV

XLSX

Density dependent mortality parameters for the four aphid species.

	Sa	k
A. fabae	0.9	0.02
A. pisum	0.9	0.012
M. persicae	0.9	0.02
S. avenae	0.9	0.012

m = Sa (1 + k * x (t − 4))⁻¹ (3)

Where:

m = Delayed density-dependent mortality (0–1).

Sa = Density independent survivorship.

k = Density dependence.

x = Aphid density in aphids per gram green biomass, at four days previously represented by (t-4)

To simulate the effect of specialist predators on the adults, a sub-model for the parasitoid Aphidius smithi is suggested. A. smithi will be modelled in a similar but simpler framework as the aphid model. The model A. smithi has two morphs; egg and adult. Fig. 5 shows the general framework of this sub-model and its interaction with the pea aphid model. For simplification, each grid cell within the simulation contains 10 A. smithi eggs as a starting point. When aphids enter a new cell, A. smithi egg development starts. A. smithi adults have a fixed lifespan of 7 days. Eggs are laid based on the ratio between parasitoids and aphids, the number of eggs in each cell is limited by the number of suitable aphids, as only 1 egg can develop in each aphid. When parasitoid eggs reach the mummification stage (73 DD) adult aphids (All morphs) aged 64+-15 DD are killed within the cell. It is assumed that each morph is equally susceptible to mortality by parasitoids. When A. smithi eggs have accumulated 180DD above 6.1 °C they become adults.

Figure 5. 

Diagram showing the interactions between the A. smithi sub-model and the aphid models.

Equation 4 was created based on data from Mackauer (1983) to simulate parasitoid fecundity at different aphid to parasitoid ratios.

F = 0.00003422F_max * x² + 0.00387022F_max * x + 0.08445296F_max (4)

Where:

F = Parasitoid fecundityx = Aphid to adult parasitoid ratio

F_max = Maximum parasitoid fecundity.

Aphid fecundity depends on age of the individual, temperature, nutrition, parasitism and size of the aphid (Hagen and Van Den Bosch 1968; Ratte 1984; Awmack and Leather 2002). The effect of parasitism and aphid size is not considered in the model, as individual aphids are not modelled. Winged aphid morphs have been shown to be less fecund then their unwinged counterparts, highlighting the trade-off between dispersal and fecundity in polymorphism in aphids (Davis 1915; Wratten 1977). Aphids reach their peak fecundity quickly after reaching adulthood. Age-specific survival first decreases after the majority of their fecundity is reached (Bieri et al. 1983).

Data on sexual morph reproduction rates is rare. An overview of the maximum fecundities of the four modelled species in presented in Table 8. Leather (1993) noted that laboratory data obtained most likely is an overestimation of the in-field fecundity and that, in effect, the average egg production is closer to half the maximum fecundity (Leather 1993).

Movement can be divided into long-distance (flying) and local movement. Local occurs if the local cell population size surpasses that of its surrounding cells, up to 2.5% of the local population being able to disperse to one of the neighbouring cells per day. The extent of dispersal is dependent on the difference of the population sizes between cells. Long distance movement is undertaken by winged adult morphs for all species. Winged morphs have a daily chance to fly at a given windspeed and have a set distance they can move (usually long distance). Wind speed is an important fact for flying since aphids are weak fliers that are not able to fly against wind direction when the wind speed is greater than 0.5 m/s (Loxdale et al. 1993).

Planned implementation

The landing of the flying winged adults will be controlled by a 2D landing mask in the subpopulation method. A brief introduction to the creation of the landing masks is included here. The detailed information can be found in the subpopulation paper (Topping and Duan 2024b).

First, two polynomial functions are used to generate a long tail landing curve along the wind direction, as shown in equation 8.

$p\left(l\right)=\left\{\begin{array}{l}-l^2+2d_{p}l,\text{ if }l\le d_p\\ \frac{d_p^2{\left(d_l-l\right)}^2}{{\left(d_l-d_p\right)}^2},\text{ if }l>d_p\text{ and }l\le d_l\end{array}\right.$ (7)

where l is the distance from the landing position to the departure position. In this way, position at d_p will have the largest landing proportion compared to other positions. The flying aphids can land at maximum distance of d_l.

After obtaining the landing curve along the wind direction, we need to convert it to a 2D mask in order to land the flying ones on a 2D landscape. Firstly, we rotate the 1D landing curve for 360 degrees through the departure location. Suppose the departure location is at (0,0), the rotation can be done by the equation 9.

$l=\sqrt{x^2+y^2}$ (8)

After the rotation, we can get the initial symmetric 2D landing mask to the original departure location using the equation 10.

$L^i(x,y)=p\left(\sqrt{x^2+y^2}\right)$ (9)

In the next step, we need to skew the mask based on the wind direction. To do this, the cosine value of the angle between the direction from the departure to the destination landing location and the wind direction is used as a weight.

Suppose the wind direction to the departure location is given by (b, −a), the direction from the destination location to the source location is given by (y, −x). Then the weight used to skew the initial landing mask can be calculated by the equation 11.

$w_c(x,y)={\left(\frac{by+ax}{\sqrt{x^2+y^2}\sqrt{a^2+b^2}}\right)}^n$ (10)

where n is a integer number to control the areas of the landing, when n is larger more flying ones will land along the wind direction. Afterwards, the 2D mask with the wind direction of (b, −a) can be calculated by the equation 12, which will be used to the calculation of the landing proportion. Examples of the resulting landing masks at different wind speeds are shown in Fig. 6.

Figure 6. 

Examples of the generated 2D landing mask with colour bar for landing proportions. Winged adults fly from (0,0). The wind direction is along the y-axis (from bottom to top) (Duan and Topping 2024).

$\mathrm{L}(\mathrm{x},\mathrm{y})=\left\{\begin{array}{cl}{L}^i(x,y)w_c(x,y),& \text{ if }{w}_c(x,y)>0\\ 0,& \text{ else }\end{array}\right.$ (11)

When a grid is with winged aphids and the wind speed is below the threshold, the winged aphids will fly to other grids based on the proportion calculated used equation 11. If the landed grid is with host plants, they will start to form a colony in the new location otherwise they will die.

Different aphid populations from the same species have been shown to specialise at particular host plants. Individuals collected from a particular host plant showed higher survival and fecundity rates on the “home” host plant (Via 1991a, b; Hufbauer and Via 1999; Ro and Long 1999; Ferrari et al. 2006). However, this intra-species specialisation is not taken into account in the simulation. Every simulated aphid population can colonise each available host plant equally well. A. pisum feeds primarily on a various leguminous plants (Fabaceae), notably peas, beans, vetches, and clovers, as detailed in Table 9. Its significant host crops include peas, beans, and fodder plants like clovers. A. fabae is host alternating between secondary hosts such as beans, tomatoes, potatoes, cucumber, beets, and their primary host Euonymus europaeus, where it overwinters. A detailed list of identified host plants is given in Table 9. M. persicae has the largest range of secondary host plants of the four modelled species. M. persicae overwinters on peach and migrates to a wide range of secondary hosts, including peas, beans, sugar beet, potatoes, cabbage, and tobacco. A detailed list of identified host plants is presented in Table 10. S. avenae spends its entire life cycle on Gramineae, colonising mainly cereal crops such wheat, oats, barley, rye and maize. A detailed list of identified host plants is presented in Table 11.

The aphid model, like any other, is a simplification of what happens in the real world. However, we do not believe that a model as detailed as this is linked to a highly realistic landscape simulation exists for aphids. As such, we hope it will be a step forward in the simulation of integrated pest management. We discuss the main exclusions in the “Framing the Model” section. These exclusions point towards areas where improvements can be made in future versions. Of these, the simplifications relating to the effect of different host plants and their distribution in the landscape are probably the most important. These simplifications have the potential to alter the emergent population patterns quite markedly. Variation within the aphids is known to occur. Not all of the modelled aphid species have the same fecundity, development thresholds, development time and longevity on the same hosts, and ecological specialisation (Frantz et al. 2006; Ferrari et al. 2008) is part of the aphid biology which we do not represent.

A notable simplification in our model relates to the assignment of specific host plants to landscape elements (tole/toc) within ALMaSS. While this method is effective for certain species and crops, like peas for A. pisum and cereals for S. avenae, it presents challenges for others, such as A. fabae and M. persicae. Designating generic winter hosts can lead to an unrealistic uniformity in host distribution across landscapes, diminishing the ecological significance of spatial host distribution—a critical factor in the long-term dynamics of A. fabae (Cammell et al. 1989).

Despite limitations, the new subpopulation model approach makes it tractable to model aphids both locally and dynamically in detail and link the population in space. Hence, we hope the trade-off with some levels of realism will lead to a valuable and applicable model. In building the formal model, some limitations to knowledge were noted. For instance, even if it were possible to include the impact of ecological specialisation on different hosts (Ferrari et al. 2008), there would be insufficient data available to parameterise this. Data on pea aphid performance on different host plants (specialised or not) is somewhat sparse. Although there is plenty of laboratory data on some of the identified host plants, most come from different studies, meaning different aphid clones, geographic locations, development thresholds, temperature optimums, and sometimes different experimental conditions. Thus, comparisons are difficult. Furthermore, the role of non-crop areas within our model, treated primarily as sources for spring migrants and overwintering sites, introduces a level of abstraction that may not fully capture the ecological nuances of these habitats. For all four species considered, this approach results in a homogenised performance of non-crop areas across landscapes, potentially oversimplifying the complex interactions and spatial dynamics essential to understanding long-term population trends, especially notable for the host alternating species A. fabae and M. persicae.

In addressing the complexities of aphid population dynamics, the new subpopulation approach offers a detailed, although abstract, perspective on spatial and temporal interactions. This approach introduces certain compromises in representing ecological nuances within the modelled species. Many other aphid models, such as those made by Ro and Long (1999), Carter (1982), improved upon by Duffy et al. (2017), adopt a more focused approach by examining population dynamics within a single host context. This narrower scope allows these models to avoid the complexities introduced by spatial dynamics. Others use simpler models but include space. For example, Damgaard Bruus, and Axelsen (2020) created a spatio-temporal stochastic model, calibrated with empirical data through a Bayesian hierarchical fitting procedure. They found low spatial variation in initial occurrence within a 10 km radius, suggesting that for agricultural decision making, localised observations could suffice for predicting outbreaks. Parry et al. (2006) created a spatially explicit agent-based model for Rhopalosiphum padi, aligning more closely to our approach, considering spatial dynamics in a 5 × 5 km landscape. Such models, including ours and Parry et al.`s, offer a framework to study complex ecological systems in detail instead of simply serving as a decision-making tool. These spatially explicit models facilitate a better understanding of the interactions and behaviours within aphid populations and their broader environmental context, such as farm management and predators, contributing to theoretical knowledge and practical pest management strategies.

This formal model paper describes the starting point for the next model development stage: implementation and calibration. The model will be calibrated using data for S. avenae from Carter (1982), Skirvin (1997), Plantegenest et al. (2001), Saussure et al. (2024) and Hasken and Poehling (1995). For the early stages of evaluation and parameterisation of A. pisum, we will rely on findings from Bommarco and Ekbom (1995) and Hutchinson and Hogg (1985). Calibration for A. fabae will largely be based on the work of Way (1966). Ro and Long’s study (1999) will be used for calibrating the M. persicae model. In early parameterisation, the focus will be put on general population dynamics, e.g., the hatching of eggs in spring, population peak, the timing of crop colonisation, production of sexual morphs, egg production in autumn and the overall spatial population dynamics of the aphid population.

However, a key challenge for calibrating this model is that data on aphid population dynamics on a large landscape scale are lacking. Most studies focus on the population growth dynamics of a single host, primarily under laboratory conditions. This lack of data means that mechanisms will be easier to calibrate than systems’ responses.

We would like to thank Jordan Chetcuti for the valuable writing advice and Oskar Alexander Thøger Rosén for the initial work on black bean aphids.