Bifurcation analysis and transient spatio-temporal dynamics for a diffusive plant-herbivore system with Dirichlet boundary conditions

1.   Introduction

Cassava Mosaic Disease (CMD) is a plant disease that reduces tuber size and starch percentage upstream in the cassava supply chain, reducing sales of the cassava crop ^[1]. This leads to downstream economic impacts since cassava is a major industrial raw material. The financial losses due to CMD in the African continent were estimated at ＄1.2-2.3 billion in 1997 ^[2], increasing to ＄1.9-2.7 billion in 2009 ^[3]. A CMD outbreak can be controlled in four ways: (1) by removing infected cassava from the plantation, (2) by promoting the use of virus-free cuttings, (3) by using resistant varieties, or (4) by killing infected whitefly in the plantation area ^[4,5,6]. These methods incur huge planning effort and expense. However, there is no clear approach to mapping the major spread factors for appropriate control planning. In this study, we analyzed the CMD spread factors and proposed an epidemic model to isolate those factors that have a major impact in an outbreak. The output may be used to formulate optimal strategies for outbreak control.

Using such an epidemic model, the authors of ^[7] showed how tomato yields in India suffered from whitefly-vectored Tomato Leaf Curl Geminivirus (TLCV). The most effective disease-control strategy was to distribute netting treated with persistent insecticide. This significantly controlled spread by decreasing vector immigration to the plantation and increasing vector mortality. The authors of ^[8] studied the outbreak of Huanglongbing disease (HLB) caused by psyllid bacteria, which seriously impacts citrus greening. They tested different intervention strategies over a short time-frame and showed that the best insecticide intervention was to spray over an optimal number of days.

Most epidemic models in the literature are based on the pioneering work in ^[9]. The models identify the major factors that contribute to the outbreak by applying differential equations. This yields the basic reproduction number ($ \mathscr{R}_0 $), used to indicate the severity of spread ^[10,11,12]. Analysis of the sensitivity of $ \mathscr{R}_0 $ determines which parameters have the greatest impact on disease control and should be targeted by prevention policies ^[13,14].

Most viral crop plant diseases radiate in vectors ^[15,16]. The models must therefore represent the dynamics of plant-virus epidemics in order to reveal relationships between vectors and infected plant numbers ^[17,18,19]. In 1994, African cassava mosaic virus (ACMV) caused widespread loss of production. The authors in ^[20] incorporated vector-population dynamics into their statistical model to empirically derive plant-virus relationships. In ^[21], an epidemiological model of ACMV was introduced to describe the dynamics of infected cassava and infective whitefly. The authors derived a strategy in which cassava yields are maximized by reducing the whitefly population. However, this approach is not cost-effective, so farmers are advised to select uninfected cuttings for planting to prevent a collapse of the healthy cassava population. This is an economical strategy that is capable of controlling an outbreak. The authors ^[22] developed a mathematical model for the spread of Cassava Brown Streak Disease (CBSD) caused by whitefly, using an approach similar to that in ^[21]. They simulated two control policies: uprooting and burning of infected cassava, and killing of whitefly in the plantation. They concluded that the optimal policy was uprooting and burning.

In previous works ^[21,22], the same authors formulated a mathematical model for vector-host dynamics using a single spread factor (whitefly). However, CMD infection requires both whitefly (Bemisia tabaci) transmission and the presence of virus-infected cuttings ^[23,24]. Outbreaks have spread beyond the African continent to countries in Asia including India and Sri Lanka following imports of virus-infected cuttings. The CMD model must incorporate this route when devising optimal strategies for outbreak control. The mathematical model developed in ^[25] for an outbreak of ACMV took account of diseased cuttings, following ^[21]. The model suggested that a strategy combining rouging and spraying performed better than those that apply a single control mechanism.

When infection levels are high, symptomatic cassava can be readily detected and removed from the plantation. However, the symptoms take five weeks to appear ^[20], and during this period the asymptomatic cassava may spread the disease. Our goal was to develop a better understanding of the relationship between asymptomatic plants and the severity of CMD spread. Our model added a latent state to that in ^[25], in an effort to analyze the comparative contribution of whitefly transmission and planting of infected cuttings. We applied sensitivity analysis to identify the most important parameters in determining the severity of an outbreak.

