- Email: [email protected]
Generalized aphid population growth models with immigration and cumulative-size dependent dynamics
Mathematical Biosciences 215 (2008) 137–143
Contents lists available at ScienceDirect
Mathematical Biosciences journal homepage: www.elsevier.com/locate/mbs
Generalized aphid population growth models with immigration and cumulative-size dependent dynamics James H. Matis a,*, Thomas R. Kiffe b, Timothy I. Matis c, C. Chattopadhyay d a
Department of Statistics, Texas A&M University, College Station, TX 77843-3143, USA Department of Mathematics, Texas A&M University, College Station, TX 77843-3368, USA c Department of Industrial Engineering, Texas Tech University, Lubbock, TX 79409-3061, USA d National Research Centre on Rapeseed-Mustard (ICAR), Sewar, Bharatpur 321303, India b
a r t i c l e
i n f o
Article history: Received 12 February 2008 Received in revised form 10 July 2008 Accepted 11 July 2008 Available online 31 July 2008 Keywords: Logistic model Mustard aphid Power-law dynamics Temperature-dependent reproduction
a b s t r a c t Mechanistic models in which the per-capita death rate of a population is proportional to cumulative past size have been shown to describe adequately the population size curves for a number of aphid species. Such previous cumulative-sized based models have not included immigration. The inclusion of immigration is suggested biologically as local aphid populations are initiated by migration of winged aphids and as reproduction is temperature-dependent. This paper investigates two models with constant immigration, one with continuous immigration and the other with restricted immigration. Cases of the latter are relatively simple to fit to data. The results from these two immigration models are compared for data sets on the mustard aphid in India. Ó 2008 Elsevier Inc. All rights reserved.
1. Introduction Aphids are a group of small, sap-sucking insects which are serious pests of agricultural crops around the world. The annual worldwide economic loss on food and feed grains due to aphids is estimated to be $5 billion [1], with every crop-producing country having aphids as occasional pests producing yield losses of between 1% and 18%. The principal economic losses are on wheat, rice, corn, and sorghum. Even aside from their huge economic impact, ‘‘Aphids have fascinated and frustrated man for a very long time . . . mainly because of their intricate life style in close association with their host plants, the polymorphism and ability to reproduce both asexually and sexually” [2, p. 1]. Some of their remarkable life cycle is briefly outlined subsequently. We have recently developed a new class of mechanistic models for aphid population growth, in which the population growth rate depends upon the ‘cumulative-size’ of past generations. The first such model [3] with ‘ordinary’ kinetics has an analytical solution which is symmetric and relatively easy to fit to data. We have shown that this solution fits diverse aphid data very well, in particular we have fitted it successfully to large data sets on the pecan, mustard, cotton and greenbug aphids [4–7], respectively. Notwith-
* Corresponding author. Tel.: +1 979 845 3187; fax: +1 979 845 3144. E-mail address: [email protected] (J.H. Matis). 0025-5564/$ - see front matter Ó 2008 Elsevier Inc. All rights reserved. doi:10.1016/j.mbs.2008.07.007
standing the practical utility of the simple symmetric solution, there is strong biological rationale to argue that aphid population growth curves should be skewed. Hence a second mechanistic model with ‘power-law’ kinetics was developed [8] which has a left-skewed population size solution. This solution tends to fit aphid data a bit better, however the solution is not available analytically and hence is far less user-friendly than the symmetric solution. Though these existing models have been successful in describing diverse aphid population data sets, it is nevertheless apparent in some other more extensive data sets that immigration is also present. This is particularly evident in some data for the mustard aphid in India, in which a relatively long period of population size initiation is manifested. This paper extends the previous cumulative-size dependent models to include immigration. Section 2 reviews the cumulative-size dependent mechanistic models without immigration, focusing on their biological basis. Section 3 discusses the case for including immigration in such models. We will emphasize in some detail the case for immigration for the mustard aphid, though immigration is arguable present for all aphid species. Sections 4 and 5 present several cumulative-size based models with immigration. The models differ in the assumed duration of immigration and also on their ease of application. The models are illustrated by fitting them to some exceptional data sets on the population growth curves of the mustard aphid. The results from these fittings are compared in Section 6, and concluding remarks are given in Section 7.
138
J.H. Matis et al. / Mathematical Biosciences 215 (2008) 137–143
2.2. Solutions to cumulative-size models
2. Review of cumulative-size dependent models without immigration In this section we review first the mechanistic basis and then the solution of these models. We let NðtÞ and N 0 ðtÞ denote the population size variable and its derivative at time t in all subsequent work. 2.1. Mechanistic basis of models The cumulative-size based models are based on three biological principles of aphid population dynamics [8]. The first two, which lead to the basic model, are: Principle 1. Aphids are prolific reproducers with a simple, efficient reproductive system. The aphids are parthenogenetic, which is apparently unique among insects, and viviparous. This asexual reproductive mechanism gives the aphid a prolific birth potential which helps explain how they can quickly appear in large numbers. Hence, we assume the straightforward linear population birth rate, kN, which is common in population modeling. The per-capita birth rate, k, is often called the ‘intrinsic rate of natural increase.’ Principle 2. Aphid population growth is constrained by the cumulative-size of the past population. Aphids suck sap out of plants, and secrete ‘honeydew’ which covers the leaves and reduces the food habitat resources of their offspring. Current aphids also degrade the environment for their offspring by, among other factors, attracting predators, increasing the likelihood of plant and aphid disease, and producing a chemical response of the plant to a loss of sap. We assume that the extent of past environmental degradation is correlated with the ‘cumulative-size’, FðtÞ, of past aphid populations, defined as
FðtÞ ¼
Z
t
NðsÞds:
ð1Þ
0
We also assume that the per-capita death rate is proportional to this measure, FðtÞ, of cumulative damage, with proportionality parameter d. Hence the population death rate is assumed to be dFN. The basic mechanistic model combines these two principles to yield model
N0 ¼ ðk dFÞN;
ð2Þ
where k, d > 0 are called the birth and death rate parameters. Principle 3. Aphids tend to migrate elsewhere (suddenly) when their local resources are depleted. Aphids exist in a comparatively ‘hostile’ environment and so the adult aphid for many species has evolved with fully developed wings. As the environment on a leaf becomes inhospitable, due to the environmental degradation from its predecessors, the adult aphids fly away seeking a more suitable environment. The aphids tend not to leave a leaf until the environment has degraded below the aphid limits of tolerance, which results qualitatively in a rather sudden, rapid population decline. This principle leads to the power-law model developed in [8], defined as 0
p
N ¼ ðk dF ÞN
ð3Þ
with p > 1. This model is shown to have a left-skewed solution. Stochastic formulations of (2) and (3) are also given in [8]. An alternative form of the model is based on CðtÞ, the cumulative aphid count, instead of FðtÞ, the cumulative density.
Kindlmann [9,10] is apparently the first to formalize models equivalent to Eqs. (2) and (3) mathematically. The analytical solution to model (2)
NðtÞ ¼ aebt ð1 þ debt Þ2 ;
ð4Þ
where a; b, and d > 0 are general regression function parameters, was proposed by Prajneshu [11]. Matis et al. [4] shows that the peak, N max , and time of peak, t max , are
tmax ¼ b
1
log d; and Nmax ¼ a=4d:
ð5Þ
These are used to reparameterize (4) as
NðtÞ ¼ 4N max ebðttmax Þ ½1 þ ebðttmax Þ 2 :
ð6Þ
Regression function (6) is user-friendly, as its parameters are interpretable biologically and as the function is also relatively easy to fit to data due to reduced statistical multicollinearity between parameters. Defining
d ¼ expðbt max Þ;
ð7Þ
the mechanistic parameters in (2) in terms of the regression parameters in (6) are
k ¼ bðd 1Þ=ðd þ 1Þ 2
d ¼ b =2Nmax
ð8Þ 2
N0 ¼ 4dNmax ð1 þ dÞ ; where N o is the initial value Nð0Þ. We noted recently [12] that solution (6) could be expressed as 2
NðtÞ ¼ Nmax sech ½bðt t max Þ=2;
ð9Þ
and hence NðtÞ could be expressed as a scaled logistic probability density function, which is symmetric with many known properties [13]. The per-capita birth rate for this function is
N0 =N ¼ b tanh½bðt t max Þ=2:
ð10Þ
As noted previously, the three-parameter regression function (6) has been fitted successfully to large data sets on a number of aphid species, and, due to its mechanistic basis, we are confident that it would describe the population growth curves of aphid species (of which there are over 4000) in general. Though some extensions, including the incorporation of immigration in this paper, have been developed which fit selected data sets better, function (6) has been shown consistently to provide a user-friendly first approximation for all aphid growth curves that we have investigated. Nevertheless, the practical value of function (6) is not merely in describing growth curves successfully, but rather in ‘explaining’ their underlying dynamics through the estimated birth and death rate parameters in (8). The use of these parameters as response variables opens up new avenues for analyzing designed experiments involving aphid management practices, (see e.g. [4,6]). The power-law model (3), as a contrast, does not have an apparent analytical solution for p > 1. The models are solved, instead, using numerical integration. Though the model has been found to provide improved fittings for many aphid data sets, it has not been used as widely in practice due to its more challenging practical implementation. 3. The case for including immigration in cumulative-size dependent models 3.1. Classic definition of immigration ‘Immigration’ in classic population dynamics usually refers to the increase (or ‘recruitment’) in a population which occurs inde-
139
J.H. Matis et al. / Mathematical Biosciences 215 (2008) 137–143
A 600 500 400 Aphids
pendently of current population size. In the simplest case, it occurs at some time-independent rate, say m, and hence at elapsed time t its cumulative effects would be to increase population size by the simple linear function mt. As a contrast, classic ‘birth’, as in Eqs. (1) and (2), is dependent upon population size and its contribution would be to increase the population size by some exponential function, e.g. by expðatÞ.
300
3.2. Theoretical case for immigration
200
We review first the case for including immigration in the model based on the underlying population dynamics. We will detail their application to the mustard aphid though, as indicated later, these principles apply to many other aphid species. The two new principles are:
100
Principle 4. For many aphid species, local populations are initiated by the migration of winged aphids from other areas.
B 600
0
10
20
30 Day
40
50
60
0
10
20
30 Day
40
50
60
500 400 Aphids
The immigration of wind-borne winged mustard aphids is a common seasonal phenomenon in parts of India, which locals sometimes term as ‘raining of aphids’. Previous studies in India ([14,15]) have indicated that the immigration of the mustard aphid could be as long as 300 km, which contributes to the build-up of its population at a particular site. Extensive nocturnal migrations at 150 m height are also evidenced [14].
300 200 100
Principle 5. Aphid reproduction (birth) is temperature dependent, and is virtually or totally absent below certain threshold temperatures. It would appear from subsequent data analysis, for example, that the temperatures often were not conducive to full reproduction until sometime in January, four weeks after aphid initial arrival. Mustard aphid appearance on sticky traps in experiments in India happened within a range of maximum temperature of 20–29 °C in the preceding week, more so when the maximum temperature was 22–25 °C [16–18]. It is due to the dependence of reproduction on temperature that, after initial appearance on sticky traps, the population did not rise rapidly till the temperature reached within the above range, when the population jumped suddenly from traces to high numbers. Such effect of temperature on the catch in traps for other aphid species have also been reported ([19–21]). In summary, the known dynamics of the mustard aphid life cycle provide strong justification for incorporating an immigration component into the model. 3.3. Empirical case for immigration from data We suggest that empirical evidence of initial immigration would consist of a simple linear increase in the population size curve, NðtÞ. This signature trend is often clearly noticeable in mustard aphid population data when there is a relatively long initial infestation period. Fig. 1 illustrates an exceptional data set in which the mean number of aphids in five sticky traps was observed daily from January 26 to March 28, 2002 at an experimental plot in Bharatpur (at the National Research Centre on Rapeseed-Mustard) in Northern India. It seems clear by inspection that the period of initial immigration lasts approximately 31 days. More generally, this daily mustard aphid data set suggests that a pronounced shift in process dynamics occurs at around day 31, and hence that neither the ordinary birth model (2) nor the power-law model with p > 1 in (3) is adequate for describing this detailed curve. Most of the aphid abundance data which we have seen has been gathered with weekly sampling intervals. Immigration is often visible in such data, though not as remarkably as with daily data. For example, Matis et al. [5] analyzes 13 individual observed data sets, which were gathered on the abundance of the mustard aphid over the period 1983–1990 on three different oilseed cultivars in North-
Fig. 1. Observed mean mustard aphid counts with daily sampling interval, and fitted curves. (A) Fitted model (6) without immigration (solid line) and model (11) with immigration (dashed line). (B) Fitted model (6) without immigration (solid line) and model (13) with immigration (dashed line).
ern India (at the Haryana Agricultural University). The data were gathered weekly, with the first aphid occurrence in late December/early January and a peak aphid count in late February. The symmetric function (6) is shown to provide a fit which is adequate for subsequent data analysis for all 13 cases. However, six of the 13 cases have a clearly visible initial simple linear increase, which leads to obvious lack-of-fit for function (6). The parameter estimates for regression function (6) and the root mean squared residuals, s, are given for each case in Table 1. These six cases are illustrated in Fig. 2. It appears that in each of these six cases there is a period of initial immigration, of four to six weeks duration, during which period the residuals are predominately positive, as the observed values tend to exceed the estimated values from regression function (6). 3.4. Summary of case for immigration Because local infestation is initiated by wind-borne (immigrant) aphids, we suggest that an immigration component is always part
Table 1 Estimates of regression parameters (N max , tmax , b) in model (6), mechanistic parameters (k; d; N 0 ), residual sum of squares and root mean squared error, ðRSS; sÞ, for 6 mustard aphid data sets Case
N max
tmax
b
k
d 103
N0
RSS
s
J8 C8 C9 C0 N7 N8
96.6 985.4 59.7 148.7 171.6 113.7
6.42 5.71 6.57 7.31 7.41 7.76
0.933 2.413 0.620 1.324 1.124 0.764
0.928 2.413 0.603 1.324 1.123 0.760
4.51 2.95 3.25 5.89 3.68 2.57
0.96 .004 3.83 .037 .166 1.20
960 142196 4382 1564 2409 2237
11.7 168.6 27.0 16.1 20.0 17.9
140
J.H. Matis et al. / Mathematical Biosciences 215 (2008) 137–143
Aphids
Aphids
C8 1400 1200 1000 800 600 400 200 6 8 Week
10
12
2
N8
200 175 150 125 100 75 50 25 2
4
6 8 Week
10
2
4
6 8 Week
10
4
10
12
6 8 Week
10
12
10
12
C0
200 175 150 125 100 75 50 25
12
6 8 Week C9
2
J8
200 175 150 125 100 75 50 25
4
200 175 150 125 100 75 50 25
12
Aphids
Aphids
4
Aphids
Aphids
2
N7
200 175 150 125 100 75 50 25
2
4
6 8 Week
Fig. 2. Curves for six mustard aphid data sets, with weekly sampling intervals, fitted to model (6) without immigration (solid line) and to model (11) with immigration (dashed line).
of the local population dynamics of the mustard aphid. However, the duration of the simple linear increase in population size attributable to immigration only, without reproduction, is dependent upon the ambient temperature. Examining all 13 mustard aphid data sets in Matis et al. [5], all four data sets for years 1988 and 1990 have simple linear phases of at least four weeks duration, indicating apparent low temperatures during this period of those years. On the other hand, none of the four data sets collected during 1983–86 have an apparent simple linear phase, perhaps due to higher temperatures during this period of those years. The data sets for 1987 and 1989 are mixed. In summary, we argue that all aphid abundance curves have an initial simple linear phase, though the duration is often so short that the symmetric function (6) may be sufficient for all practical purposes for describing weekly data with short infestation periods. At other times, however, particularly for daily data, the simple linear phase is striking. This motivates a model generalization for the mustard aphid that includes immigration. Empirical evidence for including immigration in the population models is also evident in our data for the pecan and for the cotton aphid in Texas, however the infestation period is shorter and the case is not as striking. 4. A cumulative-size dependent model with continuous immigration Though the case for immigration is arguably very conclusive, the precise ‘real-world’ description of the immigration effect is no doubt a complicated function of numerous weather and environmental variables. Our ultimate objective is, however, to fit our population models to field data, and hence we require a simple, robust mechanism, involving only a few parameters, for incorporating the immigration effect. Clearly, simulation modelers who
deal with more complex system models would not need to be as parsimonious with parameters. In this spirit, we propose two models which incorporate immigration. Both of these assume constant (time-homogeneous) immigration of size m insects/time interval. Our first model assumes that this effect lasts continuously throughout the period of observation. A second model, developed in the next section as an extreme contrast, terminates the immigration effect abruptly at some fixed time. In the first model, we assume that the cumulative-size dependent growth effect in (2) is added to the immigration effect after some time lag s. These assumptions yield model
N0 ¼
m t
ð11Þ
with m; s; k, and d P 0. This sudden start of the cumulative-size growth is also obviously an oversimplification. However, it seems appropriate as a first approximation for generalizing (2), and paraphrasing Kac [22], the model serves as a ‘caricature of reality’ which should portray some features of the ‘real-world’. This model was fitted to the observed mustard aphid data using numerical integration programs in EasyFit [23]. For simplicity, the initial value was set to 0, i.e. Nð0Þ ¼ 0, hence the model had four parameters compared to the three in model (2) without immigration. Fig. 1A illustrates the fitted curves for the daily aphid counts. The parameter estimates for the model are m ¼ 7:52 aphids/day, s ¼ 31:4 days, k ¼ 1:116/day and d ¼ 2:164 104 /aphid-day2. The sum of squared residuals, RSS, used as a measure of goodness-of-fit is RSS ¼ 51; 294. As a contrast, we also fitted model (2) without immigration. Its parameter estimates are k ¼ 0:2289/ day, d ¼ 6:764 105 =aphid-day2, and Nð0Þ ¼ 0:79 aphids, with
J.H. Matis et al. / Mathematical Biosciences 215 (2008) 137–143 Table 2 Estimates of mechanistic parameters (m; s; k; d) in model (11), estimated N max and t max , and residual sum of squares and root mean squared error, ðRSS; sÞ, for 6 mustard aphid data sets Case
m
s
k
d 103
N max
tmax
RSS
s
J8 C8 C9 C0 N7 N8
9.03 68.25 8.00 3.67 8.41 8.71
4.87 4.92 6.00 4.83 5.52 5.82
1.466 7.320 6.394 1.484 1.656 1.240
6.478 4.958 31.055 6.085 4.683 3.876
104.8 1426.5 112.8 149.1 182.5 124.9
6.5 5.6 6.7 7.3 7.5 7.8
268 29566 2514 1407 1153 780
6.7 86.0 22.4 16.8 15.2 11.4
RSS ¼ 211; 278. The addition of immigration obviously provides a far superior fit for these daily count data, as evident in Fig. 1A. The fitted curves of model (11) for the six data sets with weekly counts are illustrated in Fig. 2, and their parameters are given in Table 2. For convenience, we also calculate the residual standard error, s, as a measure of fit. Model (11) with immigration provides a noticeably better fitting curve than model (2) without immigration, for five of the six cases, with a reduction in s ranging from 16% to 48% for these cases. The largest reduction, 48%, occurred in case C8, which also has the largest s due to the high counts. The smallest reduction, 16%, occurred in case C9, for which the improvement due to immigration is obscured by an apparently outlying data point at t ¼ 9. Case C0 is instructive, as s did not decrease, but instead increased by 4% with the inclusion of immigration. It is clear that the fitted model (2) yields systematic lack-of-fit, as the residuals are positive for the first four time periods, t = 1–4. However, though model (11) has a lower sum of squared residuals, RSS, its residual standard error s, is slightly higher due to a loss of one degree of freedom from the residual degrees of freedom. As a rule, model (11) is expected to always provide a better fitting curve for mustard aphid data than model (2), even though s does not always decrease due to the limited number of data points.
As a contrast to model (11), we consider a model which recognizes that the immigration of aphids in India is seasonal and of limited duration. Instead of lasting for over the 10-week-period of the data collection, the changing climate conditions would tend to shift the wind-borne aphid immigration northward. We consider therefore a second model with immigration. In this model, rather than assuming continuous immigration, we make the other extreme assumption that immigration terminates at the end of the visible simple linear trend in the data, at the time therefore when reproduction starts. Hence the new model is
m ðk dFÞN
t
ð12Þ
The actual immigration effect obviously lies somewhere between the two assumed extremes in models (11) and (12). There is likely a gradual tapering off between time s and 1. However, through immigration plays a crucial role in initiating and sustaining the early part of the infestation, its effect even in model (11) is gradually swamped out after time s by the birth effect and the rapidly increasing death effect. Model (12) has a huge practical advantage due to its analytical solution. If one also assumes that the cumulative density, FðtÞ, applies only to the second, symmetric piece of the model, the model has piecewise solution
NðtÞ ¼
mt 4Nmax ebðttmax Þ ½1 þ ebðttmax Þ 2
t6s t > s:
The five parameters, m; s; N max ; t max , and b, may be estimated with standard nonlinear least squares software, e.g. we use SPSS [24]. Often s may be estimated by inspection from the data. In this simplifying case, the solution consists of two separate parts, each of which is easy to fit to data. The daily aphid data was fitted to regression model (13) with the fitted curve illustrated in Fig. 1B. The simple linear part fitted best for t < 32, with parameter estimate m ¼ 7:55 aphids/day. The cumulative-size model fitted best for t > 31, thus conveniently giving the bound 31 < s < 32 in (13). The cumulative-size model has parameter estimates N max ¼ 485:8 aphids, tmax ¼ 35:3 days and b ¼ 0:447/day. The corresponding mechanistic parameter estimates from (8) are k ¼ 0:447 and d ¼ 2:056 104 . The combined residual sum of squares for the two parts of model (13) is RSS ¼ 50; 527. This indicates a slightly better fit to model (13) than to model (11), which is not surprising due to its single additional parameter. Consider now the six data sets with weekly data. It appears by inspection, in each case, that the period with only immigration is at least four weeks but no longer than six weeks. Hence, for simplicity, the data points for t 6 5 were fitted to the simple linear equation in the first part of (13), and the data points for t P 5 were fitted to function (6) for the second part of (13). The parameter s is assumed to be the intersection of these two fitted curves. The fitted curves are illustrated in Fig. 3, and the parameter estimates are given in Table 3. The fitted curves of models (11) and (13) are compared in the next section. 6. Comparison of models with immigration It is obvious by inspection that both models (11) and (13) with immigration fit the daily data well, and both make a strong case for including immigration in the model. This section compares the two models for the six cases of weekly data. 6.1. Comparison of goodness-of-fit
5. A cumulative-size dependent model with restricted immigration
N0 ¼
141
ð13Þ
Both cases fit the data well qualitatively, as expected. Model (11) provides a slightly better fitting based on s for three of the six cases; whereas model (13) gives a somewhat better fitting for the other three. We are unable to make any generalizations on which model fits better based on these data. 6.2. Comparison of immigration components The estimated immigration rate m tended to be slightly higher for model (13) than for model (11). The estimated rate increased in five of the six cases, with the increase ranging from 1.0% to 13.8%, for an average increase of 6.0%. The estimated duration, s, was virtually identical in each case for the two models. Hence we conclude that the estimated immigration components of the two models are very close to one another. 6.3. Comparison of N max and tmax parameters The peak sizes, N max , of the two models are very close (within 2.5%) of one another except for cases C8 and C9, both of which are atypical. In case C8, N max ¼ 1; 426 for model (11) exceeds N max ¼ 1; 266 by 12.7%, however C8 has a smaller number of data points than the other cases which makes the determination of N max more challenging. In C9, N max ¼ 112:8 for model (11) exceeds N max ¼ 102:1 for model (13) by 10.5%, however this case has an anomalous observation at t ¼ 8 which makes any generalization questionable. In general, however, even with the cases with suspect and limited data, it appears that the estimated N max of model (13) is closer than those of model (11) to the ob-
142
J.H. Matis et al. / Mathematical Biosciences 215 (2008) 137–143
Aphids
Aphids
C8 1400 1200 1000 800 600 400 200 6 8 Week
10
12
2
N8
200 175 150 125 100 75 50 25 2
4
6 8 Week
10
2
4
6 8 Week
10
12
6 8 Week
10
12
10
12
10
12
C9
2
J8
200 175 150 125 100 75 50 25
4
200 175 150 125 100 75 50 25
12
Aphids
Aphids
4
Aphids
Aphids
2
N7
200 175 150 125 100 75 50 25
4
6 8 Week C0
200 175 150 125 100 75 50 25 2
4
6 8 Week
Fig. 3. Curves for six mustard aphid data sets, with weekly sampling intervals, fitted to model (6) without immigration (solid line) and to piecewise model (13) with immigration (dashed line).
Table 3 Estimates of regression parameters (m; N max , t max , b) in model (13), mechanistic parameters (k; d; N 0 ; s), residual sum of squares and root mean squared error, ðRSS; sÞ, for 6 mustard aphid data sets Case
m
N max
t max
b
k
d 103
N0
J8 C8 C9 C0 N7 N8
10.12 77.65 7.95 4.01 8.53 8.75
103.5 1265.5 102.1 150.3 178.2 123.9
6.57 5.65 6.79 7.32 7.45 7.85
1.129 3.414 2.174 1.360 1.234 0.960
1.128 3.414 2.174 1.360 1.234 0.956
6.16 4.61 23.15 6.15 4.27 3.72
2.48 2.12 1.59 2.85 7.25 2.64
E-1 E-5 E-4 E-2 E-2 E-2
s
RSS
s
4.9 4.9 5.9 4.8 5.3 5.7
162 31862 3214 1354 1241 619
5.2 89.2 25.4 16.5 15.8 10.2
served maximum in the data. The estimated tmax are very close to each other in each case.
bines with the higher birth rate in model (11) to yield a peak, N max , about the same size as in model (13) with lower rates.
6.4. Comparison of k and d rate parameters
6.5. Summary of comparison of models with immigration
Model (11) with continuous immigration always had a higher estimate of the birth rate, k, than model (13) with restricted immigration for the six cases. The increase ranged from 10% to 194%, i.e. almost tripling, for an average increase of 68% and a median of 32%. For the case with daily data in Fig. 1, the rate increased to k ¼ 1:116 from k ¼ 0:447, for a 150% increase, between the two models. This substantial increase seems paradoxical, as one might expect that the model in (13) would have a higher birthrate to compensate for the loss of immigration after s. However, model (11) also has a consistently higher estimate of the death rate d. The estimate of d was higher in model (11) than in model (13) for five of the six cases, the increase ranging from 4% to 34%. In the one case where the rate decreases, the decrease was only 1%. We suggest that conceptually the higher death rate com-
Our limited experience based on these mustard aphid data suggest that, as a rule: 1. Both models fit data about equally well. 2. Both models have similar shape characteristics, in particular regarding the size of immigration and the peak size, N max . In fact, the fitted curves are often visually indistinguishable, as they are for the daily data in Fig. 1A and B. 3. Model (11) with constant immigration tends to have higher birth and death rate estimates than model (13) with restricted immigration. If the accurate estimation of these rates is an important objective of the study, it seems crucial to seek expert subject-matter opinion regarding the likely duration of immigration, in order to decide between models (11) and (13).
J.H. Matis et al. / Mathematical Biosciences 215 (2008) 137–143
6.6. Comparison of models with immigration to model without immigration Though the models with immigration provide superior fittings in many cases, no doubt many experimenters will prefer to use the symmetric model in (6) for convenience. A question of interest then is how the estimated rates in (6) differ from those given by models (11) and (13). It is clear comparing the results of Tables 1–3, that in each case the estimated birth and death rates from the symmetric model (6) are lower than those estimated from either model (11) or (13). Hence, it seems crucial to use a model which includes immigration whenever immigration is clearly evident in the data, if one desires accurate estimates of these rates. 7. Discussion We suggest, based on purely theoretical grounds, that the initial increase in mustard aphid population size is due to immigration only. Often this initial period of growth due to immigration is sizeable, and in such cases a population model which includes a mechanism to incorporate immigration should be useful. Our limited data verify empirically that models with immigration often give superior fittings of data. We expect that these findings will be strengthened as more daily count data become available. This paper focuses on only two models with immigration. Both seem adequate for the data, however model (13) has a piecewise solution which is relatively simple to fit to data using widely available software. For the sake of completeness, we explored a number of other models incorporating immigration. One is the power-law model with immigration, which is model (11) with p ¼ 2. In general, this model did not fit any better, perhaps in part because the immigration already explained some of the skewness, which without immigration is attributed to the power-law. Another is a variation of model (11), in which the cumulative density FðtÞ is integrated from 0 to s over the immigration part rather than over the second, symmetric part. One cannot use the piecewise solution in (13) under this assumption, hence the numerical integration in EasyFit [23] was used again. This model also failed to improve the fittings, and hence was not investigated further due to its difficult practical implementation. One generalization which would be easy to implement would be to change the time-homogeneous immigration rate in (11) and (12) to some time varying function mðtÞ, which might be specified with few coefficients. In the absence of an alternative detailed hypothesis about mðtÞ, we did not explore such possibility. We note also that although the current investigation involves only data for the mustard aphid, similar initial population dynamics also apply no doubt to other aphid and non-aphid species. We leave this for others, probably with access to data with daily sampling intervals, to confirm as we do not have data necessary to verify this assertion.
143
Acknowledgement We express our appreciation to the associate editor and to the two reviewers for their useful suggestions which have improved this paper. References [1] R.L. Metcalf, R.A. Metcalf, Destructive and Useful Insects, fifth ed., McGrawHill, New York, 1995. [2] A.F.G. Dixon, Aphid Ecology, Chapman and Hall, New York, 1985. [3] J.H. Matis, T.R. Kiffe, T.I. Matis, D.E. Stevenson, Nonlinear stochastic modeling of aphid population growth, Math. Biosci. 198 (2005) 148. [4] J.H. Matis, T.R. Kiffe, T.I. Matis, D.E. Stevenson, Application of population growth models based on cumulative size to pecan aphids, J. Ag. Biol. Environ. Stat. 11 (2006) 1. [5] J.H. Matis, T.R. Kiffe, T.I. Matis, J.A. Jackman, H. Singh, Population size models based on cumulative size, with application to aphids, Ecol. Model. 205 (2007) 81. [6] T.I. Matis, M.N. Parajulee, J.H. Matis, R.B. Shrestra. A mechanistic model based analysis of cotton aphid population dynamics data, Ag. Forest. Ento. in press. [7] T.I. Matis, J.H. Matis, G.J. Michels Jr., M.N. Parajulee, Using the cumulative-size mechanistic model for analyzing insect data., Proc. 2007 Conf. Appl. Stat. Agr., Kansas State Dept. Stat., Manhattan, KS, 2008. [8] J.H. Matis, T.R. Kiffe, T.I. Matis, D.E. Stevenson, Stochastic modeling of aphid population growth with nonlinear, power-law dynamics, Math. Biosci. 208 (2007) 469. [9] P. Kindlemann, A model of aphid population with age structure, in: V. Capasso, E. Grosso, S.L. Paveri-Fontana (Eds.), Mathematics in Biology and Medicine, Proceedings, Bari, 1983, Lecture Notes in Biomathematics, Springer, Berlin, 1985, p. 72. [10] P. Kindlemann, R. Arditi, A.F.G. Dixon, A simple aphid population model, in: J.C. Simon, C.A. Dedryver, C. Rispe, M. Hulle (Eds.), Aphids in a New Millennium, INRA, Paris, 2004, p. 325. [11] Prajneshu, A nonlinear statistical model for aphid population growth, J. Indian Soc. Agric. Stat. 51 (1998) 73. [12] J.H. Matis, T.R. Kiffe, T.I. Matis, J.A. Jackman, W.E. Grant, H. Singh, On the use of growth rate parameters for projecting population sizes: application to aphids, Ecol. Model. 213 (2008) 133. [13] N.L. Johnson, S. Kotz, Continuous Univariate Distributions – 1, Wiley, New York, 1970. [14] J.R. Riley, D.R. Reynolds, S. Mukhopadhyay, M.R. Ghosh, T.K. Sarkar, Longdistance migration of aphids and other small insects in northeast India, Eur. J. Ento. 92 (1995) 639. [15] A. Srivastava, N.P. Kashyap, Population build-up of mustard aphid, Lipaphis erysimi (Kalt.) on brassica in relation to weather factors, J. Insect Sci. 12 (1999) 141. [16] S.K. Prasad, K.G. Phadke, K.N. Mehrotra, Forecasting of mustard aphid, Lipaphis erysimi Kaltenbach on rapeseed crop, Proc. Indian Natl. Sci. Acad. (B) 50 (1984) 146. [17] K.G. Phadke, Ecological factors affecting aphid Lipaphis erysimi incidence on mustard crop, in: K.G. Phadke (Ed.), Aphidology in India, A.R. Printers, Calcutta, 1986, p. 37. [18] S.K. Roy, K. Baral, Role of weather parameters on population build up of mustard aphid Lipaphis erysimi (Kaltenbach), J. Oilseeds Res. 19 (2002) 86. [19] A. Campbell, M. Mackauer, Reproduction and population growth of the pea aphid (Homoptera: Aphididae) under laboratory conditions, Can. Ento. 109 (1977) 277. [20] R.I.K. McVean, A.F.G. Dixon, T. Harrington, Causes of regional and yearly variation in pea aphid numbers in eastern England, J. Appl. Ento. 123 (1999) 495. [21] T.X. Liu, B. Yue, Comparison of some life history parameters between alate and apterous forma of turnip aphid (Homoptera: Aphididae) on cabbage under constant temperatures, Florida Ento. 84 (2001) 239. [22] M. Kac, Some mathematical models in science, Science 166 (1969) 695. [23] EasyFit, MathWave Technologies, (www.mathwave.com), Arlington, MA, 2005. [24] SPSS. SPSS 11.0 for Macintosh. SPSS Inc., Chicago, IL, 2002.
Contents lists available at ScienceDirect
Mathematical Biosciences journal homepage: www.elsevier.com/locate/mbs
Generalized aphid population growth models with immigration and cumulative-size dependent dynamics James H. Matis a,*, Thomas R. Kiffe b, Timothy I. Matis c, C. Chattopadhyay d a
Department of Statistics, Texas A&M University, College Station, TX 77843-3143, USA Department of Mathematics, Texas A&M University, College Station, TX 77843-3368, USA c Department of Industrial Engineering, Texas Tech University, Lubbock, TX 79409-3061, USA d National Research Centre on Rapeseed-Mustard (ICAR), Sewar, Bharatpur 321303, India b
a r t i c l e
i n f o
Article history: Received 12 February 2008 Received in revised form 10 July 2008 Accepted 11 July 2008 Available online 31 July 2008 Keywords: Logistic model Mustard aphid Power-law dynamics Temperature-dependent reproduction
a b s t r a c t Mechanistic models in which the per-capita death rate of a population is proportional to cumulative past size have been shown to describe adequately the population size curves for a number of aphid species. Such previous cumulative-sized based models have not included immigration. The inclusion of immigration is suggested biologically as local aphid populations are initiated by migration of winged aphids and as reproduction is temperature-dependent. This paper investigates two models with constant immigration, one with continuous immigration and the other with restricted immigration. Cases of the latter are relatively simple to fit to data. The results from these two immigration models are compared for data sets on the mustard aphid in India. Ó 2008 Elsevier Inc. All rights reserved.
1. Introduction Aphids are a group of small, sap-sucking insects which are serious pests of agricultural crops around the world. The annual worldwide economic loss on food and feed grains due to aphids is estimated to be $5 billion [1], with every crop-producing country having aphids as occasional pests producing yield losses of between 1% and 18%. The principal economic losses are on wheat, rice, corn, and sorghum. Even aside from their huge economic impact, ‘‘Aphids have fascinated and frustrated man for a very long time . . . mainly because of their intricate life style in close association with their host plants, the polymorphism and ability to reproduce both asexually and sexually” [2, p. 1]. Some of their remarkable life cycle is briefly outlined subsequently. We have recently developed a new class of mechanistic models for aphid population growth, in which the population growth rate depends upon the ‘cumulative-size’ of past generations. The first such model [3] with ‘ordinary’ kinetics has an analytical solution which is symmetric and relatively easy to fit to data. We have shown that this solution fits diverse aphid data very well, in particular we have fitted it successfully to large data sets on the pecan, mustard, cotton and greenbug aphids [4–7], respectively. Notwith-
* Corresponding author. Tel.: +1 979 845 3187; fax: +1 979 845 3144. E-mail address: [email protected] (J.H. Matis). 0025-5564/$ - see front matter Ó 2008 Elsevier Inc. All rights reserved. doi:10.1016/j.mbs.2008.07.007
standing the practical utility of the simple symmetric solution, there is strong biological rationale to argue that aphid population growth curves should be skewed. Hence a second mechanistic model with ‘power-law’ kinetics was developed [8] which has a left-skewed population size solution. This solution tends to fit aphid data a bit better, however the solution is not available analytically and hence is far less user-friendly than the symmetric solution. Though these existing models have been successful in describing diverse aphid population data sets, it is nevertheless apparent in some other more extensive data sets that immigration is also present. This is particularly evident in some data for the mustard aphid in India, in which a relatively long period of population size initiation is manifested. This paper extends the previous cumulative-size dependent models to include immigration. Section 2 reviews the cumulative-size dependent mechanistic models without immigration, focusing on their biological basis. Section 3 discusses the case for including immigration in such models. We will emphasize in some detail the case for immigration for the mustard aphid, though immigration is arguable present for all aphid species. Sections 4 and 5 present several cumulative-size based models with immigration. The models differ in the assumed duration of immigration and also on their ease of application. The models are illustrated by fitting them to some exceptional data sets on the population growth curves of the mustard aphid. The results from these fittings are compared in Section 6, and concluding remarks are given in Section 7.
138
J.H. Matis et al. / Mathematical Biosciences 215 (2008) 137–143
2.2. Solutions to cumulative-size models
2. Review of cumulative-size dependent models without immigration In this section we review first the mechanistic basis and then the solution of these models. We let NðtÞ and N 0 ðtÞ denote the population size variable and its derivative at time t in all subsequent work. 2.1. Mechanistic basis of models The cumulative-size based models are based on three biological principles of aphid population dynamics [8]. The first two, which lead to the basic model, are: Principle 1. Aphids are prolific reproducers with a simple, efficient reproductive system. The aphids are parthenogenetic, which is apparently unique among insects, and viviparous. This asexual reproductive mechanism gives the aphid a prolific birth potential which helps explain how they can quickly appear in large numbers. Hence, we assume the straightforward linear population birth rate, kN, which is common in population modeling. The per-capita birth rate, k, is often called the ‘intrinsic rate of natural increase.’ Principle 2. Aphid population growth is constrained by the cumulative-size of the past population. Aphids suck sap out of plants, and secrete ‘honeydew’ which covers the leaves and reduces the food habitat resources of their offspring. Current aphids also degrade the environment for their offspring by, among other factors, attracting predators, increasing the likelihood of plant and aphid disease, and producing a chemical response of the plant to a loss of sap. We assume that the extent of past environmental degradation is correlated with the ‘cumulative-size’, FðtÞ, of past aphid populations, defined as
FðtÞ ¼
Z
t
NðsÞds:
ð1Þ
0
We also assume that the per-capita death rate is proportional to this measure, FðtÞ, of cumulative damage, with proportionality parameter d. Hence the population death rate is assumed to be dFN. The basic mechanistic model combines these two principles to yield model
N0 ¼ ðk dFÞN;
ð2Þ
where k, d > 0 are called the birth and death rate parameters. Principle 3. Aphids tend to migrate elsewhere (suddenly) when their local resources are depleted. Aphids exist in a comparatively ‘hostile’ environment and so the adult aphid for many species has evolved with fully developed wings. As the environment on a leaf becomes inhospitable, due to the environmental degradation from its predecessors, the adult aphids fly away seeking a more suitable environment. The aphids tend not to leave a leaf until the environment has degraded below the aphid limits of tolerance, which results qualitatively in a rather sudden, rapid population decline. This principle leads to the power-law model developed in [8], defined as 0
p
N ¼ ðk dF ÞN
ð3Þ
with p > 1. This model is shown to have a left-skewed solution. Stochastic formulations of (2) and (3) are also given in [8]. An alternative form of the model is based on CðtÞ, the cumulative aphid count, instead of FðtÞ, the cumulative density.
Kindlmann [9,10] is apparently the first to formalize models equivalent to Eqs. (2) and (3) mathematically. The analytical solution to model (2)
NðtÞ ¼ aebt ð1 þ debt Þ2 ;
ð4Þ
where a; b, and d > 0 are general regression function parameters, was proposed by Prajneshu [11]. Matis et al. [4] shows that the peak, N max , and time of peak, t max , are
tmax ¼ b
1
log d; and Nmax ¼ a=4d:
ð5Þ
These are used to reparameterize (4) as
NðtÞ ¼ 4N max ebðttmax Þ ½1 þ ebðttmax Þ 2 :
ð6Þ
Regression function (6) is user-friendly, as its parameters are interpretable biologically and as the function is also relatively easy to fit to data due to reduced statistical multicollinearity between parameters. Defining
d ¼ expðbt max Þ;
ð7Þ
the mechanistic parameters in (2) in terms of the regression parameters in (6) are
k ¼ bðd 1Þ=ðd þ 1Þ 2
d ¼ b =2Nmax
ð8Þ 2
N0 ¼ 4dNmax ð1 þ dÞ ; where N o is the initial value Nð0Þ. We noted recently [12] that solution (6) could be expressed as 2
NðtÞ ¼ Nmax sech ½bðt t max Þ=2;
ð9Þ
and hence NðtÞ could be expressed as a scaled logistic probability density function, which is symmetric with many known properties [13]. The per-capita birth rate for this function is
N0 =N ¼ b tanh½bðt t max Þ=2:
ð10Þ
As noted previously, the three-parameter regression function (6) has been fitted successfully to large data sets on a number of aphid species, and, due to its mechanistic basis, we are confident that it would describe the population growth curves of aphid species (of which there are over 4000) in general. Though some extensions, including the incorporation of immigration in this paper, have been developed which fit selected data sets better, function (6) has been shown consistently to provide a user-friendly first approximation for all aphid growth curves that we have investigated. Nevertheless, the practical value of function (6) is not merely in describing growth curves successfully, but rather in ‘explaining’ their underlying dynamics through the estimated birth and death rate parameters in (8). The use of these parameters as response variables opens up new avenues for analyzing designed experiments involving aphid management practices, (see e.g. [4,6]). The power-law model (3), as a contrast, does not have an apparent analytical solution for p > 1. The models are solved, instead, using numerical integration. Though the model has been found to provide improved fittings for many aphid data sets, it has not been used as widely in practice due to its more challenging practical implementation. 3. The case for including immigration in cumulative-size dependent models 3.1. Classic definition of immigration ‘Immigration’ in classic population dynamics usually refers to the increase (or ‘recruitment’) in a population which occurs inde-
139
J.H. Matis et al. / Mathematical Biosciences 215 (2008) 137–143
A 600 500 400 Aphids
pendently of current population size. In the simplest case, it occurs at some time-independent rate, say m, and hence at elapsed time t its cumulative effects would be to increase population size by the simple linear function mt. As a contrast, classic ‘birth’, as in Eqs. (1) and (2), is dependent upon population size and its contribution would be to increase the population size by some exponential function, e.g. by expðatÞ.
300
3.2. Theoretical case for immigration
200
We review first the case for including immigration in the model based on the underlying population dynamics. We will detail their application to the mustard aphid though, as indicated later, these principles apply to many other aphid species. The two new principles are:
100
Principle 4. For many aphid species, local populations are initiated by the migration of winged aphids from other areas.
B 600
0
10
20
30 Day
40
50
60
0
10
20
30 Day
40
50
60
500 400 Aphids
The immigration of wind-borne winged mustard aphids is a common seasonal phenomenon in parts of India, which locals sometimes term as ‘raining of aphids’. Previous studies in India ([14,15]) have indicated that the immigration of the mustard aphid could be as long as 300 km, which contributes to the build-up of its population at a particular site. Extensive nocturnal migrations at 150 m height are also evidenced [14].
300 200 100
Principle 5. Aphid reproduction (birth) is temperature dependent, and is virtually or totally absent below certain threshold temperatures. It would appear from subsequent data analysis, for example, that the temperatures often were not conducive to full reproduction until sometime in January, four weeks after aphid initial arrival. Mustard aphid appearance on sticky traps in experiments in India happened within a range of maximum temperature of 20–29 °C in the preceding week, more so when the maximum temperature was 22–25 °C [16–18]. It is due to the dependence of reproduction on temperature that, after initial appearance on sticky traps, the population did not rise rapidly till the temperature reached within the above range, when the population jumped suddenly from traces to high numbers. Such effect of temperature on the catch in traps for other aphid species have also been reported ([19–21]). In summary, the known dynamics of the mustard aphid life cycle provide strong justification for incorporating an immigration component into the model. 3.3. Empirical case for immigration from data We suggest that empirical evidence of initial immigration would consist of a simple linear increase in the population size curve, NðtÞ. This signature trend is often clearly noticeable in mustard aphid population data when there is a relatively long initial infestation period. Fig. 1 illustrates an exceptional data set in which the mean number of aphids in five sticky traps was observed daily from January 26 to March 28, 2002 at an experimental plot in Bharatpur (at the National Research Centre on Rapeseed-Mustard) in Northern India. It seems clear by inspection that the period of initial immigration lasts approximately 31 days. More generally, this daily mustard aphid data set suggests that a pronounced shift in process dynamics occurs at around day 31, and hence that neither the ordinary birth model (2) nor the power-law model with p > 1 in (3) is adequate for describing this detailed curve. Most of the aphid abundance data which we have seen has been gathered with weekly sampling intervals. Immigration is often visible in such data, though not as remarkably as with daily data. For example, Matis et al. [5] analyzes 13 individual observed data sets, which were gathered on the abundance of the mustard aphid over the period 1983–1990 on three different oilseed cultivars in North-
Fig. 1. Observed mean mustard aphid counts with daily sampling interval, and fitted curves. (A) Fitted model (6) without immigration (solid line) and model (11) with immigration (dashed line). (B) Fitted model (6) without immigration (solid line) and model (13) with immigration (dashed line).
ern India (at the Haryana Agricultural University). The data were gathered weekly, with the first aphid occurrence in late December/early January and a peak aphid count in late February. The symmetric function (6) is shown to provide a fit which is adequate for subsequent data analysis for all 13 cases. However, six of the 13 cases have a clearly visible initial simple linear increase, which leads to obvious lack-of-fit for function (6). The parameter estimates for regression function (6) and the root mean squared residuals, s, are given for each case in Table 1. These six cases are illustrated in Fig. 2. It appears that in each of these six cases there is a period of initial immigration, of four to six weeks duration, during which period the residuals are predominately positive, as the observed values tend to exceed the estimated values from regression function (6). 3.4. Summary of case for immigration Because local infestation is initiated by wind-borne (immigrant) aphids, we suggest that an immigration component is always part
Table 1 Estimates of regression parameters (N max , tmax , b) in model (6), mechanistic parameters (k; d; N 0 ), residual sum of squares and root mean squared error, ðRSS; sÞ, for 6 mustard aphid data sets Case
N max
tmax
b
k
d 103
N0
RSS
s
J8 C8 C9 C0 N7 N8
96.6 985.4 59.7 148.7 171.6 113.7
6.42 5.71 6.57 7.31 7.41 7.76
0.933 2.413 0.620 1.324 1.124 0.764
0.928 2.413 0.603 1.324 1.123 0.760
4.51 2.95 3.25 5.89 3.68 2.57
0.96 .004 3.83 .037 .166 1.20
960 142196 4382 1564 2409 2237
11.7 168.6 27.0 16.1 20.0 17.9
140
J.H. Matis et al. / Mathematical Biosciences 215 (2008) 137–143
Aphids
Aphids
C8 1400 1200 1000 800 600 400 200 6 8 Week
10
12
2
N8
200 175 150 125 100 75 50 25 2
4
6 8 Week
10
2
4
6 8 Week
10
4
10
12
6 8 Week
10
12
10
12
C0
200 175 150 125 100 75 50 25
12
6 8 Week C9
2
J8
200 175 150 125 100 75 50 25
4
200 175 150 125 100 75 50 25
12
Aphids
Aphids
4
Aphids
Aphids
2
N7
200 175 150 125 100 75 50 25
2
4
6 8 Week
Fig. 2. Curves for six mustard aphid data sets, with weekly sampling intervals, fitted to model (6) without immigration (solid line) and to model (11) with immigration (dashed line).
of the local population dynamics of the mustard aphid. However, the duration of the simple linear increase in population size attributable to immigration only, without reproduction, is dependent upon the ambient temperature. Examining all 13 mustard aphid data sets in Matis et al. [5], all four data sets for years 1988 and 1990 have simple linear phases of at least four weeks duration, indicating apparent low temperatures during this period of those years. On the other hand, none of the four data sets collected during 1983–86 have an apparent simple linear phase, perhaps due to higher temperatures during this period of those years. The data sets for 1987 and 1989 are mixed. In summary, we argue that all aphid abundance curves have an initial simple linear phase, though the duration is often so short that the symmetric function (6) may be sufficient for all practical purposes for describing weekly data with short infestation periods. At other times, however, particularly for daily data, the simple linear phase is striking. This motivates a model generalization for the mustard aphid that includes immigration. Empirical evidence for including immigration in the population models is also evident in our data for the pecan and for the cotton aphid in Texas, however the infestation period is shorter and the case is not as striking. 4. A cumulative-size dependent model with continuous immigration Though the case for immigration is arguably very conclusive, the precise ‘real-world’ description of the immigration effect is no doubt a complicated function of numerous weather and environmental variables. Our ultimate objective is, however, to fit our population models to field data, and hence we require a simple, robust mechanism, involving only a few parameters, for incorporating the immigration effect. Clearly, simulation modelers who
deal with more complex system models would not need to be as parsimonious with parameters. In this spirit, we propose two models which incorporate immigration. Both of these assume constant (time-homogeneous) immigration of size m insects/time interval. Our first model assumes that this effect lasts continuously throughout the period of observation. A second model, developed in the next section as an extreme contrast, terminates the immigration effect abruptly at some fixed time. In the first model, we assume that the cumulative-size dependent growth effect in (2) is added to the immigration effect after some time lag s. These assumptions yield model
N0 ¼
m t
ð11Þ
with m; s; k, and d P 0. This sudden start of the cumulative-size growth is also obviously an oversimplification. However, it seems appropriate as a first approximation for generalizing (2), and paraphrasing Kac [22], the model serves as a ‘caricature of reality’ which should portray some features of the ‘real-world’. This model was fitted to the observed mustard aphid data using numerical integration programs in EasyFit [23]. For simplicity, the initial value was set to 0, i.e. Nð0Þ ¼ 0, hence the model had four parameters compared to the three in model (2) without immigration. Fig. 1A illustrates the fitted curves for the daily aphid counts. The parameter estimates for the model are m ¼ 7:52 aphids/day, s ¼ 31:4 days, k ¼ 1:116/day and d ¼ 2:164 104 /aphid-day2. The sum of squared residuals, RSS, used as a measure of goodness-of-fit is RSS ¼ 51; 294. As a contrast, we also fitted model (2) without immigration. Its parameter estimates are k ¼ 0:2289/ day, d ¼ 6:764 105 =aphid-day2, and Nð0Þ ¼ 0:79 aphids, with
J.H. Matis et al. / Mathematical Biosciences 215 (2008) 137–143 Table 2 Estimates of mechanistic parameters (m; s; k; d) in model (11), estimated N max and t max , and residual sum of squares and root mean squared error, ðRSS; sÞ, for 6 mustard aphid data sets Case
m
s
k
d 103
N max
tmax
RSS
s
J8 C8 C9 C0 N7 N8
9.03 68.25 8.00 3.67 8.41 8.71
4.87 4.92 6.00 4.83 5.52 5.82
1.466 7.320 6.394 1.484 1.656 1.240
6.478 4.958 31.055 6.085 4.683 3.876
104.8 1426.5 112.8 149.1 182.5 124.9
6.5 5.6 6.7 7.3 7.5 7.8
268 29566 2514 1407 1153 780
6.7 86.0 22.4 16.8 15.2 11.4
RSS ¼ 211; 278. The addition of immigration obviously provides a far superior fit for these daily count data, as evident in Fig. 1A. The fitted curves of model (11) for the six data sets with weekly counts are illustrated in Fig. 2, and their parameters are given in Table 2. For convenience, we also calculate the residual standard error, s, as a measure of fit. Model (11) with immigration provides a noticeably better fitting curve than model (2) without immigration, for five of the six cases, with a reduction in s ranging from 16% to 48% for these cases. The largest reduction, 48%, occurred in case C8, which also has the largest s due to the high counts. The smallest reduction, 16%, occurred in case C9, for which the improvement due to immigration is obscured by an apparently outlying data point at t ¼ 9. Case C0 is instructive, as s did not decrease, but instead increased by 4% with the inclusion of immigration. It is clear that the fitted model (2) yields systematic lack-of-fit, as the residuals are positive for the first four time periods, t = 1–4. However, though model (11) has a lower sum of squared residuals, RSS, its residual standard error s, is slightly higher due to a loss of one degree of freedom from the residual degrees of freedom. As a rule, model (11) is expected to always provide a better fitting curve for mustard aphid data than model (2), even though s does not always decrease due to the limited number of data points.
As a contrast to model (11), we consider a model which recognizes that the immigration of aphids in India is seasonal and of limited duration. Instead of lasting for over the 10-week-period of the data collection, the changing climate conditions would tend to shift the wind-borne aphid immigration northward. We consider therefore a second model with immigration. In this model, rather than assuming continuous immigration, we make the other extreme assumption that immigration terminates at the end of the visible simple linear trend in the data, at the time therefore when reproduction starts. Hence the new model is
m ðk dFÞN
t
ð12Þ
The actual immigration effect obviously lies somewhere between the two assumed extremes in models (11) and (12). There is likely a gradual tapering off between time s and 1. However, through immigration plays a crucial role in initiating and sustaining the early part of the infestation, its effect even in model (11) is gradually swamped out after time s by the birth effect and the rapidly increasing death effect. Model (12) has a huge practical advantage due to its analytical solution. If one also assumes that the cumulative density, FðtÞ, applies only to the second, symmetric piece of the model, the model has piecewise solution
NðtÞ ¼
mt 4Nmax ebðttmax Þ ½1 þ ebðttmax Þ 2
t6s t > s:
The five parameters, m; s; N max ; t max , and b, may be estimated with standard nonlinear least squares software, e.g. we use SPSS [24]. Often s may be estimated by inspection from the data. In this simplifying case, the solution consists of two separate parts, each of which is easy to fit to data. The daily aphid data was fitted to regression model (13) with the fitted curve illustrated in Fig. 1B. The simple linear part fitted best for t < 32, with parameter estimate m ¼ 7:55 aphids/day. The cumulative-size model fitted best for t > 31, thus conveniently giving the bound 31 < s < 32 in (13). The cumulative-size model has parameter estimates N max ¼ 485:8 aphids, tmax ¼ 35:3 days and b ¼ 0:447/day. The corresponding mechanistic parameter estimates from (8) are k ¼ 0:447 and d ¼ 2:056 104 . The combined residual sum of squares for the two parts of model (13) is RSS ¼ 50; 527. This indicates a slightly better fit to model (13) than to model (11), which is not surprising due to its single additional parameter. Consider now the six data sets with weekly data. It appears by inspection, in each case, that the period with only immigration is at least four weeks but no longer than six weeks. Hence, for simplicity, the data points for t 6 5 were fitted to the simple linear equation in the first part of (13), and the data points for t P 5 were fitted to function (6) for the second part of (13). The parameter s is assumed to be the intersection of these two fitted curves. The fitted curves are illustrated in Fig. 3, and the parameter estimates are given in Table 3. The fitted curves of models (11) and (13) are compared in the next section. 6. Comparison of models with immigration It is obvious by inspection that both models (11) and (13) with immigration fit the daily data well, and both make a strong case for including immigration in the model. This section compares the two models for the six cases of weekly data. 6.1. Comparison of goodness-of-fit
5. A cumulative-size dependent model with restricted immigration
N0 ¼
141
ð13Þ
Both cases fit the data well qualitatively, as expected. Model (11) provides a slightly better fitting based on s for three of the six cases; whereas model (13) gives a somewhat better fitting for the other three. We are unable to make any generalizations on which model fits better based on these data. 6.2. Comparison of immigration components The estimated immigration rate m tended to be slightly higher for model (13) than for model (11). The estimated rate increased in five of the six cases, with the increase ranging from 1.0% to 13.8%, for an average increase of 6.0%. The estimated duration, s, was virtually identical in each case for the two models. Hence we conclude that the estimated immigration components of the two models are very close to one another. 6.3. Comparison of N max and tmax parameters The peak sizes, N max , of the two models are very close (within 2.5%) of one another except for cases C8 and C9, both of which are atypical. In case C8, N max ¼ 1; 426 for model (11) exceeds N max ¼ 1; 266 by 12.7%, however C8 has a smaller number of data points than the other cases which makes the determination of N max more challenging. In C9, N max ¼ 112:8 for model (11) exceeds N max ¼ 102:1 for model (13) by 10.5%, however this case has an anomalous observation at t ¼ 8 which makes any generalization questionable. In general, however, even with the cases with suspect and limited data, it appears that the estimated N max of model (13) is closer than those of model (11) to the ob-
142
J.H. Matis et al. / Mathematical Biosciences 215 (2008) 137–143
Aphids
Aphids
C8 1400 1200 1000 800 600 400 200 6 8 Week
10
12
2
N8
200 175 150 125 100 75 50 25 2
4
6 8 Week
10
2
4
6 8 Week
10
12
6 8 Week
10
12
10
12
10
12
C9
2
J8
200 175 150 125 100 75 50 25
4
200 175 150 125 100 75 50 25
12
Aphids
Aphids
4
Aphids
Aphids
2
N7
200 175 150 125 100 75 50 25
4
6 8 Week C0
200 175 150 125 100 75 50 25 2
4
6 8 Week
Fig. 3. Curves for six mustard aphid data sets, with weekly sampling intervals, fitted to model (6) without immigration (solid line) and to piecewise model (13) with immigration (dashed line).
Table 3 Estimates of regression parameters (m; N max , t max , b) in model (13), mechanistic parameters (k; d; N 0 ; s), residual sum of squares and root mean squared error, ðRSS; sÞ, for 6 mustard aphid data sets Case
m
N max
t max
b
k
d 103
N0
J8 C8 C9 C0 N7 N8
10.12 77.65 7.95 4.01 8.53 8.75
103.5 1265.5 102.1 150.3 178.2 123.9
6.57 5.65 6.79 7.32 7.45 7.85
1.129 3.414 2.174 1.360 1.234 0.960
1.128 3.414 2.174 1.360 1.234 0.956
6.16 4.61 23.15 6.15 4.27 3.72
2.48 2.12 1.59 2.85 7.25 2.64
E-1 E-5 E-4 E-2 E-2 E-2
s
RSS
s
4.9 4.9 5.9 4.8 5.3 5.7
162 31862 3214 1354 1241 619
5.2 89.2 25.4 16.5 15.8 10.2
served maximum in the data. The estimated tmax are very close to each other in each case.
bines with the higher birth rate in model (11) to yield a peak, N max , about the same size as in model (13) with lower rates.
6.4. Comparison of k and d rate parameters
6.5. Summary of comparison of models with immigration
Model (11) with continuous immigration always had a higher estimate of the birth rate, k, than model (13) with restricted immigration for the six cases. The increase ranged from 10% to 194%, i.e. almost tripling, for an average increase of 68% and a median of 32%. For the case with daily data in Fig. 1, the rate increased to k ¼ 1:116 from k ¼ 0:447, for a 150% increase, between the two models. This substantial increase seems paradoxical, as one might expect that the model in (13) would have a higher birthrate to compensate for the loss of immigration after s. However, model (11) also has a consistently higher estimate of the death rate d. The estimate of d was higher in model (11) than in model (13) for five of the six cases, the increase ranging from 4% to 34%. In the one case where the rate decreases, the decrease was only 1%. We suggest that conceptually the higher death rate com-
Our limited experience based on these mustard aphid data suggest that, as a rule: 1. Both models fit data about equally well. 2. Both models have similar shape characteristics, in particular regarding the size of immigration and the peak size, N max . In fact, the fitted curves are often visually indistinguishable, as they are for the daily data in Fig. 1A and B. 3. Model (11) with constant immigration tends to have higher birth and death rate estimates than model (13) with restricted immigration. If the accurate estimation of these rates is an important objective of the study, it seems crucial to seek expert subject-matter opinion regarding the likely duration of immigration, in order to decide between models (11) and (13).
J.H. Matis et al. / Mathematical Biosciences 215 (2008) 137–143
6.6. Comparison of models with immigration to model without immigration Though the models with immigration provide superior fittings in many cases, no doubt many experimenters will prefer to use the symmetric model in (6) for convenience. A question of interest then is how the estimated rates in (6) differ from those given by models (11) and (13). It is clear comparing the results of Tables 1–3, that in each case the estimated birth and death rates from the symmetric model (6) are lower than those estimated from either model (11) or (13). Hence, it seems crucial to use a model which includes immigration whenever immigration is clearly evident in the data, if one desires accurate estimates of these rates. 7. Discussion We suggest, based on purely theoretical grounds, that the initial increase in mustard aphid population size is due to immigration only. Often this initial period of growth due to immigration is sizeable, and in such cases a population model which includes a mechanism to incorporate immigration should be useful. Our limited data verify empirically that models with immigration often give superior fittings of data. We expect that these findings will be strengthened as more daily count data become available. This paper focuses on only two models with immigration. Both seem adequate for the data, however model (13) has a piecewise solution which is relatively simple to fit to data using widely available software. For the sake of completeness, we explored a number of other models incorporating immigration. One is the power-law model with immigration, which is model (11) with p ¼ 2. In general, this model did not fit any better, perhaps in part because the immigration already explained some of the skewness, which without immigration is attributed to the power-law. Another is a variation of model (11), in which the cumulative density FðtÞ is integrated from 0 to s over the immigration part rather than over the second, symmetric part. One cannot use the piecewise solution in (13) under this assumption, hence the numerical integration in EasyFit [23] was used again. This model also failed to improve the fittings, and hence was not investigated further due to its difficult practical implementation. One generalization which would be easy to implement would be to change the time-homogeneous immigration rate in (11) and (12) to some time varying function mðtÞ, which might be specified with few coefficients. In the absence of an alternative detailed hypothesis about mðtÞ, we did not explore such possibility. We note also that although the current investigation involves only data for the mustard aphid, similar initial population dynamics also apply no doubt to other aphid and non-aphid species. We leave this for others, probably with access to data with daily sampling intervals, to confirm as we do not have data necessary to verify this assertion.
143
Acknowledgement We express our appreciation to the associate editor and to the two reviewers for their useful suggestions which have improved this paper. References [1] R.L. Metcalf, R.A. Metcalf, Destructive and Useful Insects, fifth ed., McGrawHill, New York, 1995. [2] A.F.G. Dixon, Aphid Ecology, Chapman and Hall, New York, 1985. [3] J.H. Matis, T.R. Kiffe, T.I. Matis, D.E. Stevenson, Nonlinear stochastic modeling of aphid population growth, Math. Biosci. 198 (2005) 148. [4] J.H. Matis, T.R. Kiffe, T.I. Matis, D.E. Stevenson, Application of population growth models based on cumulative size to pecan aphids, J. Ag. Biol. Environ. Stat. 11 (2006) 1. [5] J.H. Matis, T.R. Kiffe, T.I. Matis, J.A. Jackman, H. Singh, Population size models based on cumulative size, with application to aphids, Ecol. Model. 205 (2007) 81. [6] T.I. Matis, M.N. Parajulee, J.H. Matis, R.B. Shrestra. A mechanistic model based analysis of cotton aphid population dynamics data, Ag. Forest. Ento. in press. [7] T.I. Matis, J.H. Matis, G.J. Michels Jr., M.N. Parajulee, Using the cumulative-size mechanistic model for analyzing insect data., Proc. 2007 Conf. Appl. Stat. Agr., Kansas State Dept. Stat., Manhattan, KS, 2008. [8] J.H. Matis, T.R. Kiffe, T.I. Matis, D.E. Stevenson, Stochastic modeling of aphid population growth with nonlinear, power-law dynamics, Math. Biosci. 208 (2007) 469. [9] P. Kindlemann, A model of aphid population with age structure, in: V. Capasso, E. Grosso, S.L. Paveri-Fontana (Eds.), Mathematics in Biology and Medicine, Proceedings, Bari, 1983, Lecture Notes in Biomathematics, Springer, Berlin, 1985, p. 72. [10] P. Kindlemann, R. Arditi, A.F.G. Dixon, A simple aphid population model, in: J.C. Simon, C.A. Dedryver, C. Rispe, M. Hulle (Eds.), Aphids in a New Millennium, INRA, Paris, 2004, p. 325. [11] Prajneshu, A nonlinear statistical model for aphid population growth, J. Indian Soc. Agric. Stat. 51 (1998) 73. [12] J.H. Matis, T.R. Kiffe, T.I. Matis, J.A. Jackman, W.E. Grant, H. Singh, On the use of growth rate parameters for projecting population sizes: application to aphids, Ecol. Model. 213 (2008) 133. [13] N.L. Johnson, S. Kotz, Continuous Univariate Distributions – 1, Wiley, New York, 1970. [14] J.R. Riley, D.R. Reynolds, S. Mukhopadhyay, M.R. Ghosh, T.K. Sarkar, Longdistance migration of aphids and other small insects in northeast India, Eur. J. Ento. 92 (1995) 639. [15] A. Srivastava, N.P. Kashyap, Population build-up of mustard aphid, Lipaphis erysimi (Kalt.) on brassica in relation to weather factors, J. Insect Sci. 12 (1999) 141. [16] S.K. Prasad, K.G. Phadke, K.N. Mehrotra, Forecasting of mustard aphid, Lipaphis erysimi Kaltenbach on rapeseed crop, Proc. Indian Natl. Sci. Acad. (B) 50 (1984) 146. [17] K.G. Phadke, Ecological factors affecting aphid Lipaphis erysimi incidence on mustard crop, in: K.G. Phadke (Ed.), Aphidology in India, A.R. Printers, Calcutta, 1986, p. 37. [18] S.K. Roy, K. Baral, Role of weather parameters on population build up of mustard aphid Lipaphis erysimi (Kaltenbach), J. Oilseeds Res. 19 (2002) 86. [19] A. Campbell, M. Mackauer, Reproduction and population growth of the pea aphid (Homoptera: Aphididae) under laboratory conditions, Can. Ento. 109 (1977) 277. [20] R.I.K. McVean, A.F.G. Dixon, T. Harrington, Causes of regional and yearly variation in pea aphid numbers in eastern England, J. Appl. Ento. 123 (1999) 495. [21] T.X. Liu, B. Yue, Comparison of some life history parameters between alate and apterous forma of turnip aphid (Homoptera: Aphididae) on cabbage under constant temperatures, Florida Ento. 84 (2001) 239. [22] M. Kac, Some mathematical models in science, Science 166 (1969) 695. [23] EasyFit, MathWave Technologies, (www.mathwave.com), Arlington, MA, 2005. [24] SPSS. SPSS 11.0 for Macintosh. SPSS Inc., Chicago, IL, 2002.