Delays in Plant Virus Models and Their Stability

Abstract

Viruses infect humans and animals but also infect plants and cause great economic and ecological damage. In most cases, the virus is transmitted by a vector. After being bitten by an infected vector, the virus takes some time to replicate and spread in the plant. We present two models of the spread of viruses in plants based on ordinary differential equations, and then add either a delay or an exposed plant population. We study two ways of adding the delay. In the first one, a plant infected by a vector changes from susceptible to infective after a time equal to the delay. In the second one, immediately after the contact between a susceptible plant and infective vector, the plant is no longer susceptible, but it takes time equal to the delay for it to turn infective. To better explain the two ways of incorporating the delays, we first introduce them in a simple SIRS model. We analyze the models and study their stability numerically. We conclude by studying the interactions and the conservation of the total plant population that the first way of introducing the delay is better justified.

1. Introduction

Plants are fundamental members of all ecosystems and their health is of the utmost importance to man and animals. There is a wide variety of diseases that affect plants, from bacterial and fungal to those caused by viruses. There are more than one thousand different plant viruses [1]. Viruses may be transmitted from plant by plant by different mechanisms that can be either mechanical such as by agricultural tools or biological mainly by vectors such as insects but there is sometimes vertical transmission [2]. Transmission by insect vectors is the most common [3]. Among the vectors that transport these viruses are aphids, whiteflies and leafhoppers. Jeger et al. [4] consider that there are four types of insect transmission depending on where the virus spreads in the insect and on whether it replicates inside the insect. In this paper we will only be concerned with circulative, persistent transmission. That is, the vectors consume sap from an infected plant and the viruses circulate within the vector but do not replicate, do not cause the vector to become ill, and the vector will hold the virus for the rest of its life. When the infected vector contacts a healthy plant, some virus particles from the vector invade the plant [5,6,7]. Once inside the plant the virus propagates and the plant may use defense mechanisms to combat the virus. For example, it may use antiviral RNA silencing [8].

For a virus to replicate, it must invade a healthy plant cell and use the cell’s DNA or RNA to make copies of itself. Then the infected cell bursts and the copies of the virus exit and move to new cells to invade. See, for example, [9,10,11,12]. This infection process of the plant by the virus takes time. This process time is a delay time. Processes with delay times can be modeled using delay differential equations (DDEs). DDEs have been used to model time delays in virus propagation models, for example in [13,14]. For viruses in plants, DDEs have also been used, for example in [15,16,17]. The introduction of the delay may rise to oscillations in the solutions, but delay times can change the dynamics of the model in many other ways. For example, delays may change the solutions, cause discontinuities in the derivative, affect uniqueness, or change the stability. The global asymptotic stability of an SIR (susceptible, infective, recovered) epidemic model is described in [18] and a plant–virus model in [19]. Hopf bifurcations can happen in delayed models and have been studied in [20] for an SIRS (susceptible, infective, recovered, susceptible) model and in [17,21] for a plant–virus model. Fractional derivatives have also been introduced in [22,23]. Models have also been introduced for the within-plant interactions with the virus. See, for example, [24] and references therein. Despite the additional analytical and numerical difficulties, the results are more realistic from the biological and physical points of view. See for example [25,26].

Delays can be introduced into a plant virus model described by ordinary differential equations in different ways. Introducing the delays in different terms will produce different results. Not all of terms can be delayed. For example, death terms should not be delayed. Infection terms can be delayed as there is a time delay between contact and infection, but it can be done in different ways. The introduction of an additional exposed or latent population is another way of incorporating the infection time. The main objective of this paper is to study the different ways of incorporating the infection times and comparing the solutions and their stability for the different alternatives. For the delayed model the stability will be studied numerically. We will consider an SIRS model and two different and relatively simple models of vector-mediated virus transmission in plants. In Section 2 we present the three models. In Section 3, we study the different models and their stability, and compare them. The stability of the delayed models is done numerically. In Section 4 we discuss the results and we present conclusions in Section 5.

2. Materials and Methods

In this section we present three models. The first one is an SIRS epidemic model for which it is simpler to present the ideas on how to introduce the delay and on how to study the interactions by writing the system of ODEs in terms of biochemical reactions. Then we present a simple model of plant–virus transmission by a vector with a constant plant population. Finally, we present a second plant–virus propagation model with saturation and non-constant plant population.

