Modeling the dynamics of viral–host interaction during treatment of productively infected cells and free virus involving total immune response

Preeti Dubey, Uma S. Dubey, Balram Dubeyc,1 The Program for Experimental & Theoretical Modeling, Department of Medicine, Loyola University Chicago, Illinois, USA preetiup28@gmail.com Department of Biological Sciences, BITS Pilani, Pilani Campus, Pilani, Rajasthan, India uma@pilani.bits-pilani.ac.in Department of Mathematics, BITS Pilani, Pilani Campus, Pilani, Rajasthan, India bdubey@pilani.bits-pilani.ac.in


Onset of viral infection
Diseases caused by viral infections have had a major impact on populations (e.g. hepatitis B virus, HIV/AIDS). Viruses may be present in the body either in a free state or as an intracellular parasite. Free viruses enter a cell in a receptor specific manner and infect the cell. Viral replication is possible only within a cell. The virus replicates within the cell by utilizing the cells own machinery. Inside the cell, a virus may either exist in a latent (inactive) form for a prolonged period or it may immediately adopt the host replication machinery and start producing multiple copies of itself. Once large number of virus particles have been produced inside the cell, they come out of the cell by destroying it. Now, these viruses are free to infect the other healthy cells. Upon infection, virus produces infected cells; all of these infected cells do not contribute in virus production, and some of them remains in latent stage and can produce virus once reactivated. Within the body, virus encounters an immune response, which prevents its spread from an infected cell to adjacent uninfected cells. The present model is being proposed to understand the dynamics of interaction between uninfected cells, productively infected cells, latently infected cells, free virus, and the effect of immune response with treatment of productively infected cells and free virus.

Mathematical models for HIV infection
In the literature, numerous mathematical models have been developed to study the dynamics of HIV infection [1,4,9,17,18,[24][25][26]29,41,42]. The basic viral dynamics model has been studied in detail by researchers [1,24]: In the above model, x(t) is the number of uninfected cells, y(t) is the number of infected cells, and v(t) is the number of free virus at any time t 0. Uninfected cells grow constantly at the rate λ and die out with the rate δ 0 x. When free virus interacts with uninfected cells, it infects uninfected cells with infection rate constant α, and δ 1 y denotes the loss rate of infected cells. Virus is produced at the rate ky via lysis of infected cells, and k 0 is clearance rate constant of virus.

Inclusion of effect of immune response
Immune response plays an important role to control viral infection as it fights against the virus and reduces the virus load. Nowak and Bangham [24] extended the basic virus dynamics model involving CTL (Cytotoxic T lymphocyte) component of immune response, which suggests that an active CTL immune response may reduce virus load. Vargas-De-León and Esteva [38] studied global stability of virus dynamics model with linear immune response. Further, Vargas-De-León [37] proposed a virus dynamics model with lytic and nonlytic immune responses. The effect of nonlinear immune attack rates has been explored to understand the global dynamics of infection. It was shown that if the basic reproductive number R 0 > 1, then the virus will spread and persist within its cellular population, whereas if R 0 < 1, then the virus will not spread. They have also shown the importance of immune responses to control infection. Wang and Liu [41] proposed a class of delayed viral models with saturation infection rate and immune response. In their models, they considered the CTL response (which kills infected cells) and antibody response (which facilitates removal of viruses) separately to study the analytical behavior of systems. Roy et al. [29] studied the effect of CTL immune response on the dynamics of infected CD4 + T cells, virus producing CD4 + T cells, and virus by introducing a positive feedback parameters in the model. They also investigated an optimal control therapy using reverse transcriptase inhibitors (RTIs) that block new infection. Further, Dubey et al. [6] generalized the discussed model (1) by considering the effect of appropriate immune response on the infected cells and suggested that under some sufficient conditions, the number of uninfected cells can be maintained at an appropriate level. Hattaf et al. [13] studied the spread and transmission of HIV infection via: (i) virus-to-cell infection, and (ii) cell-to-cell transmission together with the adaptive immune response, which provides better understanding of HIV infection.