The rest of this paper is organized as follows. Section 2 discusses CMD, the construction of our model, and the most important parameters in CMD spread. In section 3 we introduce the positivity and boundedness of solutions to initial conditions. Section 4 concerns the stability of the disease-free and endemic-equilibrium points, and introduces the threshold parameter ($ \mathscr{R}_0 $), obtained using the Next-Generation method. The sensitivity analysis described in section 5 identified the most sensitive parameter with respect to $ \mathscr{R}_0 $, key to improving the design of CMD-control policies. We conducted numerical simulations to explore the behavior of the system with and without control policies. The simulation is introduced in section 6, and the results reported in section 7. Section 8 presents a summary and our conclusions.

2.   Construction of CMD outbreak model

Few mathematical models of viral cassava disease (East African Cassava Mosaic Virus (EACMV), ACMV, or CBSD) ^[21,22,26] have analyzed how the disease spreads within a plantation. To understand the dynamics of an outbreak, we used ordinary differential equations (ODEs).

The dynamics of ACMV spread were analyzed in ^[21] and ^[25], and the "SI-type epidemic model" was developed, which categorizes plants as infected or uninfected. In ^[26], an epidemic model was developed in which two viruses, ACMV and EACMV, are carried by whitefly into the plantation. The model was based on the "m-group SI-type epidemic model" of ^[27], which incorporates four plant states: uninfected, infected by ACMV, infected by EACMV, or infected by both. The authors of ^[22] developed an epidemic model based on the SEIR-type model of ^[28], and used it to explore the dynamics of CBSD. It assumed four states: healthy, latent, exposed, and rouging. The latent stage, which was missing from earlier models, increases the transmission rate from infected plants to whitefly and vice versa. Our proposed model takes account of the spread of CMD by both by planting and whitefly transmission. The design process is comprised of three steps:

(1) The dynamics are based on infection by transplanting and by whitefly transmission ^[23,24]. Defining these factors as parameters yields the ODEs, which we use to calculate equilibrium points and to establish the stability of the outbreak. The Next-Generation method is used to construct $ \mathscr{R}_0 $.

(2) The local stability of the disease-free equilibrium is proved by the Routh-Hurwitz criterion. The system is locally asymptotically stable if $ \mathscr{R}_0 < 1 $ and unstable if $ \mathscr{R}_0 > 1 $. Global stability of the disease-free equilibrium and the endemic equilibrium point are proved by Lyapunov's method and LaSalle's invariance principle. When $ \mathscr{R}_0 \leq 1 $, the disease-free equilibrium is globally asymptotically stable. If $ \mathscr{R}_0 > 1 $ then the endemic equilibrium is globally asymptotically stable.

(3) The use of sensitivity analysis to determine the parameters that most strongly affect the outbreak, and the use of these parameters to simulate outbreak control by addressing risk factors.

2.1.   Parameters

In our model, the parameters in the CMD outbreak represent disease spread factors. The values and ranges of the parameters were set with reference to previous studies or by analysis of CMD outbreaks. The parameters are listed in Table 1.

2.2.   Ordinary differential equations

The model tracks the dynamics of the cassava and whitefly populations. Cassava may be uninfected $ X $, latent $ Y $ (infected but asymptomatic), or infected and symptomatic $ Z $. Whitefly may be infected $ R $ or uninfected $ Q $. The total cassava population $ N $ is therefore $ N = X+Y+Z $ at time $ t $. Similarly, the total whitefly population $ N_V $ is $ N_V = Q+R $.

The transmission dynamics of CMD from an initial state to a second state are governed by the parameters in Table 1. Infection by whitefly is given by $ a p_1 \frac{X}{N}R $. Infection of whitefly from infected plants is given by $ a p_2 Q \frac{Y + Z}{N} $. CMD infection by planting of infected cuttings is $ r\left(1-\frac{N}{K}\right) p_3 Y $. This term was not used in previous works.

The size of the cassava population is increased by replanting into the plantation area. $ r\left(1 - \frac{N}{K}\right)X $ in Eq (2.1) and $ r \left(1 - \frac{N}{K}\right)p_3Y $ in Eq (2.2) represent the replanting terms of healthy and latent cassavas. Cassavas are removed from the system by constant harvesting ($ h $) and rouging ($ \gamma $). The change from state $ X $ (non-infected) to $ Y $ reflects the plants that become infected after planting. The change from $ Y $ to $ Z $ reflects the number of infected cuttings that begin to show CMD symptoms during the period.

