Mathematical analysis of an HIV latent infection model including both virus-to-cell infection and cell-to-cell transmission

ABSTRACT HIV can infect cells via virus-to-cell infection or cell-to-cell viral transmission. These two infection modes may occur in a synergistic way and facilitate viral spread within an infected individual. In this paper, we developed an HIV latent infection model including both modes of transmission and time delays between viral entry and integration or viral production. We analysed the model by defining the basic reproductive number, showing the existence, positivity and boundedness of the solution, and proving the local and global stability of the infection-free and infected steady states. Numerical simulations have been performed to illustrate the theoretical results and evaluate the effects of time delays and fractions of infection leading to latency on the virus dynamics. The estimates of the relative contributions to the HIV latent reservoir and the virus population from the two modes of transmission have also been provided.


Introduction
Human immunodeficiency virus (HIV) is a retrovirus that mainly infects CD4 + T cells. Upon infection, CD4 + T cells can produce new virions, which lead to more cell infection and viral production. Many mathematical models have been developed to study the within-host dynamics of HIV infection [8, 16, 31-36, 41, 42, 53-55] . Most of these models were focused on the virus-to-cell infection. However, virus can also spread by direct cell-to-cell transmission [26, 29, 30, 45-48, 50, 57]. The mechanisms underlying the cellto-cell transmission are not completely clear. Virus may be transferred from infected cells to bystander uninfected cells via the formation of virological synapses [17]. Comparing with the cell-free virus infection, cell-to-cell transmission may have some advantages for viral spread [3,45,60]. For example, virus transmitted via cell-to-cell transmission may have a lower risk of being neutralized by neutralizing antibodies or cleared by cytotoxic T lymphocytes [27]. Moreover, a number of virions can be transferred to a target cell simultaneously and the probability of being completely blocked by antiretroviral therapy is lower. CONTACT  Thus, cell-to-cell spread of HIV may explain the viral load persistence during effective combination therapy [49]. The relative contributions to the viral replication from the two modes of transmission still remain unclear [5,6,18,19,38,59]. Mathematical models have been used to study the dynamics of HIV infection including these two transmission modes [18-22, 24, 25, 39, 58]. A recent review on modelling viral spread can be found in Ref. [12]. Iwami et al. [18] and Komarova et al. [19] studied the relative contribution of the two transmission modes to the spread of HIV. Lai and Zou [24,25] developed models that incorporate two modes of viral spread with and without the logistic target cell growth. They showed that the basic reproductive number determines the infection dynamics and that the logistic growth of target cells can generate a Hopf bifurcation. A within-host viral infection model with the two transmission modes and three distributed delays was analysed in Ref. [58]. A multi-strain model including the two modes was also analysed and shown to exhibit a competitive exclusion principle [39].
Current treatment consisting of several antiretroviral drugs can suppress viral replication to a low level but cannot eradicate the virus. An important reason is that HIV provirus can reside in latently infected CD4 + T cells [7,56]. Latently infected CD4 + T cells live long, are not affected by antiretroviral drugs or immune responses, but can be activated to produce virus by relevant antigens. In the present paper, on the basis of the models by Lai and Zou [24,25] and Alshorman et al. [2], we developed an HIV latent infection model incorporating both the cell-free virus infection and cell-to-cell transmission. Two delays, one the time between viral entry and viral DNA integration into the host cell's DNA and the other the time between viral entry and viral production, have also been included in the model. We analysed the model by deriving the basic reproductive number and showing the local and global stability of the steady states. We also performed numerical simulations to explore the effects of time delays and other parameters on virus dynamics and estimate the relative contributions to the latent reservoir and viral load from the two modes of viral transmission.

