Disease outbreaks in plant-vector-virus models with vector aggregation and dispersal

While feeding on host plants, viruliferous insects serve as vectors for viruses. Successful viral transmission depends on vector behavior. Two behaviors that impact viral transmission are vector aggregation and dispersal. Vector aggregation may be due to chemical or visual cues or feeding preferences. Vector dispersal can result in widespread disease outbreaks among susceptible host plants. These two behaviors are investigated in plant-vector-virus models. Deterministic and stochastic models are formulated to account for stages of infection, vector aggregation and local dispersal between adjacent crops. First, models for a single crop are studied with aggregation included implicitly through the acquisition and inoculation rates. Second, models with aggregation and dispersal of vectors are studied when one field contains a disease-sensitive crop and another a disease-resistant crop. Analytical expressions are computed for the basic reproduction number in the deterministic models and for the probability of disease extinction in the stochastic models. These two expressions provide useful measures to assess effects of aggregation and dispersal on the rate of disease spread within and between crops and the potential for an outbreak. The modeling framework is based on cassava mosaic virus that causes significant damage in cassava crops in Africa.

1. Introduction. Plant viruses are often insect-transmitted. Begomoviruses, carried by the whitefly Bemisia tabaci, are some of the most important emerging viral pathogens of plants [39]. These viruses cause drastic reductions in crop yield, especially for cotton, soybean, tomato, and cassava crops [39]. Cassava is a staple food crop in many countries, including those in Africa, Asia, and South America. One of the most devastating diseases of cassava is cassava mosaic virus (CMV), transmitted by the whitefly [1,11] which causes severe reduction in tuber yield. Cassava mosaic viruses are found in Africa, India, and neighboring islands. Distinct virus species are known as African cassava mosaic virus, East African cassava mosaic virus, South African cassava mosaic virus, Indian cassava mosaic virus, and Sri Lankan cassava mosaic virus [1]. Root yield loss to cassava mosaic disease in sub-Saharan Africa is estimated to be 15% to 24% annually [1].
Cassava plants are vegetatively propagated. That is, new plants are propagated from stem cuttings. Therefore, infected cassava plants can be introduced if cuttings come from infected plants or if stems from infected cassavas are grafted onto healthy cassavas [1,11]. However, the major method of transmission is through the whitefly, Bemisia tabaci. Whiteflies feed on phloem sap by biting into the leaves of host plants. CMV is a persistent, nonpropagative virus in whiteflies, meaning that the virus does not replicate within the insect but the virus generally persists in the whitefly for its' entire lifetime [37].
Management of CMV in crops requires control of whiteflies or the diseased plants [45]. Control of whiteflies is difficult, especially in large crops. Pesticides are rarely used because they are ineffective against whiteflies and they may negatively impact any natural whitefly predators and the environment [1]. Other control methods include culling of infected plants and planting of cassava plants resistant to CMV [1,11]. Planting of resistant varieties prevents or slows the spread of the disease [11,44,45].
Two vector population behaviors that impact methods for control of disease spread are vector aggregation and dispersal. Vector aggregation has been observed in insect-transmitted plant viruses, which may be due either to visual or chemical cues or to feeding preferences [21,35,50]. Dispersal of small insects often depends on wind speed and direction but B. tabaci may have two morphs, those that respond to vegetative cues (visual or chemical) and those that engage in long term migration (directed by skylight and wind) [6,7]. The distribution of whiteflies within fields is often patchy [7]. Other factors may affect insect flight behavior: plant senescence, vector age, quality of the host plant, and population crowding [13].
Our goal in this investigation is to apply structured deterministic and stochastic models, ordinary differential equation (ODE) and Markov chain (MC) models, for plant, vector, and virus to predict initial growth rate and the probability of a disease outbreak when densities of infected plants or infected vectors are low. The models are structured according to healthy, infected and recovered stages for plants and healthy and infected stages for vectors. These types of models have frequently been applied to animal diseases and they have also been applied successfully to plant epidemics (e.g., [8,19,20,23,24,25,31,32,50]). Our analysis compares the dynamics of models with or without vector aggregation within a field and models with random or directed local dispersal from one field or crop to another.
We develop a modeling framework based on cassava, CMV, and the whitefly vector. Plant-vector-virus models specific to CMV include Holt et al. [20], Jeger et al. [23], and Zhang et al. [50]. We apply the ODE model formulated by Zhang et al. [50] which includes vector aggregation through the acquisition and inoculation rates. We extend this model with vector aggregation to include whitefly feeding preferences [29,35] and to account for local dispersal between adjacent crops. The basic reproduction number for the model is computed.
Only a few stochastic plant-vector-virus models have been formulated [12,33,36,43]. These latter models are spatially explicit computer simulations. Our MC models differ from these models in that they are closely related to the ODE models. An advantage of these simpler models is that they are more amenable to mathematical analysis. A branching process approximation of the MC model yields an estimate for the probability of disease extinction (when the basic reproduction number from the ODE model is greater than one). Several numerical examples of the models illustrate the effects of vector aggregation, feeding preferences, and local dispersal on the initial rate of spread and the occurrence of disease outbreaks.
2.1. ODE. We apply a plant-vector-virus model developed by Zhang et al. [50] for cassava, CMV, and the whitefly. The model includes healthy, infected, and recovered host plants and healthy and infected vectors. Let H, I, and R denote the density of the three stages of host plants and X and Z denote the two stages for the insect vector. As noted earlier, the virus is persistent, but nonpropagative [37]. For a short-lived insect, it is reasonable to assume that the insect is infected for its' entire lifetime; there is no recovery. The inoculation rate, denoted α m , is the rate at which an infected vector transmits the virus to a healthy plant and the acquisition rate, denoted β n , is the rate at which an infected plant passes the virus to an uninfected vector. The units of α m are per (host time) whereas the units of β n are per (vector time). Therefore, the rate healthy plants become infected is α m ZH and the rate healthy vectors become infected is β n IX. The total host plant population and total vector population densities are constant. Plants are continually harvested at a per capita rate r, replanted at the same rate, and all those replanted are healthy. No disease is introduced from the stem cuttings. The recovery rate of infected plants is γI. In addition, healthy and infected vectors die at the same per capita rate, cX and cZ. There is no vertical transmission of the infection among insects. With the preceding assumptions, the plant-vector-virus model takes the following form: The total population densities of plants and vectors are P and V , respectively, The plant-vector-virus model of Zhang et al. [50] is distinguished from other epidemiological models in that the model takes into account vector aggregation implicitly through the acquisition and inoculation rates. If α m and β n are treated as constants, then vectors feed at any location of the plant without preference. The aggregation of vectors on feeding locations affects inoculation and acquisition. In African cassava mosaic disease, the spatial distribution of diseased plants is clumped [10]. Therefore, Zhang et al. [50] assumed newly infected host plants are favorable to vectors. Only when there is no vacancy in the feeding sites of infected plants will vectors disperse to healthy plants. Zhang et al. [50] assumed that the inoculation term has an inverse relationship to the density of healthy plants, expressed as follows: where α is a constant, and 0 ≤ m < 1. The inoculation rate can be expressed as Similarly, vectors accumulate at ideal feeding sites such as young leaves and shoots [11]. Therefore, Zhang et al. [50] assumed interference between vectors due to aggregation decreases acquisition of the virus. As the healthy vector density increases at these sites, less space is available for new vectors to occupy. As a result, the acquisition rate will decrease as healthy vector density increases and therefore, it is expressed as an inverse relation with the abundance of healthy vectors as follows: where β > 0 is a constant and 0 ≤ n < 1 [50]. The acquisition rate can be expressed as as βIX 1−n . To distinguish between m and n, m will be called the "inoculation aggregation power" and n the "acquisition aggregation power". Applying the particular forms derived for the acquisition and inoculation rates, equations (2) and (3), and replacing R and X by R = P − H − I and X = V − Z, respectively, model (1) reduces to a system of three ODEs: Initial conditions are nonnegative with I(0) + Z(0) > 0. Therefore, solutions of model (4) are nonnegative and bounded. Applying theory from autonomous differential equations, solutions exist for t ∈ [0, ∞) [5]. The disease-free equilibrium (DFE) of model (4) is H = P , I = 0, X = V , and Z = 0. The transmission rates have the form of power functions. Such types of epidemic models have been considered by others in a more general setting (e.g., [2,18,26,30,34,40,48,49]) but for SIS or SIR models. In these latter models, the transmission rate between susceptible and infected individuals has the general form βS p I q , where p and q are strictly positive. Hochberg [18] referred to the powers p and q as the susceptible response and the infected response, respectively. For example, if p ∈ (0, 1), the response is decelerating, whereas if p ∈ (1, ∞), the response is accelerating. Applying Hochberg's description, the healthy plant and the healthy insect responses are decelerating in model (4).
For whiteflies, recent studies show that viruliferous whiteflies respond differently than non-viruliferous whiteflies to healthy and infected plants [29,35]. The virus alters vector behavior. The studies conducted with Tomato yellow leaf curl virus on tomato plants, showed that viruliferous whiteflies fed for longer periods of time and made a greater number of phloem contacts on plants than non-viruliferous whiteflies [35]. Preference assays showed that viruliferous whiteflies preferred healthy plants and non-viruliferous whiteflies preferred infected plants [29]. Also, Ingwell et al. [21] showed in a study with aphids and Barley yellow dwarf virus that viruliferous aphids preferred healthy plants over infected plants and non-viruliferous aphids preferred infected plants. Therefore, with respect to feeding preferences, the inoculation rate is increased with the density of healthy plants and acquisition rate is increased with the density of infected plants.
Based on the preceding studies, we propose an extension of the aggregation model (4 ) to include feeding preference of viruliferous vectors. If the inoculation rate were accelerating rather than decelerating with respect to the density of healthy plants, then the virus manipulates the vector to increase its spread among healthy host plants. In this extension, we allow −1 ≤ m < 1 in the inoculation power. A negative power of m implies feeding preference of viruliferous vectors on healthy plants but a positive power of m implies preference for infected plants. The power of n is on the abundance of healthy vectors (not healthy plants), so we will assume 0 ≤ n < 1.