The whitefly population is driven by two factors: birth rate $ (\Lambda) $ and death rate $ (\mu) $. Transmission within the population causes a change from state $ Q $ to state $ R $, as whitefly visit infected plants and acquire the virus. Equations (2.6) and (2.7) give the total cassava and whitefly populations:

The model makes the following assumptions:

(1) All model parameters are positive.

(2) The growth rate of the cassava population is positive, i.e., $ r-h > 0 $, where $ r $ is the replanting rate and $ h $ is the harvesting rate.

(3) The cassava population increases by one of the logistic growth equations, $ r \left(1-\frac{N}{K}\right)X $ or $ r\left(1-\frac{N}{K}\right)p_3Y $.

(4) The whitefly birth rate $ (\Lambda) $ is constant.

3.   Basic properties of the model

To confirm the biological validity of the model, we must prove that solutions to the ODEs are positive and bounded for all time values. Furthermore, the cassava and whitefly populations must remain finite since they are bounded by the plantation area. In the next section we analyze our model and demonstrate the positivity and boundedness of the system of ODEs.

3.1.   Positivity

Theorem 3.1. Let $ t \geq 0 $. In the model, if the initial conditions satisfy

then $ X(t), Y(t), Z(t), Q(t), \mathit{\text{and}}R(t) $ will remain positive in $ \mathbb{R}_+^5 $.

Proof. Let us consider $ Y(t) $ for all $ t \geq 0 $.

From Eq (2.2), we have

The integration of the inequality is

From the initial conditions (3.1), we have $ Y(t) > 0 $ for all $ t \geq 0 $.

Next, we consider $ Z(t) $ for all $ t \geq 0 $. From Eq (2.3),

A comparison argument shows that

Since $ Z(0) > 0 $, we conclude that $ Z(t) $ is always positive for all $ t \geq 0 $.

We prove that $ Q(t) > 0 $ for all $ t \geq 0 $ by contradiction. Let $ t_1 $ be the first point that satisfies $ Q(t_1) = 0 $. Equation (2.4) yields

This means that $ Q(t) < 0 $ for $ t \in (t_1-\theta, t_1) $, where $ \theta $ is an arbitrarily small positive constant. This leads to a contradiction. Hence, $ Q(t) > 0 $ for all $ t \geq 0 $.

Next, we prove that $ R(t) $ is positive for all $ t \geq 0 $. From Eq (2.5),

it then follows that

Since $ R(0) > 0 $, we know that $ R(t) $ remains positive for all $ t \geq 0 $.

Finally, we prove that $ X(t) > 0 $ for all $ t \geq 0 $. From Eq (2.1), we have

The integration of the inequality is

From the initial conditions (3.1), we have $ X(t) > 0 $ for all $ t \geq 0 $.

We conclude that the solutions of the model are positive in $ \mathbb{R}_+^5 $.

3.2.   Boundedness

We show the boundedness of the system with the initial conditions (3.1).

Theorem 3.2. All solutions of the system (2.1) to (2.5) with positive initial conditions (3.1) are ultimately bounded.

Proof. From Theorem 3.1, the solutions to the system are positive for all $ t \geq 0 $. Since $ N = X + Y + Z $. From Eq (2.6), we have

As can be observed, the solution is bounded by logistic growth

The integration of the inequality is

(assuming $ N(0) > 0 $), and hence $ \lim_{t \to \infty} N(t) \leq K $. This gives the feasible solution set of the cassava entering the region:

Next, note that $ N_V = Q + R $. From Eq (2.7),

The integration of the equation is

(assuming $ N_V(0) > 0 $), which implies that $ \lim_{t \to \infty} N_V(t) = \frac{\Lambda}{\mu} $. Thus, the feasible solution set for the CMD system is given by

The solutions of the system in $ \mathbb{R}_+^5 $ are confined to the region $ \Omega $. Here

Therefore, all solutions of the system (2.1) to (2.5) with initial conditions (3.1) remain positive invariant in the region $ \Omega $ for all $ t \geq 0 $.

4.   Equilibrium points and stability analysis

To ensure that disease eradication is independent of the initial size of the susceptible population, locally- and globally-stable disease-free equilibria and endemic equilibria are established when the region $ \Omega $ is positivity invariant.

4.1.   Equilibrium points

Our model admits two equilibrium points: a disease-free equilibrium point (DFE) and an endemic equilibrium point (EE). The notation of DFE is

The EE is denoted $ E_1 = (\bar{X}, \bar{Y}, \bar{Z}, \bar{Q}, \bar{R}) $, where

4.2.   Basic reproduction number