Model formulation
We considered both the virus-to-cell infection and cell-to-cell transmission in an HIV latent infection model [34]. Two time delays were included in our model. The first delay (τ 1 ) is the time between viral entry and latent infection, that is, the time between initial viral infection and the integration of viral DNA into the host cell's DNA. The second (τ 2 ) is the time between cell infection and viral production. It is clear that τ 1 < τ 2 according to the viral life cycle. As assumed in HIV latent infection models [2,34,43,44], a very small fraction of CD4 + T cell infection leads to HIV latency. They don't produce new virus unless activated by antigens. The model can be described by the following delayed differential equations: where T(t) denotes the concentration of uninfected CD4 + T cells at time t, L(t) represents the concentration of latently infected T cells at time t, I(t) is the concentration of productively infected T cells, and V(t) denotes the concentration of virions in plasma. The model assumes that uninfected CD4 + T cells are produced at a rate s, infected by free virus at a rate β or by direct cell-to-cell transmission at a rate k. Parameter d T is the per capita death rate of uninfected CD4 + T cells. The constants f , η ∈ (0, 1) are the proportions of infection that lead to latency. Latently infected cells die at a rate δ L per cell and productively infected cells die at a rate δ per cell. Latently infected cells can be activated by their relevant antigens to become productively infected cells at a rate α. Proliferation of latently infected cells is not considered here but can be added to the model as in [43]. N is the viral burst size, that is, the total number of virions released by one infected cell during its lifespan, and c is the viral clearance rate. The parameter δ 1 is the death rate of infected cells in which viral DNA has not integrated into the host cell's DNA. Thus, e −δ 1 τ 1 and e −δ 1 τ 2 represent the probability of an infected cell surviving to the age of τ 1 and τ 2 , respectively. We define the following Banach space [13,23]: where l is a positive constant and the norm in the space is The set C + = C((−∞, 0], R + ) is the non-negative cone of C. For a given function The initial conditions are given by By the fundamental theory of functional differential equations [14,23], we can obtain the existence and uniqueness of the solution for t > 0.

Positivity and boundedness of solution
The following result shows that a solution (T(t), L(t), I(t), V(t)) of system (1) with the initial condition (2) remains non-negative and bounded ultimately. Theorem 2.1: Let (T(t), L(t), I(t), V(t)) be a solution of system (1) with the initial condition (2). It is positive and ultimately bounded for t > 0.
Proof: We first show that T(t) > 0 for all t > 0. Assume that there exists a t 1 > 0 such that . Thus,Ṫ(t 1 ) ≤ 0. From the first equation of (1), we havė T(t 1 ) = s > 0, which is a contradiction. This implies that T(t) > 0 for all t > 0. By the last three equations of system (1), we have , then we havė L(t 2 ) ≤ 0. On the other hand, from the second equation of (1), we havė This leads to a contradiction. If I(t 2 ) = 0, I(t) > 0 for t ∈ [0, t 2 ) and L(t) > 0, V(t) > 0, t ∈ [0, t 2 ], we also havė I(t 2 ) ≤ 0. However, from the third equation, we havė This also leads to a contradiction.
Similarly, we know that Next, we prove the boundedness of the solution of system (1) with the initial condition (2). From the positivity of the solution and the first equation of (1), we obtain which yields Let We obtain This implies that Thus, By the last equation of V(t), we get Thus, Therefore, T(t), L(t), I(t) and V(t) are ultimately uniformly bounded. Moreover, if we define where a,b are given in Equation (3), then it follows from Theorem 2.1 that the region is a positive invariant with respect to system (1).

The basic reproductive number and equilibria
System (1) is an autonomous differential equation system with fixed time delays. It always has an infection-free steady-state E 0 = (T 0 , 0, 0, 0), where T 0 = s/d T . Inspired by the method in Diekmann et al. [10] and van den Driessche and Watmough [52], we consider the infection and viral production, and define matrices F and V as Thus, the basic reproductive number, R 0 , can be defined as the spectral radius of the next generation operator FV −1 , where Therefore, For a within-host virus dynamic model, the basic reproduction number is the number of virions (or infected cells) produced by one virion (or one infected cell) in its lifespan in a fully susceptible environment (i.e. all the target cells are susceptible to virus infection). In the above expression of R 0 , Nβs/(cd T ) is the basic reproduction number via the virus-to-cell infection and ks/(δd T ) is the basic reproduction number via the cell-to-cell transmission. The first term in each parenthesis is the contribution of latent infection and the second term is the contribution of productive infection. Thus, R 0 is the basic reproduction number of system (1). R 01 and R 02 represent the contribution to R 0 from the virus-to-cell infection and cell-to-cell transmission, respectively. We also notice that R 0 derived here is the same as that obtained from the threshold condition guaranteeing the existence of the infected steady state.