2.2.
Markov chain. The following Table 1 defines the continuous-time MC model based on the assumptions used to derive the ODE model (4), a birth and death process. The rates r i define infinitesimal transition probabilities, r i ∆t + o(∆t), e.g., . There are 10 different events or transitions, i = 1, . . . , 10, in Table 1.
The same notation for the random variables in the MC model is applied as for the deterministic model (4), but the variables in the MC model are discrete-valued, for n ∈ [0, 1) and m ∈ [−1, 1).  Vector Infection (X, Z) → (X − 1, Z + 1) βIX 1−n 5 Infected Vector Death Healthy Host Harvest Healthy Vector Death X → X − 1 cX 10 Vector Birth X → X + 1 cV 2.3. Reproduction numbers. The basic reproduction number R 0 and type reproduction numbers T i , i = 1, 2, are defined for the ODE model (4) at the DFE, The type reproduction number, introduced by Roberts and Heesterbeek [16,38], defines a threshold for control of type i, i = 1, 2, either the infected plants or the infected vectors. First, we define the basic reproduction number using the next generation matrix [46,47]. For the infected states, I and Z, the next generation matrix is defined as follows: The basic reproduction number is the spectral radius of K [46,47], If R 0 ≤ 1, then the DFE is globally stable and if R 0 > 1, the DFE is unstable and there exists a unique globally asymptotically stable endemic equilibrium (see Appendix). The next generation matrix K = [K ij ] is used to define the type reproduction numbers [16,38]: Since K ii = 0, this leads to the same type reproduction number for the infected plants and the infected vectors, Although both the basic and type reproduction numbers are thresholds with respect to one, the magnitude of the type reproduction numbers give a better indication of the amount of control necessary to eliminate the disease by directing the control at either infected plants or infected vectors [38]. It is easy to see that there exists a relation between T and R 0 .
For host-vector models, the type reproduction number T is used more frequently than R 0 as the threshold [17]. The formulas for R 0 and T are decreasing functions of m and n, the inoculation and acquisition aggregation factors.
2.4. Probability of disease extinction. Although the case R 0 > 1 or T > 1 predicts an outbreak in the deterministic model, in a stochastic model with discrete random variables, this is not always the case. There is a positive probability that introduction of a small number of infected plants or infected vectors does not result in a disease outbreak. Introduced infected plants or vectors may either die or recover, returning the system to a DFE. The probability of disease extinction is computed from a linear multitype branching process approximation of the MC model, where the disease-free state is an absorbing state of the Markov chain. As the name implies probability generating functions generate the probabilities for new host plants or new vectors from one infected plant or from one infected vector. For example, for two random variables, the general form of an "offspring" probability generating function for a discrete random variable Y i given Y i (0) = 1 is where P i (k 1 , k 2 ) is the probability of k 1 and k 2 offspring of types 1 and 2, respectively, from one individual of type i and u i ∈ [0, 1], i = 1, 2. The assumption is that the probabilities P i do not depend on time and that the same probability applies for each individual of type i, independent of the number of individuals. More importantly, it follows from branching process theory that the fixed point of the offspring probability generating functions can be used to calculate the probability of extinction, i.e, f i (q 1 , q 2 ) = q i , i = 1, 2 [9,14,22]. Due to the assumption of independence, the probability of disease extinction, given Y i (0) = y i , i = 1, 2 is Since the linear branching process approximates the dynamics of the MC model near the DFE, the probability of extinction calculation from the branching process is an approximation of the probability of disease extinction (P 0 ) or of a disease outbreak (1 − P 0 ) in the MC model when infected plants or infected vectors are introduced into a healthy plant population. The infinitesimal transition probabilities in Table 1 for events 1-5 and their respective changes of ±1 define the probability generating functions for the two discrete random variables I and Z [3,9,27]. If one infected plant is introduced into a disease-free system, i.e. I(0) = 1 and Z(0) = 0, then the offspring probability generating function for the infected plants is According to Table 1, one infected plant either results in one infected vector (and the plant still remains infected) with probability ) or the plant dies or recovers with probability (events 2 and 3). If one infected vector is introduced into a disease-free system, i.e. I(0) = 0 and Z(0) = 1, then the offspring probability generating function for Z is