2.1. SIRS Model

For the ease of presenting the main ideas, we first consider the well known SIRS model for epidemics [27]. It consists of three populations: susceptibles (S), infectives (I) and recovered (R)—all are functions of time t. A susceptible converts into an infective with probability $\beta$ after a contact with an infective. An infective recovers with a rate $\gamma$. A recovered loses immunity at a rate $\nu$ and turns into a susceptible. The total population is assumed to be constant, $S\left(t\right)+I\left(t\right)+R\left(t\right)=N$. Applied to plant virus transmission, it means that the transmission is by mechanical means since there are no vectors. We present it to illustrate two different ways of introducing a delay. The system of ordinary differential equations (ODEs) describing the model is:

$$\begin{array}{cc}\hfill \frac{dS}{dt}& =-\frac{\beta }{N}SI+\nu R\hfill \\ \hfill \frac{dI}{dt}& =\frac{\beta }{N}SI-\gamma I\hfill \\ \hfill \frac{dR}{dt}& =\gamma I-\nu R.\hfill \end{array}$$

It can be simplified by using the fact that the total population is constant: $R\left(t\right)=N-S\left(t\right)-I\left(t\right)$:

$$\begin{array}{cc}\hfill \frac{dS}{dt}& =-\frac{\beta }{N}SI+\nu (N-S-I)\hfill \\ \hfill \frac{dI}{dt}& =\frac{\beta }{N}SI-\gamma I.\hfill \end{array}$$

There are two common methods of introducing delays: Model 1 [16,28] given by:

$$\begin{array}{cc}\hfill \frac{dS}{dt}& =-\frac{\beta }{N}S(t-\tau )I(t-\tau )+\nu (N-S\left(t\right)-I\left(t\right))\hfill \\ \hfill \frac{dI}{dt}& =\frac{\beta }{N}S(t-\tau )I(t-\tau )-\gamma I\left(t\right),\hfill \end{array}$$

where the hypothesis is that an infective has contact with a susceptible and after a time $\tau$ the susceptible turns infective. A second model is Model 2 [29] described by:

$$\begin{array}{cc}\hfill \frac{dS}{dt}& =-\frac{\beta }{N}S\left(t\right)I(t-\tau )+\nu (N-S\left(t\right)-I\left(t\right))\hfill \\ \hfill \frac{dI}{dt}& =\frac{\beta }{N}S\left(t\right)I(t-\tau )-\gamma I\left(t\right),\hfill \end{array}$$

where the assumption is that a susceptible after contact with an infective waits a time $\tau$ before being able to infect a susceptible. So a new infective happens after a contact between the susceptible and an individual infected an amount of time ago equal to the delay.

Epidemic models can also be interpreted as biochemical reactions in which, for example, an infective reacts with a susceptible and the results are two infectives. This approach helps determine the right approach to introducing the delay. The non delayed SIRS model written in terms of biochemical reactions is:

$$\begin{array}{cc}\hfill S& \stackrel{SI\beta /N}{\to }\hfill \\ \hfill & \stackrel{SI\beta /N}{\to }I\hfill \\ \hfill & \stackrel{\nu N}{\to }S\hfill \\ \hfill S& \stackrel{\nu I}{\to }\hfill \\ \hfill S& \stackrel{\nu S}{\to }\hfill \\ \hfill I& \stackrel{\gamma I}{\to }.\hfill \end{array}$$

It is important to note that the system of reactions is not unique. This non-uniqueness has been established by [30]. In this case a simpler system is:

$$\begin{array}{cc}\hfill I+S& \stackrel{SI\beta /N}{\to }2I\hfill \\ \hfill & \stackrel{\nu N}{\to }S\hfill \\ \hfill S& \stackrel{\nu I}{\to }\hfill \\ \hfill S& \stackrel{\nu S}{\to }\hfill \\ \hfill I& \stackrel{\gamma I}{\to }.\hfill \end{array}$$

SIRS Model 1 with the delay can be represented by biochemical reactions by delaying the following reaction:

$$I+S\stackrel{SI\beta /N}{\to }2I,$$