Local stability
In this section, we study the local stability of the infection-free and the infected steady states.

Theorem 2.2:
The infection-free steady-state E 0 of system (1) is locally asymptotically stable when R 0 < 1 and unstable when R 0 > 1.

Proof:
We consider the linearization of system (1), obtain the Jacobian matrix J(E 0 ) at the infection-free equilibrium E 0 , and calculate the characteristic equation as follows: where λ is the eigenvalue and E is the identity matrix with the same dimension as J(E 0 ). Expanding the determinant, we obtain λ + d T = 0 or There is always a negative root λ = −d T . Equation (6) can be rewritten as the following form: Some coefficients of the the above characteristic polynomial are delay-dependent. Beretta and Kuang [4] used a geometric method to explore the possible stability switch in a general delay differential equation system with delay-dependent parameters. However, for our specific delay differential system, we will show the stability of steady states by direct comparison of the modulus of the characteristic equation. When R 0 < 1, we aim to show that if the eigenvalue λ = x + iy is a solution of Equation (7), then the real part x < 0. We prove by contradiction. Suppose that x ≥ 0. Equation (7) becomes the following equation by The modulus of the right-hand side of Equation (8) satisfies This leads to a contradiction because R 0 < 1. Thus, all the roots of the characteristic equation (7) have negative real parts. Therefore, E 0 is locally asymptotically stable when R 0 < 1. When R 0 > 1, we rewrite the characteristic equation (7) as where Thus, there exists at least one positive real root such that F(λ) = 0, which implies that the infection-free steady-state E 0 is unstable when R 0 > 1.
The following result addresses the local stability of the infected steady-state E * when R 0 > 1 and f = η (i.e. both modes of transmission have the same probability of generating latent infection). When f = η, the analysis is challenging and remains to be further investigated.

Theorem 2.3:
In the case of f = η, the infected steady-state E * of system (1) is locally asymptotically stable when R 0 > 1.

Proof:
Similar to what we did for the infection-free steady state, we have the Jacobian matrix at E * and calculate the characteristic equation as follows: where λ is the eigenvalue. T * , I * and V * are given in Equation (4).
Equation (10) can be rewritten as We will show that all the roots of Equation (11) have negative real parts. Suppose that λ has a non-negative real part. Using Equation (11) becomes Supposing the right side is , we have On the other hand, the modulus of the left-hand side of Equation (12) satisfies This leads to a contradiction. Thus, the infected steady-state E * is locally asymptotically stable when R 0 > 1. (1) In this section, we will show that system (1) is persistent when R 0 > 1.

Uniform persistence of system
Proof: We use the persistence theory developed by Hale and Waltman [15] for the proof. Define Let (t) be the solution semiflow of system (1) with (2) for all t > 0. It is easy to verify that Notice that the boundedness of each component does not depend on the initial condition (2). Thus, for any bounded set Y in X, the positive orbit γ In view of this property, (t) is asymptotically smooth, that is, for any non- Let A 0 be the global attractor of (t) restricted to ∂X 0 . We have {E 0 } is a compact and isolated invariant set. Thus, the covering is simply {E 0 }, which is acyclic because no orbit connects E 0 to itself in ∂X 0 .
Next, we will verify that W s (E 0 ) ∩ X 0 = ∅. That is, for the above small enough ε 1 , we have To this end, we suppose the opposite. We have lim sup t→∞ (T(t), L(t), I(t), V(t)) − E 0 < ε 1 , which implies that there exists a t 0 > 0 such that Consider the following linear comparison system for (L(t), I(t), V(t)) The spectral radius ρ(F ε 1 V −1 ) > 1 implies that the trivial solution of the linear system (16) is unstable. This, together with (15) and the comparison theorem, implies that L(t), I(t) and V(t) are unbounded, which leads to a contradiction with the boundedness of the solution. Thus, we have W s (E 0 ) ∩ X 0 = ∅. By Theorem 4.2 of Hale and Waltman [15], we know that there exists a value p > 0 such that which means that each component of the solution with the initial condition (2) satisfies

