Population model for the decline of Homalodisca vitripennis (Hemiptera: Cicadellidae) over a ten-year period

ABSTRACT The glassy-winged sharpshooter, Homalodisca vitripennis (Germar), is an invasive pest which presents a major economic threat to grape industries in California, because it spreads a disease-causing bacterium, Xylella fastidiosa. In this note we develop a time and temperature dependent mathematical model to analyze aggregate population data for H. vitripennis from a 10-year study consisting of biweekly monitoring of H. vitripennis populations on unsprayed citrus, during which H. vitripennis decreased significantly. This model was fitted to the aggregate H. vitripennis time series data using iterative reweighted weighted least squares (IRWLS) with assumed probability distributions for certain parameter values. Results indicate that the H. vitripennis model fits the phenological and temperature data reasonably well, but the observed population decrease may possibly be attributed to factors other than the abiotic effect of temperature. A key factor responsible for this decline but not analyzed here could be biotic, for example, potentially parasitism of H. vitripennis eggs by Cosmocomoidea ashmeadi. A biological control program targeting H. vitripennis utilizing the mymarid egg parasitoid Cosmocomoidea (formerly Gonatocerus) ashmeadi (Girault) is described.

H. vitripennis has exhibited high invasion potential and has become a significant pest outside of its native range. Transportation of live plants infested with H. vitripennis resulted in this pest establishing in Hawaii, the Cook Islands, Easter Island, and in French Polynesia, and its population densities reached extraordinarily high levels [15,16]. Around 1990, GWSS invaded California and represented a major economic threat to the wine, table, and raisin grape industries, because it was implicated in a significant increase of the lethal and incurable grape malady, Pierce's disease [12,25]. This disease cost approximately $104 million per year in crop damage and resources devoted to mitigating this threat [26]. The H. vitripennis population in CA has exhibited a significant decline. In this note we develop a time and temperature dependent mathematical model to analyze aggregate population data for H. vitripennis from a 10-year study (during which H. vitripennis decreased significantly) consisting of biweekly monitoring of H. vitripennis populations on unsprayed citrus. There are a number of abiotic and biotic factors that can be implicated in regulating GWSS population numbers, including temperature [1,25], rainfall [17], host plant [14,22], parasitism [24], etc. -with, for example, Cosmocomoidea ashmeadi known to be a primary parasitoid [27].
Female H. vitripennis oviposit eggs in 'mass' consisting of approximately 10 individual eggs which are deposited side-by-side under the epidermis on the underside of leaves. C. ashmeadi is a solitary endoparasitoid that has high levels of host specificity. Female parasitoids lay individual eggs inside a H. vitripennis egg. The developing parasitoid larva kills the host egg by consuming the contents. It then pupates within the host egg and emerges as a winged adult capable of flight. C. ashmeadi was accidentally introduced into California probably as parasitized H. vitripennis eggs that were infesting leaves of imported plant material that originated from the native range of this pest [24]. Since its discovery in California, significant effort has been invested in mass rearing and releasing this natural enemy for the biological control of H. vitripennis [19].
In the present initial investigation we develop a mathematical model to investigate whether abiotic density independent factors, such as temperature, could be the primary cause of the observed decrease in H. vitripennis populations over time. Commencing 13 March 2002 and ending 8 August 2012, biweekly surveys were made of ten Eureka lemon trees in an unsprayed mixed variety citrus plot in the Agricultural Operations Facility at the University of California Riverside (UCR Ag. Ops). Field surveys were conducted at approximately 1200 to 1300 h for each sample date. Experimental trees were divided into north, south, west, and east quadrants, and visually searched for 30 s each. The total number of adult H. vitripennis infesting each citrus tree was recorded. Following visual counts for adults, each tree was searched for five minutes and all H. vitripennis egg masses that were found were removed from trees, labelled by tree, and were taken to the lab where the total number of egg masses and individual eggs collected on each sampling date were recorded.