The basic reproduction number ($ \mathscr{R}_0 $) is the one of the most useful threshold parameters in epidemiology. It is defined as the expected number of secondary cases produced by a single infection in a completely susceptible population. It is used as an indicator of the stability of $ E_0 $ and $ E_1 $. The DFE $ E_0 $ is asymptotically stable if $ \mathscr{R}_0 < 1 $, as the disease cannot invade the population, but unstable if $ \mathscr{R}_0 > 1 $. $ E_1 $ is asymptotically stable if $ \mathscr{R}_0 > 1 $. We use the Next-Generation method, following ^[31,32], to determine $ \mathscr{R}_0 $.

The appearance of new infections is represented by the vector $ F $ and the transfer of existing infections by the vector $ V $. Let $ x_j $ be an infection state for $ j = 1, 2, 3 $, i.e., $ x_1 = Y, x_2 = Z, x_3 = R $. $ F $ describes new infections arising in state $ x_j $ and $ V $ represents the transfer of existing infections to state $ x_j $.

The Jacobian matrices generated by differentiating $ F $ and $ V $ with respect to the relevant subset of variables are calculated at $ E_0 $. This yields the matrices $ \mathbb{F} $ and $ \mathbb{V} $. The ($ j, k $) entry of the matrix $ \mathbb{F} $ is the rate at which infected individuals in compartment $ k $ produce new infections in compartment $ j $. The ($ j, k $) entry of the matrix $ \mathbb{V} $ represents the transfer of existing infection.

$ \mathscr{R}_0 $ is computed from the spectral radius of $ \mathbb{FV}^{-1} $ at DFE. $ \mathbb{FV}^{-1} $ is called the Next Generation Matrix ^[10] and is set as follows:

where $ \rho(\mathscr{M}) $ denotes the spectral radius of a matrix $ \mathscr{M} $. The spectral radius of $ \mathbb{FV}^{-1} $ is equal to the dominant (or maximum) eigenvalue of $ \mathbb{FV}^{-1} $.

In our model, the vectors $ F $ and $ V $ can be derived from Eqs (2.1)-(2.5):

The Jacobians of $ F $ and $ V $ with respect to the infectious classes are calculated as $ \mathbb{F} $ = $

$$\begin{bmatrix} \frac{\partial F_j(E_0)}{\partial x_k} \end{bmatrix}$$ $ and $ \mathbb{V} $ = $ $$\begin{bmatrix} \frac{\partial V_j(E_0)}{\partial x_k} \end{bmatrix}$$ $, for $ j = 1, 2, 3 $ and $ k = 1, 2, 3 $. We then have

Therefore, the Next Generation Matrix $ \mathbb{FV}^{-1} $ is

Hence, $ \mathscr{R}_0 $ is the dominant eigenvalue of Matrix (4.1).

4.3.   Local stability analysis of disease-free equilibrium point

Theorem 4.1. Under the hypotheses $ r-h > 0 $ and $ 0 \leq p_3 \leq 1 $, if $ \mathscr{R}_0 < 1 $ then $ E_0 $ is locally asymptotically stable; otherwise it will be unstable.

Proof. The local asymptotic stability is determined based on the eigenvalue ($ \lambda $) of the Jacobian at $ E_0 $. The $ E_0 $ is locally asymptotically stable if the real parts of $ \lambda $ are all negative. The Jacobian matrix at $ E_0 $ is given by

The characteristic equation is given by

with coefficients

and

Under the model assumption $ r > h $, it is clear that the first two eigenvalues of this system are negative, $ -\mu $ and $ -(r-h) $, respectively.

We now consider the cubic equation from (4.4)

As $ 0 \leq p_3 \leq 1 $, we obtain

We can observe that if $ \mathscr{R}_0 = \frac{hp_3}{2(\beta + h)} + \sqrt{\left(\frac{hp_3}{2(\beta + h)}\right)^2 + \frac{ra^2 \Lambda p_1p_2(\beta + \gamma + h)}{\mu^2K(r-h)(\gamma + h)(\beta + h)}} < 1 $, then

Therefore, we obtain

Finally, we consider

Under hypotheses $ \mathscr{R}_0 < 1 $, $ r > h $, and $ 0 \leq p_3 \leq 1 $, we obtain

with coefficients

This means that $ a_1a_2 - a_3 > 0 $ whenever $ \mathscr{R}_0 < 1 $, $ r > h $, and $ 0 \leq p_3 \leq 1 $.

