DYNAMICS OF A VIRAL INFECTIOLOGY UNDER TREATMENT

This paper deals with a nonlinear model of the viral dynamics which describes the interactions between the human immune system and the virus. The novelty of this work is the introduction of combined treatments in the dynamics to modify the model. We investigate the qualitative behavior of the model and find a threshold parameter that guarantees the asymptotic stability of the equilibrium points, this parameter is known as the basic reproduction number. We estimated the parameters of the model by least-squares minimization between the numerical solution of the system and clinical data of cell cultures. It is also demonstrated that critical drug efficacy in terms of the model parameter is greatly useful to curtail the spreading of the disease.


Introduction
Mathematical modelling has become a valuable tool for the analysis of dynamics of infectious disease and for the support of treatment strategies developed in recent years. The main convergence is to analyze the transmission patterns in vivo and vitro, and methods to assess the effectiveness of treatment strategies. Infectious disease (such as HIV [1,4], HBV, HCV [10,24], and CHIKV [5,6], etc.) models which describe within-host dynamics have been described by a system of nonlinear ordinary differential equations. These diseases are most vulnerable to the human's immune system, especially to white blood cells named T cell. Due to the recent advancements in drug development, mathematical modelling of viral kinetics under treatment continues to play an instrumental role in improving our knowledge and understanding of virus pathogenesis. uninfected) cells, infected cells, and free virus pathogens as: (1.1) A lot of considerations have been added that aim to get the best representation of the virus infection. Most notable are latent cell reservoirs which serve as a major barrier in curing virus infection. Despite the fact that the antiretroviral therapy (ART) significantly limits the level of virus in the blood, there is still a low viral load due to ongoing reactivation of latent infected cells reservoirs. Variant models have been developed to study the dynamics of the virus in the presence of latent reservoirs (see, e.g. [2,7,8]). In this paper, our target is to modify the model proposed by Nowak et al. [17] with latently infected cells compartment under combined treatment strategies. We study the qualitative behavior of the models including the existence of the steady states. We investigate the local and global stability in terms of the basic reproduction number using the Lyapunov method. Further, by using clinical data from HIV infected individuals, we determine the model parameters which best fit the data. Finally, numerical simulations are presented to support the theoretical results.

Mathematical Model
We propose a latent viral dynamics model by incorporating the efficacy of combined treatments. The model takes into account four types of compartments: uninfected target cells (T ), productively infected cells (I), latently infected cells (L) and the number of virions V (t) in plasma all at time t. We describe the infection dynamics using the following differential equations: The target cells, T , die at rate d T T and are recruited into the infection site at rate λ. When these cells become infected by virus they can be eliminated at the constant rate k, which is directly proportional to the product of the participating populations. We assume that a fraction, f , of infection, generate latently infected cells with replication-competent genomes and the remaining fraction of infection, (1 − f ), leads to productively infected cells. Latently infected cells become productively infected at a contant rate of α due to activation. We assume that each productively infected cells produce N viral particles before it dies and the virus is then killed off at a clearance rate, d V , proportional to the virus population. In some literature [17], N d I is called the virus proliferation rate p. Productively infected cells and latently infected cells die at rates d I I and d L L, respectively. In addition, the drug efficacy of combined treatment in this model involved reducing new infections and blocking virions production, which are described in fractions (1−u 1 ) and (1−u 2 ), respectively, with 0 ≤ u 1 , u 2 ≤ 1. Based on biological considerations we assume that these model parameters are positive constant and given in Table 1. Notably, it is biologically reasonable to assume that infected cells have a higher death rate than target cells, namely d I > d T . Furthermore, we assume that the total number of target cells remains approximately constant, to make each T cell is susceptible only to the virus. A schematic representation of the model (2.1) is given in figure 1.