Aggregate population data sets
Labeled egg masses were placed in ventilated vials and maintained at 25 • C in a temperature cabinet for the emergence of parasitoids. Parasitoid emergence was checked weekly for up to three weeks post collection (this allowed adequate time for parasitoids to complete development and to emerge). The number of parasitoids that emerged was recorded, and sex and species determinations were made. The average realized percent parasitism (total parasitoid emergence divided by the sum of total parasitoid emergence and total GWSS nymph emergence) over all samples was 36%. The samples taken in 2002 had the largest average percent parasitism at 82%, and 2009 had the smallest average percent parasitism at 10%. The vast majority of parasitoids that emerged were C. ashmeadi, 4714 individuals or 84% of all emerged parasitoids. The second most common parasitoid was Ufens sp. (Hymenoptera: Trichogrammatidae) (500 individuals or 9% of emerged parasitoids), followed by Cosmocomoidea (formerly Gonatocerus) morrilli Howard (384 individuals or 7% of emerged parasitoids). These population data are plotted in Figure 1. Until other factors can be ruled out, it is unclear whether lab reared parasitoids can be associated with a primary cause of the 10-year decline of H. vitripennis egg and adult populations. The first four years of this data (March 2006 to August 2006) were published in [18], but the full 10year data set has not yet been published, and none of these data have been mathematically modelled.
It is highly unlikely that the biweekly removal of egg masses and the parasitoid-infected eggs significantly affected either the GWSS population or the parasitoid population. The number of egg masses removed each sampling period was relatively low in comparison to the availablility on the trees, especially with the five minute time limit. In addition, only a small part of each experimental tree was sampled, leaving the top of and bottom ares of each tree undisturbed. The experimental Eureka lemon trees were in a large plot of over 100 unsprayed citrus trees that all had high GWSS infestations. H. vitripennis adults travel long distances to find optimal host plants [22] and readily fly from tree to tree in this experimental plot, so any removed egg masses would have been readily replaced by GWSS from neighboring trees.

Temperature data sets
The daily minimum and maximum temperatures collected from 1 January 2002 to 31 December 2012 in Riverside, CA were downloaded from an automated weather station, CIMIS 44, located at UCR Ag. Ops and used to calculate monthly averages (Figure 1(D)). For the mathematical simulations, we assume a daily sinusoidal curve in temperature fluctuation in which the maximum temperature occurred at approximately 1500 h, and the minimum temperature occured at approximately 0300 h daily (personal observation from data). Average monthly temperatures during the study period are presented with H. vitripennis egg and adult population data collected over that same period in Figure 2.
Homalodisca vitripennis are strong fliers and travel long distances to find optimal host plants [22]. Thus, it is reasonable to assume that each data point represents a different subsample of the total population, making this an aggregate [7] data set (i.e. each sample may possibly consist of differing individuals from the overall population). Consequently, H. vitripennis adults, H. vitripennis eggs, and the parasitoids were assumed to have been sampled longitudinally from the aggregate population. Therefore, parameters regulating population fluctuation are likely described by some probability distribution (see Chapter 5 in [7]). This is vital and changes the inverse problem when comparing aggregate population data to a mathematical model meant to describe an individual population.