Inclusion of awareness and treatment
The spread of infectious diseases, especially, virus mediated may be controlled by providing appropriate awareness about the infection through media [21,22,30] and by treatment with appropriate therapeutic drugs [5]. Researchers have paid attention to drug therapy or treatment (cure) of targeted cells in virus dynamics models [11, 14-16, 20, 33-35, 40, 43]. Liu et al. [20] developed an HIV pathogenesis dynamics model considering cure rate. They incorporated the full logistic proliferation term for uninfected cells as well as infected cells in the model. In this study, they obtained a critical number (the smallest virus number released by per infected CD4 + T cells) and have shown that this critical number increases with increase in cure rate and observed that the HIV infection can be controlled by increasing cure rate.
Hattaf et al. [12,14] proposed and analyzed a virus dynamics model with general incidence rate and linear cure rate of the infected cells to uninfected cells. They have shown that the virus can be cleared and the disease dies out if the basic reproduction number is less than one. The model and results in [14] were further extended by Tian and Liu [35]. Gumel and Moghadas [11] studied the role of antiretroviral therapy in controlling the HIV infection. They investigated the immunological and therapeutic control of HIV and found the optimal level of antiretroviral therapy to eradicate HIV. Recently, Hill et al. [16] summarized that viral dynamics models have provided many important insights into dynamics of infection and have been extensively used to characterize the response of HIV to antiretroviral therapy. They addressed the important questions about the treatment and immunotherapy. For example, what is the efficacy of antiretroviral drugs, and how does antiretroviral efficacy influence treatment outcomes?
It is noticed that the above models focused on the cure of infected cells, while the cure or killing (via therapeutic drugs) of free virus is equally important to reduce the count of infected cells resulting in reduction of infection.
Considering the above points in view, we propose a five-dimensional virus dynamics model with saturated infection rate and immune response with treatment of infected cells as well as free virus. The main aim of this paper is to study the control of the replication of infected cells and free virus by considering immune response (innate, cellular and humoral) and appropriate treatment for both; the infected cells and free virus.