Qualitative Study of the Model
Now we must prove that solutions to the model (2.1) exist and they are positive as well as bounded for all values of time in order to retain the biological validity of the model, . Let t 0 > 0, In the model (2.1), if the initial conditions satisfy T 0 > 0, I 0 > 0, L 0 > 0 and V 0 > 0, then for all t ∈ R the functions T (t), I(t), L(t) and V (t) will exist in a quadruple set denoted as R 4 + . Proof. The Picard-Lindelöf Theorem [14] states that for the initial value problem y ′ (t) = f (y(t)), y(t 0 ) = y 0 , t ∈ [t 0 − ϵ, t 0 + ϵ], if f is locally Lipschitz in y and continuous in t, then for some value ϵ > 0, there exists a unique solution y(t) to the initial value problem within the range [t 0 − ϵ, t 0 + ϵ]. Since the system of ODEs is autonomous, it suffices to show that the function f : is locally Lipschitz in its y argument. Note that the Jacobian matrix for some y * ∈ I. By letting |∇f (y * )| = K, we obtain |f (y 1 ) − f (y 2 )| ≤ K|y 1 − y 2 | for all y 1 , y 2 ∈ I and therefore f (y) is locally bounded for every y ∈ R 4 . Hence, f has a continuous, bounded derivative on any compact subset of R 4 and so f is locally Lipschitz in y. By the Picard-Lindelöf Theorem, there exists a unique solution, y(t), to the ordinary differential equation y ′ (t) = f (y(t)) with initial value y(0) = y 0 on [0, t 0 ] for some time t 0 > 0. The next step in analyzing our model will be to prove positivity and boundedness for the system of differential equations. Theorem 3.2 (Positivity). Let t 0 > 0, In the model (2.1), if the initial conditions satisfy T 0 > 0, I 0 > 0, L 0 > 0 and V 0 > 0, then for all t ∈ [0, t 0 ] the functions T (t), I(t), L(t) and V (t) will be remain positive in R 4 + . Proof. Since all of the parameters used in the system are positive, we can place lower bounds on each of the equations given in the model. Thus, Through basic differential equations methods we can resolve the inequalities and produce: Thus, for all t ∈ [0, t 0 ] the functions T (t), I(t), L(t) and V (t) will be positive and remain in the set R 4 + . Theorem 3.3 (Boundedness). Assume the initial conditions of (2.2) satisfy T 0 > 0, I 0 > 0, L 0 > 0, V 0 > 0 and 0 ≤ u 1 , u 2 ≤ 1. If the unique solution provided by Theorem 3.1 exists on the interval [0, t 0 ] for some t 0 > 0, then the functions T (t), I(t), L(t) and V (t) will be bounded and remain positive for all t ∈ [0, t 0 ].
Proof. The state variables we consider here represent supersolutions for given problems (2.1-2.2). From the given equations we have where R + , is the set of non-negative real numbers. The upper bound for X is also the upper bound for T, I, and L. Lastly which also leads to Since 0 ≤ u 1 , u 2 ≤ 1, then, T (t), I(t), L(t) and V (t) are bounded above with values elements of R + . Via a maximum principle [20] theory for first-order nonlinear differential equations, we obtain the solutions to the problems (2.1-2.2) bounded for all t ∈ [0, t 0 ] and lies in the compact set where the quadruple set R 4 + defines as R 4

Equilibria of the System
An equilibrium point is the constant solution of (2.1) so that the rate of change for each compartment is zero. By setting the right-hand side of (2.1) to zero, we get exactly two equilibria which are biologically meaningful. The non-infective equilibrium (viral extinction) as and the infective equilibrium (viral persistence) as

Basic Reproduction Number
To apply the next generation method [3,12] to the model (2.1), we need the compartments that spread the infection, so we need only the infected I, latent L and virions V , compartments. Let us define the model dynamics using the equations For this system, at the disease free equilibrium point Then, for the system (2.1), the next generation matrix is The dominant eigenvalue of F V −1 is the basic reproduction number for the system (2.1) and given by expression The basic reproduction number R L is the average number of secondary infections produced when one single virus cell is introduced into a host where every T cell is susceptible [12].
Remark 3.1. Using basic reproduction number R L the infected equilibrium point