Mathematical model
H. vitripennis is known to begin life as an egg and undergo five nymphal instars in which they are incapable of sexual reproduction before developing into sexually mature adults. Thus, we use the continuous-time age-structured population growth model described by Kapur [20] to model the H. vitripennis individual population. For simplicity, we model the egg stage, adult stage, and combined nymphal stage in  [25]. Developmental rates, r e (T) and r N (T), are defined by the Lactin model with T being temperature in • C and parameter values listed in Table 1. Pilkington et al. [25] uses life tables to find the net female maternity, R 0 (T), and mean generation time, T c (T), at different temperatures, and fits each of these parameters to a quadratic curve. We let is the average number of offspring produced per female per day at given temperature T.
Finally, temperature-dependent GWSS adult death rate is calculated by inverting the mean adult longevities (emergence as an adult to death) from Pilkington et al. [25] and fitting them to a quadratic curve, Since birth, developmental, and death rates are nonnegative, negative values of b(T), r e (T), r N (T), and d A (T) are set to zero. The sex ratio of H. vitripennis is assumed approximately even and unaffected by temperature [25].
H. vitripennis have a lack of effective natural enemies, abundant host plants, and no significant competitors in California [16,17,23,25], so we assume that most GWSS deaths can be attributed to unfavourable temperatures. In a controlled laboratory study, H. vitripennis spend on average 5 to 14 days as eggs and 34 to 65 days as nymphs [25]. Thus, we assume that if a GWSS egg has not developed in 20 days then it will die, so we set the GWSS egg death rate to be d e = 1/20 = 0.05. Similarly, we assume that if a GWSS nymph has not developed in 50 days then it will die, so we set the GWSS nymph death rate be d N = 1/50 = 0.02.
In Figure 2, egg and adult GWSS population data tend to increase and decrease with the temperature, and maximum yearly population counts occur during the hottest parts of the year. If we calculate the temperature-dependent birth and delelopmental rates, r e (T t ), r N (T t ), and b(T t ), at the current daily temperature T t , then the model populations x e , x N , and x A begin to increase during the hottest time of the year, and the population counts do not reach their yearly maximum until a few months later. This behaviour, which is presented alongside the data in Figure 3, is not constistent with the data. The average time between the maximum yearly model population and the next maximum yearly data population, as plotted in Figure 3, is 269 days, which we round to 270 days. Thus, temperature-dependent rates r e (T t−270 ), r N (T t−270 ), b(T t−270 ), and $d_A(T_{t-27-})$ are calculated using temperatures from 9 months or 270 days in the past. With this added delay, the maximum yearly model populations, x e , x N , and x A , occur during the hottest   Table 1. time of the year. This matches up with the behaviour of the phenology population data presented in Figure 2. We plot the observable model populations x e and x A with the 270 day temperature delay and the corresponding H. vitripennis egg and adult population data in Figure 4.
We have temperature data beginning on 1 January 2002, so with a 270 day delay we cannot model populations earlier than 28 September 2002. Thus, the first fifteen dates in the time series cannot be modelled, bringing the total number of usable H. vitripennis egg and adult data points from 544 to 514, which is still sufficient for our purposes (usually 30-40 data points [11] are suggested per parameter being estimated).