The Routh-Hurwitz criterion for the third order polynomial indicated that the remaining eigenvalues of the disease-free equilibrium system have negative real parts when $ a_1 > 0, a_3 > 0, $ and $ a_1a_2 > a_3 $.

Therefore, $ E_0 $ will be locally asymptotically stable if $ \mathscr{R}_0 < 1 $. When $ \mathscr{R}_0 > 1 $, $ E_0 $ is unstable and the disease will persist.

4.4.   Global stability analysis

The global stability of the DFE and the EE are established using Lyapunov's method and LaSalle's invariance principle to obtain the control condition under which disease can be eradicated.

Theorem 4.2. (1) If $ \mathscr{R}_0 \leq 1 $, then the disease-free equilibrium point $ E_0 $ is globally asymptotically stable in $ \Omega $. (2) If $ \mathscr{R}_0 > 1 $, then the endemic equilibrium point $ E_1 $ is globally asymptotically stable in $ \Omega $.

Proof. Theorems 4.2 (1) and 4.2 (2) can be proved by using a Lyapunov function. We adopt the Lyapunov function used in ^[33,34].

First, we prove Theorem 4.2 (1) by assuming that:

We define the Lyapunov function

by

When $ L_1 $ is $ C^1 $, a proper positive definite function, and $ E_0 $ is the global minimum of $ L_1 $ on $ \Omega $, we obtain

The time derivative of $ L_1 $ computed along solutions of Eqs (2.1)-(2.5) is

Substituting the ODEs (2.1) to (2.5) into Eq (4.8) yields

this simplifies to

Since $ Y+Z = N-X $, $ N = X^* = \frac{(r-h)K}{r} $, $ N_V = Q^* = \frac{\Lambda}{\mu} $, and observing the condition in Eq (4.6), Equation (4.9) can be rewritten as

where $ G_3 = \frac{ap_2}{N} Q^* $. Therefore, Eq (4.10) becomes

All terms in Eq (4.11) are always non-positive since $ \lim_{t \to \infty} X(t) \geq X^* $ (Eq (3.6)) and $ 0 \leq p_3 \leq 1 $ (Table 1). Note that $ \frac{dL_1}{dt} = 0 $ if and only if $ X = X^* $, $ Q = Q^* $, $ Y = 0 $, and $ Z = 0 $.

Therefore, the largest invariant compact invariant set in $ \left\lbrace(X, Y, Z, Q, R) \in \Omega \mid \frac{dL_1}{dt} = 0 \right\rbrace $ is the singleton, i.e., $ E_0 = (X^*, Y^*, Z^*, Q^*, R^*) $. Thus, $ \frac{dL_1(t)}{dt} \leq 0 $. This implies that $ E_0 $ is globally asymptotically stable in $ \Omega $.

We next prove Theorem 4.2 (2). We define the Lyapunov function

by

When $ L_2 $ is $ C^1 $, a proper positive definite function, and $ E_1 $ is the global minimum of $ L_2 $ on $ \Omega $, we obtain

The time derivative of $ L_2 $ computed along solutions of Eqs (2.1)-(2.5) is

Under Eqs (2.1)-(2.5), if the system meets $ E_1 $,

We rewrite this to obtain

The derivative of $ L_2 $ then becomes

Using the arithmetic-geometric means inequality, we can see that $ \frac{dL_2}{dt} \leq 0 $. Note that, $ \frac{dL_2}{dt} = 0 $ only if $ X = \bar{X} $, $ Y = \bar{Y} $, $ Z = \bar{Z} $, $ Q = \bar{Q} $, and $ R = \bar{R} $.

Therefore, the largest compact invariant set in $ \left\lbrace(X, Y, Z, Q, R) \in \Omega \mid \frac{dL_2}{dt} = 0 \right\rbrace $ is the singleton $ E_1 = (\bar{X}, \bar{Y}, \bar{Z}, \bar{Q}, \bar{R}) $. Thus, $ \frac{dL_2}{dt} \leq 0 $. This implies that $ E_1 $ is globally-asymptotically-stable in $ \Omega $.

We have established that the disease-free equilibrium point $ E_0 $ is globally asymptotically stable and that the disease can be controlled as long as the threshold $ \mathscr{R}_0 \leq 1 $. If $ \mathscr{R}_0 > 1 $ the endemic equilibrium point $ E_1 $ is globally-asymptotically-stable and the disease will persist. From an epidemiological point of view the goal of policy is to control CMD outbreaks by reducing $ \mathscr{R}_0 $ to below one and maximizing the uninfected population. However, an agricultural viewpoint would instead focus on maximizing the economic returns. We therefore need to identify the strategy that maximizes profits. To do this we first identify the parameter that plays the greatest role in the stability of $ \mathscr{R}_0 $.

