ARTICLE IN PRESS
JID: MBS
[m5G;April 23, 2015;15:3]
Mathematical Biosciences 000 (2015) 1–12
Contents lists available at ScienceDirect
Mathematical Biosciences journal homepage: www.elsevier.com/locate/mbs
Free-virus and cell-to-cell transmission in models of equine infectious anemia virus infection Linda J.S. Allen a,∗, Elissa J. Schwartz b a b
Department of Mathematics and Statistics, Texas Tech University, Lubbock, TX 79409-1042, United States School of Biological Sciences and Department of Mathematics, Washington State University, Pullman, WA 99164-3113, United States
a r t i c l e
i n f o
Article history: Available online xxx Keywords: Branching process Markov chain Within-host model Virus dynamics Equine infectious anemia virus
a b s t r a c t Equine infectious anemia virus (EIAV) is a lentivirus in the retrovirus family that infects horses and ponies. Two strains, referred to as the sensitive strain and the resistant strain, have been isolated from an experimentallyinfected pony. The sensitive strain is vulnerable to neutralization by antibodies whereas the resistant strain is neutralization-insensitive. The sensitive strain mutates to the resistant strain. EIAV may infect healthy target cells via free virus or alternatively, directly from an infected target cell through cell-to-cell transfer. The proportion of transmission from free-virus or from cell-to-cell transmission is unknown. A system of ordinary differential equations (ODEs) is formulated for the virus-cell dynamics of EIAV. In addition, a Markov chain model and a branching process approximation near the infection-free equilibrium (IFE) are formulated. The basic reproduction number R0 is defined as the maximum of two reproduction numbers, R0s and R0r , one for the sensitive strain and one for the resistant strain. The IFE is shown to be globally asymptotically stable for the ODE model in a special case when the basic reproduction number is less than one. In addition, two endemic equilibria exist, a coexistence equilibrium and a resistant strain equilibrium. It is shown that if R0 > 1, the infection persists with at least one of the two strains. However, for small infectious doses, the sensitive strain and the resistant strain may not persist in the Markov chain model. Parameter values applicable to EIAV are used to illustrate the dynamics of the ODE and the Markov chain models. The examples highlight the importance of the proportion of cell-to-cell versus free-virus transmission that either leads to infection clearance or to infection persistence with either coexistence of both strains or to dominance by the resistant strain. © 2015 Elsevier Inc. All rights reserved.
1. Introduction The lentiviruses human immunodeficiency virus type 1 (HIV-1) and equine infectious anemia virus (EIAV) establish persistent, lifelong infections in their hosts. EIAV infects cells of the monocyte/macrophage lineage in its host, which are equids including horses and ponies. While HIV-1 infection results in CD4 T cell depletion and progression to AIDS, immune responses including antibodies and cytotoxic T lymphocytes are able to control EIAV infection. EIAV infected animals live, for the most part, without disease, and viral loads are maintained at low levels. Thus the study of EIAV infection may provide critical clues to the pathogenesis and control of lentiviral infection in other species, notably humans. Since mathematical modeling has been a useful tool to advance our understanding of virus dynamics, we investigate a model of
∗
Corresponding author. Tel.: +1 8067422580. E-mail address: [email protected] (L.J.S. Allen).
viral infection at the within-host level that includes two modes of viral transmission, free-virus and direct cell-to-cell transmission. Our model is applicable to EIAV, as it displays both modes of transmission [34]. Cell-to-cell transmission has been observed in enveloped viruses, viruses surrounded by a lipid envelope such as herpesviruses, paramyxoviruses, pox viruses, rhabdoviruses and retroviruses [23]. In cell-to-cell transmission, viruses are transferred from infected cells to uninfected cells via cell adhesions called virological synapses. One advantage of cell-to-cell transmission over free-virus transmission is that the virus is not exposed to the extracellular environment. This allows the virus to evade immune responses, such as the neutralizing antibody response [6,7,23]. Furthermore, cell-tocell transmission permits viral spread despite antiretroviral therapy [28]. In EIAV, a neutralization-resistant variant emerged over time in an experimentally-infected pony [29]. This resistant variant shows lower infectiousness but more efficient transmission than the sensitive strain [34]. In fact, the efficiency of EIAV cell-to-cell spread for this strain was found to be 1000-fold greater than free virus transmission in vitro [34]. To better understand the role of transmission
http://dx.doi.org/10.1016/j.mbs.2015.04.001 0025-5564/© 2015 Elsevier Inc. All rights reserved.
Please cite this article as: L.J.S. Allen, E.J. Schwartz, Free-virus and cell-to-cell transmission in models of equine infectious anemia virus infection, Mathematical Biosciences (2015), http://dx.doi.org/10.1016/j.mbs.2015.04.001
JID: MBS 2
ARTICLE IN PRESS
[m5G;April 23, 2015;15:3]
L.J.S. Allen, E.J. Schwartz / Mathematical Biosciences 000 (2015) 1–12
modes among neutralization sensitive and resistant strains, here we incorporate both modes of transmission in our EIAV model, which includes two viral strains. Others have modeled cell-to-cell viral transmission, but they have concentrated on models with a single viral strain. Komarova et al. [15,16] and Sigal et al. [28] each looked at sensitivity to antiretroviral therapy (ART), using deterministic and stochastic infection models, respectively. Komarova et al. considered cells multiply infected with a single HIV strain and reported how the relative probability of individual virions to infect a cell during free virus and cell-to-cell transmission affects ART susceptibility [16]. Komarova et al. also fit a mathematical model to previously published and newly generated in vitro data and determined that free virus and cell-to-cell transmission contribute approximately equally to the growth of the virus population [15]. Pourbashash et al. [22] presented a global analysis for a single strain model and for a special case of a multiple strain competition model. Both models included cell-to-cell and free-virus transmission. Analysis of their multi-strain model showed that only one strain dominated, the one with the largest strain reproduction number [22]. The study by Pourbashash et al. [22] is one of the first to include multiple viral strains and two modes of transmission. Lentiviruses such as EIAV and HIV are characterized by coinfection by multiple strains. Therefore it is critical to consider multiple strain competition, and emergence of resistance, particularly in studies of virus escape from control. The current work expands upon previous models that include multiple viral strains (see e.g. [3,9,22]) and models with two modes of transmission [15,16,22,28] by investigating how cell-to-cell transmission plays a role in multiple strain coexistence versus dominance by a resistant strain. The current work also advances the study of EIAV viral dynamics by considering the effect of cell-to-cell transmission. Previous work in EIAV modeling has investigated a single strain model with immune responses [24], the role of antibody neutralization on competition and escape with two viral strains [8], and the correlates of antibody protection and escape in EIAV infection [25], but these studies did not address how the mode of infection (i.e., cell-to-cell versus free virus transmission) affects virus escape and the development of resistance, as in the current study. Our goals in this investigation are to provide a preliminary mathematical analysis of within-host viral infection for EIAV using both modes of transmission to determine conditions that allow one or both strains to persist and to provide an analysis for the probability of infection clearance. In Section 2, the new within-host model with two viral strains and with free-virus and cell-to-cell transmission is formulated. The model is a system of ordinary differential equations (ODEs), an extension of the model recently studied by Pourbashash et al. [22] to account for the mutation from the sensitive strain to the resistant strain and for the proportion of cell-to-cell versus freevirus transmission. Five variables in the model include healthy target cells, sensitive and resistant free viruses, and cells infected by the sensitive or the resistant viral strain. In Section 3, the basic reproduction number for the model is computed and conditions for existence of three equilibria are verified, the infection-free equilibrium (IFE) and two endemic equilibria. Global stability of the IFE is verified in a special case when the basic reproduction number is less than one and uniform persistence of the infection when the basic reproduction number exceeds one. In Section 4, a continuous-time Markov chain model is formulated based on the ODE model. The probability that an infection is cleared is estimated from multitype branching process theory. In Section 5, parameter values representing two variants of EIAV are used to illustrate the persistence of one or both viral strains in the ODE model and infection clearance in the Markov chain model. 2. ODE model The within-host ODE model includes two viral variants transmitted via free viral entry into a cell or direct cell-to-cell transmission.
The five variables in the ODE model are denoted T = target cells, Is = target cells infected with the sensitive strain, Ir = target cells infected with the resistant strain, Vs = free virions of the sensitive strain, and Vr = free virions of the resistant strain. The terms β fi VT and β ci IT represent transmission rates of free-virus and cell-to-cell, respectively, from either the sensitive i = s or the resistant viral strain, i = r. The parameter pi represents the proportion of transmission from the free virus and 1 − pi the proportion that is cell-to-cell for either the sensitive (i = s) or the resistant (i = r) strain, 0 pi 1. The particular cases pi = 0 or pi = 1 represent models with exclusive cell-to-cell transmission or with exclusive free-virus transmission, respectively. Mutation occurs within a cell infected with the sensitive strain. The proportion of cells infected with the sensitive virus that mutates to the resistant strain is μ, 0 < μ < 1. We assume there is no back mutation from the resistant strain to the sensitive strain. With these assumptions, the ODE model takes the following form:
T˙ = − dT T − [ps βf s Vs + pr βf r Vr + (1 − ps )βcs Is + (1 − pr )βcr Ir ]T I˙s = [ps (1 − μ)βf s Vs + (1 − ps )(1 − μ)βcs Is ]T − ds Is I˙r = [pr βf r Vr + (1 − pr )βcr Ir + ps μβf s Vs + (1 − ps )μβcs Is ]T − dr Ir V˙ s = ps Ns ds Is − cs Vs − ps βf s Vs T V˙ r = pr Nr dr Ir − cr Vr − pr βf r Vr T.
(1)
Healthy target cells are produced at a constant rate . Parameters dT , ds and dr are the death rates of healthy or infected target cells, either sensitive or resistant. It is assumed that the death rates for infected cells are greater than or equal to that of the healthy cells, that is,
dT ≤ max{dr , ds }. Parameters cr and cs are the clearance rates of the resistant and sensitive viral strains. The parameter Ni , i = s, r, is referred to as the burst number, the number of free viruses produced by either a sensitive or a resistant infected cell. The probability pi , i = s, r, is multiplied by the burst number, representing the proportion of virions released for free viral transmission. The terms pi β fi Vi T in the differential equations for Vi account for loss of free virus via entry into target cells. Often these terms are neglected in mathematical models since they are small compared to the viral clearance rates. All parameters in model (1) are positive with the exception of pi (0 pi 1). Initial values for model (1) are nonnegative, T(0) = T0 > 0, Ii (0) = Ii0 0, and Vi (0) = Vi0 0, i = r, s with Vs0 > 0 and T0 + Is0 + Ir0 /dT . Also, if ps = 0, then Is (0) > 0. 3. Analysis of the ODE model 3.1. Bounded, nonnegative solutions Solutions of model (1) are unique, nonnegative and bounded. The healthy target cell population T¯0 = /dT is an upper bound for the total cell population and in addition, an upper bound for the viral population depends on T¯0 . Specifically, the set defined in Theorem 3.1 is forward invariant for model (1). Solutions of model (1) exist for t [0, ). Verification of these results is relegated to the Appendix. Theorem 3.1. Solutions of model (1) are unique, nonnegative and bounded. The set ⊂ R5+ is forward invariant for model (1), where
= {(T, Is , Ir , Vs , Vr ) ∈ R5+ |0 ≤ T + Ir + Is ≤ T¯0 , 0 ≤ Vs + Vr ≤ V¯ }, (2) T¯0 = /dT , V¯ = max{Vs (0) + Vr (0), K T¯0 /c}, K = ps Ns ds + pr Nr dr and c = min {cr , cs }.
Please cite this article as: L.J.S. Allen, E.J. Schwartz, Free-virus and cell-to-cell transmission in models of equine infectious anemia virus infection, Mathematical Biosciences (2015), http://dx.doi.org/10.1016/j.mbs.2015.04.001
ARTICLE IN PRESS
JID: MBS
[m5G;April 23, 2015;15:3]
L.J.S. Allen, E.J. Schwartz / Mathematical Biosciences 000 (2015) 1–12
values of T¯r , I¯r and V¯ r are positive. Let the right side of the differential equation for T˙ + I˙r = 0 and let Is = 0 = Vs . These assumptions yield an expression for cells infected with the resistant strain which depends on T:
3.2. Infection-free equilibrium The infection-free equilibrium (IFE) for model (1) is
E0 =
dT
, 0, 0, 0, 0 = (T¯0 , 0, 0, 0, 0).
(3)
Ir =
Linearization of model (1) about the IFE yields the following Jacobian matrix, J, separated into a matrix representing sources of infection F and another matrix for other transitions V, where J = F − V with
⎤
⎡
(1 − ps )(1 − μ)βcs T¯0 0 ps (1 − μ)βf s T¯0 0 ⎥ ⎢ ⎢ (1 − ps )μβcs T¯0 ( 1 − pr )βcr T¯0 ps μβf s T¯0 pr βf r T¯0 ⎥ ⎥ ⎢
F =⎢
⎣
0
0
0
0
0
0
0
0
⎥ ⎦
(4) and
⎡ ⎢ ⎢
V=⎢ ⎢
ds
0
0
0
0
dr
0
0
0
cs + ps βf s T¯0
0
−pr Nr dr
0
cr + pr βf r T¯0
⎣ −ps Ns ds 0
(5)
The spectral radius of the next generation matrix F V −1 is the basic reproduction number, R0 = ρ(F V −1 ) [33]. Specifically, R0 = max{R0s , R0r },
(6)
where R0s = (1 − ps )
R0r = (1 − pr )
(1 − μ)βcs T¯0 ds
βcr T¯0 dr
+ ps
(1 − μ)βf s T¯0 ps Ns cs + ps βf s T¯0
βf r T¯0 pr Nr + pr . cr + pr βf r T¯0
− dT T dr
.
(9)
The differential equation for Vr yields an expression for the equilibrium Vr in terms of Ir and T. Let V˙ r = 0, then
Vr =
pr Nr dr I¯r . cr + pr βf r T
(10)
Thus, it is only necessary to show there exists a unique positive value for T. Substituting the preceding expressions for Ir and Vr into the right side of T˙ and setting it equal to zero, leads to an identity for T, G(T) = 1, where
G(T ) = (1 − pr )
⎤ ⎥ ⎥ ⎥. ⎥ ⎦
3
βcr T dr
+ pr
βf r Tpr Nr . cr + pr βf r T
(11)
Evaluating G at the IFE value yields G(T¯0 ) = R0r . Since G(0) = 0 and G is a strictly increasing function of T on [0, ), it follows that if R0r > 1, then there exists a unique positive solution T¯r such that G(T¯r ) = 1. In addition, T¯r < T¯0 . This fact and the expressions for Ir in (9) and for Vr in (10) show that a unique endemic equilibrium Er exists (with Is and Vs equal to zero) provided R0r > 1. Second, conditions will be derived for existence and uniqueness of a feasible equilibrium Ers = (Tˆrs , Iˆs , Iˆr , Vˆs , Vˆr ). Identities for Vs and Vr can be obtained in terms of T, Is and Ir :
ps Ns ds Is pr Nr dr Ir and Vr = . cs + ps βf s T cr + pr βf r T
(7)
Vs =
(8)
As in the preceding analysis, an identity for T can be derived from the differential equation for Is set equal to zero and the expression for Vs . This yields H(T) = 1, where
In this calculation, the infected cells are the sources of infections in matrix F and not the viruses, which leads to the commonly-defined expression for the basic reproduction number given in (6), (7) and (8). If the viruses are included in matrix F as sources of infections through the burst number, then an alternate form for the basic reproduction number is obtained which is the square root of the expressions in (7) and (8). This alternate form is equivalent to (6) in the sense that it is also a threshold with respect to one. It follows from [33] that if R0 < 1, then the IFE is locally asymptotically stable and unstable if R0 > 1. The reproduction numbers for the sensitive and resistant strains in (7) and (8) consist of a sum of two terms representing the two modes of transmission. The first term in the sum represents transmission due to direct cell-to-cell transmission and the second term due to freevirus transmission. The factor 1 − μ in R0s accounts for the fraction of infected cells with virus that does not mutate to the resistant strain. The mutated resistant cells are not included in the expression for the basic reproduction number, due to the fact that the cells infected with the sensitive strain change their reproductive strategy after they become infected. Mutation is important for establishment of endemic equilibria containing the resistant strain. 3.3. Endemic equilibria For model (1), there exist two endemic equilibria, an equilibrium with only the resistant strain, denoted Er , and a second equilibrium with both the sensitive and resistant strains, denoted Ers . Given that μ > 0 with no back mutation from resistant to sensitive strain, there does not exist an endemic equilibrium with only the sensitive strain. First, conditions will be derived for existence and uniqueness of a feasible equilibrium Er = (T¯r , 0, I¯r , 0, V¯ r ), that is, conditions so that the
H(T ) = (1 − ps )
(1 − μ)βcs T ds
+ ps
(1 − μ)βf s ps Ns T . cs + ps βf s T
(12)
(13)
Note that H(T¯0 ) = R0s and H is nonnegative and a strictly increasing function of T on [0, ). Therefore, it follows that if R0s > 1, then a unique positive solution for Tˆrs exists satisfying H(Tˆrs ) = 1. Also, Tˆrs < T¯0 . The value of Tˆrs depends on the mutation proportion μ. Let f(T) = − dT T. Then f (Tˆrs ) > 0. Applying the differential equation for T, the identity H(Tˆrs ) = 1, the relation f(T) − dr Ir − ds Is = 0, and either the expression for Vr or for Vs , then the following two identities must be satisfied by Ir and Is at the endemic equilibrium Ers :
f (Tˆrs )μ 1 − G(Tˆrs ) Ir = 1−μ dr (1 − μ)
(14)
1 f (Tˆrs )[1 − G(Tˆrs )] − G(Tˆrs ) Is = , 1−μ ds
(15)
where the function G is defined in (11). To ensure there exist positive solutions Iˆr and Iˆs to (14) and (15), respectively, requires that G(Tˆrs ) < 1. Existence and uniqueness of a feasible Ers requires two conditions: R0s > 1 and G(Tˆrs ) < 1. Either one or two endemic equilibria exist (in R5+ ) if R0r > 1 and R0s > 1. The relation between the functions G(T) and H(T) at Tˆrs determines whether there exist one or two endemic equilibria. If G(Tˆrs ) < H(Tˆrs ) = 1, there exist two endemic equilibria, Ers and Er , and the healthy cells satisfy
0 < Tˆrs < T¯r < T¯0 . The Jacobian matrix of model (1) evaluated at Er shows that Er is unstable when H(T¯r ) > 1. Since H is an increasing function of T,
Please cite this article as: L.J.S. Allen, E.J. Schwartz, Free-virus and cell-to-cell transmission in models of equine infectious anemia virus infection, Mathematical Biosciences (2015), http://dx.doi.org/10.1016/j.mbs.2015.04.001
ARTICLE IN PRESS
JID: MBS 4
[m5G;April 23, 2015;15:3]
L.J.S. Allen, E.J. Schwartz / Mathematical Biosciences 000 (2015) 1–12
H(Tˆrs ) = 1 < H(T¯r ) implies if both Ers and Er exist, then Er is unstable. On the other hand, if H(Tˆrs ) = 1 ≤ G(Tˆrs ), only the resistant strain Er exists. Note that the parameters in H(T) depend on the sensitive strain and the mutation parameter μ, whereas the parameters in G(T) depend on the resistant strain. If, for example, in the sensitive strain, the mutation is large or transmission rates are smaller than in the resistant strain so that H(T) < G(T), then the resistant strain Er will dominate; the coexistence equilibrium Ers will not exist. Only one endemic equilibrium exists if either R0s < 1 < R0r or R0r < 1 < R0s . If R0s < 1 < R0r , then the endemic equilibrium Er ∈ / R5+ . If R0r < 1 < R0s , then Ers ∈ R5+ and Er ∈ / R5+ . In either R5+ and Ers ∈ case, the healthy cell population at Er or Ers is bounded above by T¯0 . The preceding results concerning E0 , Er and Ers are summarized in the following theorem.
(ii) If R0r > 1, then the resistant strain in model (1) is uniformly persistent. That is, there exists a positive constant κ > 0 such that lim inft→∞ X (t) > κ for X equal to Ir or Vr . In the case that R0 > 1, there exist either one or two endemic equilibria in R5+ . Numerical results indicate that if Ers ∈ R5+ , it is stable but if Ers ∈ / R5+ and Er ∈ R5+ , then Er is stable. When the endemic equilibria are feasible, the numerical results and the results of Theorems 3.2–3.4 suggest that either Er or Ers is globally asymptotically stable. If stochasticity in the birth, death and transmission process is included, then even if R0 > 1, the sensitive and the resistant strains may not persist. These results are shown for the Markov chain model in the next section. 4. Markov chain model
Theorem 3.2. There exists at most three feasible equilibria for model (1) which depend on the values of R0 , R0s and R0r , defined in (6), (7) and (8), respectively. The IFE E0 = (T0 , 0, 0, 0, 0) always exists and E0 is the only feasible equilibrium when R0 ≤ 1. If R0 < 1, then E0 is locally asymptotically stable and if R0 > 1, then E0 is unstable and in addition, either one or two endemic equilibria exist in R5+ . An endemic equilibrium Er = (T¯r , 0, I¯r , 0, V¯ r ) exists in R5+ if R0r > 1, where the value of T¯r is defined implicitly by G(T¯r ) = 1 in (11) and I¯r and V¯ r are defined in (9) and (10), respectively. An endemic equilibrium Ers = (Tˆrs , Iˆs , Iˆr , Vˆs , Vˆr ) exists in R5+ if R0s > 1, where Tˆrs is defined implicitly by H(Tˆrs ) = 1 in (13) and G(Tˆrs ) < 1, and the remaining endemic values Iˆr , Iˆs , Vˆr , and Vˆs are defined by the expressions in (12), (14) and (15), respectively. If both Ers and Er exist in R5+ , then Er is unstable. Although we assumed μ > 0 if mutation μ = 0, then the coexistence equilibrium Ers simplifies to a single strain equilibrium with only the sensitive strain, Es . It is interesting to note that Pourbashash et al. [22] results on multi-strain competition can be applied to model (1) if μ = 0, 0 < pi < 1, i = r, s and the terms pi β fi Vi T are removed from the differential equations for V˙ i , i = r, s. In their model, if R0s > max{1, R0r }, then Es is globally asymptotically stable whereas if R0r > max{1, R0s }, then Er is globally asymptotically stable (Theorem 5.2, p. 3352 [22]). 3.4. Global stability of the IFE and uniform persistence
Introduction of virus or infected cells into a healthy individual may not result in infection, especially if the infectious dose (virus and infected cells) is sufficiently small. To study whether an infection can be cleared, we formulate a Markov chain model, where the variables are discrete-valued and random. In addition, a branching process approximation is used to estimate the probability an infection is cleared. The continuous-time Markov chain model is based on the parameters in the ODE model. Fourteen distinct events and their corresponding infinitesimal transition rates, ri t + o(t), are defined in Table 1. The same notation is used as in the ODE model but the random variables are discrete, T, Is , Ir , Vs , Vr {0, 1, 2, . . . }, and time is continuous, t [0, ). The time between events is exponentially distributed with parameter 14 i=1 ri . The changes in the random variables are ±1, as in a birth and death process. For example, sensitive virus entry into target cells with no mutation, event 3, results in the transition (T, Is , Vs ) → (T − 1, Is + 1, Vs − 1) with probability ps (1 − μ)β fs Vs Tt + o(t). 4.1. Probability of infection clearance A Galton–Watson multitype branching process approximation of the Markov chain model is defined at the IFE. We estimate the probability that the infection is cleared when either the resistant strain or the sensitive strain is introduced. Offspring probability generating functions (pgfs) are defined for the infected cells and the viruses. The
It is shown that if
ˆ 0s = (1 − ps ) R
(1 − μ)βcs T¯0 ds
+ ps
(1 − μ)βf s T¯0 ps Ns cs
≤1
and
ˆ 0r = (1 − pr ) R
βcr T¯0 dr
+ pr
βf r T¯0 pr Nr cr
≤ 1,
Table 1 Discrete events in a Markov chain model for free-virus and cell-to-cell transmission with either sensitive and resistant viral strains based on the ODE model (1). FVTi = free-virus transmission of either sensitive strain or resistant strain, CCTi = cell-to-cell transmission of either sensitive strain or resistant strain, i = s, r. Event i
Description
Transition
Rate ri
T→T+1 T→T−1 (T, Is , Vs ) → (T − 1, Is + 1, Vs − 1) (T, Is ) → (T − 1, Is + 1) (T, Ir , Vs ) → (T − 1, Ir + 1, Vs − 1) (T, Ir ) → (T − 1, Ir + 1) (T, Ir , Vr ) → (T − 1, Ir + 1, Vr − 1) (T, Ir ) → (T − 1, Ir + 1) Is → Is − 1 Ir → Ir − 1 Vs → Vs + 1 Vr → Vr + 1 Vs → Vs − 1 Vr → Vr − 1
then the IFE is globally asymptotically stable. In this case R0 < 1. A Lyapunov function is constructed [27]. The proof of this result is ˆ 0i correspond to relegated to the Appendix. Note that the values of R model (1) assuming the terms pi β fi Vi T, i = r, s are removed.
1 2 3
Birth T Death T FVTs
4
CCTs
ˆ 0s } ≤ 1 and 0 < pi < 1, i = r, s, then ˆ 0r , R ˆ 0 = max{R Theorem 3.3. If R the IFE E0 of model (1) is globally asymptotically stable.
5
7
FVTs & mutate CCTs & mutate FVTr
8
CCTr
In addition, if R0 > 1, it is shown that the infection persists. Mathematically, it is shown that either the sensitive strain or the resistant strain is uniformly persistent. The theorem is stated below and the proof is given in the Appendix. Theorem 3.4. Assume R0 > 1 and 0 < pi < 1, i = r, s. (i) If R0s > 1 ≥ R0r , then the sensitive strain in model (1) is uniformly persistent. That is, there exists a positive constant κ > 0 such that lim inft→∞ X (t) > κ for X equal to Is or Vs .
6
9 10 11 12 13 14
Death Is Death Ir Birth Vs Birth Vr Clearance Vs Clearance Vr
dT T ps (1 − μ)β fs Vs T (1 − ps )(1 − μ)β cs Is T ps μβ fs Vs T (1 − ps )μβ cs Is T pr β fr Vr T (1 − pr )β cr Ir T ds Is dr Ir ps Ns ds Is pr Nr dr Ir cs Vs cr Vr
Please cite this article as: L.J.S. Allen, E.J. Schwartz, Free-virus and cell-to-cell transmission in models of equine infectious anemia virus infection, Mathematical Biosciences (2015), http://dx.doi.org/10.1016/j.mbs.2015.04.001
ARTICLE IN PRESS
JID: MBS
[m5G;April 23, 2015;15:3]
L.J.S. Allen, E.J. Schwartz / Mathematical Biosciences 000 (2015) 1–12
general form of an offspring pgf for a random variable Ri from a total of n random variables is
fi (u1 , . . . , un ) =
∞
∞
···
j1 =0
Pi (j1 , . . . , jn )u11 · · · unn , j
j
jn =0
where Pi (j1 , . . . jn ) is the probability that one type i gives “birth” to jk type k, k = 1, . . . , n. The values of ui [0, 1]. For example, fi /uj evaluated at u1 = 1, . . . , un = 1 is the expected number of offspring of type j from one of type i. For our model, Pi is the probability new infected cells or new virions (offspring) are generated from introduction of either one infected cell or one virion. It follows from branching process theory that the fixed points of the pgfs on [0, 1]n approximate the probability of extinction [1,2,10,12,14]. The minimal fixed point of fi (q1 , . . . , qn ) = qi , where qi [0, 1] determines the probability of infection extinction (clearance). In the branching process approximation, sample paths either approach zero or infinity. An important assumption is that zero is an absorbing state for the process. The probability that sample paths for the infected states hit zero yields an estimate for the probability of infection clearance. Given n infected states, Ri , i = 1, . . . , n, with associated pgfs fi , minimal fixed point (q1 , . . . , qn ) and relatively small initial values Ri (0) = ri0, then the branching process estimate for the probability of extinction (Pext ) of the Ri is r20 rn0 Pext = qr10 1 q2 · · · qn .
It is important to note that Pext is an asymptotic estimate based on the linear branching approximation. The time until infection clearance varies and depends on the initial concentrations and the magnitude of R0 . Let T = T¯0 . Then the events 3–14, listed in Table 1, have linear transition rates ri with respect to the random variables Ii and Vi , i = s, r. These transition rates are used to define offspring pgfs [1,2,10]. Let fi (u1 , u2 , u3 , u4 ) for ui [0, 1], i = 1, 2, 3, 4, denote the pgf for each of the four random variables Is , Ir , Vs , and Vr , respectively. Assume (Is (0), Ir (0), Vs (0), Vr (0)) = (δ i1 , δ i2 , δ i3 , δ i4 ), where
δik =
0, i = k 1, i = k
is the Kronecker’s delta symbol (introduction of either one infected cell or one virion). The pgfs fi fi (u1 , u2 , u3 , u4 ) for i = 1, 2, 3, 4 are defined as follows:
f1 =
ds +(1−ps )(1−μ)βcs T¯0 u21 + (1−ps )μβcs T¯0 u1 u2 + ps Ns ds u1 u3 ds + (1 − ps )βcs T¯0 + ps Ns ds
f2 =
dr + (1 − pr )βcr T¯0 u22 + pr Nr dr u2 u4 dr + (1 − pr )βcr T¯0 + pr Nr dr
f3 =
cs + ps (1 − μ)βf s T¯0 u1 + ps μβf s T¯0 u2 cs + ps βf s T¯0
f4 =
(16)
Note fi (1, 1, 1, 1) = 1. For example,
cr cr + pr βf r T¯0
, P4 (0, 1, 0, 0) =
(i) If R0 < 1, then Pext = 1. vs0 (ii) If R0r < 1 < R0s , then Pext = qis0 11 q31 , where q11 and q31 are defined as cs + ps μβf s T¯0 ps (1 − μ)βf s T¯0 1 1 + q11 = and q31 = R0s R0s cs + pss βf s T¯0 cs + ps βf s T¯0 (Eq. (19) in the Appendix). ir0 vs0 vr0 (iii) If min{R0r , R0s } > 1, then Pext = qis0 12 q22 q32 q42 , where qi2 , i = 2, 3, 4 are defined as
q22 =
1 R0r
, q42 =
pr βf r T¯0
cr + pr βf r T¯0
1 R0r
+
cr cr + pr βf r T¯0
and q32 = f3 (u1 , q22 , q32 , q42 ) =
cs + ps (1 − μ)βf s T¯0 u1 + ps μβf s T¯0 q22 cs + ps βf s T¯0
(Eqs. (18) and (20) in the Appendix). The value of q12 is the solution u1 (0, 1) of the quadratic equation f1 (u1 , q22 , q32 , q42 ) = u1 .
5. Numerical examples Parameter values for EIAV for model (1) are summarized in Table 2. For the parameter values in Table 2 the value of R0 ranges from a minimum of 21.6 (ps = pr = 1) to a maximum of 457 (pr = ps = 0). The basic reproduction number R0 > 1 for all parameter values. Large values of R0 correspond to greater cell-to-cell transmission when the values of pr and ps are small. Note that dT = ds = dr and cs = cr . The terms representing viral loss due to entry into target cells have a negligible effect on the value of the basic reproduction number due to the large viral clearance rates, cr and cs . FVTi = free-virus transmission rate of either sensitive strain or resistant strain, CCTi = cell-to-cell transmission rate of either sensitive strain or resistant strain, i = s, r, D = day. The reproduction rate of healthy target cells, , was calculated using = dT T(0), assuming that dT/dt = 0 before infection. The death rate of healthy target cells, dT , was calculated using the estimate for a macrophage lifespan of 21 days [32]. The free virus transmission rate for sensitive virus (β fs ), death rates of infected cells (ds and dr ), burst numbers (Ns and Nr ), and virus clearance rates (cs and cr ) were based upon the estimated values for EIAV infection in SCID horses [26]. Values here have been expressed using virions instead of vRNA copies (given that one virion contains two vRNA copies). Transmission rate β fr was determined given that resistant free virus EIAVPND-5 is 4.5-fold less infectious than sensitive free virus EIAVPND-1 [34]. Transmission rate β cs was determined given that cell-to-cell transmission Table 2 Parameter values for EIAV corresponding to model (1).
cr + pr βf r T¯0 u2 . cr + pr βf r T¯0
P4 (0, 0, 0, 0) =
5
pr βf r T¯0 cr + pr βf r T¯0
with all other probabilities P4 (j1 , j2 , j3 , j4 ) = 0. Details of the calculation for the fixed points of (16) are provided in the Appendix and the results are summarized in the following theorem. Theorem 4.1. Assume T (0) = T¯0 with relatively small initial values for the infected states, (Is (0), Ir (0), Vs (0), Vr (0)) = (is0, ir0, vs0, vr0). Then for the Markov chain model, there are three different cases for the estimate of the probability of infection clearance.
Parameter
Units
Definition
Value
Reference
ds dr Ns
cells/(ml D) D−1 ml/(virions D) ml/(virions D) ml/(cells D) ml/(cells D) (base cycle)−1 D−1 D−1 virions/cell
2019 1/21 6.50 × 10−7 1.44 × 10−7 5.13 × 10−4 5.13 × 10−4 3 × 10−5 1/21 1/21 5302.5
Calculated [32] [26] [34] [34] [34] [19] [26] [26] [26]
Nr
virions/cell
5302.5
[26]
cs
D−1
6.73
[26]
cr
−1
Reprod rate T Death rate T FVTs FVTr CCTs CCTr Mutation rate Death rate Is Death rate Ir Burst number Sens virus Burst number Res virus Clearance rate Sens virus Clearance rate Res virus
6.73
[26]
dT
β fs β fr β cs β cr μ
D
Please cite this article as: L.J.S. Allen, E.J. Schwartz, Free-virus and cell-to-cell transmission in models of equine infectious anemia virus infection, Mathematical Biosciences (2015), http://dx.doi.org/10.1016/j.mbs.2015.04.001
ARTICLE IN PRESS
JID: MBS 6
[m5G;April 23, 2015;15:3]
L.J.S. Allen, E.J. Schwartz / Mathematical Biosciences 000 (2015) 1–12 Table 3 Initial values for EIAV corresponding to model (1). Initial value
Units
Value
T(0) Is (0) Ir (0) Vs (0) Vr (0)
cells/ml cells/ml cells/ml virions/ml virions/ml
42,390 0 0 233 Varies
Existence and stability of the endemic equilibria for the ODE model (1) depend on the relation between G(T) and H(T) (Theorem 3.2). If G(T) < H(T), then the coexistence equilibrium exists and the resistant strain equilibrium is unstable, and if H(T) < G(T), then the coexistence equilibrium does not exist; only the resistant strain equilibrium exists. In Figs. 3 and 4 are graphs of the contours of H and G as a function of T and ps . If ps = pr , then G(T, ps ) < H(T, ps ) (Figs. 1 and 3), but if pr = 0.8ps , then H(T, ps ) < G(T, ps ) (Figs. 2 and 4). If the ratio of cell-to-cell to free-virus transmission is greater in the resistant strain as compared to the sensitive strain such that H(T, ps ) < G(T, ps ), e.g.,
Reference [13,20]
[5,30,31]
of EIAV in vitro is 789 times more efficient than sensitive cell-free virus transmission [34]. Cell-to-cell transmission of sensitive and resistant viruses has similar efficiency [7,34], so we set β cs = β cr . Initial conditions for model (1) are summarized in Table 3. The target cells for EIAV infection are tissue macrophages. This cell population arises from monocytes in the circulation that enter the tissue and differentiate into macrophages. The initial value for healthy target cells T(0) was estimated by multiplying the average monocyte count in the peripheral blood of four horses (measured as 353.25 cells/μl [20]) by the rate of monocyte turnover in one day (due to dying or differentiating into tissue macrophages), which is 0.12 [13]. The value for Vs (0) was estimated from the initial inoculum of infectious virus for an experimental infection of 106 TCID50 [20,30,31] and the horse plasma volume (3 L, since the horses weighed 60 kg at the time of the experiment and plasma volume is 5% of body weight) using the conversion factor between TCID50 and plaque forming units (PFU) of 0.7 [5]. Figs. 1 and 2 illustrate the dynamics of the ODE model in two cases. In Fig. 1, ps = pr = 0.9, the coexistence equilibrium Ers is stable and is reached quickly, after only 5 days. In Fig. 2, ps = 0.9, pr = 0.8ps and H(T) < G(T). Initially both strains coexist, but eventually the sensitive strain is replaced by the resistant strain equilibrium due to mutation. We note that the equilibrium values for virus and cells are at physiological levels in the simulations in Figs. 1 and 2.
1 − ps 1 − 0.8ps > , 0.8ps ps then the resistant strain dominates (Figs. 2 and 4). Graphs of the reproduction number for the sensitive strain R0s are illustrated as a function of ps (pr = ps ) and the parameters β fs , β cs , ds and cs in Figs. 5–8. If ps = pr , then R0 = R0s . In Figs. 5–8, the graph of R0s and the contours are on the left and right, respectively. To eradicate EIAV, R0 should be reduced below one (Theorem 3.3). The value of R0 can be reduced by decreasing the transmission rates β fs or β cs or by increasing the death rate of infected cells ds or the clearance rate of cs . The contours show that the greatest reduction in R0s can be achieved by decreasing β cs (cell-to-cell transmission) and increasing ds (death rate of sensitive cells). These two parameters β cs and ds have a greater effect on R0s than free-virus transmission β fs or the clearance rate of the virus cs . Similar results hold for R0r and parameters β fr , β cr , dr and cr . In Fig. 9, one sample path is plotted for the Markov chain model for initial concentrations of
(T (0), Is (0), Ir (0), Vs (0), Vr (0)) = (T¯0 , 0, 0, 233, 0) with either ps = pr = 0.9 or with ps = 0.9, pr = 0.8ps = 0.72 and parameter values from Table 2. Cell and virion concentrations are per milliliter. Therefore clearance is with respect to 1 ml, as in, for
EIAV Sensitive Strain
EIAV Resistant Strain 5 Log10 of T, Ir, Vr
Log10 of T, Is, Vs
6
4 T I
2
s
V
T Ir
4
Vr
3 2 1
s
0 0
5 Days
0 0
10
5 Days
10
5 Log10 of T, Ir, Vr
Log10 of T, Is, Vs
6
4 T Is
2
Vs 0 0
100
200 300 Days
400
500
T Ir
4
Vr
3 2 1 0 0
100
200 300 Days
400
500
Fig. 1. Graph of a solution of the ODE model (1) is plotted for parameter values in Table 2 and initial values (T (0), Is (0), Ir (0), Vs (0), Vr (0)) = (T¯0 , 0, 0, 233, 0), where R0s = 63.2, R0r = 49.6, ps = pr = 0.9. The coexistence strain equilibrium is stable, Ers = (Tˆrs , Iˆs , Iˆr , Vˆs , Vˆr ) = (670, 4.17 × 104 , 5.8, 1.41 × 106 , 195), G(Tˆrs ) < H(Tˆrs ) = 1. The figures on the top show the dynamics of the sensitive and resistant strains for the time interval [0, 10], whereas the figures on the bottom are for the time interval [0, 500]. The coexistence equilibrium is reached after 5 days.
Please cite this article as: L.J.S. Allen, E.J. Schwartz, Free-virus and cell-to-cell transmission in models of equine infectious anemia virus infection, Mathematical Biosciences (2015), http://dx.doi.org/10.1016/j.mbs.2015.04.001
ARTICLE IN PRESS
JID: MBS
[m5G;April 23, 2015;15:3]
L.J.S. Allen, E.J. Schwartz / Mathematical Biosciences 000 (2015) 1–12
EIAV Sensitive Strain
7
EIAV Resistant Strain
6 Log10 of T, Ir, Vr
Log10 of T, Is, Vs
6 4
T I
2
s
4 T I
2
r
V
V
s
0 0
r
5 Days
0 0
10
5 Days
10
T Is
6 Log10 of T, Ir, Vr
Log10 of T, Is, Vs
6
V
4
s
2
4 T Ir
2
V
r
0 0
100
200 300 Days
400
500
0 0
100
200 300 Days
400
500
Fig. 2. Graph of the solution of the ODE model (1) is plotted for parameter values in Table 2 and initial values (T (0), Is (0), Ir (0), Vs (0), Vr (0)) = (T¯0 , 0, 0, 233, 0), where R0s = 63.2, R0r = 130.4, ps = 0.9, pr = 0.8ps . The resistant strain equilibrium is stable, Er = (T¯r , I¯r , V¯ r ) = (325, 4.21 × 104 , 1.14 × 106 ), H(T¯r ) < G(T¯r ) = 1. The figures on the top show the dynamics of the sensitive and resistant strains for the time interval [0, 10], whereas the figures on the bottom are for the time interval [0, 500]. The resistant strain appears after 5 days but the sensitive strain is present until 500 days.
H(T,ps)
G(T,ps)
1
1
50
p =p r
0.8
50
s
0.8 150
0.6
150
ps
p
s
0.6 250
0.4
250
0.4
100
100
0.2
350
200 300
0
1
2
0.2
300
400
3
0
4
T
350
200
1
2
3 T
4
x 10
400 4 4
x 10
Fig. 3. Contours of G and H, Eqs. (11) and (13), are plotted as functions of T and ps when pr = ps ; G(T, ps ) < H(T, ps ). Coexistence of both strains.
H(T,p )
G(T,p )
s
s
1
1
50
p =0.8p r
0.8
s
50 150
0.8 150
0.6
250
p
s
ps
0.6 250
0.4
0.4
100
100 0.2
350
200 300
0
1
2
3 T
350
200
0.2
300
400
400 0
4 4
x 10
1
2
3 T
4 4
x 10
Fig. 4. Contours of G and H, Eqs. (11) and (13), are plotted as functions of T and ps when pr = 0.8ps ; H(T, ps ) < G(T, ps ). Resistant strain dominates.
Please cite this article as: L.J.S. Allen, E.J. Schwartz, Free-virus and cell-to-cell transmission in models of equine infectious anemia virus infection, Mathematical Biosciences (2015), http://dx.doi.org/10.1016/j.mbs.2015.04.001
ARTICLE IN PRESS
JID: MBS 8
[m5G;April 23, 2015;15:3]
L.J.S. Allen, E.J. Schwartz / Mathematical Biosciences 000 (2015) 1–12
1
500
0.8
300
0.6
150 200
250
R
ps
0s
400
100
50
200
0.4
300 350
100 0.2 0 5
400
0 −6
x 10
βfs
1
0
450
0 0
0.5 ps
1
2
βfs
3
4 −6
x 10
Fig. 5. Graph of R0s and contours as a function of ps and β fs .
1 50 500
0.8 150
0.6 s
300
p
R0s
400
200
250
0.4 100
100 0.2
0 1
350
200 300
5 0.5
0 0
−4
2
x 10 ps
0 0
4
β
cs
βcs
400 −4
x 10
Fig. 6. Graph of R0s and contours as a function of ps and β cs .
1 50 600
100
0.6
200 250
R0s
ps
400
0.8
150
0.4
200
350 0.2 450
0 0.1
0 0.5 0.05 d s
1
500
0 0.04
400 0.06
p
300
ds
0.08
0.1
s
Fig. 7. Graph of R0s and contours as a function of ps and ds .
example, separate chambers in experimentation in vitro. The probability an infection is not cleared is 1 − Pext = 0.57 whereas the probability of infection clearance is Pext = q233 32 = 0.43 (case (iii) in Theorem 4.1). The estimate for Pext is in good agreement with the Markov chain model. The proportion of the sample paths out of 10,000 for which Is + Ir + Vs + Vr hit zero before reaching a size of 10,000 in the Markov chain model is close to 0.43 (e.g., 0.4241 for ps = pr = 0.9 and 0.4316 for pr = 0.9, ps = 0.8pr ). The branching process approximation
accurately predicts clearance but does not predict which infected state is approached when the infection persists. As can be seen in the sample paths in Fig. 9, the sensitive strain is present in both sample paths. After five days, mutation to the resistant strain has only occurred in the sample path in Fig. 9 (bottom) when ps = 0.8pr . The top and bottom graphs in Fig. 9 should be compared with the top graphs in Figs. 1 and 2, respectively, where the corresponding solutions of the ODE model are plotted.
Please cite this article as: L.J.S. Allen, E.J. Schwartz, Free-virus and cell-to-cell transmission in models of equine infectious anemia virus infection, Mathematical Biosciences (2015), http://dx.doi.org/10.1016/j.mbs.2015.04.001
ARTICLE IN PRESS
JID: MBS
[m5G;April 23, 2015;15:3]
L.J.S. Allen, E.J. Schwartz / Mathematical Biosciences 000 (2015) 1–12
9
1 500
50
100 0.8
400
150 200
0.6
250
p
R0s
s
300 200
0.4
300 350
100
0.2
0 1
450
0
20
0.5
400
5
10 0
p
s
0
10 c
15
20
s
cs
Fig. 8. Graph of R0s and contours as a function of ps and cs .
EIAV Sensitive Strain
EIAV Resistant Strain 5 of T, Ir, Vr
ps=pr=0.9
T I
2
Log
Log
10
4
10
of T, Is, Vs
6
s
Vs 0 0
1
2
3
4
4 ps=pr=0.9
3 2 1 0 0
5
T I r
Vr 1
2
Days
3
4
5
3
4
5
Days 5 of T, Ir, Vr
ps=0.9, pr=0.72
T I
2
Log
Log
10
4
10
of T, Is, Vs
6
s
Vs 0 0
1
2
3
4
5
Days
4 ps=0.9, pr=0.72
3 2 1 0 0
T I r
Vr 1
2 Days
Fig. 9. Graphs of two sample paths of the Markov chain model are plotted for parameter values in Table 2 and initial values (T (0), Is (0), Ir (0), Vs (0), Vr (0)) = (T¯0 , 0, 0, 233, 0). In the top two figures, sample paths corresponding to the sensitive strain and the resistant strain are plotted for R0s = 63.2, R0r = 49.6, ps = pr = 0.9 (compare with Fig. 1). In the bottom two figures, sample paths corresponding to the sensitive strain and the resistant strain are plotted for R0s = 63.2, R0r = 130.4, ps = 0.9, pr = 0.8ps (compare with Fig. 2). The probability of infection clearance is approximately Pext = q233 32 = 0.4302.
Some estimates for the probability of extinction Pext from the branching process approximation are plotted in Fig. 10. The value ir0 vs0 vr0 of Pext = qis0 12 q22 q32 q42 depends on the parameter values (Table 2) and on the initial conditions. An expression for the probabilities qij is given in Theorem 4.1 part (iii). A set of four initial conditions is assumed: (is0, ir0, vs0, vr0) = (1, 0, 0, 0), (0, 1, 0, 0), (0, 0, 233, 0) and (0, 0, 0, 233) with either one infected cell/ml or 233 virions/ml of the sensitive strain or the resistant strain, respectively. Here, we assume pr = ps . Therefore, Pext is plotted as a function of ps . As ps increases, cell-to-cell transmission decreases and free-virus transmission increases. For ps = 0, there is exclusive cell-to-cell transmission and no free virus is produced. For ps = 0, the graphs illustrate that exposure of healthy tissue to infected cells, as low as one infected
cell/ml, is almost certain to result in infection with Pext = qi2 ≈ 0, i = 1, 2 whereas exposure to a virion concentration of 233/ml (or any convr0 centration), results in clearance of infection with Pext = qvs0 32 q42 = 1. On the other hand, for ps = 1, there is exclusive free-virus transmission and no cell-to-cell transmission. In this case, both exposure to infected cells or virions at concentrations of 1/ml or 233/ml, respectively, may either result in infection clearance with Pext > 0, or viral persistence with 1 − Pext > 0. As the value of ps increases, the probability of infection clearance decreases with exposure to a fixed concentration of virions (233 virions/ml, no infected cells) whereas the probability of infection clearance increases with exposure to a fixed concentration of infected cells (1 cell/ml, no virions) but at a much lower rate.
Please cite this article as: L.J.S. Allen, E.J. Schwartz, Free-virus and cell-to-cell transmission in models of equine infectious anemia virus infection, Mathematical Biosciences (2015), http://dx.doi.org/10.1016/j.mbs.2015.04.001
ARTICLE IN PRESS
JID: MBS 10
[m5G;April 23, 2015;15:3]
L.J.S. Allen, E.J. Schwartz / Mathematical Biosciences 000 (2015) 1–12
1
0.8
Pext
0.6
q12 q22
q233 0.4
32 233
q42
0.2
0 0
0.2
0.4
ps
0.6
0.8
1
Fig. 10. Probability of extinction Pext is plotted as a function of ps (pr = ps ) and initial values (T (0), Is (0), Ir (0), Vs (0), Vr (0)) = (T¯0 , is0, ir0, vs0, vr0), where Pext = ir0 vs0 vr0 qis0 12 q22 q32 q42 . The four initial conditions are (is0, ir0, vs0, vr0) = (1, 0, 0, 0), (0, 1, 0, 0), (0, 0, 233, 0) and (0, 0, 0, 233) representing either one infected cell/ml or 233 virions/ml from either the sensitive strain or the resistant strain. With exclusive cell-to-cell transmission, ps = pr = 0, healthy tissue exposed to infected cells at concentrations as low as 1 cell/ml results in almost no clearance of infection with Pext = qi2 ≈ 0, i = 1, 2, whereas exposure to virions at any concentration results in vr0 complete clearance with Pext = qvs0 32 q42 = 1. With exclusive free-virus transmission, ps = pr = 1, exposure to infected cells or virions results either in infection clearance with Pext > 0 or viral persistence with 1 − Pext > 0.
6. Conclusions In this study, we presented novel mathematical models of EIAV infection dynamics to evaluate the role of cell-to-cell transmission among two competing viral strains. One strain is sensitive to antibody neutralization and the other strain is neutralization resistant. These models, which compare the effects of the proportion of cell-tocell transmission versus free virus transmission, are among the first to consider multiple viral strains as well as multiple transmission modes together (see e.g., Pourbashash et al. [22]). Our models also include mutation of sensitive virus to the resistant strain. In our analyses of an ODE model and a Markov chain model with a branching process approximation, we calculated R0 and determined the probabilities of infection clearance. We also provided analyses of the IFE and two endemic equilibria, one with the coexistence of both strains and another with existence of only the resistant strain. Our results indicate that whether an EIAV infection evolves to the resistant virus steady state (Er ) or the coexistence steady state (Ers ) depends crucially on the proportion of cell-to-cell transmission versus free virus transmission (i.e., the values of pr and ps ) (Figs. 1–4). The same proportion of cell-to-cell transmission for the sensitive and resistant viral strains (i.e., pr = ps ) leads to coexistence of both strains. If the ratio of cell-to-cell to free virus transmission is greater in the resistant strain as compared to the sensitive strain under specified conditions, however, the resistant strain will eventually dominate. For eradication of both EIAV strains, we found that R0 is more sensitive to changes in cell-to-cell transmission (that is, the cell-to-cell transmission rates β cs and β cr and the death rates of infected cells ds and dr ) than free-virus transmission parameters (that is, the free virus transmission rates β fs and β fr and the virus clearance rates cs and cr ) (Figs. 5–8). Clearance of infection depends on cell-to-cell versus free-virus transmission, i.e. whether pr and ps are small or large, that is, whether cell-to-cell or free-virus transmission dominates, respectively (Figs. 9 and 10). This investigation adds to previous work in within-host virus dynamics by meeting the call by Perelson and Ribeiro for research on models of the role that cell-to-cell transmission may play in virus
evasion of immune responses [21]. Since only some broadly neutralizing antibodies effectively block cell-to-cell transmission of HIV-1 [18], cell-to-cell transmission is a mechanism by which HIV-1 can evade antibody control. As such, this evasion mechanism will need to be taken into account in vaccine development. A few limitations of this work should be mentioned. Our models assume that back mutation from the resistant strain to the sensitive strain is negligible. This is a reasonable assumption, however, because reversion to the sensitive strain would be rare in the continued presence of antibody pressure. Back mutation could easily be added to the model. Another simplifying model assumption is that the cells and viruses are homogeneously mixed. Also, the results shown in our simulations are specific to the values chosen for pr and ps . We also note that some of the parameter values we use are specific for EIAVinfected SCID horses, which lack a competent immune system. Consequently, these parameters are representative of the rates of these processes in the absence of immune responses; rates may be different in the presence of immune responses. In this study we have focused on parameters relevant to EIAV infection, and it will be interesting to see if similar behaviors arise with related viruses, such as HIV-1, which also has multiple strains and uses both modes of transmission. Our results may also be relevant for understanding strain competition in other viruses like HTLV-1 that use only cell-to-cell transmission [4]. Importantly, our results have implications for potential vaccines for EIAV infection. This study predicts that treatment modalities that only block free virus infection may lead to resistance, whereas treatments that can block transmission by cell-to-cell could reduce R0 and the development of resistance. Furthermore, our results suggest that an optimal vaccine would induce immune responses that remove infected cells over responses that improve virus clearance. Future work in this area that would be worthwhile with this model includes using the model to predict the proportion of transmission via free virus infection for the sensitive and resistance strains (i.e., the values of ps and pr ) given known viral loads for sensitive and resistant virus (i.e., Vs and Vr ), as well as using the model to predict transient benefits in the short term or long term. Future work in this area that would be worthwhile with new models include modeling efforts that examine coinfection of cells by both strains, or models of within host infection dynamics of two strains with two modes of transmission with ART or vaccination. Acknowledgments We acknowledge the support of the Mathematical Biosciences Institute as long-term visitors in spring 2012 at which time this research was initiated. We thank an anonymous reviewer for suggestions that improved the original manuscript. Appendix Proof of Theorem 3.1. To verify solutions of (1) are nonnegative, it is only necessary to show that solutions to do not cross the boundary of R5+ . Each of the differential equations evaluated on their respective zero boundaries, T˙ |T=0 , I˙i |Ii =0 and V˙ i |Vi =0 , i = s, r is nonnegative; hence solutions of model (1) are nonnegative. To show solutions are bounded above let Ttot = T + Ir + Is and Vtot = Vr + Vs . Then
T˙ tot ≤ − dT Ttot and
V˙ tot ≤ K T¯0 − cVtot , where T¯0 = /dT , K = ps Ns ds + pr Nr dr and c = min {cr , cs }. By comparison of Ttot with the solution of x˙ = − dT x, x(0) = Ttot (0) ≤ T¯0 and the solution Vtot with the solution of y˙ = K T¯0 − cy, y(0) = Vtot (0) leads to Ttot (t) ≤ x(t) ≤ T¯0 and Vtot (t) ≤ y(t) ≤ max{Vtot (0), K T¯0 /c}.
Please cite this article as: L.J.S. Allen, E.J. Schwartz, Free-virus and cell-to-cell transmission in models of equine infectious anemia virus infection, Mathematical Biosciences (2015), http://dx.doi.org/10.1016/j.mbs.2015.04.001
ARTICLE IN PRESS
JID: MBS
[m5G;April 23, 2015;15:3]
L.J.S. Allen, E.J. Schwartz / Mathematical Biosciences 000 (2015) 1–12
To verify Theorem 3.3, we apply mathematical results of Shuai and van den Driessche to construct a Lyapunov function (Theorem 2.1, p. 515 [27]). Shuai and van den Driessche also state sufficient conditions for global asymptotic stability of the IFE (Theorem 2.2, p. 516 [27]). Unfortunately, these conditions do not apply directly to model (1) but the method of their proof can be used to verify global asymptotic stability. Proof of Theorem 3.3.. Let F = [Fij ] be the matrix defined in (4), x = (Is , Ir , Vs , Vr ) and y = T. Define the vector f(x, y) as
⎤
⎡
(F11 Is + F13 Vs )(T¯0 − T ) ⎥ ⎢ ⎢ (F21 Is + F22 Ir + F23 Vs + F24 Vr )(T¯0 − T ) ⎥ ⎥ f (x, y) = ⎢ ⎥ ⎢ ps βf s Vs T ⎦ ⎣ pr βf r Vr T and matrix Vˆ as
⎡ ⎢ ⎢
0
0
0
dr
0
0
cs
⎥ 0⎥ ⎥. ⎥ 0⎦
−pr Nr dr
0
cr
Vˆ = ⎢ ⎢
⎣ −ps Ns ds 0
0
⎤
ds
Since T (t) ≤ T¯0 for t [0, ) and solutions to model (1) lie in , it follows that f(x, y) 0 for (x, y) ∈ ⊂ R5+ . There are two cases ˆ 0r > R ˆ 0s and (ii) R ˆ0 = R ˆ 0s ≥ R ˆ 0r . In case (i), ˆ0 = R to consider, (i) R ˆ 0s , denote the left eigenvector of Vˆ −1 F corresponding to the ˆ 0r > R R ˆ 0r as ωrT . An eigenvector ωr is ˆ0 = R eigenvalue R
⎡
ωr
⎤ ˆ 0r F21 R ⎢ ˆ 0s ) ⎥ ˆ 0r − R ⎢ F22 (R ⎥ ⎥ =⎢ ⎢ ⎥ ˆ F R ⎣ ⎦ 23 0r ˆ 0s ) ˆ 0r − R F24 (R
ˆ 0s ≥ (or any nonzero constant multiple). Thus, ωr > 0. In case (ii), R ˆ 0r , a left eigenvector ωsT of Vˆ −1 F corresponding to the eigenvalue R ˆ 0s is ˆ0 = R R
⎡
F11
⎤
⎥ ⎢ ⎢ 0 ⎥ ⎥ ⎥ ⎣ F13 ⎦
ωs = ⎢ ⎢
0 (or any nonzero constant multiple), ωs 0. Define a Lyapunov function for model (1) as
11
To verify part (i), assume R0s > 1 ≥ R0r . Let S0 = {(T, Is , Ir , Vs , Vr ) ∈ R5+ |0 ≤ T + Ir + Is ≤ T¯0 , Is > 0, Vs > 0} and S0 = S\S0 . In S0 , either Is = 0 or Vs = 0. Let φ (t, X0 ) be the solution of model (1) with initial condition X0 = (T(0), Is (0), Ir (0), Vs (0), Vr (0)) S0 . That is, if X0 S0 , one of Is (0) or Vs (0) equals zero. The only invariant set in S0 is {E0 }. Equilibrium Er is not feasible since R0r ≤ 1. Let Ws (E0 ) be the stable manifold of E0 in R5+ . To show uniform persistence of the sensitive strain, it suffices to show Ws (E0 ) S0 =
(Theorem 4.3 [11]). Let ε > 0 and define
M =
(1 − ps )(1 − μ)βcs (T¯0 − ) − ds ps (1 − μ)βf s (T¯0 − ) . ps Ns ds −cs − ps βf s T¯0
Choose α sufficiently small such that if ε [0, α ], then Mε has a dominant positive eigenvalue and associated positive eigenvector. This is possible since the spectral abscissa s(M0 ) > 0 iff R0s > 1. In addition, choose α sufficiently small such that T (0) ≥ T¯0 − α > 0 and X0 − E0 α . Suppose φ (t, X0 ) − E0 α for all t 0, that is, φ (t, X0 ) Ws (E0 ). Then
I˙s ≥ [ps (1 − μ)βf s Vs + (1 − ps )(1 − μ)βcs Is ](T¯0 − α) − ds Is V˙ s ≥ ps Ns ds Is − cs Vs − ps βf s Vs T¯0 . The solution Is (t) and Vs (t) can be compared to the solution of the linear system Y = Mε Y, with Y(0) = (Is (0), Vs (0))T . It follows that the solution Is (t) and Vs (t) → as t → , a contradiction to boundedness of solutions. Hence, Ws (E0 ) S0 = . By Theorem 4.3 [11], the conclusion of Theorem 3.4 part (i) follows. To verify part (ii), assume R0r > 1. Let S0 = {(T, Is , Ir , Vs , Vr ) ∈ R5+ |0 ≤ T + Ir + Is ≤ T¯0 , Ir > 0, Vr > 0} and S0 = S\S0 . In S0 , either Ir = 0 or Vr = 0. Let φ (t, X0 ) be the solution of model (1) with initial condition X0 = (T(0), Is (0), Ir (0), Vs (0), Vr (0)) S0 . The only invariant set in S0 is {E0 }. As in the proof of part (i), to show uniform persistence of the resistant strain, it suffices to show Ws (E0 ) S0 = (Theorem 4.3 [11]). Let ε > 0 and define
M =
(1 − pr )βcr (T¯0 − ) − dr
pr βf r (T¯0 − )
pr Nr dr
−cr − pr βf r T¯0
.
Choose α sufficiently small such that if ε [0, α ], then Mε has a positive dominant eigenvalue s(Mε ) > 0 and associated positive eigenvector (since s(M0 ) > 0 iff R0r > 1). In addition, choose α sufficiently small such that T (0) ≥ T¯0 − α and X0 − E0 α . Suppose φ (t, X0 ) − E0 α for all t 0. Then
I˙r ≥ [pr βf r Vr + (1 − pr )βcr Ir ](T¯0 − α) − dr Ir
Li (x, y) = ωiT Vˆ −1 x,
V˙ r ≥ pr Nr dr Ir − cr Vr − pr βf r Vr T¯0 .
where i = r for case (i) or i = s for case (ii). Then as in the proof of Theorem 2.1 in [27], differentiation of Li along solutions of (1) leads to
By comparison with Y = Mε Y, Y(0) = (Ir (0), Vr (0))T , it follows that solutions Ir (t) and Vr (t) → as t → , a contradiction. Hence, Ws (E0 ) S0 = . By Theorem 4.3 [11], the conclusion of Theorem 3.4 part (ii) follows.
ˆ 0 − 1)ωiT x − ωiT Vˆ −1 f (x, y) L i (x, y) = (R ˆ 0 ≤ 1, then Li is a Lyasince dx/dt = (F − Vˆ )x − f (x, y). Therefore, if R punov function for model (1) [27]. To show global asymptotic stability of the IFE, E0 , we need to show that the only invariant set of L i = 0 is ˆ 0r ≤ 1 and L r = 0, the the IFE [17]. In the case i = r, ωr > 0, so that if R ˆ 0r ≤ R ˆ 0s ≤ 1, only invariant set is the IFE, {E0 }. In the case i = s, if R ˆ 0r ≤ 1 then L s = 0 if Is = 0 = Vs or if Vs = 0 and T = T¯0 . Since R (R0r < 1) the only equilibrium solution to model (1) with this restriction (Theorem 3.1) is the IFE. That is, the only invariant set of L s = 0 is the IFE, {E0 }. It follows that the IFE E0 is globally asymptotically stable in [17]. Proof of Theorem 3.4.. Solutions of (1) lie in a compact set in R5+ (Theorem 3.1). Therefore, model (1) is point-dissipative. Let S = {(T, Is , Ir , Vs , Vr ) ∈ R5+ |0 ≤ T + Ir + Is ≤ T¯0 }.
We compute the fixed points for the probability of extinction in the Galton–Watson branching process approximation of the Markov chain model (Table 1). The extinction results for the branching process approximation are summarized in Theorem 4.1. The expectation matrix of the pgfs f1 , . . . , f4 in (16) is M = [fi /uj ] evaluated at ui = 1, i = 1, 2, 3, 4 is ⎡ ⎢ ⎢ ⎢ W −1 ⎢ ⎢ ⎣
(1−ps )βcs T¯0 (2−μ) + ps Ns ds
(1 − ps )μβcs T¯0
ps Ns ds
0
2(1−pr )βcr T¯0 + pr Nr dr
0
ps (1 − μ)βf s T¯0
ps μβf s T¯0
0
0
pr βf r T¯0
0
0
⎤
⎥ pr Nr dr ⎥ ⎥ ⎥ 0 ⎥ ⎦ 0
where W = diag[w1 , w2 , w3 , w4 ] and wi is the denominator of each of the fi , i = 1, 2, 3, 4. It can be easily verified that [W (M − I)]T = F − V
Please cite this article as: L.J.S. Allen, E.J. Schwartz, Free-virus and cell-to-cell transmission in models of equine infectious anemia virus infection, Mathematical Biosciences (2015), http://dx.doi.org/10.1016/j.mbs.2015.04.001
ARTICLE IN PRESS
JID: MBS 12
L.J.S. Allen, E.J. Schwartz / Mathematical Biosciences 000 (2015) 1–12
and I is the 4 × 4 identity matrix [2]. In particular, ρ (M) > 1 (ρ (M) < 1) iff R0 > 1 (R0 < 1). The probability of extinction is one if ρ (M) < 1. If matrix M is irreducible and ρ (M) > 1, then there is a unique fixed point of the pgfs, (q1 , q2 , q2 , q4 ) (0, 1)4 [2]. As can be seen matrix M is reducible; there is not a unique fixed point of the pgfs in (0, 1)4 . In particular, we will show that there are two fixed points, other than (1,1,1,1), in the region (0, 1]4 when ρ (M) > 1 (i.e., R0 > 1). These two fixed points approximate infection clearance when R0 > 1 in the two cases: R0r ≤ 1 < R0s and R0r > 1. Given (Is (0), Ir (0), Vs (0), Vr (0)) = (is0 , ir0 , vs0 , vr0 ) are initially small values, then the probability of infection clearance is approximately (see e.g. [1,2,10,14]): 0 vr0 Pext = qis1 0 qir2 0 qvs 3 q4 .
(17)
Alternately, with probability 1 − Pext , an infection becomes established, which for our model may be either the resistant strain or both sensitive and resistant strains. If R0 ≤ 1, the pgfs (16) have only one fixed point, namely (1, 1, 1, 1). But if R0 > 1, there exist two more fixed points which lie in (0, 1]4 . Denote these two fixed points as (q1i , q2i , q3i , q4i ), i = 1, 2. The pgfs f2 and f4 are functions of u2 and u4 ; therefore, they can be solved independently of f1 and f3 . The functions f2 and f4 have two fixed points: (q21 , q41 ) = (1, 1) if R0 ≤ 1 and a unique fixed point (q22 , q42 ) (0, 1)2 if R0r > 1 given by
q22 =
1
and q42 =
R0r
pr βf r T¯0 cr + pr βf r T¯0
1
+
R0r
cr cr + pr βf r T¯0
.
(18)
These two fixed points of f2 and f4 are used to find the corresponding fixed points of f1 and f3 . For the case, R0r ≤ 1 with fixed point (q2 , q4 ) = (1, 1) of f2 and f4 , the values of the fixed point of f1 and f3 is (q1 , q3 ) = (1, 1) if R0s ≤ 1. However, if R0s > 1, then the values of the fixed point (q1 , q3 ) of f1 and f3 are
q11 =
[m5G;April 23, 2015;15:3]
1
and q31 =
R0s
ps (1 − μ)βf s T¯0 cs + pss βf s T¯0
1 R0s
+
cs + ps μβf s T¯0 . cs + ps βf s T¯0 (19)
For the case R0r > 1 with fixed point values (q22 , q42 ) of f2 and f4 given by (18), the value of q32 depends only on u1 :
q32 = f3 (u1 , q22 , q32 , q42 ) =
cs + ps (1 − μ)βf s T¯0 u1 + ps μβf s T¯0 q22 . cs + ps βf s T¯0 (20)
Note that if u1 = q12 [0, 1], then q32 (0, 1). Substitution of q32 into f1 leads to a quadratic function in u1 : f12 (u1 ) = f1 (u1 , q22 , q32 , q42 ). Since f12 (0) > 0 and f12 (1) < 1, there is only one fixed point of f12 in [0, 1]. In particular, the fixed point q12 (0, 1). Since f11 (u1 ) = f1 (u1 , 1, f3 (u1 , 1), 1) > f12 (u1 ) for u1 (0, 1) and f11 (q11 ) = q11 = 1/R0s , it follows that
0 < q12
JID: MBS
[m5G;April 23, 2015;15:3]
Mathematical Biosciences 000 (2015) 1–12
Contents lists available at ScienceDirect
Mathematical Biosciences journal homepage: www.elsevier.com/locate/mbs
Free-virus and cell-to-cell transmission in models of equine infectious anemia virus infection Linda J.S. Allen a,∗, Elissa J. Schwartz b a b
Department of Mathematics and Statistics, Texas Tech University, Lubbock, TX 79409-1042, United States School of Biological Sciences and Department of Mathematics, Washington State University, Pullman, WA 99164-3113, United States
a r t i c l e
i n f o
Article history: Available online xxx Keywords: Branching process Markov chain Within-host model Virus dynamics Equine infectious anemia virus
a b s t r a c t Equine infectious anemia virus (EIAV) is a lentivirus in the retrovirus family that infects horses and ponies. Two strains, referred to as the sensitive strain and the resistant strain, have been isolated from an experimentallyinfected pony. The sensitive strain is vulnerable to neutralization by antibodies whereas the resistant strain is neutralization-insensitive. The sensitive strain mutates to the resistant strain. EIAV may infect healthy target cells via free virus or alternatively, directly from an infected target cell through cell-to-cell transfer. The proportion of transmission from free-virus or from cell-to-cell transmission is unknown. A system of ordinary differential equations (ODEs) is formulated for the virus-cell dynamics of EIAV. In addition, a Markov chain model and a branching process approximation near the infection-free equilibrium (IFE) are formulated. The basic reproduction number R0 is defined as the maximum of two reproduction numbers, R0s and R0r , one for the sensitive strain and one for the resistant strain. The IFE is shown to be globally asymptotically stable for the ODE model in a special case when the basic reproduction number is less than one. In addition, two endemic equilibria exist, a coexistence equilibrium and a resistant strain equilibrium. It is shown that if R0 > 1, the infection persists with at least one of the two strains. However, for small infectious doses, the sensitive strain and the resistant strain may not persist in the Markov chain model. Parameter values applicable to EIAV are used to illustrate the dynamics of the ODE and the Markov chain models. The examples highlight the importance of the proportion of cell-to-cell versus free-virus transmission that either leads to infection clearance or to infection persistence with either coexistence of both strains or to dominance by the resistant strain. © 2015 Elsevier Inc. All rights reserved.
1. Introduction The lentiviruses human immunodeficiency virus type 1 (HIV-1) and equine infectious anemia virus (EIAV) establish persistent, lifelong infections in their hosts. EIAV infects cells of the monocyte/macrophage lineage in its host, which are equids including horses and ponies. While HIV-1 infection results in CD4 T cell depletion and progression to AIDS, immune responses including antibodies and cytotoxic T lymphocytes are able to control EIAV infection. EIAV infected animals live, for the most part, without disease, and viral loads are maintained at low levels. Thus the study of EIAV infection may provide critical clues to the pathogenesis and control of lentiviral infection in other species, notably humans. Since mathematical modeling has been a useful tool to advance our understanding of virus dynamics, we investigate a model of
∗
Corresponding author. Tel.: +1 8067422580. E-mail address: [email protected] (L.J.S. Allen).
viral infection at the within-host level that includes two modes of viral transmission, free-virus and direct cell-to-cell transmission. Our model is applicable to EIAV, as it displays both modes of transmission [34]. Cell-to-cell transmission has been observed in enveloped viruses, viruses surrounded by a lipid envelope such as herpesviruses, paramyxoviruses, pox viruses, rhabdoviruses and retroviruses [23]. In cell-to-cell transmission, viruses are transferred from infected cells to uninfected cells via cell adhesions called virological synapses. One advantage of cell-to-cell transmission over free-virus transmission is that the virus is not exposed to the extracellular environment. This allows the virus to evade immune responses, such as the neutralizing antibody response [6,7,23]. Furthermore, cell-tocell transmission permits viral spread despite antiretroviral therapy [28]. In EIAV, a neutralization-resistant variant emerged over time in an experimentally-infected pony [29]. This resistant variant shows lower infectiousness but more efficient transmission than the sensitive strain [34]. In fact, the efficiency of EIAV cell-to-cell spread for this strain was found to be 1000-fold greater than free virus transmission in vitro [34]. To better understand the role of transmission
http://dx.doi.org/10.1016/j.mbs.2015.04.001 0025-5564/© 2015 Elsevier Inc. All rights reserved.
Please cite this article as: L.J.S. Allen, E.J. Schwartz, Free-virus and cell-to-cell transmission in models of equine infectious anemia virus infection, Mathematical Biosciences (2015), http://dx.doi.org/10.1016/j.mbs.2015.04.001
JID: MBS 2
ARTICLE IN PRESS
[m5G;April 23, 2015;15:3]
L.J.S. Allen, E.J. Schwartz / Mathematical Biosciences 000 (2015) 1–12
modes among neutralization sensitive and resistant strains, here we incorporate both modes of transmission in our EIAV model, which includes two viral strains. Others have modeled cell-to-cell viral transmission, but they have concentrated on models with a single viral strain. Komarova et al. [15,16] and Sigal et al. [28] each looked at sensitivity to antiretroviral therapy (ART), using deterministic and stochastic infection models, respectively. Komarova et al. considered cells multiply infected with a single HIV strain and reported how the relative probability of individual virions to infect a cell during free virus and cell-to-cell transmission affects ART susceptibility [16]. Komarova et al. also fit a mathematical model to previously published and newly generated in vitro data and determined that free virus and cell-to-cell transmission contribute approximately equally to the growth of the virus population [15]. Pourbashash et al. [22] presented a global analysis for a single strain model and for a special case of a multiple strain competition model. Both models included cell-to-cell and free-virus transmission. Analysis of their multi-strain model showed that only one strain dominated, the one with the largest strain reproduction number [22]. The study by Pourbashash et al. [22] is one of the first to include multiple viral strains and two modes of transmission. Lentiviruses such as EIAV and HIV are characterized by coinfection by multiple strains. Therefore it is critical to consider multiple strain competition, and emergence of resistance, particularly in studies of virus escape from control. The current work expands upon previous models that include multiple viral strains (see e.g. [3,9,22]) and models with two modes of transmission [15,16,22,28] by investigating how cell-to-cell transmission plays a role in multiple strain coexistence versus dominance by a resistant strain. The current work also advances the study of EIAV viral dynamics by considering the effect of cell-to-cell transmission. Previous work in EIAV modeling has investigated a single strain model with immune responses [24], the role of antibody neutralization on competition and escape with two viral strains [8], and the correlates of antibody protection and escape in EIAV infection [25], but these studies did not address how the mode of infection (i.e., cell-to-cell versus free virus transmission) affects virus escape and the development of resistance, as in the current study. Our goals in this investigation are to provide a preliminary mathematical analysis of within-host viral infection for EIAV using both modes of transmission to determine conditions that allow one or both strains to persist and to provide an analysis for the probability of infection clearance. In Section 2, the new within-host model with two viral strains and with free-virus and cell-to-cell transmission is formulated. The model is a system of ordinary differential equations (ODEs), an extension of the model recently studied by Pourbashash et al. [22] to account for the mutation from the sensitive strain to the resistant strain and for the proportion of cell-to-cell versus freevirus transmission. Five variables in the model include healthy target cells, sensitive and resistant free viruses, and cells infected by the sensitive or the resistant viral strain. In Section 3, the basic reproduction number for the model is computed and conditions for existence of three equilibria are verified, the infection-free equilibrium (IFE) and two endemic equilibria. Global stability of the IFE is verified in a special case when the basic reproduction number is less than one and uniform persistence of the infection when the basic reproduction number exceeds one. In Section 4, a continuous-time Markov chain model is formulated based on the ODE model. The probability that an infection is cleared is estimated from multitype branching process theory. In Section 5, parameter values representing two variants of EIAV are used to illustrate the persistence of one or both viral strains in the ODE model and infection clearance in the Markov chain model. 2. ODE model The within-host ODE model includes two viral variants transmitted via free viral entry into a cell or direct cell-to-cell transmission.
The five variables in the ODE model are denoted T = target cells, Is = target cells infected with the sensitive strain, Ir = target cells infected with the resistant strain, Vs = free virions of the sensitive strain, and Vr = free virions of the resistant strain. The terms β fi VT and β ci IT represent transmission rates of free-virus and cell-to-cell, respectively, from either the sensitive i = s or the resistant viral strain, i = r. The parameter pi represents the proportion of transmission from the free virus and 1 − pi the proportion that is cell-to-cell for either the sensitive (i = s) or the resistant (i = r) strain, 0 pi 1. The particular cases pi = 0 or pi = 1 represent models with exclusive cell-to-cell transmission or with exclusive free-virus transmission, respectively. Mutation occurs within a cell infected with the sensitive strain. The proportion of cells infected with the sensitive virus that mutates to the resistant strain is μ, 0 < μ < 1. We assume there is no back mutation from the resistant strain to the sensitive strain. With these assumptions, the ODE model takes the following form:
T˙ = − dT T − [ps βf s Vs + pr βf r Vr + (1 − ps )βcs Is + (1 − pr )βcr Ir ]T I˙s = [ps (1 − μ)βf s Vs + (1 − ps )(1 − μ)βcs Is ]T − ds Is I˙r = [pr βf r Vr + (1 − pr )βcr Ir + ps μβf s Vs + (1 − ps )μβcs Is ]T − dr Ir V˙ s = ps Ns ds Is − cs Vs − ps βf s Vs T V˙ r = pr Nr dr Ir − cr Vr − pr βf r Vr T.
(1)
Healthy target cells are produced at a constant rate . Parameters dT , ds and dr are the death rates of healthy or infected target cells, either sensitive or resistant. It is assumed that the death rates for infected cells are greater than or equal to that of the healthy cells, that is,
dT ≤ max{dr , ds }. Parameters cr and cs are the clearance rates of the resistant and sensitive viral strains. The parameter Ni , i = s, r, is referred to as the burst number, the number of free viruses produced by either a sensitive or a resistant infected cell. The probability pi , i = s, r, is multiplied by the burst number, representing the proportion of virions released for free viral transmission. The terms pi β fi Vi T in the differential equations for Vi account for loss of free virus via entry into target cells. Often these terms are neglected in mathematical models since they are small compared to the viral clearance rates. All parameters in model (1) are positive with the exception of pi (0 pi 1). Initial values for model (1) are nonnegative, T(0) = T0 > 0, Ii (0) = Ii0 0, and Vi (0) = Vi0 0, i = r, s with Vs0 > 0 and T0 + Is0 + Ir0 /dT . Also, if ps = 0, then Is (0) > 0. 3. Analysis of the ODE model 3.1. Bounded, nonnegative solutions Solutions of model (1) are unique, nonnegative and bounded. The healthy target cell population T¯0 = /dT is an upper bound for the total cell population and in addition, an upper bound for the viral population depends on T¯0 . Specifically, the set defined in Theorem 3.1 is forward invariant for model (1). Solutions of model (1) exist for t [0, ). Verification of these results is relegated to the Appendix. Theorem 3.1. Solutions of model (1) are unique, nonnegative and bounded. The set ⊂ R5+ is forward invariant for model (1), where
= {(T, Is , Ir , Vs , Vr ) ∈ R5+ |0 ≤ T + Ir + Is ≤ T¯0 , 0 ≤ Vs + Vr ≤ V¯ }, (2) T¯0 = /dT , V¯ = max{Vs (0) + Vr (0), K T¯0 /c}, K = ps Ns ds + pr Nr dr and c = min {cr , cs }.
Please cite this article as: L.J.S. Allen, E.J. Schwartz, Free-virus and cell-to-cell transmission in models of equine infectious anemia virus infection, Mathematical Biosciences (2015), http://dx.doi.org/10.1016/j.mbs.2015.04.001
ARTICLE IN PRESS
JID: MBS
[m5G;April 23, 2015;15:3]
L.J.S. Allen, E.J. Schwartz / Mathematical Biosciences 000 (2015) 1–12
values of T¯r , I¯r and V¯ r are positive. Let the right side of the differential equation for T˙ + I˙r = 0 and let Is = 0 = Vs . These assumptions yield an expression for cells infected with the resistant strain which depends on T:
3.2. Infection-free equilibrium The infection-free equilibrium (IFE) for model (1) is
E0 =
dT
, 0, 0, 0, 0 = (T¯0 , 0, 0, 0, 0).
(3)
Ir =
Linearization of model (1) about the IFE yields the following Jacobian matrix, J, separated into a matrix representing sources of infection F and another matrix for other transitions V, where J = F − V with
⎤
⎡
(1 − ps )(1 − μ)βcs T¯0 0 ps (1 − μ)βf s T¯0 0 ⎥ ⎢ ⎢ (1 − ps )μβcs T¯0 ( 1 − pr )βcr T¯0 ps μβf s T¯0 pr βf r T¯0 ⎥ ⎥ ⎢
F =⎢
⎣
0
0
0
0
0
0
0
0
⎥ ⎦
(4) and
⎡ ⎢ ⎢
V=⎢ ⎢
ds
0
0
0
0
dr
0
0
0
cs + ps βf s T¯0
0
−pr Nr dr
0
cr + pr βf r T¯0
⎣ −ps Ns ds 0
(5)
The spectral radius of the next generation matrix F V −1 is the basic reproduction number, R0 = ρ(F V −1 ) [33]. Specifically, R0 = max{R0s , R0r },
(6)
where R0s = (1 − ps )
R0r = (1 − pr )
(1 − μ)βcs T¯0 ds
βcr T¯0 dr
+ ps
(1 − μ)βf s T¯0 ps Ns cs + ps βf s T¯0
βf r T¯0 pr Nr + pr . cr + pr βf r T¯0
− dT T dr
.
(9)
The differential equation for Vr yields an expression for the equilibrium Vr in terms of Ir and T. Let V˙ r = 0, then
Vr =
pr Nr dr I¯r . cr + pr βf r T
(10)
Thus, it is only necessary to show there exists a unique positive value for T. Substituting the preceding expressions for Ir and Vr into the right side of T˙ and setting it equal to zero, leads to an identity for T, G(T) = 1, where
G(T ) = (1 − pr )
⎤ ⎥ ⎥ ⎥. ⎥ ⎦
3
βcr T dr
+ pr
βf r Tpr Nr . cr + pr βf r T
(11)
Evaluating G at the IFE value yields G(T¯0 ) = R0r . Since G(0) = 0 and G is a strictly increasing function of T on [0, ), it follows that if R0r > 1, then there exists a unique positive solution T¯r such that G(T¯r ) = 1. In addition, T¯r < T¯0 . This fact and the expressions for Ir in (9) and for Vr in (10) show that a unique endemic equilibrium Er exists (with Is and Vs equal to zero) provided R0r > 1. Second, conditions will be derived for existence and uniqueness of a feasible equilibrium Ers = (Tˆrs , Iˆs , Iˆr , Vˆs , Vˆr ). Identities for Vs and Vr can be obtained in terms of T, Is and Ir :
ps Ns ds Is pr Nr dr Ir and Vr = . cs + ps βf s T cr + pr βf r T
(7)
Vs =
(8)
As in the preceding analysis, an identity for T can be derived from the differential equation for Is set equal to zero and the expression for Vs . This yields H(T) = 1, where
In this calculation, the infected cells are the sources of infections in matrix F and not the viruses, which leads to the commonly-defined expression for the basic reproduction number given in (6), (7) and (8). If the viruses are included in matrix F as sources of infections through the burst number, then an alternate form for the basic reproduction number is obtained which is the square root of the expressions in (7) and (8). This alternate form is equivalent to (6) in the sense that it is also a threshold with respect to one. It follows from [33] that if R0 < 1, then the IFE is locally asymptotically stable and unstable if R0 > 1. The reproduction numbers for the sensitive and resistant strains in (7) and (8) consist of a sum of two terms representing the two modes of transmission. The first term in the sum represents transmission due to direct cell-to-cell transmission and the second term due to freevirus transmission. The factor 1 − μ in R0s accounts for the fraction of infected cells with virus that does not mutate to the resistant strain. The mutated resistant cells are not included in the expression for the basic reproduction number, due to the fact that the cells infected with the sensitive strain change their reproductive strategy after they become infected. Mutation is important for establishment of endemic equilibria containing the resistant strain. 3.3. Endemic equilibria For model (1), there exist two endemic equilibria, an equilibrium with only the resistant strain, denoted Er , and a second equilibrium with both the sensitive and resistant strains, denoted Ers . Given that μ > 0 with no back mutation from resistant to sensitive strain, there does not exist an endemic equilibrium with only the sensitive strain. First, conditions will be derived for existence and uniqueness of a feasible equilibrium Er = (T¯r , 0, I¯r , 0, V¯ r ), that is, conditions so that the
H(T ) = (1 − ps )
(1 − μ)βcs T ds
+ ps
(1 − μ)βf s ps Ns T . cs + ps βf s T
(12)
(13)
Note that H(T¯0 ) = R0s and H is nonnegative and a strictly increasing function of T on [0, ). Therefore, it follows that if R0s > 1, then a unique positive solution for Tˆrs exists satisfying H(Tˆrs ) = 1. Also, Tˆrs < T¯0 . The value of Tˆrs depends on the mutation proportion μ. Let f(T) = − dT T. Then f (Tˆrs ) > 0. Applying the differential equation for T, the identity H(Tˆrs ) = 1, the relation f(T) − dr Ir − ds Is = 0, and either the expression for Vr or for Vs , then the following two identities must be satisfied by Ir and Is at the endemic equilibrium Ers :
f (Tˆrs )μ 1 − G(Tˆrs ) Ir = 1−μ dr (1 − μ)
(14)
1 f (Tˆrs )[1 − G(Tˆrs )] − G(Tˆrs ) Is = , 1−μ ds
(15)
where the function G is defined in (11). To ensure there exist positive solutions Iˆr and Iˆs to (14) and (15), respectively, requires that G(Tˆrs ) < 1. Existence and uniqueness of a feasible Ers requires two conditions: R0s > 1 and G(Tˆrs ) < 1. Either one or two endemic equilibria exist (in R5+ ) if R0r > 1 and R0s > 1. The relation between the functions G(T) and H(T) at Tˆrs determines whether there exist one or two endemic equilibria. If G(Tˆrs ) < H(Tˆrs ) = 1, there exist two endemic equilibria, Ers and Er , and the healthy cells satisfy
0 < Tˆrs < T¯r < T¯0 . The Jacobian matrix of model (1) evaluated at Er shows that Er is unstable when H(T¯r ) > 1. Since H is an increasing function of T,
Please cite this article as: L.J.S. Allen, E.J. Schwartz, Free-virus and cell-to-cell transmission in models of equine infectious anemia virus infection, Mathematical Biosciences (2015), http://dx.doi.org/10.1016/j.mbs.2015.04.001
ARTICLE IN PRESS
JID: MBS 4
[m5G;April 23, 2015;15:3]
L.J.S. Allen, E.J. Schwartz / Mathematical Biosciences 000 (2015) 1–12
H(Tˆrs ) = 1 < H(T¯r ) implies if both Ers and Er exist, then Er is unstable. On the other hand, if H(Tˆrs ) = 1 ≤ G(Tˆrs ), only the resistant strain Er exists. Note that the parameters in H(T) depend on the sensitive strain and the mutation parameter μ, whereas the parameters in G(T) depend on the resistant strain. If, for example, in the sensitive strain, the mutation is large or transmission rates are smaller than in the resistant strain so that H(T) < G(T), then the resistant strain Er will dominate; the coexistence equilibrium Ers will not exist. Only one endemic equilibrium exists if either R0s < 1 < R0r or R0r < 1 < R0s . If R0s < 1 < R0r , then the endemic equilibrium Er ∈ / R5+ . If R0r < 1 < R0s , then Ers ∈ R5+ and Er ∈ / R5+ . In either R5+ and Ers ∈ case, the healthy cell population at Er or Ers is bounded above by T¯0 . The preceding results concerning E0 , Er and Ers are summarized in the following theorem.
(ii) If R0r > 1, then the resistant strain in model (1) is uniformly persistent. That is, there exists a positive constant κ > 0 such that lim inft→∞ X (t) > κ for X equal to Ir or Vr . In the case that R0 > 1, there exist either one or two endemic equilibria in R5+ . Numerical results indicate that if Ers ∈ R5+ , it is stable but if Ers ∈ / R5+ and Er ∈ R5+ , then Er is stable. When the endemic equilibria are feasible, the numerical results and the results of Theorems 3.2–3.4 suggest that either Er or Ers is globally asymptotically stable. If stochasticity in the birth, death and transmission process is included, then even if R0 > 1, the sensitive and the resistant strains may not persist. These results are shown for the Markov chain model in the next section. 4. Markov chain model
Theorem 3.2. There exists at most three feasible equilibria for model (1) which depend on the values of R0 , R0s and R0r , defined in (6), (7) and (8), respectively. The IFE E0 = (T0 , 0, 0, 0, 0) always exists and E0 is the only feasible equilibrium when R0 ≤ 1. If R0 < 1, then E0 is locally asymptotically stable and if R0 > 1, then E0 is unstable and in addition, either one or two endemic equilibria exist in R5+ . An endemic equilibrium Er = (T¯r , 0, I¯r , 0, V¯ r ) exists in R5+ if R0r > 1, where the value of T¯r is defined implicitly by G(T¯r ) = 1 in (11) and I¯r and V¯ r are defined in (9) and (10), respectively. An endemic equilibrium Ers = (Tˆrs , Iˆs , Iˆr , Vˆs , Vˆr ) exists in R5+ if R0s > 1, where Tˆrs is defined implicitly by H(Tˆrs ) = 1 in (13) and G(Tˆrs ) < 1, and the remaining endemic values Iˆr , Iˆs , Vˆr , and Vˆs are defined by the expressions in (12), (14) and (15), respectively. If both Ers and Er exist in R5+ , then Er is unstable. Although we assumed μ > 0 if mutation μ = 0, then the coexistence equilibrium Ers simplifies to a single strain equilibrium with only the sensitive strain, Es . It is interesting to note that Pourbashash et al. [22] results on multi-strain competition can be applied to model (1) if μ = 0, 0 < pi < 1, i = r, s and the terms pi β fi Vi T are removed from the differential equations for V˙ i , i = r, s. In their model, if R0s > max{1, R0r }, then Es is globally asymptotically stable whereas if R0r > max{1, R0s }, then Er is globally asymptotically stable (Theorem 5.2, p. 3352 [22]). 3.4. Global stability of the IFE and uniform persistence
Introduction of virus or infected cells into a healthy individual may not result in infection, especially if the infectious dose (virus and infected cells) is sufficiently small. To study whether an infection can be cleared, we formulate a Markov chain model, where the variables are discrete-valued and random. In addition, a branching process approximation is used to estimate the probability an infection is cleared. The continuous-time Markov chain model is based on the parameters in the ODE model. Fourteen distinct events and their corresponding infinitesimal transition rates, ri t + o(t), are defined in Table 1. The same notation is used as in the ODE model but the random variables are discrete, T, Is , Ir , Vs , Vr {0, 1, 2, . . . }, and time is continuous, t [0, ). The time between events is exponentially distributed with parameter 14 i=1 ri . The changes in the random variables are ±1, as in a birth and death process. For example, sensitive virus entry into target cells with no mutation, event 3, results in the transition (T, Is , Vs ) → (T − 1, Is + 1, Vs − 1) with probability ps (1 − μ)β fs Vs Tt + o(t). 4.1. Probability of infection clearance A Galton–Watson multitype branching process approximation of the Markov chain model is defined at the IFE. We estimate the probability that the infection is cleared when either the resistant strain or the sensitive strain is introduced. Offspring probability generating functions (pgfs) are defined for the infected cells and the viruses. The
It is shown that if
ˆ 0s = (1 − ps ) R
(1 − μ)βcs T¯0 ds
+ ps
(1 − μ)βf s T¯0 ps Ns cs
≤1
and
ˆ 0r = (1 − pr ) R
βcr T¯0 dr
+ pr
βf r T¯0 pr Nr cr
≤ 1,
Table 1 Discrete events in a Markov chain model for free-virus and cell-to-cell transmission with either sensitive and resistant viral strains based on the ODE model (1). FVTi = free-virus transmission of either sensitive strain or resistant strain, CCTi = cell-to-cell transmission of either sensitive strain or resistant strain, i = s, r. Event i
Description
Transition
Rate ri
T→T+1 T→T−1 (T, Is , Vs ) → (T − 1, Is + 1, Vs − 1) (T, Is ) → (T − 1, Is + 1) (T, Ir , Vs ) → (T − 1, Ir + 1, Vs − 1) (T, Ir ) → (T − 1, Ir + 1) (T, Ir , Vr ) → (T − 1, Ir + 1, Vr − 1) (T, Ir ) → (T − 1, Ir + 1) Is → Is − 1 Ir → Ir − 1 Vs → Vs + 1 Vr → Vr + 1 Vs → Vs − 1 Vr → Vr − 1
then the IFE is globally asymptotically stable. In this case R0 < 1. A Lyapunov function is constructed [27]. The proof of this result is ˆ 0i correspond to relegated to the Appendix. Note that the values of R model (1) assuming the terms pi β fi Vi T, i = r, s are removed.
1 2 3
Birth T Death T FVTs
4
CCTs
ˆ 0s } ≤ 1 and 0 < pi < 1, i = r, s, then ˆ 0r , R ˆ 0 = max{R Theorem 3.3. If R the IFE E0 of model (1) is globally asymptotically stable.
5
7
FVTs & mutate CCTs & mutate FVTr
8
CCTr
In addition, if R0 > 1, it is shown that the infection persists. Mathematically, it is shown that either the sensitive strain or the resistant strain is uniformly persistent. The theorem is stated below and the proof is given in the Appendix. Theorem 3.4. Assume R0 > 1 and 0 < pi < 1, i = r, s. (i) If R0s > 1 ≥ R0r , then the sensitive strain in model (1) is uniformly persistent. That is, there exists a positive constant κ > 0 such that lim inft→∞ X (t) > κ for X equal to Is or Vs .
6
9 10 11 12 13 14
Death Is Death Ir Birth Vs Birth Vr Clearance Vs Clearance Vr
dT T ps (1 − μ)β fs Vs T (1 − ps )(1 − μ)β cs Is T ps μβ fs Vs T (1 − ps )μβ cs Is T pr β fr Vr T (1 − pr )β cr Ir T ds Is dr Ir ps Ns ds Is pr Nr dr Ir cs Vs cr Vr
Please cite this article as: L.J.S. Allen, E.J. Schwartz, Free-virus and cell-to-cell transmission in models of equine infectious anemia virus infection, Mathematical Biosciences (2015), http://dx.doi.org/10.1016/j.mbs.2015.04.001
ARTICLE IN PRESS
JID: MBS
[m5G;April 23, 2015;15:3]
L.J.S. Allen, E.J. Schwartz / Mathematical Biosciences 000 (2015) 1–12
general form of an offspring pgf for a random variable Ri from a total of n random variables is
fi (u1 , . . . , un ) =
∞
∞
···
j1 =0
Pi (j1 , . . . , jn )u11 · · · unn , j
j
jn =0
where Pi (j1 , . . . jn ) is the probability that one type i gives “birth” to jk type k, k = 1, . . . , n. The values of ui [0, 1]. For example, fi /uj evaluated at u1 = 1, . . . , un = 1 is the expected number of offspring of type j from one of type i. For our model, Pi is the probability new infected cells or new virions (offspring) are generated from introduction of either one infected cell or one virion. It follows from branching process theory that the fixed points of the pgfs on [0, 1]n approximate the probability of extinction [1,2,10,12,14]. The minimal fixed point of fi (q1 , . . . , qn ) = qi , where qi [0, 1] determines the probability of infection extinction (clearance). In the branching process approximation, sample paths either approach zero or infinity. An important assumption is that zero is an absorbing state for the process. The probability that sample paths for the infected states hit zero yields an estimate for the probability of infection clearance. Given n infected states, Ri , i = 1, . . . , n, with associated pgfs fi , minimal fixed point (q1 , . . . , qn ) and relatively small initial values Ri (0) = ri0, then the branching process estimate for the probability of extinction (Pext ) of the Ri is r20 rn0 Pext = qr10 1 q2 · · · qn .
It is important to note that Pext is an asymptotic estimate based on the linear branching approximation. The time until infection clearance varies and depends on the initial concentrations and the magnitude of R0 . Let T = T¯0 . Then the events 3–14, listed in Table 1, have linear transition rates ri with respect to the random variables Ii and Vi , i = s, r. These transition rates are used to define offspring pgfs [1,2,10]. Let fi (u1 , u2 , u3 , u4 ) for ui [0, 1], i = 1, 2, 3, 4, denote the pgf for each of the four random variables Is , Ir , Vs , and Vr , respectively. Assume (Is (0), Ir (0), Vs (0), Vr (0)) = (δ i1 , δ i2 , δ i3 , δ i4 ), where
δik =
0, i = k 1, i = k
is the Kronecker’s delta symbol (introduction of either one infected cell or one virion). The pgfs fi fi (u1 , u2 , u3 , u4 ) for i = 1, 2, 3, 4 are defined as follows:
f1 =
ds +(1−ps )(1−μ)βcs T¯0 u21 + (1−ps )μβcs T¯0 u1 u2 + ps Ns ds u1 u3 ds + (1 − ps )βcs T¯0 + ps Ns ds
f2 =
dr + (1 − pr )βcr T¯0 u22 + pr Nr dr u2 u4 dr + (1 − pr )βcr T¯0 + pr Nr dr
f3 =
cs + ps (1 − μ)βf s T¯0 u1 + ps μβf s T¯0 u2 cs + ps βf s T¯0
f4 =
(16)
Note fi (1, 1, 1, 1) = 1. For example,
cr cr + pr βf r T¯0
, P4 (0, 1, 0, 0) =
(i) If R0 < 1, then Pext = 1. vs0 (ii) If R0r < 1 < R0s , then Pext = qis0 11 q31 , where q11 and q31 are defined as cs + ps μβf s T¯0 ps (1 − μ)βf s T¯0 1 1 + q11 = and q31 = R0s R0s cs + pss βf s T¯0 cs + ps βf s T¯0 (Eq. (19) in the Appendix). ir0 vs0 vr0 (iii) If min{R0r , R0s } > 1, then Pext = qis0 12 q22 q32 q42 , where qi2 , i = 2, 3, 4 are defined as
q22 =
1 R0r
, q42 =
pr βf r T¯0
cr + pr βf r T¯0
1 R0r
+
cr cr + pr βf r T¯0
and q32 = f3 (u1 , q22 , q32 , q42 ) =
cs + ps (1 − μ)βf s T¯0 u1 + ps μβf s T¯0 q22 cs + ps βf s T¯0
(Eqs. (18) and (20) in the Appendix). The value of q12 is the solution u1 (0, 1) of the quadratic equation f1 (u1 , q22 , q32 , q42 ) = u1 .
5. Numerical examples Parameter values for EIAV for model (1) are summarized in Table 2. For the parameter values in Table 2 the value of R0 ranges from a minimum of 21.6 (ps = pr = 1) to a maximum of 457 (pr = ps = 0). The basic reproduction number R0 > 1 for all parameter values. Large values of R0 correspond to greater cell-to-cell transmission when the values of pr and ps are small. Note that dT = ds = dr and cs = cr . The terms representing viral loss due to entry into target cells have a negligible effect on the value of the basic reproduction number due to the large viral clearance rates, cr and cs . FVTi = free-virus transmission rate of either sensitive strain or resistant strain, CCTi = cell-to-cell transmission rate of either sensitive strain or resistant strain, i = s, r, D = day. The reproduction rate of healthy target cells, , was calculated using = dT T(0), assuming that dT/dt = 0 before infection. The death rate of healthy target cells, dT , was calculated using the estimate for a macrophage lifespan of 21 days [32]. The free virus transmission rate for sensitive virus (β fs ), death rates of infected cells (ds and dr ), burst numbers (Ns and Nr ), and virus clearance rates (cs and cr ) were based upon the estimated values for EIAV infection in SCID horses [26]. Values here have been expressed using virions instead of vRNA copies (given that one virion contains two vRNA copies). Transmission rate β fr was determined given that resistant free virus EIAVPND-5 is 4.5-fold less infectious than sensitive free virus EIAVPND-1 [34]. Transmission rate β cs was determined given that cell-to-cell transmission Table 2 Parameter values for EIAV corresponding to model (1).
cr + pr βf r T¯0 u2 . cr + pr βf r T¯0
P4 (0, 0, 0, 0) =
5
pr βf r T¯0 cr + pr βf r T¯0
with all other probabilities P4 (j1 , j2 , j3 , j4 ) = 0. Details of the calculation for the fixed points of (16) are provided in the Appendix and the results are summarized in the following theorem. Theorem 4.1. Assume T (0) = T¯0 with relatively small initial values for the infected states, (Is (0), Ir (0), Vs (0), Vr (0)) = (is0, ir0, vs0, vr0). Then for the Markov chain model, there are three different cases for the estimate of the probability of infection clearance.
Parameter
Units
Definition
Value
Reference
ds dr Ns
cells/(ml D) D−1 ml/(virions D) ml/(virions D) ml/(cells D) ml/(cells D) (base cycle)−1 D−1 D−1 virions/cell
2019 1/21 6.50 × 10−7 1.44 × 10−7 5.13 × 10−4 5.13 × 10−4 3 × 10−5 1/21 1/21 5302.5
Calculated [32] [26] [34] [34] [34] [19] [26] [26] [26]
Nr
virions/cell
5302.5
[26]
cs
D−1
6.73
[26]
cr
−1
Reprod rate T Death rate T FVTs FVTr CCTs CCTr Mutation rate Death rate Is Death rate Ir Burst number Sens virus Burst number Res virus Clearance rate Sens virus Clearance rate Res virus
6.73
[26]
dT
β fs β fr β cs β cr μ
D
Please cite this article as: L.J.S. Allen, E.J. Schwartz, Free-virus and cell-to-cell transmission in models of equine infectious anemia virus infection, Mathematical Biosciences (2015), http://dx.doi.org/10.1016/j.mbs.2015.04.001
ARTICLE IN PRESS
JID: MBS 6
[m5G;April 23, 2015;15:3]
L.J.S. Allen, E.J. Schwartz / Mathematical Biosciences 000 (2015) 1–12 Table 3 Initial values for EIAV corresponding to model (1). Initial value
Units
Value
T(0) Is (0) Ir (0) Vs (0) Vr (0)
cells/ml cells/ml cells/ml virions/ml virions/ml
42,390 0 0 233 Varies
Existence and stability of the endemic equilibria for the ODE model (1) depend on the relation between G(T) and H(T) (Theorem 3.2). If G(T) < H(T), then the coexistence equilibrium exists and the resistant strain equilibrium is unstable, and if H(T) < G(T), then the coexistence equilibrium does not exist; only the resistant strain equilibrium exists. In Figs. 3 and 4 are graphs of the contours of H and G as a function of T and ps . If ps = pr , then G(T, ps ) < H(T, ps ) (Figs. 1 and 3), but if pr = 0.8ps , then H(T, ps ) < G(T, ps ) (Figs. 2 and 4). If the ratio of cell-to-cell to free-virus transmission is greater in the resistant strain as compared to the sensitive strain such that H(T, ps ) < G(T, ps ), e.g.,
Reference [13,20]
[5,30,31]
of EIAV in vitro is 789 times more efficient than sensitive cell-free virus transmission [34]. Cell-to-cell transmission of sensitive and resistant viruses has similar efficiency [7,34], so we set β cs = β cr . Initial conditions for model (1) are summarized in Table 3. The target cells for EIAV infection are tissue macrophages. This cell population arises from monocytes in the circulation that enter the tissue and differentiate into macrophages. The initial value for healthy target cells T(0) was estimated by multiplying the average monocyte count in the peripheral blood of four horses (measured as 353.25 cells/μl [20]) by the rate of monocyte turnover in one day (due to dying or differentiating into tissue macrophages), which is 0.12 [13]. The value for Vs (0) was estimated from the initial inoculum of infectious virus for an experimental infection of 106 TCID50 [20,30,31] and the horse plasma volume (3 L, since the horses weighed 60 kg at the time of the experiment and plasma volume is 5% of body weight) using the conversion factor between TCID50 and plaque forming units (PFU) of 0.7 [5]. Figs. 1 and 2 illustrate the dynamics of the ODE model in two cases. In Fig. 1, ps = pr = 0.9, the coexistence equilibrium Ers is stable and is reached quickly, after only 5 days. In Fig. 2, ps = 0.9, pr = 0.8ps and H(T) < G(T). Initially both strains coexist, but eventually the sensitive strain is replaced by the resistant strain equilibrium due to mutation. We note that the equilibrium values for virus and cells are at physiological levels in the simulations in Figs. 1 and 2.
1 − ps 1 − 0.8ps > , 0.8ps ps then the resistant strain dominates (Figs. 2 and 4). Graphs of the reproduction number for the sensitive strain R0s are illustrated as a function of ps (pr = ps ) and the parameters β fs , β cs , ds and cs in Figs. 5–8. If ps = pr , then R0 = R0s . In Figs. 5–8, the graph of R0s and the contours are on the left and right, respectively. To eradicate EIAV, R0 should be reduced below one (Theorem 3.3). The value of R0 can be reduced by decreasing the transmission rates β fs or β cs or by increasing the death rate of infected cells ds or the clearance rate of cs . The contours show that the greatest reduction in R0s can be achieved by decreasing β cs (cell-to-cell transmission) and increasing ds (death rate of sensitive cells). These two parameters β cs and ds have a greater effect on R0s than free-virus transmission β fs or the clearance rate of the virus cs . Similar results hold for R0r and parameters β fr , β cr , dr and cr . In Fig. 9, one sample path is plotted for the Markov chain model for initial concentrations of
(T (0), Is (0), Ir (0), Vs (0), Vr (0)) = (T¯0 , 0, 0, 233, 0) with either ps = pr = 0.9 or with ps = 0.9, pr = 0.8ps = 0.72 and parameter values from Table 2. Cell and virion concentrations are per milliliter. Therefore clearance is with respect to 1 ml, as in, for
EIAV Sensitive Strain
EIAV Resistant Strain 5 Log10 of T, Ir, Vr
Log10 of T, Is, Vs
6
4 T I
2
s
V
T Ir
4
Vr
3 2 1
s
0 0
5 Days
0 0
10
5 Days
10
5 Log10 of T, Ir, Vr
Log10 of T, Is, Vs
6
4 T Is
2
Vs 0 0
100
200 300 Days
400
500
T Ir
4
Vr
3 2 1 0 0
100
200 300 Days
400
500
Fig. 1. Graph of a solution of the ODE model (1) is plotted for parameter values in Table 2 and initial values (T (0), Is (0), Ir (0), Vs (0), Vr (0)) = (T¯0 , 0, 0, 233, 0), where R0s = 63.2, R0r = 49.6, ps = pr = 0.9. The coexistence strain equilibrium is stable, Ers = (Tˆrs , Iˆs , Iˆr , Vˆs , Vˆr ) = (670, 4.17 × 104 , 5.8, 1.41 × 106 , 195), G(Tˆrs ) < H(Tˆrs ) = 1. The figures on the top show the dynamics of the sensitive and resistant strains for the time interval [0, 10], whereas the figures on the bottom are for the time interval [0, 500]. The coexistence equilibrium is reached after 5 days.
Please cite this article as: L.J.S. Allen, E.J. Schwartz, Free-virus and cell-to-cell transmission in models of equine infectious anemia virus infection, Mathematical Biosciences (2015), http://dx.doi.org/10.1016/j.mbs.2015.04.001
ARTICLE IN PRESS
JID: MBS
[m5G;April 23, 2015;15:3]
L.J.S. Allen, E.J. Schwartz / Mathematical Biosciences 000 (2015) 1–12
EIAV Sensitive Strain
7
EIAV Resistant Strain
6 Log10 of T, Ir, Vr
Log10 of T, Is, Vs
6 4
T I
2
s
4 T I
2
r
V
V
s
0 0
r
5 Days
0 0
10
5 Days
10
T Is
6 Log10 of T, Ir, Vr
Log10 of T, Is, Vs
6
V
4
s
2
4 T Ir
2
V
r
0 0
100
200 300 Days
400
500
0 0
100
200 300 Days
400
500
Fig. 2. Graph of the solution of the ODE model (1) is plotted for parameter values in Table 2 and initial values (T (0), Is (0), Ir (0), Vs (0), Vr (0)) = (T¯0 , 0, 0, 233, 0), where R0s = 63.2, R0r = 130.4, ps = 0.9, pr = 0.8ps . The resistant strain equilibrium is stable, Er = (T¯r , I¯r , V¯ r ) = (325, 4.21 × 104 , 1.14 × 106 ), H(T¯r ) < G(T¯r ) = 1. The figures on the top show the dynamics of the sensitive and resistant strains for the time interval [0, 10], whereas the figures on the bottom are for the time interval [0, 500]. The resistant strain appears after 5 days but the sensitive strain is present until 500 days.
H(T,ps)
G(T,ps)
1
1
50
p =p r
0.8
50
s
0.8 150
0.6
150
ps
p
s
0.6 250
0.4
250
0.4
100
100
0.2
350
200 300
0
1
2
0.2
300
400
3
0
4
T
350
200
1
2
3 T
4
x 10
400 4 4
x 10
Fig. 3. Contours of G and H, Eqs. (11) and (13), are plotted as functions of T and ps when pr = ps ; G(T, ps ) < H(T, ps ). Coexistence of both strains.
H(T,p )
G(T,p )
s
s
1
1
50
p =0.8p r
0.8
s
50 150
0.8 150
0.6
250
p
s
ps
0.6 250
0.4
0.4
100
100 0.2
350
200 300
0
1
2
3 T
350
200
0.2
300
400
400 0
4 4
x 10
1
2
3 T
4 4
x 10
Fig. 4. Contours of G and H, Eqs. (11) and (13), are plotted as functions of T and ps when pr = 0.8ps ; H(T, ps ) < G(T, ps ). Resistant strain dominates.
Please cite this article as: L.J.S. Allen, E.J. Schwartz, Free-virus and cell-to-cell transmission in models of equine infectious anemia virus infection, Mathematical Biosciences (2015), http://dx.doi.org/10.1016/j.mbs.2015.04.001
ARTICLE IN PRESS
JID: MBS 8
[m5G;April 23, 2015;15:3]
L.J.S. Allen, E.J. Schwartz / Mathematical Biosciences 000 (2015) 1–12
1
500
0.8
300
0.6
150 200
250
R
ps
0s
400
100
50
200
0.4
300 350
100 0.2 0 5
400
0 −6
x 10
βfs
1
0
450
0 0
0.5 ps
1
2
βfs
3
4 −6
x 10
Fig. 5. Graph of R0s and contours as a function of ps and β fs .
1 50 500
0.8 150
0.6 s
300
p
R0s
400
200
250
0.4 100
100 0.2
0 1
350
200 300
5 0.5
0 0
−4
2
x 10 ps
0 0
4
β
cs
βcs
400 −4
x 10
Fig. 6. Graph of R0s and contours as a function of ps and β cs .
1 50 600
100
0.6
200 250
R0s
ps
400
0.8
150
0.4
200
350 0.2 450
0 0.1
0 0.5 0.05 d s
1
500
0 0.04
400 0.06
p
300
ds
0.08
0.1
s
Fig. 7. Graph of R0s and contours as a function of ps and ds .
example, separate chambers in experimentation in vitro. The probability an infection is not cleared is 1 − Pext = 0.57 whereas the probability of infection clearance is Pext = q233 32 = 0.43 (case (iii) in Theorem 4.1). The estimate for Pext is in good agreement with the Markov chain model. The proportion of the sample paths out of 10,000 for which Is + Ir + Vs + Vr hit zero before reaching a size of 10,000 in the Markov chain model is close to 0.43 (e.g., 0.4241 for ps = pr = 0.9 and 0.4316 for pr = 0.9, ps = 0.8pr ). The branching process approximation
accurately predicts clearance but does not predict which infected state is approached when the infection persists. As can be seen in the sample paths in Fig. 9, the sensitive strain is present in both sample paths. After five days, mutation to the resistant strain has only occurred in the sample path in Fig. 9 (bottom) when ps = 0.8pr . The top and bottom graphs in Fig. 9 should be compared with the top graphs in Figs. 1 and 2, respectively, where the corresponding solutions of the ODE model are plotted.
Please cite this article as: L.J.S. Allen, E.J. Schwartz, Free-virus and cell-to-cell transmission in models of equine infectious anemia virus infection, Mathematical Biosciences (2015), http://dx.doi.org/10.1016/j.mbs.2015.04.001
ARTICLE IN PRESS
JID: MBS
[m5G;April 23, 2015;15:3]
L.J.S. Allen, E.J. Schwartz / Mathematical Biosciences 000 (2015) 1–12
9
1 500
50
100 0.8
400
150 200
0.6
250
p
R0s
s
300 200
0.4
300 350
100
0.2
0 1
450
0
20
0.5
400
5
10 0
p
s
0
10 c
15
20
s
cs
Fig. 8. Graph of R0s and contours as a function of ps and cs .
EIAV Sensitive Strain
EIAV Resistant Strain 5 of T, Ir, Vr
ps=pr=0.9
T I
2
Log
Log
10
4
10
of T, Is, Vs
6
s
Vs 0 0
1
2
3
4
4 ps=pr=0.9
3 2 1 0 0
5
T I r
Vr 1
2
Days
3
4
5
3
4
5
Days 5 of T, Ir, Vr
ps=0.9, pr=0.72
T I
2
Log
Log
10
4
10
of T, Is, Vs
6
s
Vs 0 0
1
2
3
4
5
Days
4 ps=0.9, pr=0.72
3 2 1 0 0
T I r
Vr 1
2 Days
Fig. 9. Graphs of two sample paths of the Markov chain model are plotted for parameter values in Table 2 and initial values (T (0), Is (0), Ir (0), Vs (0), Vr (0)) = (T¯0 , 0, 0, 233, 0). In the top two figures, sample paths corresponding to the sensitive strain and the resistant strain are plotted for R0s = 63.2, R0r = 49.6, ps = pr = 0.9 (compare with Fig. 1). In the bottom two figures, sample paths corresponding to the sensitive strain and the resistant strain are plotted for R0s = 63.2, R0r = 130.4, ps = 0.9, pr = 0.8ps (compare with Fig. 2). The probability of infection clearance is approximately Pext = q233 32 = 0.4302.
Some estimates for the probability of extinction Pext from the branching process approximation are plotted in Fig. 10. The value ir0 vs0 vr0 of Pext = qis0 12 q22 q32 q42 depends on the parameter values (Table 2) and on the initial conditions. An expression for the probabilities qij is given in Theorem 4.1 part (iii). A set of four initial conditions is assumed: (is0, ir0, vs0, vr0) = (1, 0, 0, 0), (0, 1, 0, 0), (0, 0, 233, 0) and (0, 0, 0, 233) with either one infected cell/ml or 233 virions/ml of the sensitive strain or the resistant strain, respectively. Here, we assume pr = ps . Therefore, Pext is plotted as a function of ps . As ps increases, cell-to-cell transmission decreases and free-virus transmission increases. For ps = 0, there is exclusive cell-to-cell transmission and no free virus is produced. For ps = 0, the graphs illustrate that exposure of healthy tissue to infected cells, as low as one infected
cell/ml, is almost certain to result in infection with Pext = qi2 ≈ 0, i = 1, 2 whereas exposure to a virion concentration of 233/ml (or any convr0 centration), results in clearance of infection with Pext = qvs0 32 q42 = 1. On the other hand, for ps = 1, there is exclusive free-virus transmission and no cell-to-cell transmission. In this case, both exposure to infected cells or virions at concentrations of 1/ml or 233/ml, respectively, may either result in infection clearance with Pext > 0, or viral persistence with 1 − Pext > 0. As the value of ps increases, the probability of infection clearance decreases with exposure to a fixed concentration of virions (233 virions/ml, no infected cells) whereas the probability of infection clearance increases with exposure to a fixed concentration of infected cells (1 cell/ml, no virions) but at a much lower rate.
Please cite this article as: L.J.S. Allen, E.J. Schwartz, Free-virus and cell-to-cell transmission in models of equine infectious anemia virus infection, Mathematical Biosciences (2015), http://dx.doi.org/10.1016/j.mbs.2015.04.001
ARTICLE IN PRESS
JID: MBS 10
[m5G;April 23, 2015;15:3]
L.J.S. Allen, E.J. Schwartz / Mathematical Biosciences 000 (2015) 1–12
1
0.8
Pext
0.6
q12 q22
q233 0.4
32 233
q42
0.2
0 0
0.2
0.4
ps
0.6
0.8
1
Fig. 10. Probability of extinction Pext is plotted as a function of ps (pr = ps ) and initial values (T (0), Is (0), Ir (0), Vs (0), Vr (0)) = (T¯0 , is0, ir0, vs0, vr0), where Pext = ir0 vs0 vr0 qis0 12 q22 q32 q42 . The four initial conditions are (is0, ir0, vs0, vr0) = (1, 0, 0, 0), (0, 1, 0, 0), (0, 0, 233, 0) and (0, 0, 0, 233) representing either one infected cell/ml or 233 virions/ml from either the sensitive strain or the resistant strain. With exclusive cell-to-cell transmission, ps = pr = 0, healthy tissue exposed to infected cells at concentrations as low as 1 cell/ml results in almost no clearance of infection with Pext = qi2 ≈ 0, i = 1, 2, whereas exposure to virions at any concentration results in vr0 complete clearance with Pext = qvs0 32 q42 = 1. With exclusive free-virus transmission, ps = pr = 1, exposure to infected cells or virions results either in infection clearance with Pext > 0 or viral persistence with 1 − Pext > 0.
6. Conclusions In this study, we presented novel mathematical models of EIAV infection dynamics to evaluate the role of cell-to-cell transmission among two competing viral strains. One strain is sensitive to antibody neutralization and the other strain is neutralization resistant. These models, which compare the effects of the proportion of cell-tocell transmission versus free virus transmission, are among the first to consider multiple viral strains as well as multiple transmission modes together (see e.g., Pourbashash et al. [22]). Our models also include mutation of sensitive virus to the resistant strain. In our analyses of an ODE model and a Markov chain model with a branching process approximation, we calculated R0 and determined the probabilities of infection clearance. We also provided analyses of the IFE and two endemic equilibria, one with the coexistence of both strains and another with existence of only the resistant strain. Our results indicate that whether an EIAV infection evolves to the resistant virus steady state (Er ) or the coexistence steady state (Ers ) depends crucially on the proportion of cell-to-cell transmission versus free virus transmission (i.e., the values of pr and ps ) (Figs. 1–4). The same proportion of cell-to-cell transmission for the sensitive and resistant viral strains (i.e., pr = ps ) leads to coexistence of both strains. If the ratio of cell-to-cell to free virus transmission is greater in the resistant strain as compared to the sensitive strain under specified conditions, however, the resistant strain will eventually dominate. For eradication of both EIAV strains, we found that R0 is more sensitive to changes in cell-to-cell transmission (that is, the cell-to-cell transmission rates β cs and β cr and the death rates of infected cells ds and dr ) than free-virus transmission parameters (that is, the free virus transmission rates β fs and β fr and the virus clearance rates cs and cr ) (Figs. 5–8). Clearance of infection depends on cell-to-cell versus free-virus transmission, i.e. whether pr and ps are small or large, that is, whether cell-to-cell or free-virus transmission dominates, respectively (Figs. 9 and 10). This investigation adds to previous work in within-host virus dynamics by meeting the call by Perelson and Ribeiro for research on models of the role that cell-to-cell transmission may play in virus
evasion of immune responses [21]. Since only some broadly neutralizing antibodies effectively block cell-to-cell transmission of HIV-1 [18], cell-to-cell transmission is a mechanism by which HIV-1 can evade antibody control. As such, this evasion mechanism will need to be taken into account in vaccine development. A few limitations of this work should be mentioned. Our models assume that back mutation from the resistant strain to the sensitive strain is negligible. This is a reasonable assumption, however, because reversion to the sensitive strain would be rare in the continued presence of antibody pressure. Back mutation could easily be added to the model. Another simplifying model assumption is that the cells and viruses are homogeneously mixed. Also, the results shown in our simulations are specific to the values chosen for pr and ps . We also note that some of the parameter values we use are specific for EIAVinfected SCID horses, which lack a competent immune system. Consequently, these parameters are representative of the rates of these processes in the absence of immune responses; rates may be different in the presence of immune responses. In this study we have focused on parameters relevant to EIAV infection, and it will be interesting to see if similar behaviors arise with related viruses, such as HIV-1, which also has multiple strains and uses both modes of transmission. Our results may also be relevant for understanding strain competition in other viruses like HTLV-1 that use only cell-to-cell transmission [4]. Importantly, our results have implications for potential vaccines for EIAV infection. This study predicts that treatment modalities that only block free virus infection may lead to resistance, whereas treatments that can block transmission by cell-to-cell could reduce R0 and the development of resistance. Furthermore, our results suggest that an optimal vaccine would induce immune responses that remove infected cells over responses that improve virus clearance. Future work in this area that would be worthwhile with this model includes using the model to predict the proportion of transmission via free virus infection for the sensitive and resistance strains (i.e., the values of ps and pr ) given known viral loads for sensitive and resistant virus (i.e., Vs and Vr ), as well as using the model to predict transient benefits in the short term or long term. Future work in this area that would be worthwhile with new models include modeling efforts that examine coinfection of cells by both strains, or models of within host infection dynamics of two strains with two modes of transmission with ART or vaccination. Acknowledgments We acknowledge the support of the Mathematical Biosciences Institute as long-term visitors in spring 2012 at which time this research was initiated. We thank an anonymous reviewer for suggestions that improved the original manuscript. Appendix Proof of Theorem 3.1. To verify solutions of (1) are nonnegative, it is only necessary to show that solutions to do not cross the boundary of R5+ . Each of the differential equations evaluated on their respective zero boundaries, T˙ |T=0 , I˙i |Ii =0 and V˙ i |Vi =0 , i = s, r is nonnegative; hence solutions of model (1) are nonnegative. To show solutions are bounded above let Ttot = T + Ir + Is and Vtot = Vr + Vs . Then
T˙ tot ≤ − dT Ttot and
V˙ tot ≤ K T¯0 − cVtot , where T¯0 = /dT , K = ps Ns ds + pr Nr dr and c = min {cr , cs }. By comparison of Ttot with the solution of x˙ = − dT x, x(0) = Ttot (0) ≤ T¯0 and the solution Vtot with the solution of y˙ = K T¯0 − cy, y(0) = Vtot (0) leads to Ttot (t) ≤ x(t) ≤ T¯0 and Vtot (t) ≤ y(t) ≤ max{Vtot (0), K T¯0 /c}.
Please cite this article as: L.J.S. Allen, E.J. Schwartz, Free-virus and cell-to-cell transmission in models of equine infectious anemia virus infection, Mathematical Biosciences (2015), http://dx.doi.org/10.1016/j.mbs.2015.04.001
ARTICLE IN PRESS
JID: MBS
[m5G;April 23, 2015;15:3]
L.J.S. Allen, E.J. Schwartz / Mathematical Biosciences 000 (2015) 1–12
To verify Theorem 3.3, we apply mathematical results of Shuai and van den Driessche to construct a Lyapunov function (Theorem 2.1, p. 515 [27]). Shuai and van den Driessche also state sufficient conditions for global asymptotic stability of the IFE (Theorem 2.2, p. 516 [27]). Unfortunately, these conditions do not apply directly to model (1) but the method of their proof can be used to verify global asymptotic stability. Proof of Theorem 3.3.. Let F = [Fij ] be the matrix defined in (4), x = (Is , Ir , Vs , Vr ) and y = T. Define the vector f(x, y) as
⎤
⎡
(F11 Is + F13 Vs )(T¯0 − T ) ⎥ ⎢ ⎢ (F21 Is + F22 Ir + F23 Vs + F24 Vr )(T¯0 − T ) ⎥ ⎥ f (x, y) = ⎢ ⎥ ⎢ ps βf s Vs T ⎦ ⎣ pr βf r Vr T and matrix Vˆ as
⎡ ⎢ ⎢
0
0
0
dr
0
0
cs
⎥ 0⎥ ⎥. ⎥ 0⎦
−pr Nr dr
0
cr
Vˆ = ⎢ ⎢
⎣ −ps Ns ds 0
0
⎤
ds
Since T (t) ≤ T¯0 for t [0, ) and solutions to model (1) lie in , it follows that f(x, y) 0 for (x, y) ∈ ⊂ R5+ . There are two cases ˆ 0r > R ˆ 0s and (ii) R ˆ0 = R ˆ 0s ≥ R ˆ 0r . In case (i), ˆ0 = R to consider, (i) R ˆ 0s , denote the left eigenvector of Vˆ −1 F corresponding to the ˆ 0r > R R ˆ 0r as ωrT . An eigenvector ωr is ˆ0 = R eigenvalue R
⎡
ωr
⎤ ˆ 0r F21 R ⎢ ˆ 0s ) ⎥ ˆ 0r − R ⎢ F22 (R ⎥ ⎥ =⎢ ⎢ ⎥ ˆ F R ⎣ ⎦ 23 0r ˆ 0s ) ˆ 0r − R F24 (R
ˆ 0s ≥ (or any nonzero constant multiple). Thus, ωr > 0. In case (ii), R ˆ 0r , a left eigenvector ωsT of Vˆ −1 F corresponding to the eigenvalue R ˆ 0s is ˆ0 = R R
⎡
F11
⎤
⎥ ⎢ ⎢ 0 ⎥ ⎥ ⎥ ⎣ F13 ⎦
ωs = ⎢ ⎢
0 (or any nonzero constant multiple), ωs 0. Define a Lyapunov function for model (1) as
11
To verify part (i), assume R0s > 1 ≥ R0r . Let S0 = {(T, Is , Ir , Vs , Vr ) ∈ R5+ |0 ≤ T + Ir + Is ≤ T¯0 , Is > 0, Vs > 0} and S0 = S\S0 . In S0 , either Is = 0 or Vs = 0. Let φ (t, X0 ) be the solution of model (1) with initial condition X0 = (T(0), Is (0), Ir (0), Vs (0), Vr (0)) S0 . That is, if X0 S0 , one of Is (0) or Vs (0) equals zero. The only invariant set in S0 is {E0 }. Equilibrium Er is not feasible since R0r ≤ 1. Let Ws (E0 ) be the stable manifold of E0 in R5+ . To show uniform persistence of the sensitive strain, it suffices to show Ws (E0 ) S0 =
(Theorem 4.3 [11]). Let ε > 0 and define
M =
(1 − ps )(1 − μ)βcs (T¯0 − ) − ds ps (1 − μ)βf s (T¯0 − ) . ps Ns ds −cs − ps βf s T¯0
Choose α sufficiently small such that if ε [0, α ], then Mε has a dominant positive eigenvalue and associated positive eigenvector. This is possible since the spectral abscissa s(M0 ) > 0 iff R0s > 1. In addition, choose α sufficiently small such that T (0) ≥ T¯0 − α > 0 and X0 − E0 α . Suppose φ (t, X0 ) − E0 α for all t 0, that is, φ (t, X0 ) Ws (E0 ). Then
I˙s ≥ [ps (1 − μ)βf s Vs + (1 − ps )(1 − μ)βcs Is ](T¯0 − α) − ds Is V˙ s ≥ ps Ns ds Is − cs Vs − ps βf s Vs T¯0 . The solution Is (t) and Vs (t) can be compared to the solution of the linear system Y = Mε Y, with Y(0) = (Is (0), Vs (0))T . It follows that the solution Is (t) and Vs (t) → as t → , a contradiction to boundedness of solutions. Hence, Ws (E0 ) S0 = . By Theorem 4.3 [11], the conclusion of Theorem 3.4 part (i) follows. To verify part (ii), assume R0r > 1. Let S0 = {(T, Is , Ir , Vs , Vr ) ∈ R5+ |0 ≤ T + Ir + Is ≤ T¯0 , Ir > 0, Vr > 0} and S0 = S\S0 . In S0 , either Ir = 0 or Vr = 0. Let φ (t, X0 ) be the solution of model (1) with initial condition X0 = (T(0), Is (0), Ir (0), Vs (0), Vr (0)) S0 . The only invariant set in S0 is {E0 }. As in the proof of part (i), to show uniform persistence of the resistant strain, it suffices to show Ws (E0 ) S0 = (Theorem 4.3 [11]). Let ε > 0 and define
M =
(1 − pr )βcr (T¯0 − ) − dr
pr βf r (T¯0 − )
pr Nr dr
−cr − pr βf r T¯0
.
Choose α sufficiently small such that if ε [0, α ], then Mε has a positive dominant eigenvalue s(Mε ) > 0 and associated positive eigenvector (since s(M0 ) > 0 iff R0r > 1). In addition, choose α sufficiently small such that T (0) ≥ T¯0 − α and X0 − E0 α . Suppose φ (t, X0 ) − E0 α for all t 0. Then
I˙r ≥ [pr βf r Vr + (1 − pr )βcr Ir ](T¯0 − α) − dr Ir
Li (x, y) = ωiT Vˆ −1 x,
V˙ r ≥ pr Nr dr Ir − cr Vr − pr βf r Vr T¯0 .
where i = r for case (i) or i = s for case (ii). Then as in the proof of Theorem 2.1 in [27], differentiation of Li along solutions of (1) leads to
By comparison with Y = Mε Y, Y(0) = (Ir (0), Vr (0))T , it follows that solutions Ir (t) and Vr (t) → as t → , a contradiction. Hence, Ws (E0 ) S0 = . By Theorem 4.3 [11], the conclusion of Theorem 3.4 part (ii) follows.
ˆ 0 − 1)ωiT x − ωiT Vˆ −1 f (x, y) L i (x, y) = (R ˆ 0 ≤ 1, then Li is a Lyasince dx/dt = (F − Vˆ )x − f (x, y). Therefore, if R punov function for model (1) [27]. To show global asymptotic stability of the IFE, E0 , we need to show that the only invariant set of L i = 0 is ˆ 0r ≤ 1 and L r = 0, the the IFE [17]. In the case i = r, ωr > 0, so that if R ˆ 0r ≤ R ˆ 0s ≤ 1, only invariant set is the IFE, {E0 }. In the case i = s, if R ˆ 0r ≤ 1 then L s = 0 if Is = 0 = Vs or if Vs = 0 and T = T¯0 . Since R (R0r < 1) the only equilibrium solution to model (1) with this restriction (Theorem 3.1) is the IFE. That is, the only invariant set of L s = 0 is the IFE, {E0 }. It follows that the IFE E0 is globally asymptotically stable in [17]. Proof of Theorem 3.4.. Solutions of (1) lie in a compact set in R5+ (Theorem 3.1). Therefore, model (1) is point-dissipative. Let S = {(T, Is , Ir , Vs , Vr ) ∈ R5+ |0 ≤ T + Ir + Is ≤ T¯0 }.
We compute the fixed points for the probability of extinction in the Galton–Watson branching process approximation of the Markov chain model (Table 1). The extinction results for the branching process approximation are summarized in Theorem 4.1. The expectation matrix of the pgfs f1 , . . . , f4 in (16) is M = [fi /uj ] evaluated at ui = 1, i = 1, 2, 3, 4 is ⎡ ⎢ ⎢ ⎢ W −1 ⎢ ⎢ ⎣
(1−ps )βcs T¯0 (2−μ) + ps Ns ds
(1 − ps )μβcs T¯0
ps Ns ds
0
2(1−pr )βcr T¯0 + pr Nr dr
0
ps (1 − μ)βf s T¯0
ps μβf s T¯0
0
0
pr βf r T¯0
0
0
⎤
⎥ pr Nr dr ⎥ ⎥ ⎥ 0 ⎥ ⎦ 0
where W = diag[w1 , w2 , w3 , w4 ] and wi is the denominator of each of the fi , i = 1, 2, 3, 4. It can be easily verified that [W (M − I)]T = F − V
Please cite this article as: L.J.S. Allen, E.J. Schwartz, Free-virus and cell-to-cell transmission in models of equine infectious anemia virus infection, Mathematical Biosciences (2015), http://dx.doi.org/10.1016/j.mbs.2015.04.001
ARTICLE IN PRESS
JID: MBS 12
L.J.S. Allen, E.J. Schwartz / Mathematical Biosciences 000 (2015) 1–12
and I is the 4 × 4 identity matrix [2]. In particular, ρ (M) > 1 (ρ (M) < 1) iff R0 > 1 (R0 < 1). The probability of extinction is one if ρ (M) < 1. If matrix M is irreducible and ρ (M) > 1, then there is a unique fixed point of the pgfs, (q1 , q2 , q2 , q4 ) (0, 1)4 [2]. As can be seen matrix M is reducible; there is not a unique fixed point of the pgfs in (0, 1)4 . In particular, we will show that there are two fixed points, other than (1,1,1,1), in the region (0, 1]4 when ρ (M) > 1 (i.e., R0 > 1). These two fixed points approximate infection clearance when R0 > 1 in the two cases: R0r ≤ 1 < R0s and R0r > 1. Given (Is (0), Ir (0), Vs (0), Vr (0)) = (is0 , ir0 , vs0 , vr0 ) are initially small values, then the probability of infection clearance is approximately (see e.g. [1,2,10,14]): 0 vr0 Pext = qis1 0 qir2 0 qvs 3 q4 .
(17)
Alternately, with probability 1 − Pext , an infection becomes established, which for our model may be either the resistant strain or both sensitive and resistant strains. If R0 ≤ 1, the pgfs (16) have only one fixed point, namely (1, 1, 1, 1). But if R0 > 1, there exist two more fixed points which lie in (0, 1]4 . Denote these two fixed points as (q1i , q2i , q3i , q4i ), i = 1, 2. The pgfs f2 and f4 are functions of u2 and u4 ; therefore, they can be solved independently of f1 and f3 . The functions f2 and f4 have two fixed points: (q21 , q41 ) = (1, 1) if R0 ≤ 1 and a unique fixed point (q22 , q42 ) (0, 1)2 if R0r > 1 given by
q22 =
1
and q42 =
R0r
pr βf r T¯0 cr + pr βf r T¯0
1
+
R0r
cr cr + pr βf r T¯0
.
(18)
These two fixed points of f2 and f4 are used to find the corresponding fixed points of f1 and f3 . For the case, R0r ≤ 1 with fixed point (q2 , q4 ) = (1, 1) of f2 and f4 , the values of the fixed point of f1 and f3 is (q1 , q3 ) = (1, 1) if R0s ≤ 1. However, if R0s > 1, then the values of the fixed point (q1 , q3 ) of f1 and f3 are
q11 =
[m5G;April 23, 2015;15:3]
1
and q31 =
R0s
ps (1 − μ)βf s T¯0 cs + pss βf s T¯0
1 R0s
+
cs + ps μβf s T¯0 . cs + ps βf s T¯0 (19)
For the case R0r > 1 with fixed point values (q22 , q42 ) of f2 and f4 given by (18), the value of q32 depends only on u1 :
q32 = f3 (u1 , q22 , q32 , q42 ) =
cs + ps (1 − μ)βf s T¯0 u1 + ps μβf s T¯0 q22 . cs + ps βf s T¯0 (20)
Note that if u1 = q12 [0, 1], then q32 (0, 1). Substitution of q32 into f1 leads to a quadratic function in u1 : f12 (u1 ) = f1 (u1 , q22 , q32 , q42 ). Since f12 (0) > 0 and f12 (1) < 1, there is only one fixed point of f12 in [0, 1]. In particular, the fixed point q12 (0, 1). Since f11 (u1 ) = f1 (u1 , 1, f3 (u1 , 1), 1) > f12 (u1 ) for u1 (0, 1) and f11 (q11 ) = q11 = 1/R0s , it follows that
0 < q12