Aggregate model and parameter estimation methodology
Equations (1)-(3) model individual populations, where all individual glassy-winged sharpshooters are assumed to have the same growth, birth, and death rates, i.e. the same parameter values from Table 1 that determine these rates. However, this assumption does not apply to the aggregate data which consists of frequent sampling of the changing (over time) populations with possibly different individuals at each time point having different parameter values (e.g. ρ e or λ e or a 0 , etc., varying over a range of values θ ∈ G). Thus, using the individual model we must formulate an aggregate model with which we can compare the aggregate data (see Chapter 5 of [7] for a more complete discussion).
Consider the 3-dimensional mathematical model or dynamical system, (1), (2), and (3), and θ is a parameter chosen from Table 1 (e.g. θ = ρ e or λ e or a 0 , etc.). There is an aggregate population vector u = [u e , u N , u A ] T corresponding to the individual population vector x and given by where P = P(θ ) is now a random variable, x(t; θ) is the solution to (8), G is the collection of admissible parameter values for θ , and P is a probability measure on G. Note that u is the expected value of x, which is also a random vector. Under the assumption that the probability distribution, P, possesses a density, p, the population count is given by where the density P = dP/dθ = p(θ ). Now that we have an aggregate model which can be compared to the data, we follow techniques from Banks, Hu, and Thompson [7] and Banks, Bekele-Maxwell, Everett, Stephenson, Shao, and Morgenstern [10] to estimate parameters in our mathematical model. For the inverse problem, consider the 3-dimensional dynamical system, defined in (9) to be estimated using the data. We are interested in determining the probability density P = p(θ ) which gives the best fit of the underlying model to the aggregate data. However, this parameter estimation problem involves an infinite dimensional parameter space (the space P(G) of probability measures defined on the set G). Instead of using a specific probability density function in the aggregate model, we use a family of finite approximations P M (G). Based on [3,6,11], we are guaranteed convergence in the Prohorov metric. Thus, the finite dimensional approximation P M (G) to the probability measure space P(G) is defined using linear splines given by where the chosen and known piecewise linear splines are represented by l M k (θ ), and the M-dependent (unknown) spline coefficients are M k for k = 0, . . . , M. Thus, in order to estimate P in equation (9), we only need to estimate the set of M+1 spline coefficients, where P M is the finite approximation of the probability distribution P and { M k } M k=0 are the M+1 spline coefficients of the probability density function which uniquely defines P M . If we assume G = [θ l , θ u ] is a closed bounded interval of possible values for the parameters, then is the new approximate 3-dimensional dynamical system to consider in the inverse problem.
Define the m-dimensional aggregate observation process where is an m × 3 matrix with m = 2, since our data set consists of only egg and adult aggregate populations. Let u j = (u 1 j , u 2 j ) represent the H. vitripennis egg and adult population data or observations collected at time t j for j = 1, . . . , n. We note that there is some uncertainty between the actual phenomenon, which is represented through the data, and in the above observation process. This uncertainty is accounted for in the statistical model (again, see Chapter 3 of [7] where it is explained that both a mathematical model and a statistical model are required to carry out an inverse problem properly with uncertainty in observations) where P M 0 is the nominal probability density approximation. The values γ = (γ 1 , γ 2 ), γ i ≥ 0 in the error term are weighting factors representing the dependency of error on the dynamics, and the model itself corresponds to the choices of data (see Section 4.3 below for further comments on the choice of γ ).
Note that°is the Hammond or Schur product of component wise multiplication of two vectors. The n × 1 random error vectors E 1 and E 2 respectively are assumed to be independent and identically distributed (i.i.d.) with mean zero, Var(E 1 j ) = σ 2 01 , and Var(E 2 j ) = σ 2 02 . The corresponding realizations (for the random variables [U 1 This multiplicative structure of the observational error in the above statistical model exists, because often in biological models the size of the observation error is proportional to the size of the observations. A rather thorough discussion of these issues along with concrete examples is given in Chapter 3 of [7]. For γ ≥ 0, a generalized least squares method or an iterative reweighted weighted least squares (IRWLS) method [10] as used below is appropriate to perform the inverse problem. In order to estimateP M ≈ P M 0 , we want to minimize the distance between the collected data and aggregate mathematical model, where the observables are weighted according to their variability and, for each observable, the observations over time are weighted unequally (once again, we refer the reader to [7,10] for a more detailed discussion and relevant examples).
The iterative reweighted weighted least squares estimate,P M , is numerically determined by iteratively solving the following system: for n data points. The matrixŴ j = diag{f 2γ 1 1 (t j ;P M ), f 2γ 2 2 (t j ;P M )} is made up of squared error weights, andV j is the estimated covariance matrix at data point j = 1, . . . , n. We use the following iterative procedure [5,7]: (1) EstimateP M(0) using (10) withV j = I. Set l = 0.
Note that this is not the same as taking the derivative of the argument in the right side of (14) and setting it equal to zero. We iteritavely want to minimize The estimated variances for each observable σ 2 01 and σ 2 02 are approximated by the following: If we assume γ = (γ 1 , γ 2 ) = (0, 0), then our statistical model is called an absolute error model and an ordinary least squares method is appropriate for parameter estimation. However, we believe it is more biologically realistic to assume the observation error is proportional to the size of the observed quantity i.e. a relative error model for our data sets and models investigated here.

Uncertainty quantification: standard errors
In order to quantify uncertainty in estimating the linear spline coefficients M k for the aggregate model in equation (11), the standard errors and confidence intervals can be computed using standard asymptotic theory for the IRWLS spline coordinate estimators M k [4,5,7]. For each time point t j , j = 1, . . . , n, the m × (M + 1) sensitivity matrix for the i = 1, . . . , m observations is where x 1 = x e and x 2 = x A , and f 1 = u e and f 2 = u A is defined in equations (12). The is approximated using the sensitivity matrices D j (P M ), the covariance matricesV j from equation (15), and the matrices of error weightsŴ j = diag{f Thus, if the condition number of the Fisher Information Matrix, κ(F ), is large then the probability density coefficient estimates are more uncertain. Before we perform the inverse problem, we must determine a correct statistical model by selecting γ . One method for doing this is described in the next section.

Difference-based methods and modified residuals
We use a second order difference-based method [9] to determine the correct statistical model ( γ value). Another method often implemented consists of performing an inverse problem with some γ value and computing the modified residuals, for each observable i at time t j , j = 1, . . . , n. The plots of M i j vs. t j should be randomly scattered around the horizontal axis. If an undesired non-random shape is present (e.g. a fan or inverted fan shape), then a different γ i value is chosen and the process is repeated until a γ i value produces the desired random scatter plot. However, this method does not consider both the mathematical model and statistical model misspecifications; it determines the correct statistical model under the tacit assumption that one has a correct mathematical model. It is also time consuming as it might take several attempts of performing an inverse problem (each of which may be computationally expensive) and plotting the modified residuals until a good statistical model is chosen.
We follow [9] and first apply the second order difference-based method directly to the data to determine the correct γ i value, which is both computationally economical as well as time efficient and independent of any assumed correct mathematical model. We first calculate the following pseudo measurement errorŝ for each observable i at time t j , j = 1, . . . , n. Next, we calculate the modified pseudo errors for each observable i at time t j , j = 1, . . . , n for different values of γ i ranging between 0 and 2 and plot these modified pseudo errors versus time. By inspecting these plots, we find that γ = (γ 1 , γ 2 ) = (0.3, 0.6) produce the most visually randomized scatter plots (again, see [7] for methodical discussions). The unweighted modified pseudo errors and the modified pseudo errors with the chosen weights γ = (γ 1 , γ 2 ) = (0.3, 0.6) are plotted in Figure 5. We use these values in the statistical model to perform the inverse problem and compute modified residuals in (21). If the modified residual plots are not randomly distributed around the horizontal axis, then the error must be due to mathematical model misspecification, implying another iteration of the modelling process is needed. Before we peform the inverse problem, we must select a parameter from Table 1 over which we assume the populations are distributed.

Local sensitivity analysis for parameter selection
Usually we choose to estimate parameters to which the model observations are most sensitive. However, the aggregate observations, f (t; P M ) = [u e , u A ] T from (12) and (13), depend explicitly on the parameter chosen to estimate. Thus, we perform sensitivity analysis on the individual population model, x e and x A defined in equations (1) and (3). The local sensitivity of f to a parameter θ =θ is where t is time, and f = x e or f = x A . Many of the parameters have different orders of magnitude, so in order to better compare the sensitivities of these parameters, it is useful to observe the normalized sensitivities, The sensitivities depend on our initial parameter value choice,θ . The individual popuations x e and x A die out almost immediately with parameter values from Table 1, making  Table 1. The individual populations x e and x A are temperature-dependent, so we use real temperature data downloaded from an automated weather station, CIMIS 44, located at UCR Ag. Ops in Riverside, CA. We use a fourth order Runge Kutta method and the complex step method [8] to numerically evaluate the sensitivities at the times and temperatures from the data set. The sensitivites of x e and x A to parameters T ue , T uN , ρ e , ρ N , λ e , λ N , a 0 , a 1 , and a 2 are positive, because as these parameters increase either the birth rate or one of the developmental rates increases, so the GWSS populations also increase. As these rates increase, the egg and adult populations should also increase. Similarly, the sensitivities of x e and x A to parameters e , N , b 0 , b 1 , and b 2 are negative, because as these parameters increase either the birth rate or one of the developmental rates decreases, so the GWSS populations also decrease. The sensitivities of x e and x A to the death rate parameters, d e , d N , μ 0 , μ 1 , and μ 2 are negative, because if these parameters increase, the GWSS populations tend to decrease. The normalized sensitivites are either positive and linearly increasing with time, or negative and linearly decreasing with time on a log scale. We plot the sensitivities and normalized sensitivities of x e to λ e and b 2 in Figure 6 to illustrate this behaviour. For each parameter, we find the sensitivity and normalized sensitivity of x e and x A with the largest magnitude over time, take the absolute value of each of these values, and plot them in Figure 7. According to these plots, x e and x A are most sensitive to the parameter ρ N , ρ e , and λ N . The observations have the largest normalized sensitivity to the parameters a 1 , λ N , and a 2 , followed by the parameters a 0 , b 1 , b 2 , b 0 , and λ e . Thus, we choose the egg developmental rate parameters ρ e and λ e , the nymph developmental rate parameters ρ N and λ N , and the birth rate parameters a 1 and a 2 to estimate in the inverse problem. Each of the birth rate parameters (a 0 , a 1 , a 2 , b 0 , b 1 , b 2 ) have approximately the same effect on the birth rate equation within an appropriately chosen range of values, so the results when estimating a 0 , b 0 , b 1 , and b 2 are redundant and not shown.

Estimating the probability density of egg developmental rate parameters
First, we assume a probability distribution over parameter ρ e , which is part of the egg developmental rate equation. At different closed intervals the parameter ρ e still converges and exhibits the same behaviour, so for these results we choose to let ρ e vary over the interval [0.001, 0.060] which permits convergence to a centred value of the distribution for ρ e . Next, we hold ρ e fixed and assume a probability distribution over parameter λ e , which is    (16), the estimated variances of the eggs and adults from equation (17), and the condition number of the Fisher Information Matrix from equation (19) are listed in Table 2.
The estimated probability distributions of ρ e are clustered around ρ e ≈ 0.045, and the estimated probability distributions of λ e are centred around λ e ≈ 0.4. These estimates are much larger than the parameter values originally chosen and listed in Table 1 of ρ e = 0.0073 and λ e = −1.087. In both cases, the probability distributions are multimodal, and the estimated variances of the GWSS eggs listed in Table 2 are much larger and different from the results of the other inverse problems (see Tables 3 and 4). Thus, the results of these inverse problems are more unstable. Furthermore, it usually takes a GWSS egg at least 5 days to develop [25], but at the estimated values, the egg developmental time is as short as 12 h at optimal temperatures. Table 2. Optimal cost value (J * ), estimated variances for eggs and adults respectively (σ 2 01 andσ 2 01 ), and the condition number of the Fisher Information Matrix, F when computing the inverse problem for ρ e and when computing the inverse problem for λ e .  Table 3. Optimal cost value (J * ), estimated variances for eggs and adults respectively (σ 2 01 andσ 2 01 ), and the condition number of the Fisher Information Matrix, F when computing the inverse problem for ρ N and when computing the inverse problem for λ N .  Table 4. Optimal cost value (J * ), estimated variances for eggs and adults respectively (σ 2 01 andσ 2 01 ), and the condition number of the Fisher Information Matrix, F when computing the inverse problem for a 1 and when computing the inverse problem for a 1 and a 2 , respectively.

Estimating the probability density of nymph developmental rate parameters
We assume a probability distribution over parameter ρ N , which is part of the nymph developmental rate equation. At different closed intervals the parameter ρ N still converges and exibhits the same behaviour, so we arbitrarily let ρ N vary over the interval [0.002, 0.005]. Next, we hold ρ N fixed and assume a probability distribution over parameter λ N , which is also part of the nymph developmental rate equation. At different closed intervals the parameter λ N still converges and exhibits the same behaviour, so we arbitrarily let λ N vary over the interval [−1.04, −0.96]. The temperature-dependent nymph developmental rate is plotted in Figure 11 for different values of ρ N and λ N , respectively. We run the two separate inverse problems with M = 20, 24, 28, and 32 spline functions to estimate the probability distributions of ρ N and λ N . The resulting probability densities   M = 32. The optimal cost value J * from equation (16), the estimated variances of the eggs and adults from equation (17), and the condition number of the Fisher Information Matrix from equation (19) are listed in Table 3.
The estimated probability distributions of ρ N are clustered around ρ N ≈ 0.0035, and the estimated probability distributions of λ N are centred around λ N ≈ −1. The probability distributions in both cases appear to be approaching a single value for both parameters. These estimates are much larger than the parameter values originally chosen and listed in Table 1 of ρ N = 0.0022 and λ N = −1.03. GWSS usually spend at least 50 days as nymphs [25], but at the estimated values, GWSS only spend 20 days as nymphs at optimal temperatures.

Estimating the probability density of birth rate parameters
Next, we assume that the glassy-winged sharpshooter eggs and adults have a probability distribution over the parameters a 1 and a 2 , which are part of the birth rate equation. We let a 1 vary over the interval [57.4, 58.2] while fixing a 2 , and then we let a 2 vary over the interval [−1.14, −1.08] while fixing a 1 . At different closed intervals the parameters a 1 and a 2 still converge and exhibit the same behaviour, so we choose these intervals arbitrarily. The temperature-dependent birth rate is plotted in Figure 14 for different values of a 1 and a 2 , respectively. Since these two parameters will be estimated separately and independently, Figure 14. Temperature-dependent birth rate of H. vitripennis from equation (6) for different values of a 1 (A) and a 2 (B). Note that a 1 and a 2 are parameters in the birth rate equation. Parameter values not stated in the plot legends are taken from Table 1. whenever a 1 is varied, a 2 is held fixed at its value from Table 1, and whenever a 2 is varied, a 1 is held fixed at its value from Table 1.
We run the two separate inverse problems with M = 20, 24, 28, and 32 spline functions to estimate the probability distributions of a 1 and a 2 . The resulting probability densities   (16), the estimated variances of the eggs and adults from equation (17), and the condition number of the Fisher Information Matrix from equation (19) are listed in Table 4. The estimated probability distributions of a 1 are clustered around a 1 ≈ 58.05, and the estimated probability distributions of a 2 are centred around a 2 ≈ −1.09. The probability distributions in both cases appear to be approaching a single value for both parameters. Based on data from [25], GWSS adults lay, on average, one egg per five days, but at the estimated values GWSS lay as many as 5 eggs per day at optimal temperatures.

Discussion
The goal of this study is to understand what led to the H. vitripennis population decrease during a 10-year study. Using only the biotic effect of temperature, our results are stable and appear to follow the trends in the data reasonably well. Although assuming a probability distribution over the estimated parameters was necessary due to the aggregate nature of the data, most of the parameters appear to approach a single value as M, the number of linear spline functions used to estimate the probability density, approaches infinity. While we do not form a 'best value' of M to use in our approximations of the probability measures, we found consistently that M = 32 leads to decreasing values of the weighted least squares measure J * while at the same time leads to only slight increases in the ill-conditioning of the Fisher Information Matrix κ(F ), suggesting a 'sweet spot' for the reported values of the parameters.
We obtain inferior results when we estimate the egg developmental rate parameters ρ e and λ e , chiefly when fitting the GWSS egg population model to the data. The aggregate egg populations are estimated to be zero (see Figures 9(B-C)), and the egg observation variance is estimated to be σ 2 1 ≈ 2 × 10 7 which is much larger compared to the other inverse problems (σ 2 1 ≈ 8 × 10 6 ). Furthermore, while other estimated probability distributions are centred around a single value, the probability distributions of the egg developmental rate parameters are spread over a large range of values (0.03 < ρ e < 0.06 and −0.1 < λ e < 1.0) without any signs of convergence as M increases. Thus, the results when estimate the egg developmental rate appear to be more uncertain than when estimating the nymphal developmental rate and the birth rate.
We obtain much better results when estimating birth rate parameters, and we obtain the best results when estimating the nymphal developmental rates. The observation variances (σ 2 1 ≈ 7.8 × 10 6 for eggs and σ 2 1 ≈ 2.5 × 10 6 for adults) when estimating the nymphal developmental rate parameters are smaller than when estimating birth rate parameters (σ 2 1 ≈ 8 × 10 6 for eggs and σ 2 1 ≈ 3 × 10 6 for adults), because more variance is captured in the model. Furthermore, the probability distributions of the nymphal developmental rate parameters λ N and ρ N exhibit consistent unimodal behaviour as M increases (see In all cases, the estimated parameters are larger than the parameter values originally estimated in controlled lab experiments [25] (see Table 1), indicating that this population of H. vitripennis have exceptionally high developmental and birth rates. In the non-aggregate model described in (1)-(3), using all of the chosen parameter values from Table 1 and real temperature data, the GWSS popoulations die out almost immediately. Thus, in order to fit the nonzero data, it follows that the inverse problem estimates this population to have larger parameter values and more robust birth and developmental rates. In other words, this population of H. vitripennis is laying more eggs, developing more quickly from eggs to nymphs and nymphs to adults, or some combination of these than we anticipated. This does not necessarily rule out parasitism as the cause of the GWSS decline, but there is some disconnect between GWSS populations reared in the lab at constant temperatures and GWSS populations in the wild. In a study conducted by Bahar, Soroka, and Dosdall [2], populations of the insect Plutella xylostella had shorted developmental times and higher adult longevity compared with populations reared under constant temperatures. Thus, there may be similar difference between lab-reared GWSS populations from which initial parameter values were gathered, and GWSS in the wild.

Conclusion
This study takes a first look at aggregate population phenology data for H. vitripennis collected over the course of 10 years at the Agricultural Operations Facility at the University of California Riverside. The decline in H. vitripennis population observed in this data could potentially be attributed to many factors. The abiotic effect of temperature is investigated here and does a reasonable job describing the data. The phenology data fit the model best when the birth and developmental rates are much larger than originally estimated in lab studies.
For our future work, we intend to add the effect of parasitism to the mathematical model, since this population of H. vitripennis are known to have been attacked by C. ashmeadi. We intend to analyze the parasitoid and GWSS data together, and test for statistically significantly better fits to our data.

Disclosure statement
No potential conflict of interest was reported by the authors.