The mathematical model
It is reported in the virus dynamics models (see [10,15,16,23,25,27,28,39]) that upon infection, a fraction of infected cells may be in latent state and thus do not take part in producing virus. However, latently infected cells become active upon stimulation and then join the other group of actively virus producing infected cells. In order to understand the interaction between uninfected cells (x(t)), productively infected cells (y 1 (t)), latently infected cells (y 2 (t)), virus (v(t)), and total immune response (z(t)) during treatment of both (productively infected cells and virus), we propose a mathematical model consisting of 5 ordinary differential equations x(0) > 0, y 1 (0) 0, y 2 (0) 0, v(0) 0, z(0) > 0. In most of the previous work, it has been observed that immune response increases either due to infected cells (CTL immune response) [24,37,38] or virus (humoral immunity [18]). In this study, we assume that the immune response increases due to infected cells as well as virus (as considered in our previous study [4,7,8]). Keeping the above aspects in view, the virus dynamics can be governed by system (2) of five ordinary differential equations. The interaction of subpopulations can be visualized from schematic diagram of model (2) (Fig. 1). Here x(t), y 1 (t), y 2 (t), v(t), and z(t) represent the uninfected cells, productively infected cells, latently infected cells, free virus, and total immune response, respectively. Uninfected cells grow from the sources inside the body (thymus) at rate constant λ and die out at rate δ 0 x. As shown in our previous work [4], the infection rate f (t) = αx(t)v(t) increases linearly with x or with v, which may not be realistic in real life. If we take the saturated incidence rate then it gives a rich dynamics of the system [4,18,32,39]. Thus, we considered that free virus infects uninfected cells at saturated infection rate, which is given by Holling type II  Figure 1. Schematic diagram of model (2). A mathematical model of viral dynamics augmented to include the total immune response (IR) during treatment of productively infected cells and free virus. Uninfected cells x, which proliferate at influx rate λ and die at rate δ 0 , get infected at rate α (dotted curved arrow) and divided into two types of infected cells: productively infected cells y 1 with fraction q and latently infected cells y 2 with fraction 1 − q. Both type of infected cells loss at rate δ 1 , and latently infected cells get activated at rate β. Productively infected cells produce virions v at rate k. The total immune response z proliferates with innate immune response at rate µ and stimulates with CTL-mediated immune response (dashed elbow connector) at rate µ 1 via productively infected cells and with humoral immune response (dashed elbow connector) at rate µ 2 via virus and depletes at rate µ 0 . Productively infected cells and virus are cured or killed via drug at rates θ 1 and θ 2 , respectively, and are killed by CTL-mediated and humoral immune response at the rates γ 1 and γ 2 (dotted curved lines), respectively. Virus is cleared at rate k 0 . functional response (αxv/(1 + α 1 v), where α and α 1 are positive constants) [4,18]. We also considered two states of infected cells [23,39]: virus producing state (y 1 ), in this state, infected cells produce new virions; latently infected state (y 2 ), in this state, infected cells do not contribute in producing new virions, but can be reactivated to produce new virions. The probability that upon infection, a cell enters to the virus producing state is given by q, and 1 − q is the probability that upon infection, a cell enters to the latently infected state, and δ 1 notates the loss rate constant of these infected cells. The latently infected cells are getting reactivated at the rate constant β and join the virus producing infected cells to contribute in production of virus. Free virions are produced via lysis of productively infected cells at the rate constant k and get cleared from the body with the rate constant k 0 . We assumed that µ is the innate immune response of the body. When virus enters into the body and targets the uninfected cells to infect them, then the infected cell-specific lymphocytes (cellular immunity) proliferate at rate µ 1 y 1 z, and the virusspecific lymphocytes (humoral immunity) proliferate at rate µ 2 vz. The corresponding decrease in the number of productively infected cells and virus is given by rates γ 1 y 1 z and γ 2 vz, respectively. The immune response decays at the rate µ 0 z due to natural factors. It is assumed that the productively infected cells and free virus are killed or removed with the rates θ 1 y 1 and θ 2 v due to therapeutic drugs. In modeling, it has been assumed that drug acts on the productively infected cells by killing these cells or inducing these https://www.journals.vu.lt/nonlinear-analysis cells to suicide by programmed cell death (PCD) or apoptosis. They do not re-enter the susceptible category as they are eliminated from the system. Total removal of virus from within cell to make it normal is not biologically possible, thus can not enter the susceptible class.
In the next section, we show that all the solutions of system (2) are positively invariant and bounded.
3 Boundedness and positivity of solutions of model (2) (2) (2) In order to see the well-posedness of the model, from the first equation of model (2) we Similarly, from other four equations we obtain Thus, we state the following Lemma.
Now, in order to see the boundedness and invariance of system (2), from the first equation of model (2) we haveẋ λ − δ 0 x. Using elementary calculus, we get that lim sup t→∞ Differentiating L along all the solutions of system (2), we getL This shows that solutions of system (2) point towards the region Ω defined in Lemma (1). Hence Ω is positively invariant, and solutions of system (2) are bounded. Thus, we can state the following lemma.
This proves that the model is biologically well behaved. In the next section, first we find equilibrium points of system (2). Then determine the basic reproduction number of system (2). Afterwards, the local and global stability of the equilibrium points has been studied.
4 Stability analysis of model (2) (2) (2) From model (2) it is clear that the system has two nonnegative equilibria: VFE is trivial and given by E 0 (x 0 , 0, 0, 0, z 0 ) = E 0 (λ/δ 0 , 0, 0, 0, µ/µ 0 ). https://www.journals.vu.lt/nonlinear-analysis We use the next generation operator method [3,36] to determine the basic reproduction number. We are considering y 1 , y 2 , and v to be the infection compartments. Let X = (y 1 , y 2 , v) T , F 1 (X) is the nonnegative infection matrix (gain in infection terms), and V 1 (X) (the matrix for transfer terms between compartments) associated with system (2). Then we can rewrite system (2) asẊ The Jacobian of matrices F 1 (X) and V 1 (X) evaluated at VFE; E 0 (λ/δ 0 , 0, 0, 0, µ/µ 0 ) is given by Then ρ(F 1 V −1 1 ) gives the spectral radius (largest eigenvalue) of the next generation matrix (F 1 V −1 1 ) as defined in [36]. Thus, where R I is the basic reproduction number in the presence of immune response. This can be defined as the average number of newly infected cells produced by single infected cell when introduced into a completely susceptible cells. It is apparent from the expression of the basic reproductive ratio (R I ) that it decreases with the increase in therapeutic drug (for both infected cells and virus) as well as with the improvement of immune response. In the present model, we consider immune response, so the given basic reproduction number is determined in the presence of immune response. Further, we can state the following theorem using the above results and Theorem 2 of [36].
Local stability of virus-free equilibrium shows that the infection will die out if basic reproduction number is less than one and infection will further go on if basic reproduction number is more than one. To ensure that the virus-free equilibrium is independent of initial concentration of cells, there is a need to study the global stability, which is shown in the next theorem.
Theorem 2. The virus-free equilibrium E 0 is globally asymptotically stable if R I < 1.
Proof. Let us define the following Lyapunov function: Taking the derivative of L along all the solutions of model (2), we geṫ [19] this implies that E 0 is globally asymptotically stable.
Biologically, this theorem depicts that infection is cleared from the body if the average number of newly infected cells is less than one, which is independent of initial concentration of cells.
(ii) Now, when z = 0, in equation (11), The above analysis shows that the two isoclines (10) and (11) intersects each other in the positive quadrant if z 2 > z 1 , i.e. R I > 1. Thus, the interior equilibrium exists if R I > 1. This shows that if the number of newly infected cells produced by a single infected cell when all other cells are healthy is more than one, then infection persists and interior equilibrium exists.