(events 1 and 5).
A unique fixed point of f 1 and f 2 exists in (0, 1) × (0, 1) if R 0 > 1 (or T > 1). The fixed point can be expressed in terms of the type reproduction number T , as defined in (6), The expression for q 1 can be interpreted as the probability an acquisition event does not result in a new infected plant plus the probability the infected plant either dies or recovers. Either events result in disease extinction. The expression for q 2 can be interpreted similarly, as the probability an inoculation event does not result in a new infected vector plus the probability the infected vector dies. If R 0 < 1, then the probability of disease extinction is one but if R 0 > 1 in a healthy crop, there is a likelihood that introduction of a few infected plants or vectors will die or recover before the disease spreads. The dependence of q i on T imply as T increases, the probability of disease extinction decreases (greater likelihood of a disease outbreak). The values of q i , i = 1, 2 decrease and the probability of an outbreak increases, as the value of m decreases such as with feeding preference of viruliferous vectors for healthy plants. Table 2. Most of the parameter values in set (i) have been applied to cassava, CMV, and the whitefly [50]. Parameter set (ii) differs from parameter set (i) in the inoculation and acquisition rates (chosen for illustration purposes). For parameter set (ii), the inoculation rate is accelerating due to feeding preferences of viruliferous vectors for healthy plants.  The endemic equilibrium (H,Ī,Z) = (5151, 970, 1878) is globally asymptotically stable (see Appendix). Given three sample paths, the disease persists in all three paths as seen in Figure 1, which describes the prevalence of the disease over 80 days. Disease extinction (absorption into the zero state) is likely to occur near the moment of introduction of the disease, during the first few days as in Figure 1. 2.6. With and without vector aggregation. We compare the dynamics of model (4) with a model that does not have aggregation. With no aggregation, model (4) simplifies to a model with bilinear incidence, n = m = 0, α 0 = α, and β 0 = β. If we let m = n = 0 and all other parameter values remain the same, then the model with no aggregation has a larger value for the basic reproduction number than when m, n > 0 and a smaller value than compared to m < 0. Hence, with aggregation as defined by Zhang et al. [50], m, n > 0, the probability of an outbreak decreases and the initial rate of growth is less than for the model with no aggregation. However, for m < 0, such as with feeding preferences, the probability of an outbreak increases and the initial rate of growth is greater than for the model with no aggregation. However, alternate methods to compare the two models with and without vector aggregation are to assume that either the incidence rate or the prevalence in the models are in agreement.