and Model 2 by delaying

$$\stackrel{SI\beta /N}{\to }I.$$

The concern with the reaction for Model 2 is that the propensity $S\left(t\right)I(t-\tau )\beta /N$ has terms evaluated at different times. It can be fixed by adding a reaction that takes an infective and produces an infective and delays it, but for a model with no delays one would never add a reaction, saying that an infective converts to an infective.

2.2. Plant Virus with Vector Transmission, Model A

We consider two models of the plant virus transmission due to vectors. Both are simple and the object is to see what the effect of different ways of introducing the delay is. In the first one, Model A, there are two populations of plants; susceptible, $S\left(t\right)$, healthy and subject to be infected, and infective, $I\left(t\right)$, already infected. We consider that plants do not recover so there is no recovered class. There are also two populations of vectors; susceptible, $X\left(t\right)$ and infective, $Y\left(t\right)$. This model is a simplified version of the models presented, for example, in [5,16].

The assumptions of Model A are: new plants and vectors are susceptible, the total population of plants is a constant K since a farmer can plant new healthy plants to replace any dead ones, the interaction between vector and plant is of mass action type, the viruses kill plants but not vectors that do not get the disease, and neither plants or vectors recover from the disease. The parameters of the model are: $\beta$ infection rate of plants due to vectors, ${\beta }_1$ infection rate of vectors due to plants, $\mu$ natural death rate of plants, d extra death rate of plants due to disease, m natural death rate of vectors, and $\Lambda$ replenishing rate of vectors (due to birth and/or migration).

The system of ODEs for Model A is:

$$\begin{array}{cc}\hfill \frac{dS}{dt}& =\mu (K-S)+dI-\beta YS\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill \frac{dI}{dt}& =\beta YS-(d+\mu )I\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill \frac{dX}{dt}& =\Lambda -{\beta }_{1}IX-mX\hfill \end{array}$$

(3)

$$\begin{array}{cc}\hfill \frac{dY}{dt}& ={\beta }_{1}IX-mY.\hfill \end{array}$$

(4)

In virus transmission via a vector there are two delays. The first one is the time it takes the virus to spread in the plant after contagion. The second is the time it takes the virus to spread in the vector after contagion. This second one is much smaller than the first since the virus is not replicating in the vector. For simplicity we take this second delay to be zero.

As was done for the SIRS model, we will introduce the delay in two forms; the first one using the assumption that a susceptible after contact with an infective takes the delay time to become infective itself [16,28]. This will be Model A1:

$$\begin{array}{cc}\hfill \frac{dS}{dt}& =\mu (K-S(t\left)\right)+dI\left(t\right)-\beta Y(t-\tau )S(t-\tau )\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill \frac{dI}{dt}& =\beta Y(t-\tau )S(t-\tau )-(d+\mu )I\left(t\right)\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill \frac{dX}{dt}& =\Lambda -{\beta }_{1}I\left(t\right)X\left(t\right)-mX\left(t\right)\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill \frac{dY}{dt}& ={\beta }_{1}I\left(t\right)X\left(t\right)-mY\left(t\right).\hfill \end{array}$$

(8)

The second model with delay uses the assumption that after a contagion the susceptible immediately stops being susceptible but it takes the delay time to become infective. It also takes into account that the newly infected plant may die before becoming infective by introducing the term $exp(-\mu \tau )$. This is Model A2:

$$\begin{array}{cc}\hfill \frac{dS}{dt}& =\mu (K-S(t\left)\right)+dI\left(t\right)-\beta Y\left(t\right)S\left(t\right)\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill \frac{dI}{dt}& =exp(-\mu \tau )\beta Y(t-\tau )S(t-\tau )-(d+\mu )I\left(t\right)\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill \frac{dX}{dt}& =\Lambda -{\beta }_{1}I\left(t\right)X\left(t\right)-mX\left(t\right)\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill \frac{dY}{dt}& ={\beta }_{1}I\left(t\right)X\left(t\right)-mY\left(t\right).\hfill \end{array}$$

(12)

Again, writing the models in terms of biochemical reactions helps clarify the validity of each approach. For Model A the reactions are:

$$\begin{array}{cc}\hfill S& \stackrel{SY\beta }{\to }I\hfill \\ \hfill I& \stackrel{(\mu +d)I}{\to }S\hfill \\ \hfill X& \stackrel{{\beta }_{1}IX}{\to }Y\hfill \\ \hfill & \stackrel{\Lambda }{\to }X\hfill \\ \hfill X& \stackrel{mX}{\to }\hfill \\ \hfill Y& \stackrel{mY}{\to }.\hfill \end{array}$$

An alternative for accounting for the delay is to introduce instead an exposed population (E). After a contact with an infective a susceptible is converted into an exposed or latent, that cannot yet infect. The exposed turns into an infective at a rate $ϵ=1/\tau$:

$$\begin{array}{cc}\hfill \frac{dS}{dt}& =\mu (K-S(t\left)\right)+dI\left(t\right)-\beta Y\left(t\right)S\left(t\right)+\mu E\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill \frac{dE}{dt}& =\beta Y\left(t\right)S\left(t\right)-\mu E-ϵE\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill \frac{dI}{dt}& =ϵE-(d+\mu )I\left(t\right)\hfill \end{array}$$

(15)

$$\begin{array}{cc}\hfill \frac{dX}{dt}& =\Lambda -{\beta }_{1}I\left(t\right)X\left(t\right)-mX\left(t\right)\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill \frac{dY}{dt}& ={\beta }_{1}I\left(t\right)X\left(t\right)--mY\left(t\right).\hfill \end{array}$$

(17)

Epidemic models with an exposed class are very common. For plant virus propagation see, for example, [31,32]. The advantages of using models with exposed populations is that they do not require the initial conditions to be given in an interval of length equal to the delay as do DDEs. Methods for determining the stability are also better known and simpler to use in the case of ODEs. However, there is an extra equation for the exposed population, and, from the biological point of view, the fact that it takes on average only the delay time for an exposed to become infective, the change can be instantaneous.

2.3. Plant Virus with Vector Transmission, Model B

As a second plant virus transmission model we construct a model based on the one presented in [15,33], but modified to use mass action interaction for the contagion, and including healthy vectors. It considers four populations; susceptible plants, S, infective plants, I, susceptible vectors, X, and infective vectors, Y. The plants have logistic growth, so the total plant population is not constant. All new vectors are susceptible and their growth rate due to births and migration is constant. As in model A, plants do not recover and insects do not get the disease.

Without delays plant virus propagation Model B is:

$$\begin{array}{cc}\hfill \frac{dS}{dt}& =rS(1-\frac{S+I}{k})-\frac{\beta SY}{1+\alpha S+cY}\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill \frac{dI}{dt}& =\beta SY-(m+d)I\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill \frac{dX}{dt}& =\Lambda -{\beta }_{1}XI-mX\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill \frac{dY}{dt}& ={\beta }_{1}XI-mY.\hfill \end{array}$$

(21)

Here r is the growth rate of the plants, k their carrying capacity and $\beta$ the infection rate from an infective vector to a susceptible plant, $\mu$ is the natural death rate of the plants and d the additional death rate due to the virus. $\Lambda$ is the recruitment rate of susceptible vectors, ${\beta }_1$ the infection rate from infective plant to susceptible vector, and m the natural death rate of the vectors. and d the additional death rate of the.

Similarly to what was done for Model A, we will introduce the delay in two forms. The first one using the assumption that a susceptible after a contact with an infective takes the delay time to become infective itself [16,28]. This will be Model B1:

$$\begin{array}{cc}\hfill \frac{dS}{dt}& =rS\left(t\right)(1-\frac{S\left(t\right)+I\left(t\right)}{k})-\frac{\beta S(t-\tau )Y(t-\tau )}{1+\alpha S(t-\tau )+cY(t-\tau )}\hfill \end{array}$$

(22)

$$\begin{array}{cc}\hfill \frac{dI}{dt}& =\frac{\beta S(t-\tau )Y(t-\tau )}{1+\alpha S(t-\tau )+cY(t-\tau )}-(m+d)I\left(t\right)\hfill \end{array}$$

(23)

$$\begin{array}{cc}\hfill \frac{dY}{dt}& =bI\left(t\right)-mY\left(t\right).\hfill \end{array}$$