Global stability of system (1) with f = η
In this section, we construct Lyapunov functionals to study the global dynamics of the steady states of system (1). The following results show that the basic reproductive number again provides a threshold value determining the global stability. Similar to the local stability of the infected steady state, we assume f = η. The case of f = η remains to be further studied.

Theorem 4.1:
In the case of f = η, the infection-free steady state E 0 of system (1) is globally asymptotically stable when R 0 < 1.

Proof: We define a Lyapunov functional
and The Lyapunov functional W(·) is non-negative. At the infection-free steady-state E 0 , we have W(E 0 ) = 0. Calculating the derivative of W 1 (t) and W 2 (t) along the solution of system (1), we obtain and Here Thus, Recalling s = d T R 0 and R 0 = R 01 + R 02 < 1, it follows from Equations (19) and (21) that We also have Therefore, from Equation (22), we obtain Thus, dW(t)/dt < 0 when R 0 < 1. Moreover, dW(t)/dt = 0 if and only if V(t) = 0 and T(t) = T 0 . The largest invariant set in {(T, L, I, V) : dW(t)/dt = 0} is the singleton set {E 0 }. Therefore, by Lyapunov-LaSalle invariance principle [14,23] and Theorem 2.2, the infection-free steady-state E 0 is globally asymptotically stable. The next theorem addresses the global stability of the infected steady state when R 0 > 1.

Theorem 4.2:
In the case of f = η, the infected steady-state E * of system (1) is globally asymptotically stable when R 0 > 1.

Proof: We define another Lyapunov functional
and The Lyapunov functional U(·) is non-negative and at the infected steady-state E * it vanishes.
We calculate the time derivatives of U 1 (t), U 2 (t) and U 3 (t) along the solution of system (1).
From the steady state, we have The derivative of U 1 (t) becomes In view of d dt we have the derivatives of U 2 (t) and U 3 (t) as and Using Equations (29)-(31), we obtain From the infected equilibrium E * and the definition of R 0 , we have The function is non-negative for any x > 0 and g(x) = 0 if and only if x = 1. Using this property, we obtain The largest invariant set in {(T, L, I, V) : dU/dt = 0} is a singleton set {E * }. Thus, by the same arguments as the infection-free steady state, the infected steady-state E * is globally asymptotically stable when R 0 > 1.

Numerical results
In this section, we performed some numerical simulations to illustrate the stability results, examine the effect of time delays and the fractions of infection that lead to latency on the virus dynamics, and also explore the relative contributions from the two modes of viral transmission to the reservoir of latently infected cells and the viral load.