Parameter values. Two sets of parameter values are listed in
Suppose the initial incidence rate in both models are equal. That is, the value of α m ZH and β n IH are the same in both models near the DFE. The parameters  Table 2 when there is one infected plant initially. The graph on the left is over a time interval of 80 days, whereas a blow-up of the graph on the left shows disease progression over 5 days. The initial conditions are (H, I, R, X, Z)| t=0 = (9999, 1, 0, 2000, 0). The probability of disease extinction in the MC model is P 0 = q 1 = 0.095 and R 0 = 3.38.
α and β in model (4) are replaced by α 0 and β 0 in the model with no aggregation so that α 0 = αP −m and β 0 = βV −n . Then both models have the same basic reproduction number. In this case, the two models have the same initial growth rate, the same basic reproduction numbers, but the final prevalence differs in these models. The model with aggregation, with parameter values (i) in Table 2, has a greater disease prevalence in the plant and the vector population (Figure 2(a)). Aggregation is advantageous to virus persistence in this instance.   Table 2 (i), m, n > 0 and without aggregation. (a) Incidence rates are equal, α m ZH and β n IX, with and without aggregation but the endemic equilibrium values are greater with aggregation, R 0 = 3.38. (b) Disease prevalence at the endemic equilibrium is the same with and without aggregation but the rate of growth is greater without aggregation, R 0 = 3.38 and (q 1 , q 2 ) = (0.095, 0.92) versus R 0,no = 5.63 and (q 1 , q 2 ) no = (0.037, 0.86).
Next, suppose the prevalence of infection is the same in both models. That is, the two models have the same endemic equilibrium values. In this case, if the values of α and β in model (4) are replaced by α 0 = αH m and β 0 = βX n in the model with no aggregation, then the two models have the same prevalence at the endemic equilibrium, but the initial rate of growth differs in these two models. The model with no aggregation has a greater growth rate than the model with aggregation (Figure 2(b)). In this instance, aggregation is not advantageous to virus spread.
We consider one more example with parameters taken from set (ii) in Table 2, feeding preference of viruliferous vectors on healthy plants, m < 0 and n = 0 to compare with no feeding preference. Assume that the prevalence of infection at the endemic equilibrium is the same without feeding preference. In Figure 3 is the reverse of the situation in Figure 2(b). The rate of growth with vector feeding preference is greater than without feeding preference. Hence, in this case, feeding preference of healthy plants is advantageous to virus spread [21].  . Model (4) with viruliferous vectors preference for healthy plants, m < 0 and n = 0, is compared with no feeding preference, m = 0 and n = 0, where disease prevalence at the endemic equilibrium is the same with and without aggregation. The parameter set (ii) in Table 2 is applied. With aggregation, R 0 = 6.32 and (q 1 , q 2 ) = (0.033, 0.76) versus R 0,no = 5.18 and (q 1 , q 2 ) no = (0.045, 0.82) with no aggregation.