Local Stability of the Equilibria
In this section, we investigate the local stability properties of the equilibrium points by approximating the nonlinear system of differential equations (2.1) with a linear system at the points E 0 and E * . Then, we locally perturb the system from equilibrium and examine the resulting long time behavior. This is done by linearizing the system about each equilibria, using the Jacobian for (2.1): Theorem 4.1. If R L < 1, then the non-infective equilibrium is locally asymptotically stable. If R L > 1 then the non-infective equilibrium is an unstable saddle point, and the endemic equilibrium is locally asymptotically stable.
Proof. We proceed by linearizing the system and using the Routh-Hurwitz criterion [9] to determine conditions under which the linear system possesses only negative eigenvalues. Then, as a consequence of the Hartman Grobman theorem [11], the local behavior of the linearized system is equivalent to that of the nonlinear system. First, we compute the Jacobian evaluated at the non-infective equilibrium E 0 = From this, we compute the associated characteristic polynomial for eigenvalues η as Since η < −d T < 0 is the one negative eigenvalue of the system, After expanding the remaining terms and ordering by powers of η, this equation ultimately simplifies to

According to the Routh-Hurwitz criteria, all roots of this cubic equation possess negative real part if and only if
Using the inequality and the previous condition and the Routh-Hurwitz criteria are satisfied. Thus, R L < 1 implies that all eigenvalues of the linearized system are negative, and hence the local asymptotic stability of E 0 follows. Conversely, if R L > 1, then the linearized system possesses at least one positive eigenvalue, and the equilibrium is unstable.
The analysis for E * is similar to that of E 0 . Linearizing (2.1) about E * , we find the Jacobian and this results in the characteristic equation After expanding the terms and ordering by powers of η, this equation ultimately simplifies to a quartic polynomial According to the Routh-Hurwitz criteria, all roots of this quartic equation possess negative real part if and only if As for the E * analysis, the positivity of A 1 follows directly from the positivity of the coefficients, and after rewriting A 4 , we find Hence, it is necessary that R L > 1 in order to satisfy the criteria. Similarly, we rewrite A 3 as In this inequality we have canceled the third term with the last term and utilized the inequality (4.1) to bound the fourth term. The only nonpositive term in A 2 can be rewritten as By the definition of A 1 , we have A 1 > d I + d V and using the above inequality for A 2 , we find Finally, we verify the last inequality Hence, we obtain With this, all of the criteria have been satisfied and E * is stable if R L > 1. Conversely, if R L < 1, then the Jacobian possesses at least one positive eigenvalue, and the endemic state is unstable. Our analysis reveals if R L < 1 and cell values begin within a sufficiently close distance of E 0 , then they will tend to E 0 as t → ∞. Contrastingly, if R L > 1 and initial populations are sufficiently close to E * , they will tend to E * in the long run. Theorem 4.1 also emphasizes the crucial feature that equilibria are not stable simultaneously, that is, bistability of E 0 and E * does not occur.