(24)

In the second form, we follow [15,33] in that after a contact between a susceptible plant and an infective insect, the plant immediately stops being susceptible, but it takes the delay time to become infective. The plant may also die before, so the survival rate is proportional to $exp(-m\tau )$, where m is the death rate of the plant and $\tau$ the delay.

$$\begin{array}{cc}\hfill \frac{dS}{dt}& =rS\left(t\right)(1-\frac{S\left(t\right)+I\left(t\right)}{k})-\frac{\beta SY}{1+\alpha S\left(t\right)+cY\left(t\right)}\hfill \end{array}$$

(25)

$$\begin{array}{cc}\hfill \frac{dI}{dt}& =exp(-m\tau )\frac{\beta S(t-\tau )Y(t-\tau )}{1+\alpha S(t-\tau )+cY(t-\tau )}-(m+d)I\left(t\right)\hfill \end{array}$$

(26)

$$\begin{array}{cc}\hfill \frac{dX}{dt}& =\Lambda -{\beta }_{1}I\left(t\right)X\left(t\right)-mX\left(t\right)\hfill \end{array}$$

(27)

$$\begin{array}{cc}\hfill \frac{dY}{dt}& ={\beta }_{1}I\left(t\right)X\left(t\right)--mY\left(t\right).\hfill \end{array}$$

(28)

As before, writing the model in terms of biochemical reactions can clarify the interactions. Model B in terms of biochemical reactions is:

$$\begin{array}{cc}\hfill & \stackrel{rS}{\to }S\hfill \\ \hfill S\stackrel{r{S}^2/k}{\to }\\ \hfill S\stackrel{rSI/k}{\to }\\ \hfill S& \stackrel{k_1}{\to }I,\mathrm{where}{k}_1=\frac{\beta SY}{1+\alpha S+cY}I\hfill \\ \hfill I& \stackrel{(\mu +d)I}{\to }S\hfill \\ \hfill X& \stackrel{{\beta }_{1}IX}{\to }Y\hfill \\ \hfill & \stackrel{\Lambda }{\to }X\hfill \\ \hfill X& \stackrel{mX}{\to }\hfill \\ \hfill Y& \stackrel{mY}{\to }.\hfill \end{array}$$

Note that in Model B1 an infected vector infects a plant at time $(t-\tau )$ and the susceptible plant turns infected at time t. In Model B2, an infected vector infects a susceptible plant that takes a time $\tau$ to become infective and the term $exp(-\tau m$) is the average number of infected susceptibles that survive in time $\tau$. Only Model B1 has the same steady state as the non-delayed model B but both models are correct from biochemical reactions point of view since the propensities have terms evaluated at the same time.

As before, an alternative, with the same pros and cons, is to added an exposed or latent class E. Model B with the exposed class is:

$$\begin{array}{cc}\hfill \frac{dS}{dt}& =rS(1-\frac{S+I}{k})-\frac{\beta SY}{1+\alpha S+cY}\hfill \end{array}$$

(29)

$$\begin{array}{cc}\hfill \frac{dE}{dt}& =\beta Y\left(t\right)S\left(t\right)-\mu E-ϵE\hfill \end{array}$$

(30)

$$\begin{array}{cc}\hfill \frac{dI}{dt}& =ϵE-(d+\mu )I\left(t\right)\hfill \end{array}$$

(31)

$$\begin{array}{cc}\hfill \frac{dX}{dt}& =\Lambda -{\beta }_{1}I\left(t\right)X\left(t\right)-mX\left(t\right)\hfill \end{array}$$

(32)

$$\begin{array}{cc}\hfill \frac{dY}{dt}& ={\beta }_{1}I\left(t\right)X\left(t\right)--mY\left(t\right).\hfill \end{array}$$

(33)

2.4. Numerical Methods

To determine the steady solutions of any of the models, the right-hand sides of the systems of equations is equaled to zero, and the algebraic system solved, either directly or by using a package such as Matlab [34], Mathematica [35] or GNU maxima [36].

The time simulations of the models given in terms of ODEs was done using the routine ode45 from Matlab. Most numerical methods for solving ODEs can be modified for DDEs. For example, Runge–Kutta and Linear Multistep methods have been used combined with the step by step method. We use the routine dde23 [37].