Analysis at R I = 1
To study the behaviour of the equilibrium points at R I = 1, we use the center manifold theorem [2] as, at R I = 1, linearization has a simple zero eigenvalue, hence linearization is inconclusive [31]. We made following assumptions to apply center manifold theorem to system (2).
Theorem 3. At R I = 1, virus-free equilibrium changes its stability from stable to unstable, which shows that system (2) exhibits a transcritical bifurcation at R I = 1 (Fig. 2).
In Fig. 2, we plotted the V component of the steady states against R I . When R I < 1, only E 0 exists and is stable, and as R I crosses 1, i.e. when R I becomes greater than 1, a transcritical bifurcation of steady state E 0 occurs, and E 0 becomes unstable, while E 1 comes into existence.
In the next theorem, we show that E 1 is locally asymptotically stable under certain conditions. Theorem 4. The interior equilibrium E 1 (x * , y * 1 , y * 2 , v * , z * ) is locally asymptotically stable in R 5 + if the following inequalities hold true: where c 4 < min 1 3 be the small perturbations about the interior equilibrium E 1 . Using the above new variables, we linearize system (2) around the interior equilibrium E 1 . Then, in the linear model, we consider the following positive definite function: where c 1 , c 2 , c 3 , and c 4 are positive constants to be chosen suitably. Now differentiating W 1 with respect to time t along the solutions of linearized version of model (2), a little algebraic manipulation yieldṡ Sufficient conditions forẆ 1 to be negative definite are given as follows: then conditions (17) 1 and (18) 1 are correspondingly satisfied. Further, we note that (14) ⇒ (17) 2 and (15) ⇒ (18) 2 . Hence the theorem follows.
In the next theorem, we are able to find sufficient conditions for E 1 to be globally asymptotically stable.
Theorem 5. Let the following inequalities hold in the interior of the positive orthant Ω: where m 4 < min 1 3 Then E 1 (x * , y * 1 , y * 2 , v * , z * ) is globally asymptotically stable with respect to all solutions in the interior of the positive orthant Ω.
Proof. We consider the following positive definite function about E 1 : where m 1 , m 2 , m 3 , and m 4 are positive constants to be chosen suitably. Now differentiating W 2 with respect to time t along the solutions of model (2), we geṫ Sufficient conditions forẆ 2 to be negative definite are given as follows: Let us choose Theorem 5 shows that by decreasing α and by increasing θ 1 and θ 2 , appropriately, conditions (20)-(24) may be satisfied. This shows that stability of the system increases by decreasing the infection rate of uninfected cells and by increasing the treatment rates of infected cells and virus at some appropriate levels.
Under the similar analysis as done in the previous sections, one can prove the following theorem easily.
The stability behaviour of the interior equilibrium e 1 can be studied in a similar manner as done in the previous section. By comparing R 0 and R I we note that R I < R 0 . This shows that rate of infection will be slow in the presence of immune response.