Numerical illustration of stability results
Parameter values were chosen from several modelling papers [2,39,[42][43][44]. 0.01 day −1 , δ L = 0.004 day −1 , N = 2000 per cell per day, and c = 23 day −1 . The death rate of infected cells that have not started to produce virus is between d T and δ, and was assumed to be δ 1 = 0.05 day −1 [2]. The time delay τ 1 between viral entry and viral DNA integration was chosen to be 0.25 days and the delay τ 2 between viral entry and viral production was chosen to be 0.5 days [40]. Different values of τ 1 and τ 2 were used to evaluate the effect of time delays on the virus dynamics. A very small fraction of infection events lead to latent infection. Thus, we assumed that f = η = 0.001 as in [43]. Different values of f and η and the case of f = η have also been used in simulation to show the sensitivity of the model prediction to these parameters. Using these parameters, the basic reproductive number was calculated to be R 0 = 3.0 > 1. By Theorems 2.3 and 4.2, the infected steady-state E * = (3.3 × 10 5 , 471, 6.5 × 10 3 , 5.7 × 10 5 ) is locally and also globally asymptotically stable (Figure 1). Under treatment that reduces the infection rates, we assumed the infection rates to be β = 2.4 × 10 −10 cells ml −1 day −1 and k = 8 × 10 −7 cells ml −1 day −1 . The effective reproductive number under treatment becomes R 0 = 0.8 < 1. It follows from Theorems 2.2 and 4.1 that the infection-free steady-state E 0 = (1.0 × 10 6 , 0, 0, 0) is locally and globally asymptotically stable. Numerical simulation using the steady state obtained from Figure 1 as the initial condition was shown in Figure 2.

Effects of time delays and fractions of infection leading to latency
To study the influence of time delays τ 1 , τ 2 on the dynamics of the latent reservoir and the viral load, we increased τ 1 from 0.1 days to 0.3 days but fixed τ 2 = 0.5 days (notice that τ 2 > τ 1 ) and other parameter values in the upper panel of Figure 3. As τ 1 increases, the time to reach the peak levels of the latent reservoir and the viral load also increases but the magnitude of the peak levels decreases. There is almost no change in the steady states. In the lower panel of Figure 3, we evaluated the influence of the fraction of infection that leads to latency on the dynamics. We fixed the fraction leading to latency via the cell-to-cell transmission (η = 0.005) but increased the latency fraction via the virus-to-cell infection (f ) from 0.001 to 0.005. As f increases, the latent reservoir size also increases. However, there is almost no change in the viral load dynamics. This shows that the contribution to the viral replication from latent infection is relatively small.
In Figure 4, we performed the simulation using different values of τ 2 and η. In the upper panel, we fixed τ 1 = 0.25 days but varied τ 2 between 0.5 and 1 days. We found no  difference in the latent reservoir and the viral load as τ 2 changes. In the lower panel of Figure 4, we increased the fraction that leads to latency via the cell-to-cell transmission (η = 0.001, 0.003, and 0.005) but fixed f = 0.001. Similar to the prediction in Figure 3, the latent reservoir increases as η increases but the vial load change is not affected. Figures 3 and 4 showed the sensitivity of modelling prediction to the time delay and latency fraction parameters when the basic reproductive number is greater than 1. In parallel, we conducted simulations with R 0 < 1 using the reduced infection rates β and k as in Figure 2. Figures 5 and 6 showed the simulations with R 0 < 1 using different time delays (τ 1 and τ 2 ) and different fractions that lead to latent infection (f and η). We found that there is almost no difference in the dynamics of the latent reservoir and the viral load. Therefore, although the time delay and the fraction leading to latency may affect the dynamics during primary infection before treatment, they have no influence on the decay dynamics of the latent reservoir and the viral load under suppressive therapy.

Relative contribution from the two modes of viral transmission
A critical question is the relative contribution to the latent reservoir and the viral load from the two modes of viral transmission. We used the following fraction C L to quantify  Figure 2. The other parameters are the same as those in Figure 3.
the relative contribution to the latent reservoir from the virus-to-cell and the cell-to-cell viral transmission, and C V to quantify the relative contribution to the viral load from the two transmission modes: and In the expression of C V , the second terms in the numerator and denominator represent the contribution to viral production from the activation of latently infected cells via the virus-to-cell and cell-to-cell transmission, respectively. Because of the very small fractions (f and η) of infection that lead to latency and the quasi-steady-state assumption (the viral load is proportional to the level of productively infected cells, that is, V(t) = NδI(t)/c), we simplified C V as Under the quasi-steady-state assumption, we also had Thus, Recall that are the contribution to the basic reproductive number R 0 from the virus-to-cell and cellto-cell transmission, respectively. We called R 01 the basic reproductive number from the virus-to-cell infection and R 02 the basic reproductive number from the cell-to-cell viral transmission. Because the fractions f and η are very small, we had Therefore, we obtained the estimates of the relative contributions to the viral load and the latent reservoir from the two transmission modes In Figure 7, we plotted the relative contribution to the latent reservoir and the viral load according to Equations (34) and (35), respectively. We found that the relative contributions agree with the estimates in Equation (36) after a very short time. From the estimates (36), the relative contribution to the viral load from the two transmission modes is the same as the ratio of the two basic reproductive numbers R 01 and R 02 . The fractions of infection that lead to latency have a negligible effect on the relative contribution to the viral load. However, the relative contribution to the latent reservoir is highly dependent on the fractions that lead to HIV latency via these two modes of viral transmission.

Conclusion and discussion
Both in vitro and in vivo experiments have showed that HIV can infect cells via two types of strategies: one is the virus-to-cell infection and the other is the cell-to-cell viral transmission [28,45]. Many existing mathematical models of HIV infection studied virus dynamics with only the virus-to-cell infection. In this paper, we developed a model including both two modes of virus transmission and studied their contributions to the latent reservoir and the viral load. Time delays between viral entry and the establishment of latent infection or virus production were also included in the model. We obtained the basic reproductive number for the model, which consists of two parts: one is the contribution from the virusto-cell infection and the other is the contribution from the cell-to-cell transmission. Similar to many within-host or between-host virus dynamic models, the basic reproductive number provides a threshold value for determining when the virus is predicted to die out and when the chronic infection is established.
Numerical simulations showed that time delays between viral entry and DNA integration or viral production don't have a significant impact on the virus and the latent reservoir dynamics. Although the fractions of infection that lead to viral latency via the two transmission modes can affect the latent reservoir size during primary and chronic infection, they have a very minor effect on the decay dynamics of the virus and the latent reservoir in patients under suppressive therapy. We also obtained the estimates for the relative contributions from the two transmission modes to the latent reservoir and the viral load. The relative contribution to the viral load is approximately the ratio of the basic reproductive numbers via the two transmission modes (i.e. R 01 /R 02 ). The fractions of infection leading to latency have a minor effect on this relative contribution. However, the relative contribution to the latent reservoir depends heavily on these fractions of latency, as well as the ratio of the two basic reproductive numbers.
Our analysis provides a theoretical estimate of the relative contribution to the virus population from the two transmission modes. However, the relative contribution remains unknown because of unknown parameter values. Different estimates of the relative contribution from the two modes to the viral replication have been obtained in the references [18,19]. Komarova et al. [19] estimated that the two transmission modes contributed approximately equally to the viral population . Sourisseau et al. [50] found that mildly shaking the cell culture can prevent the cell-to-cell infection. Iwami et al. [18] used this experimental method, combined with mathematical model to fit the data, to estimate that the cell-to-cell infection may account for approximately 60% of viral infection. The cell-to-cell transmission was also predicted to shorten the generation time of virus. Whether the shaking condition of cell culture can completely prevent cell-to-cell infection and whether the prediction from the in vitro cell culture can be applied to in vivo studies are unclear. More comparison of these two transmission modes remain to be further investigated.
Antiretroviral therapy was assumed to reduce the infection rates of both transmission strategies. However, the model is sensitive to the treatment effectiveness. Once the drug efficacy is beyond a threshold value (such that the basic reproductive number R 0 is less than 1), the infection is predicted to die out. The cell-to-cell transmission can transfer multiple virions to the target cell simultaneously. Thus, antiretroviral drugs may not be effective in blocking the infection of all virions and thus are less effective in inhibiting the cell-to-cell transmission than the virus-to-cell infection [49]. However, some other studies have shown that the cell-to-cell spread of HIV between cells can be effectively blocked by highly active antiretroviral drugs [1,28,37,51]. Extending our model by including multiple infection and considering how the cell-to-cell transmission affects the efficacy of drugs in interfering with the viral replication cycle [11,20] will improve the understanding of the dynamics of the latent reservoir and the viral load in patients receiving long-term combination therapy.

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