5.   Sensitivity analysis

Sensitivity analysis is used to assess the relative impact of different factors on the stability of a model under data uncertainty. The analysis is also able to determine which parameters play critical roles. We perform the analysis by calculating the sensitivity indices of the basic reproduction number, $ \mathscr{R}_0 $ to the parameters in the model using both local and global methods.

5.1.   Local sensitivity analysis

The local sensitivity analysis is based on the normalized forward sensitivity index $ \mathscr{R}_0 $. The sensitivity index of $ \mathscr{R}_0 $ with respect to the parameters in our model are derived as follows ^[35]:

Here, $ P $ is a parameter value from Table 2, while the value of $ \mathscr{R}_0 $ is computed from Eq (4.2). We analyze the sensitivity of $ \mathscr{R}_0 $ by substituting the parameter values into Eq (5.1). For example, the sensitivity index $ \mathscr{R}_0 $ with respect to $ h $ is

Table 2 lists the sensitivity indices of $ \mathscr{R}_0 $ obtained using Maple.

We first consider the pathogen parameters $ p_1, p_2 $ and $ p_3 $. Our model limits the severity of an outbreak by controlling three infection terms:

(1) Transmission by whitefly, $ a p_1 \frac{X}{N} R $.

(2) Infection of whitefly by CMD virus from infected cassava, $ a p_2 Q \frac{Y + Z}{N} $.

(3) CMD spread by replanting of infected cuttings, $ r \left(1 - \frac{N}{K}\right) p_3 Y $.

Whitefly-cassava transmission $ (p_1) $ and cassava-whitefly transmission $ (p_2) $ were given the same probability, and these parameters had the highest sensitivity index magnitude of +0.4965 (Table 2). This suggested that the most important determinants of the severity of a CMD outbreak are transmission by whitefly to cassava and transmission to whitefly from infected plants.

The model also assigned a probability to CMD spread by planting of infected cuttings $ (p_3) $. As can be seen from Table 2, this was the parameter with the lowest sensitivity index magnitude of + 0.0071. Replanting of infected cuttings does increase the number of plants exposed to the virus, but contributes less to the spread than does whitefly transmission. Control of the whitefly population is therefore the key policy goal.

We next analyzed the parameters that determine the whitefly population. The whitefly death rate $ (\mu = -0.9929) $ and the average number of cassava plants visited $ (a = +0.9929) $ turn out to be significant. The birth rate of whitefly also affects the stability of the system $ (\Lambda = +0.4965) $. As $ \mu $ is negative, the whitefly population will decrease as $ \mu $ increases. However, the sensitivity index values of $ a $ and $ \Lambda $ are positive, so that $ \mathscr{R}_0 $ will increase as these parameters increase. The significance of $ \mu $, $ a $ and $ \Lambda $ in determining $ \mathscr{R}_0 $ (Table 2) suggests that whitefly elimination is the optimal approach.

5.2.   Uncertainty and global sensitivity analysis

Global sensitivity analysis is used to examine the model which responses to parameter variation within a wider range in the parameter space. The parameter values and ranges are listed in Table 1. Partial rank correlation coefficient (PRCC) between the $ \mathscr{R}_0 $ and each parameter for the model are derived as in references ^[36,37,38]. We computed the PRCC by setting of input parameter values sampled using the Latin Hypercube Sampling (LHS) method. Table 3 shows the results of PRCC and P-value from 1000 independent simulations.

Parameters with absolute maximum PRCC values as well as corresponding lowest P-values are the most critical factor in disease spread. As shown in Table 3, we noticed that the parameter $ a $ has the most effect on $ \mathscr{R}_0 $, followed in decreasing order by $ p_2 $ (cassava-whitefly transmission), $ \mu $ (death rate of whitefly), and $ p_1 $ (whitefly-cassava transmission). The most important parameters for $ \mathscr{R}_0 $ from the global sensitivity analysis match those from the local sensitivity analysis, i.e., $ a $, $ p_2 $, $ \mu $, and $ p_1 $. These results allow us to considerably reduce $ \mathscr{R}_0 $ by decreasing the whitefly population.