Global Stability of the Equilibria
We also establish global asymptotic stability of the equilibria using a Lyapunov function which demonstrates the stronger result that initial values of cells have no effect on their long term (t → ∞) limiting values. Before proceeding with the global stability analysis for the model (2.1), we present some inequalities developed in [13], which will be used in the proofs. To begin with, we consider the function x n be positive numbers. Then, Summing over i = 1 to n, from above equation we obtain If p 1 , p 2 , · · · , p n = q 1 , q 2 , · · · , q n , then Proof. To investigate the global stability of E 0 , consider the following Lyapunov function Notice that U is nonnegative, and U is identically zero if and only if it is evaluated at the non-infective equilibrium point (T 0 , I 0 , L 0 , V 0 ) = λ d T , 0, 0, 0 . We compute the derivative along trajectories and find After using the definition of T 0 , we are left with Thus, under the assumption that R L ≤ 1, we see that dU dt ≤ 0 for all positive values of T, I, L, and V , and the global asymptotic stability follows by LaSalles Invariance Principle [15].
Turning to the endemic equilibrium, none of the end values are zero, so we denote this steady state by (T * , I * , L * , V * ) and define a Lyapunov function as This function is nonnegative and identically zero only when evaluated at the endemic equilibrium E * = (T * , I * , L * , V * ). Computing the derivative along trajectories yields For U 1 , we factor out a d T T * term and use the form of T * = λ d T R L to find For U 2 , we factor out a L * term and use the following identities to find Thus, combining the rearrangements of U 1 and U 2 and using the following relation Since Since, each of the resulting terms above are nonpositive because the arithmetic mean is greater than the geometric mean, using the inequality ( [15], the endemic equilibrium E * is globally asymptotically stable if R L > 1.
This analysis reveals one very important fact about the overall system: the end states of the disease dynamics are only dependent on the value of R L , and not any other parameter or initial value. If R L > 1, then the system tends to E * , an end state with non-zero infected cells and virions, but if R L < 1, then the final equilibrium is E 0 , which contains neither virions nor infected T -cells.

Parameter Estimation for Model (2.1) without Treatment
The mathematical analysis of models is very useful for understanding asymptotic behaviors and longtime qualitative outcomes, while the outcomes of a model critically depend on the values of the model parameters. Since models are confronted with disease data, an accurate estimation of parameter values is essential for reliable quantitative predictions within a finite time interval. For estimation of multiple parameters, a systematic approach for the fitting is desirable. Different techniques were used for estimating the parameters in [21][22][23]. In this section, we have used a straight forward method to calculate the parameters, which is known as nonlinear least square (NLS) method.

Nonlinear Least Square Method
In least-squares approach, we assume that the time coordinates of the data are exact, but their corresponding y-coordinates (virions) may be noisy or distorted. We fit the solution curve through the data so that the sum of the squares of the vertical distances from the data points to the point on the curve is as small as possible. This distance is commonly known as least squares error. Next we illustrate how to use NLS method to estimate unknown parameters

Step 1. Data Collection
In particular, suppose we are fitting the virions V (t), with the given data Step 2. NLS fitting So the basic problem is to identify the set parameters θ such that the following sum-of-squares error (SSE) is as small as possible: where V (t i , θ) represents the virus concentration at time t i with parameter θ and V (t i ) represents the data of patient's virus concentration at time t i . Such a problem is clearly a nonlinear least-squares problem, since the dependence of a solution on the parameter θ is through a highly nonlinear system of differential equations.

Step 3. Solve the NLS problem numerically
We use a Matlab functions fminsearch which takes the least-squares error function SSE(θ) and an initial guess of the parameter value θ 0 , and uses a direct search routine to find a minimum value of least-squares error.
Certain parameters such as production rate λ and natural death rates d T of T cells can be estimated directly from population data as given in Table 1. The rest of the parameters θ = (k, f, d I , α, d L , N, d V ) are estimated from the set of data gathered from plasma donor samples obtained in [17] at primary stage of HIV infection. We have taken most of our intial parameters from previous literature [19] except the fraction of latent infection f . Since f ∈ (0, 1), we take the intial guess as θ 0 = (2 × 10 −7 , 0.1, 0.5, 0.4, 0.004, 50, 5) with intial conditions (T 0 , I 0 , L 0 , V 0 ) = (10 6 , 0, 0, 15.8), under no treatments (u 1 = u 2 = 0) we obtained estimated parameters of model (2.1) in the Table 1.

Critical Drug Efficacy
One of the important feature from the mathematical analysis reveals that long time disease dynamics depends on the infected steady state which explicitly depends on the basic reproductive number (R L ). So a larger basic reproduction number retains disease progression for larger period of time. By choosing the new parameter values from Table 1, under no treatments (u 1 = u 2 = 0) the value of R L turns out to be In the system (2.1), the efficacies of drugs are incorporated through the terms (1 − u 1 ) and (1 − u 2 ) respectively. The values, u i = 0 and u i = 1 (i = 1, 2), reflect completely ineffective and perfectly effective treatment respectively. For brevity, the efficacies of drugs are combined to obtain a new term to reflect the overall efficacy for this combination treatment and is given by , this rearrangement indicates that the drugs act independently of one another. Note that u = u 1 + u 2 − u 1 u 2 represents the total combined drug efficacy. This choice is motivated by the condition for stability of E 0 and E * . Recalling that the stability criterion (see theorem 5.1) for E 0 is R L < 1, which equivalent to .
Similarly the condition R L > 1, for E * to be stable is equivalent to Thus, there is a transcritical bifurcation point given by .
Motivated by this we define the critical efficacy, u c by .
Thus, in order to achieve a successful treatment by way of elimination of virus persistence, i.e., the uninfected steady state E 0 being stable we need u > u c (≡ R L < 1). On the other hand, whenever u < u c (≡ R L > 1), the infected steady state E * remains stable and the infection persists. With the base-case parameters given in Table 1, in absence of treatments, the basic reproductive number is R L = 18.37, by equation (6.1). This shows that to avoid infection the combined drug efficacy u c should be maintained at a constant greater than i.e., maintaining constant drug effectiveness of at least 95% should theoretically avoid infection. The goal is to choose u 1 and u 2 so that u c > 0.95(≡ R L < 1) hereby resulting in a stable uninfected steady state.

Numerical Results
In this section, we explore the dynamics of virus infection model (2.1) to study the impact of different treatments strategies on the proliferation of the viral and infected cells within the host. Using various combinations of the two drugs, we investigate and compare the numerical results from simulations. In doing so, we are able to numerically illustrate how the efficacy of the drugs affects the level of infection in order to achieve viral clearance.

Without Treatment Strategy
The results obtained for the stability of the uninfected and the infected steady states are also numerically illustrated in this section. For this purpose, we take into account two sets of parameters corresponding to the cases of stability of the infected steady state R L > 1 and uninfected steady state R L < 1. All the models are numerically solved using the Runge-Kutta 4th order scheme. First, figure 2 and 3 illustrate the system dynamic interaction between the cells T (t), I(t), L(t), and V (t). We can see that upon initiation of infection, the population of the infected cells (I, L) and virus V -cells increases significantly until it reaches the peak. After achieving a peak, these cells decay until it reaches a steady state. As we see that during the increase of the virus-cell population, the population of the target T -cells decreases (from 10 6 cells ml −1 d −1 to 5 × 10 4 cells ml −1 d −1 ). However, after reaching the minimum, the target cell population begins to increases until it ultimately reaches a steady state. In this case, the steady state (5.4 × 10 4 cells ml −1 d −1 ), which is approximately 5.4% of the original population of T -cells.  Using the parameter values from Table 1, under no treatments (u 1 = u 2 = 0) the value of R L turns out to be R L = 18.37 > 1 and thereby indicating that the infected steady state is asymptotically stable. For this purpose, we choose three different initial conditions of (T 0 , I 0 , L 0 , V 0 ) as IC1 = (10 6 , 0, 0, 15.8), IC2 = (10 4 , 10, 10, 158), and IC3 = (10 5 , 100, 100, 1580). The evolution of the dynamics of the modified model for this scenario was observed for a duration of 60 days and we found the states of the system eventually converges to the infected steady state E * = (5.44 × 10 4 , 1.17 × 10 5 , 1.78 × 10 4 , 5.4 × 10 6 ) for all the three initial conditions. This is illustrated in Figure 4 which supports the result that the infected steady state, E * is asymptotically stable whenever R L > 1 and eventually patient does not recover without treatments. In order to study the case when R L < 1, we now choose a different value of k, namely k = 1 × 10 −8 ml d −1 , while retaining the other parameter values. Then the value of R L = 0.57 < 1. Consequently, for this scenario, the uninfected steady state E 0 would have to be asymptotically stable. To illustrate we again choose three different initial conditions as IC1 = (10 6 , 100, 100, 15.8), IC2 = (10 4 , 10, 10, 158), and IC3 = (10 5 , 0, 0, 1580) and ran the simulation for a duration of 120 days. It can be observed, from Figure 5, that all the state variables of the system eventually approach to the uninfected steady state E 0 = (10 6 , 0, 0, 0) indicating the asymptotic stability of the uninfected steady state.

Constant Treatment Strategy
Now, we numerically examine the impact of the constant efficacies u 1 and u 2 on the basic reproduction number R L . Recall that the infection clears out or persists whenever R L < 1 or R L > 1, which is equivalent to u > u c or u < u c . In section 6.2, we show that our goal is to choose u 1 and u 2 such that R L is driven to a value less than 1 and combined drug efficacy We illustrate this by a surface plot and a contour plot in Figure 6. We can easily  observe that for u 1 = 0 and u 2 = 0 the value of R L attains its maximum value of 18.37. We increase u 1 and u 2 from 0 to 1 and observe that the value of R L gradually decreases and eventually tends towards 0 (corresponding to u 1 = 1, u 2 = 1). This clearly reflects the impact of the efficacies in terms of clearance of the infection. Now, we simulated the dynamics of the system as a result of administration of the combined constant treatments. Recall that the parameter values (without treatment) chosen are for R L > 1. In Figure 7, we fix the efficacy of u 1 = 0.75 and consider three different efficacies of u 2 = 0.7, 0.86, 0.9. For u 2 = 0.7, the combination efficacy, u = 0.925, which is less than the critical efficacy u c = 0.95 with R L = 1.38 > 1. In this case, the levels of infected I-cell, Latent L-cell and virions V show some signs of decline over a period of 60 days as can be seen in Figure  7. But if the simulation is run for a longer time period, then we can see that despite the initial signs of patient recovery, the levels of all three cells will rebound and eventually move towards the infected steady state. Further, for u 2 = 0.86, 0.9, the combination efficacy u is always greater than the critical efficacy u c i.e., R L < 1. For these cases, the levels of I, L and V show a gradual decline over the period of 60 days and simulation for a longer period also confirms that the populations tend towards the levels for the uninfected steady state E 0 . We observe that this decline is biphasic in nature in case of V with a more rapid decline in the first phase of a couple of days followed by a slower decline, which is consistent with clinical results [18] for an HIV patient. We observe similar results by fixing u 2 = 0.75 and varying the values of u 1 = 0.7, 0.86, 0.9. These results are presented in Figure 8. We note that in this case also there is a biphasic decline that is observed earlier.
Interruptions in treatment can happen due to a variety of reasons such as side effects and financial constraints for a continued long term treatment [16]. To illus- trate one such scenario, we consider three sets of combination treatment (u 1 , u 2 ) as (0.88, 0.6), (0.83, 0.75), and (0.83, 0.75) for a period of 60 days. For these three pairs of (u 1 , u 2 ), u > u c . We see in Figure 9, as to how the viral load declines in a biphasic manner if full treatment is adminstrated. The decline is significant (approximately 10 2 folds). The discontinuation of treatment after 45 days results in the rebound of the levels of HIV virions. Once the treatment period of 45 days is over, we observed the dynamics of the system for another 15 days starting from the levels at the cessation of treatment after 45 days. It can be seen that the peak viral load on an average is lower with this on-off therapeutic protocol as compared to the scenario when no treatment is administered over the entire period of 60 days.

Conclusion
In this study, we sought to learn more about virus infection dynamics by introducing and analyzing mathematical models of the immune system in the presence of constant treatments. We proved existence, uniqueness, positivity, and boundedness for the model and derived the conditions on basic reproduction number that guarantees the asymptotic stability of the equilibria. Further, by using clinical data from HIV infected individuals, we determine the model parameters which best fit the data for long time dynamics. Our model exhibits a wider variety of parameter values that lead to long term viral persistence as t → ∞ due to the appearance of latent T cells.
In addition, we also examined how treatments impact the proliferation of the virus.
In doing so, we used asymptotic stability analyses to define treatment thresholds in order to eliminate the virus and clear the infection. Additionally, we were able to estimate the necessary drug efficacy of treatments for infected patients. All the results are further illustrated by a number of numerical simulations. To speculate long term latency of a viral infection, model (2.1) can be extended further by using delay in the different cell compartments or by imposing effect of immune cells like B cells, killer cells, Dendritic cells, Myeloblasts, etc. Moreover, to intensify our discussion about parameter estimation, an alternative method like Monte Carlo Markov Chain (MCMC) can be used to compare the accuracy.