Model II: Vector aggregation and dispersal.
3.1. ODE. We include vector population dynamics that are limited by the population density of the particular host crop [20] and dispersal that is proportional to the density of the vector population. Vector movement does not depend explicitly on aggregation. We assume there are k different crops or patches, where the spatial scale is representative of the cropping areas. The total host population density in each patch remains constant, H i (t) + I i (t) + R i (t) = P i , but the total vector population, X i (t) + Z i (t) = V i (t) may change with time, i = 1, . . . , k. The vector population is bounded by η i P i , where η i is the maximum number of vectors per host plant in patch i [20]. The vector density may be directed toward or away from a particular patch, d ij , from patch i to patch j. The vector population in patch i is a solution of the following ODE: In the remaining analysis, we consider two crops or patches. In patch 1, the plants are disease-sensitive, whereas in patch 2 the plants are disease-resistant. The preceding vector population dynamics are expressed in terms of the healthy and infected vectors in two patches, k = 2. The vector population with dispersal is: for i, j = 1, 2 and i = j. In formulating the corresponding MC model, the term represents births of healthy vectors which are limited by the density of plants, In addition, we assume there is no immigration into or out of the two patches, only local dispersal between these two patches. These two simplifications, two patches and neglect of immigration or emigration, are unrealistic for the long-term dynamics but they are useful in studying the initiation of an outbreak when a small number of infected plants or infected vectors are introduced into a healthy plant population. The plant-vector-virus model with two patches is given by for i, j = 1, 2 and i = j. The initial conditions are nonnegative with i (I i (0) + Z i (0)) > 0 and V i (0) ≤η i P i , i = 1, 2. Solutions to (11) are nonnegative and bounded provided the vector population size in patch i satisfies V i (t) ≤η i P i . Under these conditions, solutions exist on [0, ∞) [5]. There exists a unique DFE, 3.2. Markov chain. The transition rates for the MC model with k patches or crops, corresponding to rates given in model (11), are summarized in Table 3. The same notation is used for each of the discrete random variables in this MC model as for the ODE model (11). The random variables Inf. Host Recovery Vector Infection Healthy Host Harvest Recovered Host Harvest Replanting Vector Death Reproduction numbers. For infected states, I 1 , I 2 , Z 1 , and Z 2 , the next generation matrix is and The basic reproduction number for model (11) is and the four type reproduction numbers equal Local stability of the DFE follows from conditions (A1)-(A5) and Theorem 1 in [47] (see Appendix). If R 0 < 1, the DFE is locally asymptotically stable and if R 0 > 1, the DFE is unstable. With no dispersal, d ij = 0, then T 1 = K 13 K 31 and T 2 = K 24 K 42 which are equivalent to the type reproduction numbers for the individual patches with  Table 3. Events 1-5 and 11 are applicable. Given an introduction of one infected plant or one infected vector into a disease-free system, (I 1 (0), I 2 (0), Z 1 (0), Z 2 (0)) = (δ 1j , δ 2j , δ 3j , δ 4j ), j = 1, 2, 3, 4, then there are four probability generating functions corresponding to each of the four random variables as follows: An explicit formula for the fixed point of the four preceding probability generating functions is not possible for this system. However, it follows from the theory of branching processes that there exists a unique minimal fixed point (q 1 , q 2 , q 3 , q 4 ) that lies in (0, 1) 4 if R 0 > 1 [4,14,15,22]. The minimal fixed point determines the probability of extinction. For the numerical examples in the next sections, the fixed point is calculated numerically. Table 2 to include a disease-resistant crop, where R 0 < 1. The parameters for patch 1, containing disease-sensitive plants, are similar to the parameters defined in Table 2. Patch 2 contains disease-resistant plants that provide a buffer from disease from the diseasesensitive crops in patch 1. Disease-resistant plants are less likely to become infected, so the inoculation rate for patch 2 is selected so that it is a factor of 0.05 the value in patch 1. We assume the inoculation rate is lower in a disease-resistant plant than in a disease-sensitive plant.