For systems of ODEs, the localized stability is given by the eigenvalues of the jacobian matrix of the linearized system. An alternative is to calculate the basic reproduction number $R_0$ by using, for example, the next generation matrix method [38]. Both approaches require the numerical calculation of the eigenvalues that can be done with any of the packages mentioned above. An alternative is to use xppauto [39], which not only calculates the eigenvalues for given values of the model’s parameters but also allows the search, by doing a continuation, in one parameter, of the points where the stability changes.

To determine the stability of DDEs, either by finding the eigenvalues of the jacobian of the linearized system, or by calculating the basic reproductive number, requires solving transcendental equations [25,26,40]. These equations in their majority require numerical methods, as is in our case. Consider the following system of non-linear delay differential equations:

$$\begin{array}{c}\hfill \frac{dx}{dt}=f(x\left(t\right),x(t-\tau )),\end{array}$$

(34)

where $x\in {\mathbb{R}}^n$, $\tau \ge 0$, is a constant discrete time. The linearization of (34) at $(x^{\ast },x^{\ast })$ has the form [41]:

$$\begin{array}{c}\hfill \frac{dx}{dt}=A_{0}x\left(t\right)+A_{1}x(t-\tau ),\end{array}$$

(35)

where

$$x\in {\mathbb{R}}^n,f:{\mathbb{R}}^n\times C\to {\mathbb{R}}^n,A_j=\frac{\partial f}{\partial x_j}{|}_{(x^{\ast },x^{\ast })},j=0,1.$$

Its transcendental characteristic equation is:

$$\begin{array}{c}\hfill det\left[\begin{array}{c}\lambda I-A_0-A_1{e}^{-\lambda \tau }\end{array}\right]=0.\end{array}$$

(36)

In [42] it is shown that if all of the eigenvalues have a negative real part, then solutions to the delay differential equation decay exponentially to 0, exactly as is the case for ordinary differential equations. Therefore, to determine the local stability of the DDE it is necessary to calculate the solutions of the characteristic equation. However, since it is a transcendental equation it has to be done numerically. We use the program dde-biftool [43] for Matlab can calculate the eigenvalues with the largest real parts and can also be used to find the values of the parameters for which there is a change of stability.

3. Results

In this section we present the results for Model A with a constant plant population and Model B with a non-constant plant population. Results for the SIRS model are well known and not presented.

We first present the results for Model A. By setting the right-hand side of the system of equations as equal to zero and solving for S, I, X, and Y, we get two equilibrium solutions. The first one is the disease free: $S^{\ast }=K,I^{\ast }=0,X^{\ast }=\frac{\Lambda }{m},Y^{\ast }=0$. For the non delayed ODE model this equilibrium is easy to analyze, either by calculating the eigenvalues of the jacobian matrix of the linearized system and checking the sign of their real parts, or by calculating the basic reproduction number $R_0$. This number can be calculated using the next generation matrix [38]. This equilibrium point is locally stable when

$$R_0^2=\frac{m^2(d+\mu )}{K\Lambda }>1$$

and locally unstable when $\frac{m^2(d+\mu )}{K\Lambda }<1$. The other equilibrium state, the endemic, is:

$$\begin{array}{cc}\hfill I^{\ast }& =\frac{K\beta {\beta }_1\Lambda -m^2(\mu +d)}{{\beta }_{1}m(\mu +d)+\beta {\beta }_1\Lambda }\hfill \end{array}$$

(37)

$$\begin{array}{cc}\hfill S^{\ast }& =K-I^{\ast }\hfill \end{array}$$

(38)

$$\begin{array}{cc}\hfill Y^{\ast }& =\frac{K\beta {\beta }_1\Lambda -m^2(\mu +d)}{\beta m^2+K\beta {\beta }_{1}m}\hfill \end{array}$$

(39)

$$\begin{array}{cc}\hfill X^{\ast }& =\frac{\Lambda }{{\beta }_1{I}^{\ast }+m},\hfill \end{array}$$

(40)

which is only non-negative for $R_0>1$.