Simulation results
In this section, we present the numerical simulation results to elucidate the analytical findings. MatLab 7.10 and Mathematica 7.0 are used for simulation purposes.

Global stability and phase plane analysis of the subsystem
We note that the interior equilibrium E 1 (x * , y * 1 , y * 2 , v * , z * ) exists for the set of values of parameters given in Table 1, and (x * , y * 1 , y * 2 , v * , z * ) are given as follows: x * = 948.6799, y * 1 = 1.2014, y * 2 = 0.5565, v * = 0.2357, and z * = 0.1260. We observe that all conditions given in Theorems 4 and 5 are satisfied. This shows that E 1 is locally as well as globally asymptotically stable for the set of values of parameters in Table 1.  (2) and their values.

Parameters (units) Description
Values Source Death rate constant of x(t) 0.01 [4,33,34] Loss rate constant of y 1 (t) and y 2 (t) 0.015 [33,34] Inhibition rate constant to infection 5 - Killing rate constant of y 1 (t) by CTL IR Activation rate constant of CTL IR 0.0035 - Probability of infected cell joining y 1 0.55 [25]  Since the system is globally asymptotically stable, therefore it is independent of initial status of the subpopulations.

Effect of saturated transmission rate on uninfected cells and productively infected cells
Figures 4(a) and 4(b) represent the effect of α on uninfected cells and productively infected cells, respectively. We observe that the concentration of uninfected cells is high ( Fig. 4(a)) and that of productively infected cells is low (at zero level) for α = 0.001 and, in this case, R I = 0.4598 < 1. This shows that infection is no more at this stage. Further, the number of uninfected cells decreases with increase in α, which is obvious for real world scenarios when infection increases the number of uninfected or healthy cells declines. In Fig. 4(b), the number of productively infected cells increases with increase in infection rate α, which is also usual in real situations. Thus, infection can be controlled by controlling the spread of infection within the cells with the help of immune response and appropriate drugs. Figures 4(c) and 4(d) show the effect of α 1 , whereas α 1 changes reversibly with infection rate. In Fig. 4(c), we observe oscillations for small value of α 1 = 0.005. In the expression for infection rate, when α 1 approaches zero, then this infection rate term αxv/(1 + α 1 v) becomes bilinear, which is not realistic in case of large number of productively infected cells. It may be noted that behavior of productively infected cells is complementary to the behavior of uninfected cells, so similar explanation can be given for the oscillations in Fig. 4(d).