Local and global sensitivity analysis are used to examine the impact of parameters on $ \mathscr{R}_0 $. We found that killing whitefly is the optimal policy. In practice, whitefly density in plantation area is also a factor that affects to the severity of CMD outbreaks. We therefore investigated the sensitivity of $ \mathscr{R}_0 $ in Eq (4.2) to vary the different sizes of the plantation area $ (K) $.

Figure 1a shows that the sensitivity of $ \mathscr{R}_0 $ to $ a $ increases as $ a $ increases and $ K $ decreases. Figure 1b shows that the sensitivity of $ \mathscr{R}_0 $ to $ \Lambda $ increases as $ \Lambda $ increases and $ K $ decreases, which means that $ \mathscr{R}_0 $ is more sensitive to the whitefly birth rate if the plantation area is small. In contrast, the sensitivity of $ \mathscr{R}_0 $ to $ \mu $ increases as $ \mu $ and $ K $ decrease. This means that CMD outbreaks are more severe when $ K $ is small, and a whitefly elimination policy will fail to eliminate the disease. A cost-effective control strategy must therefore take account of the plantation area. To determine the most cost-effective policy, we conducted numerical simulations.

6.   Numerical simulations

To determine the effectiveness of a control policy, it is necessary to assess the number of whitefly in the plantation. Equations (2.4) and (2.5) are modified by introducing two control parameters: $ \epsilon\in[0, 1] $, the efficacy of pesticide spray in controlling whitefly, and $ g $, the resulting whitefly death rate. Adding terms $ (\mu + g\epsilon) $ to Eqs (2.4) and (2.5) produces Eqs (6.1) and (6.2). These take account of deaths of uninfected and infected whitefly from natural causes and from pesticide spraying.

6.1.   Simulation of CMD outbreak dynamics without control

We used Eqs (2.1)-(2.3), (6.1) and (6.2) to represent the dynamics of a CMD outbreak when no control is applied $ (g\epsilon = 0) $ and to determine the resulting cassava yield. We used two plantation capacity scenarios ($ K $) to simulate the dynamics of such a CMD outbreak: (1) $ K = 15, 000 $ and (2) $ K = 10, 000 $. The other parameters were given the values in Table 1. Figures 2-5 assume a $ 300 $ d growing season, no control policy, and initial conditions $ X(0) = 9, 000 $, $ Y(0) = 1, 000 $, $ Z(0) = 0 $, $ Q(0) = 900 $, and $ R(0) = 100 $.

6.1.1.   Scenario 1: maximum plantation capacity 15, 000 $ m^2 $

Figure 2 shows the number of infected whitefly increasing to a maximum at D 17, then decreasing. Figure 3 shows the numbers of healthy, latent, and symptomatic cassava. Latent and asymptomatic plants outnumbered healthy plants at D 19, or 2 d after the whitefly infection peak. Comparing Figure 2 with 3, the number of infected plants tracked the whitefly population. The impact was mainly on the rate of latent cassava, which exceeded the number of healthy cassava. The harvest of 1, 173 tubers from an initial planting of 9, 000 gave a daily rate of $ h = 0.003 $.

6.1.2.   Scenario 2: maximum plantation capacity 10, 000 $ m^2 $

The numerical results in Figure 5 show that latent and symptomatic plants outnumbered healthy plants at D 17, lagging the trend in infected whitefly. Infected whitefly reached a maximum at D 17 (Figure 4), then decreased. By the time of harvest at D 300, 714 tubers out of 9, 000 remained healthy. In the smaller plantation, this increased whitefly density and CMD spread.

6.2.   Simulations and cost-effectiveness analysis

Our sensitivity analysis showed that the rapidity of CMD spread was driven primarily by the whitefly population. A simple approach to CMD control is to spray pesticide. However, this may not be economical as more efficient pesticide is also more costly. It is therefore necessary to analyze the cost-effectiveness of spraying as well as its effectiveness.

We simulated the three control strategies shown in Table 4, representing pesticide sprays with effectiveness of $ 5\% $, $ 10\% $ and $ 15\% $. A whitefly control target $ (g) $ of $ 5\% $ per day was set. We then determined the cost-effectiveness of using each spray. Table 4 shows the operating costs. Two kilograms of pesticide are sprayed once a week, requiring two laborers at a cost of $ ＄10 each per 7, 500 $ m^2 $. The operating costs are shown in the final column. The Cost-Effectiveness Ratio (CER) of ^[39] is then applied. CER is the ratio of total control costs to yield loss averted, and is used to measure the economic value of an intervention (Eq (6.3)).