Figure 1 has the simulation results for Model A. The ODE models were done using ode45 from Matlab and the DDE using dde23 also of Matlab. The parameter values used are $\tau =2, K=100, \beta =0.01, {\beta }_1=0.01, d=0.1, \mu =0.1, m=0.2, \Lambda =10$, which were taken from [5,16] and are not intended to be the values for a certain species of plant and virus. They correspond to solutions that, as time increases, go to the endemic equilibrium. From the simulations it can be seen that Model A2 does not conserve the number of plants, even though it is assumed to be constant.

Figure 2 has the simulations of Model B. It includes the ODE model, the two different delay models and the model with the exposed population. The parameter used are $\tau =2,$ $k=100, \beta =0.01, m=0.2, \mu =0.1, r=5, \Lambda =10$ which are only assumed to be typical values. They correspond to solutions that as time increase go to the endemic equilibrium. For this model the total population of plants is not conserved.

Figure 3 has the bifurcation diagrams for Model A. The parameters are the same as before, and the continuation of the branches was done by changing m which is the death rate of the vector. The objective was to see if, by increasing the dead rate of the insects, it was possible to change the stability and make the disease free equilibrium stable. As can be seen from the figures at $m=0.707$ the real part of one eigenvalue for the endemic solution goes through zero. EP− in red is the stable branch and EP+ in black is the unstable branch. The graph for the infective plant population I is not shown since $I=K-S$. As the parameter changes the stability changes and the disease free equilibrium changes stability when the real part of one of the eigenvalues changes sign. It coincides when the endemic solution has zero infectives.

Figure 4 has the same plots as Figure 3 but for the Model A with the exposed population. There is now a change of stability of the branches at $m=0.636$. EP− in red is the stable branch and EP+ in black is the unstable branch. The graph for I is not shown since it is equal to $K-I-E$.

Figure 5 has the bifurcation diagrams for Model B with no delays and no exposed. Again the continuation is in m with a bifurcation point at $m=0.894$. EP− in red is the stable branch and EP+ in black is the unstable branch. The bifurcation diagrams for Model B with exposed are not presented since they are very similar to the ones with no exposed, but the bifurcation point is now $m=0.852$.

Figure 6 is the bifurcation diagram for Model A1 with a delay. The top graphs are the eigenvalues with the largest real parts for the disease free branch (left) and the endemic equilibrium branch (right) for $m=0.2$. The bottom graphs are the diagram of the real parts of the largest eigenvalues as m, the death rate of the vectors is varied, to find the bifurcation points. The endemic branch is stable from $m=0.09$ to $m=0.68$, and the disease free branch is stable for $m\ge 0.68$.

Figure 7 is the bifurcation diagram for Model A2 with a delay. The top graphs are the eigenvalues with the largest real parts for the disease free branch (left) and the endemic equilibrium branch (right) for $m=0.2$. The bottom graphs are the diagram of the real parts of the largest eigenvalues as m, the death rate of the vectors is varied, to find the bifurcation points. The endemic equilibrium branch is stable for $m\le 0.67$, and the disease free branch for $m\ge 0.67$.

Figure 8 is the bifurcation diagram for Model B1 with a delay. The top graphs are the eigenvalues with the largest real parts for the disease free branch (left) and the endemic equilibrium branch (right) for $m=0.2$. The bottom graphs are the diagram of the real parts of the largest eigenvalues as m, the death rate of the vectors is varied, to find the bifurcation points. The endemic equilibrium branch is stable between $m=0.084$ and $m=0.90$, and the disease free branch is stable for $m\ge 0.90$.

Figure 9 is the bifurcation diagram for Model B2 with a delay. The top graphs are the eigenvalues with the largest real parts for the disease free branch (left) and the endemic equilibrium branch (right) for $m=0.2$. The bottom graphs are the diagram of the real parts of the largest eigenvalues as m, the death rate of the vectors is varied, to find the bifurcation points. The endemic equilibrium branch is stable for values of m between 0.093 and 0.3, and the endemic branch is stable for $m\ge 0.73$.

There are also some Hopf bifurcation points where a pair of complex conjugate eigenvalues have real parts that move through zero. All these Hopf bifurcation points happen when the branch is unstable. That is when another eigenvalue is already positive. For example, for Model A1, on the endemic branch there is a Hopf bifurcation point at $m=0.036$. Figure 10 has the eigenvalues with the largest real parts corresponding to this value of m.