Effect of treatment and immune response on productively infected cells and free virus
The effect of therapeutic drugs θ 1 and θ 2 given to productively infected cells in the absence of treatment (θ 1 = θ 2 = 0), the number of productively infected cells and free virus are high. As efficacy of therapeutic drugs θ 1 and θ 2 increases, the number of productively infected cells and free virus decreases and settles down at their respective equilibrium levels. This highlights the efficacy of the drug in controlling the concentration of productively infected cells and free virus, respectively. concentration of the productively infected cells and free virus decreases with increase in immune response (µ) and approaches to zero for higher values of µ. In Figs. 6(c) and 6(d), the trajectories settle at high level for µ 1 = 0, and for higher values of µ 1 , these trajectories are settled down at lower concentration of productively infected cells and free virus, but not approaching to zero as µ 1 is the increase in immune response due to stimulation from productively infected cells. Similar explanation can be given for Figs. 6(e) and 6(f), where µ 2 is the increase in immune response due to stimulation from free virus. It is interesting to note that if immune responses are at high levels, then productively infected cells and free virus can be brought back to zero level with adequate drugs, and thus, the infection may lead to reduction in viral load.

Variation of reproduction number
In the last figure (Fig. 7), we have shown the variation of R I with respect to therapeutic drugs given to productively infected cells (θ 1 ) and virus (θ 2 ), respectively. We observe that for the less values of the therapeutic drug given to the productively infected cells (θ 1 = 0.15) and that of viruses (θ 2 = 0.3), the corresponding value of R I is 5. It is apparent from Fig. 7 that if we increase the amount of therapeutic drug given to the productively infected cells (θ 1 = 1) and virus (θ 2 = 1.1), the corresponding value of R I decreases, i.e. R I = 0.5. This emphasize that combination of therapy is useful in reducing the viral load and, ultimately, to lower the infection.

Conclusions
In this paper, we constructed a virus dynamics model to understand the viral-host interaction and elucidate the effect of total immune response (innate, CTL, and humoral) together with the treatment given to the free virus as well as the productively infected cells. The model is comprised of 5-state variables, that is, uninfected cells, productively infected cells, latently infected cells, free virus, and total immune response. We performed stability analysis of the model and observed that there exists two nonnegative equilibria: virus-free equilibrium (VFE, no persistent of virus) and interior equilibrium (IE, persistent of virus). From the global stability analysis of the virus-free equilibrium (VFE) it is depicted that the clearance of virus from the body is independent of initial status of subpopulations (Theorem 2). From the analysis of the model at R I = 1 it is observed that VFE loses its stability from the stable state to unstable state, which gives birth to the transcritical bifurcation (Fig. 2), and the interior equilibrium exists for R I > 1.
Simulation results show that the saturated infection rate gives realistic dynamics of transmission as bilinear infection rate leads to oscillatory behavior in subpopulations (Figs. 4(c) and 4(d)). The reduction in infection rate and increment in the effectiveness of therapy given to the productively infected cells and virus appropriately will lead to lower the count of infection and viral load. It is also shown that the number of new infection can be lowered via providing adequate treatment to both free virus and productively infected cells (Fig. 7). However, virus dynamics models suggest that the infection will extinct once the therapeutic drugs can target the latently infected cells, which lead to new infection upon activation.
The model also provides the insights into the dynamics of viral infection in absence of immune response. It is evident from the analysis of model (28) that the virus-free equilibrium (VFE) is globally asymptotically stable if R 0 < 1. The basic reproductive number in absence of immune response R 0 is greater than basic reproductive number in the presence of immune response R I , i.e. R 0 > R I . This implies that in the presence of immune response, the number of secondary infections will be less, which suggests that infection may be eradicated if R I < 1. From the expression of R I it is observed that the number of secondary infections decreases with the enhancement of immune response and drug efficacy. This shows that R I may be made less than one by increasing drug efficacy and improving the immune conditions. Thus, increase in treatment is effective in controlling the number of productively infected cells and free viruses. In addition, action of immune response also reduces the virus load. It has also been shown in Fig. 7 that the reproduction number can be reduced by applying adequate combination of therapeutic drugs, which helps in reduction of viral load.