Parameter values. We extend the parameter values from
The values of the basic and type reproduction numbers, probability of disease extinction, and the DFE and the endemic equilibrium (EE) are summarized in Table 5 for d ij = 0.1. In the absence of dispersal, the disease-sensitive crop has   Figure 4 are graphs of the ODE solution and two sample paths of the twopatch MC model (11) applying parameter values in Table 4. Patch 1 contains the disease-sensitive plants and patch 2 disease-resistant. One infected sensitive plant is introduced into patch 1. The value of R 0 = 2.17 with probability of disease outbreak equal to 1 − P 0 = 1 − q 1 = 0.78.
3.6. Probability of disease extinction. With heterogeneity in type and location of crops, the direction of insect dispersal will impact whether an outbreak occurs in adjacent crops. The extinction probabilities q i are calculated as a function of the dispersal direction, d ij , in the case of the two crops, disease-sensitive and disease-resistant. In practice, the dispersal direction depends on wind and other environmental or vegetative cues. The parameters from Table 4 are applied but the dispersal rates d ij are varied.
Three cases of dispersal are considered. In the first case, insect dispersal is dominated by movement from patch 2 to patch 1 (toward the disease-sensitive  Parameter values are given in Table 4. Patch 1, disease-sensitive plants, is on the left and patch 2, disease-resistant plants, is on the right. Solutions are graphed on the interval [0, 100] in the top and on the interval [0,20] in the bottom. In both sample paths, the disease persists resulting in an outbreak. The probability of an outbreak is 1 − P 0 = 1 − q 1 ≈ 0.78. patch), with d 12 = 0.1 but d 21 an increasing function from 0.1 to 0.6. In the second case, insect dispersal is dominated by movement from patch 1 to patch 2 (toward the disease-resistant patch). In the third case, insect dispersal is equal in both directions, d 12 = d 21 . The three cases are summarized below.
(1) d 12 = 0.  Figure 5 are plots of the equilibrium densities of the healthy vector populations at the DFE for both patches, (X 1 ,X 2 ), and the corresponding values of the basic reproduction number R 0 in each of the three cases. As movement is directed toward patch j, then the healthy vector population density decreases in patch i, i = j. With increasing dispersal rates, the basic reproduction number decreases with random dispersal or movement directed toward the disease-resistant patch. The decrease is much greater when movement is directed toward the disease-resistant patch and in fact R 0 becomes less than one for values d 12 ≥ 0.5. There is no outbreak in either patch if dispersal is of sufficient magnitude and directed toward the disease-resistant patch (when R 0 < 1).  Figure 6 is a plot of the extinction probabilities q 1 , q 2 , q 3 , and q 4 as a function of the dispersal rate d ij for the three cases (1)- (3). The probability of disease extinction is P 0 = q i1 1 q i2 2 q z1 3 q z2 4 , where I j (0) = i j and Z j (0) = z j are the initial number of infected plants or infected insects. In all cases the value of q i increases most rapidly for insect dispersal in case (2), directed toward the disease-resistant patch. When R 0 < 1, the probability of disease extinction is one, that is, for d 12 ≥ 0.5 in case (2).
The values of q 1 , q 3 and q 4 have similar relations among the three cases of dispersal. When dispersal is directed toward the disease-sensitive crop, case (1), there is a decrease in probability of extinction and this probability is less than for equal dispersal, case (3). However, with dispersal directed toward the disease-sensitive crop, there is an initial decrease, but then a slight increase in the probability of extinction in the disease-resistant crop (value of q 2 ).
This last example emphasizes the effect of dispersal on outbreaks. Aggregation and feeding preferences determine within-crop spread, while dispersal determines between-crop spread. A similar effect for between-crop spread is observed for the model with feeding preferences m ∈ [−1, 0). 4. Discussion. A major concern in sustainable agriculture is the control of pests and disease. Vector aggregation and insect dispersal are two important behaviors that impact methods of control. The preferred methods of control for CMV are phytosanitation and planting of resistant varieties [11,45]. The models considered here, show that inclusion of aggregation in models implicitly through the acquisition and inoculation rates can have distinctly different effects on the initial rate of growth when inoculation rate is driven by aggregation or by preference for healthy plants (Figures 2 and 3). A better understanding of the factors affecting aggregation behavior and feeding preferences through empirical studies and models are needed to design effective control measures [21,29,33,43]. The effects of insect dispersal on disease spread also needs further studies and modeling, to determine whether environmental or vegetative cues promote disease spread and to quantify the effects of local versus long distance spread [6,7]. As illustrated in the two-patch model in Figure 6, the dispersal direction can have unexpected and undesirable consequences if the dispersal direction is from infected crops to healthy resistant crops. However, there may be more subtle effects, not considered here, that relate insect aggregation and dispersal. McElhany et al. [33] found that an increase in spatial patchiness may decrease the rate of disease spread if vector dispersal is limited. In addition, vector aggregation may affect timing of dispersal.
In this investigation, we assumed a low frequency of diseased plants (near the DFE) and a persistent virus. We compared models with and without feeding preferences or aggregation. When prevalence at the endemic equilibrium agree in models with and without aggregation, viruliferous vectors preference for healthy plants has a faster rate of spread than no preference and no preference has a faster rate of spread than preference for diseased plants (Figure 3 versus Figure 2 (b)). Sisterson [43] also found feeding preferences for healthy plants increased rates of disease spread compared with no feeding preference and McElhany et al. [33] found that persistence increased the rate of spread when vectors preferred healthy plants more than if vectors preferred diseased plants.
The reproduction numbers and extinction probabilities discussed in this investigation provide useful analytical measures for assessing the consequences of vector behavior on disease outbreaks during the early stages of infection with few diseased plants. These analytical measures, applied in conjunction with more detailed individual-based and spatially explicit models [33,43] will provide a better understanding of the relations among vector preferences, aggregation, and dispersal on disease spread and the impact of control measures in plant-vector-virus models.
Acknowledgments. This work was assisted through participation in the Vector Transmission of Plant Viruses Investigative Workshop in April 2014 at the National Institute for Mathematical and Biological Synthesis, sponsored by the National Science Foundation, the U.S. Department of Homeland Security, and the U.S. Department of Agriculture through NSF Awards #EF-0832858 and #DBI-1300426, with additional support from The University of Tennessee, Knoxville. Summer support of Hebert was funded by the Paul W. Horn Professorship account. Discussions on plant-vector-virus systems with V. Bokil contributed to model formulation. We thank the anonymous referees for many helpful comments about plant-vector-virus associations and for suggesting the Lyapunov function for global asymptotic stability.