4. Discussion

4.1. SIRS

The main objective of this paper is to study how to incorporate delays in epidemic models. While there is agreement that growth terms can be delayed due to maturation times and that dead terms should not be delayed, there are options for how to introduce delays into the infection terms. In models where the infection occurs by contact from an infective individual with a susceptible individual, like in the SIRS models, two common methods of introducing delays are based on the following assumptions: method 1 in which an infective has contact with a susceptible and after a time $\tau$ the susceptible turns infective; and method 2 in which an infective waits a time $\tau$ before being able to infect a susceptible that immediately turns infective. Usually a term like $exp(-b\tau )$, which is the fraction of new infected susceptibles that do not die before becoming infectives, is also included. Both methods conserve the total population, both methods have the same steady state as the non delayed model and both methods make biological sense. However, looking from a biochemical reaction point of view, the propensity should involve terms evaluated at the same time, so method 1 makes more sense, even though adding a reaction that converts an infective to an infective can fix this issue. In an SIRS epidemic model with demographics, births and deaths, and taking the total population constant, it is easy to see that method 2 does not conserve the total population, while method 1 does.

4.2. Plant Virus with Vector Transmission, Models A and B

For the plant-vector models, we present two models—Model A with constant plant total population, and Model B with logistic growth for the plants, which implies that the total plant population is not constant. For both Models A and B, the delay is introduced in two alternative ways. In Model A1 an infected vector infects a plant at time $(t-\tau )$ and the susceptible plant turns infected at time t. In Model A2, an infected vector infects a susceptible plant that takes time $\tau$ to become infective and the term $exp(-\tau m$) is the average number of infected susceptibles that survive in time $\tau$. However, only Model A1 conserves the total population and only Model A1 has the same steady state as the deterministic model. Finally, both models are fine from the biochemical reactions point of view since both the propensities have all the terms evaluated at the same time.

Similarly, for Model B, in Model B1 an infected vector infects a plant at time $(t-\tau )$ and the susceptible plant turns infected at time t. In Model B2, an infected vector infects a susceptible plant that takes a time $\tau$ to become infective and the term $exp(-\tau m$) is the average number of infected susceptibles that survive in time $\tau$. In Model B, the total plant population is not conserved, so the fact that Model B2 does not conserve the total plant population does not seem important, but the second method of introducing the delay has clear limitations. Finally, both models are good from the biochemical reactions point of view since both the propensities have all the terms evaluated at the same time.

A third way of introducing a delay is to instead add an exposed or latent population. In a model with a conserved plant population, introducing the exposed population also satisfies this conservation. Since there are no delay times, the corresponding biochemical reactions make sense. Additionally, since they are still described in terms of ODEs, there is no need to have the initial conditions in an interval as is the case for DDEs. The downside is that the infection can happen at any time from the contact including zero.

Looking at the local stability of the steady states for all the models, the change of stability from stable to unstable as the death rate of the vectors m increases, happens at values that are close. The models with delays also have Hopf bifurcation points, but they happen when the branch is already unstable. So they would not be observed in real life. Therefore, neglecting the Hopf bifurcation points, the solutions produced by all models have similar stability properties.

In summary, the way that delays are introduced should preserve the total population when the model has a constant population. Looking at the corresponding biochemical reactions, the delayed model should also have propensities evaluated at the same time, since a reaction involves the contact of two individuals.

5. Conclusions

Delayed models are more realistic because they take into account the time between the contact and the infection. The delay affects the processes and, mathematically, the delay affects the dynamics, including the stability. The delay can be introduced in different forms, including adding an exposed or latent population. On models that conserve the total population, the delays should be added to still satisfy this constraint. It was shown that there are ways of introducing the delay that satisfy the conservation of the total population but not all ways do. Adding delays should not change the total population. Therefore, adding the delay in the first way—after a contact between an infective and a susceptible takes time equal to the delay for the susceptible to turn into an infective—is more realistic. It is also desired that adding the delay does not change the non-delayed equilibria; again, some ways of introducing the delay do and some do not. Finally, writing the model in terms of biochemical reactions helps to decide the correct form. Future work will include studying the case of multiple delays, and analyzing the local and global stabilities.