6.2.1.   Scenario 1: maximum plantation capacity 15, 000 $ m^2 $

The simulation used the parameter values given in Table 1. The numerical results in Figures 6 and 7 show that, when the three control strategies were applied, the increase in the healthy cassava population and decrease in the infected whitefly population tracked the efficiency of the pesticide spray. This gave a yield of 1173 tubers when no control strategy was applied. Infection of 518, 1221 and 1851 tubers were averted when Type I, II and III sprays were used. This represents a significant increase in profit. We then used CER to compare the cost-effectiveness of the three control strategies.

The CERs were calculated from Eq (6.3):

Table 5 compares the cost-effectiveness. The sixth column gives the CER of each approach for $ K = 15, 000 $, where a smaller number represents a more desirable outcome. As can be seen, the Type II spray yielded the greatest overall economic benefit. This result provides policymakers with a useful tool for optimizing their control strategies.

6.2.2.   Scenario 2: maximum plantation capacity 10, 000 $ m^2 $

We next simulated the strategy with plantation capacity $ K = 10, 000 $, using the same parameters. The results are shown in Figures 8 and 9. The harvests of healthy tubers were 714 with no control, 868 with Type I spray, 1163 with Type II, and 1491 with Type III. Equation (6.3) was used to identify the most cost-effective policy.

The final column of Table 5 confirms that Type II spraying yielded the greatest economic benefit. However, this strategy was less cost-effective when applied to the smaller plantation as the whitefly density was higher.

7.   Discussion

In this study, we modeled control of a CMD outbreak driven both by whitefly transmission and planting of infected cuttings. From an epidemiological perspective, $ \mathscr{R}_0 < 1 $ is a sufficient condition for eradication of the disease from a cassava population. The Routh-Hurwitz criterion confirmed that $ E_0 $ was locally-asymptotically-stable when $ \mathscr{R}_0 < 1 $, becoming unstable when $ \mathscr{R}_0 > 1 $. $ E_0 $ was globally-asymptotically-stable if $ \mathscr{R}_0 \leq 1 $, while $ E_1 $ was globally-asymptotically-stable if $ \mathscr{R}_0 > 1 $, as confirmed by application of Lyapunov's method and LaSalle's invariance principle.

Our goal was to clarify the vectors that determine the severity of a CMD outbreak. To achieve this, we modeled increases in latent cassava due to whitefly transmission and due to planting of infected cuttings. By identifying the parameters that most affect the stability of the system, we hoped to optimize policy design. A sensitivity analysis determined that the most critical factor in disease spread was the presence of whitefly. We then sought to identify the most cost-effective pesticide for whitefly elimination. Numerical simulations showed that an increase in the death rate of whitefly through spraying was associated with a reduction in the CMD infection rate, and that the efficacy of the pesticide was key to control of the whitefly population. However, the increase in yield had to be offset against the cost of spraying.

We compared the cost-effectiveness of three pesticide sprays to determine the optimal balance between input cost and yield. The CER showed that Type II spraying was the most cost-effective in both scenarios, and particularly in Scenario 2. Whitefly population control was more costly in a smaller plantation, where the whitefly density was higher and CMD spread more rapidly.

The effects of rouging of infected plants should be investigated further, as should the effect of combining pesticide treatment with replanting. This may make spraying less efficient, thereby encouraging spread of CMD. Development of resistant strains is an alternative approach. However, it is typically costly and time consuming, and may also reduce the population of beneficial insects. Mathematical models may be used to compare such strategies to identify those that are most cost-effective in the real world.

8.   Conclusions

We developed a mathematical model to clarify the dynamics of CMD spread. The normalized forward sensitivity index of the reproduction number was used to compare the relative contribution of different factors to a CMD outbreak. We found that the two most important were the death rate of whitefly and the number of cassava plants visited. We also identified the most cost-effective policy for reducing the whitefly population using numerical simulations. CER was used to compare the economic effectiveness of three spraying strategies in farms of different sizes. The results may help stakeholders, including cassava farmers and government agencies, develop optimal policies for control of CMD outbreaks. This should increase yields, income, and profits.

Acknowledgments

We would like to thank Mr. John Winward and Dr. Sopone Wongkaew for comments and suggestions.

Conflict of interest

The authors declare that there is no conflict of interests regarding the publication of this paper. \newpage