Appendix.
Global stability of the DFE. Applying techniques in [42] (Theorems 2.1 and 2.2), it can be easily shown that the DFE of model (4) is globally asymptotically stable. Let x = (I, Z) T and y = H. Then dx/dt = (F − V)x − f (x, y), where the next generation matrix K = FV −1 , Let w = ( βV 1−n /c, αP 1−m /(γ + r)), the left eigenvector of V −1 F corresponding to eigenvalue R 0 . The following Lyapunov function L(x, y) = wV −1 x can be used to show global asymptotic stability of the DFE (Theorem 2.2, [42]).
Existence and global stability of the EE. For model (4), the EE is a solution of Expressing the second equation as a function of H, The functions g 1 and g 2 have the following properties for H ∈ (0, P ) and Z ∈ (0, V ): g 1 (P ) = 0, g 1 (H) > 0, g 1 (H) < 0, g 1 (H) > 0, g 2 (H) > 0, g 2 (H) < 0 (see Figure 7). There exists a unique point of intersection of these two curves in the interior of the rectangular region [0, P ] × [0, V ] of the H-Z plane if g 2 (H 0 ) = 0 for H 0 < P . This latter inequality ensures that g 2 crosses the Z = 0 axis to the left of H = P and it is equivalent to the inequality R 0 > 1. Hence, a unique EE exists to model (4), if R 0 > 1. To show global stability of the EE for model (4) when R 0 > 1, we construct a Lyapunov function. The positive EE is denoted as (H * , I * , X * , Z * ). The equilibrium values for Z and I can be expressed as Z * = c 1 c and I * = c 2 γ + r , where c 1 = βI * (X It is easy to see that the first two expressions inL are negative on R 4 + if H = H * and X = X * . By the arithmetic-geometric mean inequality, the third expression is negative on R 4 , except at the EE, whereL = 0. The function L satisfies the criteria for a Lyapunov function and hence, the EE is globally asymptotically stable in R 4 + when m = 0 and n = 0 [28]. To prove global stability of the EE in the case m = 0, the term −(H * ) 1−m (H m − (H * ) m )/m in L is replaced by −H * ln(H/H * ). ThenL is given by the expression in (15) with m = 0. Similarly if n = 0, the term −(X * ) 1−n (X n − (X * ) n )/n in L is replaced by −X * ln(X/X * ). Two-patch model. With no immigration into or emigration out of the k patches, then Let T (t) = k i=1 V i (t) be the total vector population. Then The total vector population satisfies where a = min{a i } and b = max ai ηiPi . The plant population size is bounded below. From comparison theory for scalar differential equations [5], it follows that Sufficient conditions for the vector population in the two-patch model (11) to be bounded above byη i P i are that each of dV i /dt evaluated atη i P i are nonpositive for i = 1, 2. These latter conditions are satisfied for two patches if For the case of two patches, a unique DFE exists for the vector population. The equilibrium values X i =X i > 0, i = 1, 2 are solutions of dX i /dt = 0 in model (11): where a i and η i are defined in (10). The equations in (17) can be expressed as follows: Equations (18) and (19) are parabolas and continuous in the X 1 -X 2 plane. The DFE is the unique point (X 1 ,X 2 ) in the interior of R 2 + where the two curves intersect [15].
In the absence of disease, Z i = 0 and I i = 0, we show the DFE is globally asymptotically stable. From the differential equations for H i it follows that H i → P i . From the differential equations for X 1 and X 2 , there are only two equilibria, the origin and (X 1 ,X 2 ). The origin is a saddle (stable manifold outside the first quadrant), solutions are bounded, and there does not exist any periodic solutions (Dulac's criteria ∂(f 1 /(X 1 X 2 ))/∂X 1 + ∂(f 2 /(X 1 X 2 ))/∂X 2 < 0). Poincaré-Bendixson theory implies global stability of the DFE in the absence of disease (condition (A5) in [47]).