NONLINEAR EFFECTS IN THE DYNAMICS OF HIV-1 INFECTION PREDICTED BY MATHEMATICAL MODEL WITH MULTIPLE DELAYS

. We formulate a novel mathematical model to describe the development of acute HIV-1 infection and the antiviral immune response. The model is formulated using a system of delay diﬀerential- and integral equations. The model is applied to study the possibility of eradication of HIV infection during the primary acute phase of the disease. To this end a combination of analytical examination and computational treatment is used. The model belongs to the Wazewski diﬀerential equation systems with delay. The conditions of asymptotic stability of a trivial steady state solution are expressed in terms of the algebraic Sevastyanov–Kotelyanskii criterion. The results of the computational experiments with the model calibrated according to the available estimates of parameters suggest that a complete elimination of HIV-1 infection after acute phase of infection is feasible.


1.
Introduction. Mathematical modelling in immunology represents a rapidly expanding field at the interface of applied mathematics and life sciences. The formulation of mathematical models of the Human Immunodeficiency Virus type 1 (HIV-1) infection dynamics is among the most actively developing branches of mathematical immunology (see [6], [16], [19], [27], [21], [2] and the references therein). The proposed mathematical models aim to improve the understanding of the following aspects of HIV infection: • the regularities of HIV infection dynamics, including the effects of antiretroviral therapy (ART); • the mechanistic analysis of various hypothesis about the pathogenesis of HIV infection determined by the interaction of virus, target cells and immune defence; • the data analysis in the context of the phenotypic diversity of HIV infection dynamics.
Many mathematical models of HIV infection take the form of the systems of delay differential equations, see the examples in [17], [1], [20], [25], [30], [4]. The justification of the structure and the incorporation of time-delays into mathematical models of immune response can be found in monographs [14], [15]. The presence of single-or multiple delays in conjunction with nonlinear terms on the right-hand side of model equations can dramatically change the behavior of the solutions concerning the following: • asymptotic stability versus instability of steady states; • the transition time to a neighborhood of given steady state; • the features of transient dynamics; • the emergence and properties of oscillatory behavior, etc.
In this paper we formulate the mathematical model to describe the development of HIV-1 infection and the antiviral immune response taking into account the following processes [28], [29], [13]: The set of statements [F1-F3] considers some essential stages of virus infection and replication cycle and the induction of antiviral immune response. The statement [F3] suggests a positive feedback of the virus population on the proliferation and differential of CTLs. Factor [F4] is needed in order to address the role of latent reservoir which is recognized as a major barrier to curing HIV-1 infection [28]. The destruction of productively infected target cells due to the cytopathic effects of the virus and killing by CTLs is taken into account by implementing [F5,F6]. The reduction of the life-span of effector CTLs because of a number of processes including but not limited to the killing of target cells, expression of programmed death-1 receptor is parameterised in the model under assumption [F7]. As the number of uninfected target cells is quite large, we take into account the reduction of availability of freely circulating viruses which are spent on the infection of the cells as well as the latently infected cells and the cells undergoing the activation before the production of the virus [F8]. This loss of the free virus is assumed to be negligible in vast majority of simple models of HIV infection. The set of above empirical features of HIV infection suggests the formulation of mathematical model which is novel as compared to the existing models mentioned above.
In this study we develop and analyse both qualitatively and computationally the mathematical model of primary HIV-1 infection, built according to [F1-F8]. In sections 2, 3 the set of the model delay differential equations (DDEs) is formulated. In Section 4 we consider the existence and uniqueness of solutions of the initial value problem for the DDE system, and analyze the stability of its steady states. Section 5 is devoted to the formulation of the conditions for the eradication of HIV-1 infection. Computational experiments with the model are presented in Section 6. These are used to gain additional insight into necessary conditions for a complete elimination of HIV-1 infection shortly after primary infection. Section 7 is devoted to the discussion of the results and the need to study the problem of eradication of HIV-1 infection on the basis of hybrid mathematical models.

DDE model of HIV infection (description and assumptions).
We consider the following time-dependent characteristics of HIV-1 infection as the state space variables of the mathematical model: 1. U (t) -immature budded virus particles; 2. V (t) -infectious virus particles; 3. A(t) -uninfected target cells (Dendritic cells, Macrophages, CD4+ T lymphocytes); 4. L(t) -latently infected cells (cells with a stably integrated but transcriptionally silent form of the virus); 5. C(t) -cells productively-infected before the virus release; 6. I(t) -infected cells releasing the virus particles; 7. E 0 (t) -naive precursor CD8+ T lymphocytes; 8. E 1 (t) -activated proliferating CD8+ T lymphocytes; 9. E(t) -effector cytotoxic CD8+ T lymphocytes; 10. Q(t) -central memory CD8+ T lymphocytes; 11. E 2 (t) -activated central memory CD8+ T cells; Let S denote the cell source in bone marrow and thymus generating the target cells and precursor CD8+ T cells; D -standing for the dead cells and degraded virus particles, and W -the viruses spent on infection of the target cells A, L, C. The scheme of the model specifying the populations of cells and viruses, the elementary processes and the transitions between the specified states is shown in Figure 1. Now we proceed with formulating the assumptions about elementary processes which determine the dynamics of the cell and virus populations.

P1:
The target and precursor cells A and E 0 are produced at given rates r A , r E0 , respectively P2: The cells A, L, C, I, E 0 , E, Q, and virus particles U , V decline via natural death with the following per capita death rate constants so that the following transitions take place P3: The target cells A interact with the infectious viruses V at a rate γ A,V > 0 per pair (A, V ); the free viruses are removed from circulation due to infection; the fraction of target cells that become latently infected is p L and the rest represents the productively infected cells 1 − p L : P4: The target cells A interact with the productively infected cells I at a rate γ A,I > 0 per pair (A, I) to become infected cells:

P7
: The cell C, infected by the virus at time t − ω C , ω C > 0 and survived for the time period (t − ω C , t), becomes a virus producing infected cell at time t:

P8
: The virus particles U are produced by infected cells at a per capita rate ν I > 0 (the budding of U wherein virion crosses the plasma membrane and obtains its lipid envelope [29]): P9: The virus producing infected cells I die at a per capita rate σ I ν I with the fraction of cells dying due to virion budding, specified by parameter 0 < σ I < 1: P10: The virion U which is released at time t − ω U , ω U > 0 and survived for the time period (t − ω U , t), maturates at time t maturation (the virion changes structure and becomes infectious [29]): P11: Infectious viruses V can interact with the cells L and C at rates γ L,V > 0, γ C,V > 0 per pair (L, V ) and (C, V ), respectively; the viruses leave the circulation; however, the cells L, C do not change their phenotype: P12: HIV-specific cytotoxic T lymphocytes E interact with the productively infected cells I at a per pair (E, I) rate γ E,I > 0; this killing results in elimination of virus producing cells and the fraction p E of CTLs mediating the killing is assumed to die while the rest 1 − p E remains alive and continues to function, with 0 ≤ p E ≤ 1: P13: HIV-specific CD8+ T cells E 1 become activated due to contacts of the precursor cells E 0 with target cells A providing that the viruses V are around at a per triple (E 0 , A, V ) rate constant γ E0,A,V > 0: P14: The clonal expansion of activated HIV-specific CTLs E 1 requires some time ω E1 > 0. Due to division and differentiation of cell E 1 activated at time t − ω E1 the clone of cells is generated at time t consisting of n E1 > 0 cytotoxic T cells and n Q1 > 0 central memory T cells: P15: Memory CD8 T cells E 2 start proliferation being activated after contacts of memory cells Q with target cells A providing that the viruses V are around at a per triple (Q, A, V ) rate constant γ Q,A,V > 0: P16: The activated memory T cells E 2 implements an automomous proliferation and differential programm over the period ω E2 > 0. The cell E 2 activated at time t−ω E2 gives rise to a clone of n E2 > 0 effector CTLs and n Q2 > 0 central memory CTLs: 3. Equations of the model. Let us consider the cell and virus population densities at time t introduced in Section 2: Note that the activated proliferating cells E 1 , E 2 do not explicitly enter the model equations and were used only to explain the clonal expansion of naive and memory T cells. Introduce the following notation where ρ C (t) is the total rate of appearance of infected cells C, ρ U (t) -denotes the rate of secretion of premature virus particles U , ρ E1 (t) and ρ E2 (t) represent the rates of activation of naive precursor (E 1 ) and central memory (E 2 ) CD8+ T cells at time t.
Using the assumptions P1-P16, the set of the model equations can be formulated as follows: Consider the following initial data for system (1)-(9) where ω max = max{ω C , ω U , ω E1 , ω E2 }. All functions appearing in (10), are assumed to be non-negative and continuous. The time derivatives at t = 0 in equations (1)-(7) are defined as right-hand derivatives of the respective variables.

4.
Well-posedness of the model and steady state analysis. Equations (8), (9) can be considered as latent enabling the quantification of the cell and virion numbers at any given time t. Substituting the expression for C(t) into equation (4), we arrive at the following equivalent system of delay differential equations to be studied instead of the set (1)- (10): supplemented with initial data (10). Again, at t = 0 the time derivatives of (11)- (17) are considered to be the right-hand derivatives of the respective variables.
Let us denote by Z the following vector  17) for all t ∈ [0; ∞) and the initial data (10).
Using the results of study [22], the following statement can be proved. Proof of Theorem 4.1. We note that from definition (10) and equation (11) , I(t) are non-negative and continuous functions, then Let us apply the method of steps [7] to the initial value problem specified by equations (10)- (17). Let 0 < τ < ∞ be fixed. We consider the problem on consecutive closed intervals [a; b], with the set of bounds a, b, including the time point a = 0, the points, which are equal to the multiples of time-delays ω C , ω U , ω E1 , ω E2 , and the final time point b = τ . For each [a; b] we consider the initial value problem for equations (11)-(17) and an equivalent system of integral equations on [a; b], which is obtained by a variation of constants formula. The initial and the integrand functions, appearing in the non-linear integral equations satisfy the assumptions (H0) and (H3) in [22]. Specifically, the functions are non-negative and bounded by some constants which depend on the solution history. If functions A(t), V (t) are non-negative and continuous, then inequality The integrand functions are locally Lipschitz-continuous. It follows that the system of nonlinear integral equations can be bounded by a majorant system of linear integral equations for all non-negative continuous functions. Hence the integral operator, determined by the system of nonlinear integral equations and the initial data from the previous integration time interval has an invariant subset of non-negative continuous functions defined on the closed interval [a; b]. Using Lemmas 1, 5 and Theorem 1 from [22], we obtain that the problem (10)-(17) has a unique solution on fixed interval [0; τ ], with nonnegative components which continuously depend on the initial data (10) given on the closed interval [−ω max ; 0]. As the time τ is taken arbitrary, Theorem 4.1 is proved.
Theorem 4.1 and the expressions (8), (9) suggest that the functions C(t), U (t) are defined, continuous and non-negative on the whole interval [0; ∞). Note that, the formulas (8) and (9) can be considered as a compact rewriting of the following differential equations for the variables C(t) and U (t) in integral form: The derivatives of the variables C(t), U (t) for t = −ω C ; 0, t = −ω U ; 0 are considered to be the right-hand derivatives of the functions. For t ≤ 0 the functions ρ C (t), ρ U (t) are specified by the initial data (10). The solutions to equations (18) and (20) define the initial data for delay differential equations (19) and (21).

Steady states.
The equations of the model (1)-(9) have the following steady state The above equilibrium (22) will be called a trivial steady state. The trivial steady state (22) represents the solution of the model which corresponds to a healthy human organism with no HIV-1 infection: We derive conditions of asymptotic stability (unstability) of the trivial steady state (22). For the clarity of presentation let us introduce the following parameters: Following the adopted terminology, (e.g., [27], [21]), the parameter R 0 will be called the basic reproductive ratio. The parameters R 0,1 , R 0,2 determining the overall values of R 0 admit the following interpretation:  7) and (19). Consider the linearization of the system vector field about a trivial equilibrium point, i.e., the linear variational equation for the variables . The corresponding linear system of equations we subdivide into blocks, with the first block considering only variables x 2 (t), x 3 (t), x 4 (t): The terms being neglected after linearization of the model system meet the necessary smallness requirement [7], [5], [12]. The block structure of system (23)- (29) and the positivity condition on the parameters µ jj > 0, j = 1, 5, 6, 7, 8, reduce the stability analysis of the trivial steady state of the full system (23)- (30) to the study of stability of a smaller system (23)-(25) (see further details in [24]).
The system of differential equations (23)-(25) possesses a special structure and belongs to the Wazewski differential equation systems with delay. Necessary and sufficient conditions of asymptotic stability of a trivial steady state solution of the Wazewski system are established in [18]. The conditions are expressed in terms of the algebraic Sevastyanov-Kotelyanskii criterion [8]. We apply the criterion to the matrix with the elements of the matrix row corresponding to the parameters of the model equations (23), (24), (25). By applying the Sevastyanov-Kotelyanskii criterion, one needs to check the fulfilment of the inequalities M 1 < 0, M 2 > 0, M 3 < 0, where M k is a diagonal minor of matrix H * of order k = 1, 2, 3. Thus, we obtain that where M 3 is the determinant of H * : Using the above expressions for M 1 , M 2 , M 3 , we are led to the following statement.
Theorem 4.2. If the basic reproductive ratio R 0 is R 0 < 1, then the trivial steady state (22) is asymptotically (locally) stable.
We consider now the case R 0 ≥ 1. Following a standard procedure elaborated in [7], [12], we search the solutions of system (23)- (25) in the form x(t) = c exp(λt), where c ∈ R 3 is non-zero and λ is a complex number being the root of the characteristic equation with diag(1, 1, 1) being an identity matrix. Equation (32) contains the matrix which appears after substitution of x(t) = c exp(λt) into equations (23)- (25). Consider equation (32), assuming that λ ∈ R (real root). Note that H(0) = H * and according to (31) P (0) = det(H * ) ≥ 0 for R 0 ≥ 1. If R 0 = 1, then λ = 0 is a zero root of equation (32). According to [7] the case when R 0 = 1 does not allow one to conclude about the stability or instability of the trivial steady state (22) and additional examination is needed. If we assume that R 0 > 1 then P (0) > 0. Using the explicit expression for P (λ), one can see that for the large enough λ > 0 the following estimate holds true P (λ) ≤ −λ 3 + const. Taking into account this upper estimate and the continuity of P (λ), we conclude that there exists the root of equations (32) such that λ 1 > 0. Using the results of [7], we can formulate the following statement.  The rank of matrix H(λ 1 ) − λ 1 diag(1, 1, 1) is 2, and c (1) has the following form where β 1 = 0 -arbitrary real constant. The function is a particular solution of equations (23)- (25). Note, secondly, that the state vector function component x 8 (t) does not enter (23)- (30). Hence, the equations (30) can be transformed to an equivalent integral form similar to (8).

5.
Specifying the conditions representing the elimination of HIV-1 infection. In the previous section it was shown that the initial value problem for the model (1)-(10) admits the existence of solutions with some components being equal to the respective equilibrium components (22). This type of solution can be interpreted as representing the infection free state of the host organism. If the initial data (10) are such that at least one of the initial functions is not equal to zero, then the host organism is HIV infected. Then the question arises of whether the eradication of HIV-1 infection can take place over a finite time interval.
It follows from Theorems 4.1 and 4.2 that the elimination of HIV-1 infection can be achieved over an infinite time interval when R 0 < 1 and for small enough initial doses of infection V 0 (t). Indeed, these ensure that there exists asymptotic limits of the solution components as follows This conclusion follows from the instability of the trivial steady state (22) (see Theorem 4.3). Indeed, when the steady state (22) in unstable, then for any initial dose of infection V 0 (t) > 0 the asymptotic behavior of solution components (33) will not take place. Therefore, for any T > 0 there exist t 1 > T and at least one variable in the set (33) such that one of the following conditions are met: V (t 1 ) > ε, U (t 1 ) > ε, L(t 1 ) > ε, I(t 1 ) > ε, C(t 1 ) > ε for some ε > 0. In turn, this implies that for R 0 > 1 the eradication of HIV-1 infection is not possible.
Following the approach presented in [14], [15], consider the initial data (34) These data set (34) describes the primary infection of a healthy individual at time t = 0 with V 0 viruses.
In contrast to general initial data set (10), the initial data (34) consider the virus dose function V 0 (t), that has a jump discontinuity at t = 0. The other initial functions are constants. This allows one to implement the method of steps for studying the model behaviour [7]. This method can be consecutively applied to the system (1)- (7), (19) with the initial data (34) over the time interval of interest [0; T mod ]. The equation (1)- (7), (19) are considered on sequential time intervals [a, b], composing the set of interval limits which include a = 0, the instants a, b being multiples of the delays ω C , ω U , ω E1 , ω E2 , and the final time b = T mod .
Note that the initial functions (34) which are zero simplify the structure of the right-hand terms in model equations for first time intervals [a, b]. The solution of the model on [0; T mod ] will be the set of functions obtained by the method of steps on time intervals [a, b]. Similar to the Theorem 4.1, it can be shown that the resulting solution is a non-negative continuous vector-valued function over the whole interval [0; T mod ] and continuously differentiable on intervals [a, b] with the derivative at a taken as a the right-hand derivative.
In reality, one need to consider the discrete nature of the virus population and the subsets of infected cells including the productively-and latently infected cells. Let us assume that the variable V (t) is such that V (t) < 1 for t in [t 1 ; t 2 ), 0 < t 1 < t 2 . Then, the inequality 0 ≤ V (t) < 1, t ∈ [t 1 ; t 2 ) can be interpreted as representing the state of the infected host with the size of the virus population being zero over the time interval [t 1 ; t 2 ).
Consider now the case when the solution of model (1)- (7), (19), (34) is characterized by the following features: there exist times 0 < τ 0 < τ , such that We will regard the fulfilment of inequalities (35) as representing in the model the situation of a complete eradication of HIV-1 infection at time τ 0 . The nonlinearities of the equations and the presence of multiple delays admits the existence of the model solutions which satisfy the conditions (35) for some τ 0 , τ . This conjecture takes into account that for R 0 < 1 the solution can enter the attraction domain of the asymptotically stable trivial steady state. Besides, the expression for R 0 does not depend on the parameters of the model which define the dynamics of immune response components E 0 , E, Q and their interaction with target cells A, I. Therefore, it follows that for some values of R 0 including the case of R 0 ≥ 1, the system of model equations can have other (nontrivial) steady states. The asymptotic stability (unstability) of these complementary steady states essentially depends on the parameters which underly the functioning of the immune system. Various possible behaviours of the model solutions in the neighborhood of these nontrivial steady states can result in the fulfilment of inequalities (35) for some values of τ 0 , τ . 6. Computational experiments. As the conditions (35) underlying the eradication of HIV-1 infection are impossible to derive analytically, we examined the issue by performing a series of computational experiments aimed to find out the possibility of infection elimination over some fixed interval [0; T mod ]. To this end the initial value problem with data in the form of (34) was considered. To solve the intial value problem numerically, a semi-implicit Euler method with a constant time step h was used. The value of the integration step h was taken to be a multiple of the delays ω C , ω U , ω E1 , ω E2 . The model solution was advanced using the Euler scheme at the mesh points t j = jh ∈ [0; T mod ), whereas the values of the delayed variables at times were defined either by the initial data (34) or by the already computed values at earlier mesh points with respect to t j .
The computational experiments were performed to study in particular the dynamics of the following variables of the model y 1 (t) = ln(V (t) + U (t) + 1), y 2 (t) = ln(I(t) + C(t) + 1), y 3 (t) = ln(L(t) + 1). (36) The value of every component in (36) was compared to ln 2. It is straightforward to see that if the inequalities y i (t) < ln 2, i = 1, 2, 3 are satisfied then The examination of the dynamics of functions indicated in (36) provides the mean to check the fulfilment of inequalities (35) and estimate the respective value of τ 0 .
The duration of the integration interval was T mod = 150 days and the integration stepsize was h = 2 · 10 −4 day. In computational experiments, the dependence of the solution on the dose of infection V 0 and the parameters α L , γ A,V , γ A,I , γ Q,A,V was examined. Table 1 presents the times τ 0 (day) at which the inequalities (35) hold, i.e. the eradication of HIV-1 infection takes place, for various infection doses, or does not happen (-), i.e. the τ 0 does not exist. The behavior of the model variables y 1 (t), Table 1  place for lower virus doses V 0 as compared to the scenario of Experiment 3, e.g., for V 0 ≥ 10 2 .
6.5. Experiment 5. The simulated scenario is similar to the one of Experiment 3, except for the values of parameters γ A,V and γ A,I for which larger values are used: γ A,V = 0.7 · 10 −6 , γ A,I = 1.5 · 10 −6 . In this case we have R 0 = 13.2135. Table 1 and Figure 6 show that the eradication of infection takes place only for the high virus dose V 0 = 10 6 .
6.6. Experiment 6. The parameters of the model take the reference values, except for γ A,V , γ A,I , α L . We consider γ A,V = 0.7 · 10 −6 , γ A,I = 1.5 · 10 −6 to be the same as in Experiment 5 except for α L = 0.2. This set of parameters results in R 0 = 13.3789. The increase of α L with respect to the reference value can be interpreted as an engagement at a higher rate of the latently infected cells in the spreading of HIV infection, as a consequence, the infection dynamics is different compared to the previous experiment 5. The results presented in Table 1 and Figure 7 show that an increase of α L can lead to the dynamics reminiscent of infection eradication in a relatively short time after initial infection in contrast to the behavior observed in Experiment 5. The presented model is formulated in the form of a system of delay differential and integral equations which take into account the existence of distinct stages in the intracellular virus replication cycle, the infection fates of target cells, and the clonal immune response. Theorem 4.1 proofs that the mathematical model is wellposed. The model admits the existence of a trivial virus free steady state with its stability being determined by the basic reproduction ratio as stated in Theorems 4.2 and 4.3. In combination with the behavior of the model solutions revealed by computational experiments which resembles some basic well-known features of HIV infection dynamics, the model can be considered as an adequate analytical tool to gain additional insight into the pathogenesis of HIV-1 infection.
The suggested formalization of the conditions for the eradication of HIV infection takes into account that when conditions (35) are fulfilled then the values of the model variables representing the infection components become close to zero. These small values are preserved over some time interval. Therefore, it is reasonable to guess that the real populations of viruses and infected cells would extinct during this time interval. To justify this view one needs to reformulate the model in a stochastic form with discrete variables. With the stochastic models, the states with zero populations of viruses and infected cells can be treated as absorbing states for which the probability of entering these states can be computed. If the probability of getting in absorbing state for some final time is one, then the respective solution of the stochastic model can be considered as HIV-1 infection eradication. Examples of the stochastic models of HIV-1 infection dynamics are given in [26], [23]. The use of the deterministic and stochastic formulations of the HIV infection model sets a step towards the development of hybrid mathematical models.
The results of our computational experiments with the model calibrated according to the available estimates of parameters suggest that the complete elimination of HIV-1 infection in a relatively short time post infection of a healthy individual is feasible. The respective values of the model parameters for which the eradication takes place, suggest that a strong CTL response is a necessary condition for eradication. In addition, the reversal of the latency in cells harboring the virus resulting in their transformation into productively infected cells is another key requirement. When both the clonal expansion of HIV-specific CTL and the activation of latently infected cells are below some threshold level then the eradication of HIV-1 infection after primary acute phase is not possible.
In conclusion, the limitations of the modelling results need to be discussed. The eradication of HIV-1 infection requires consideration of the following additional factors and processes: the virus mutations leading to escape variants, the heterogeneity of the target cells, the spatial organization of the immune and lymphatic system in humans, the application of antiretroviral drugs, etc. The role of these factors in the pathogenesis of HIV-1 infection is highly specific and context dependent and one needs to consider them in extended equations of the mathematical model. However, as far as the primary initial phase of infection is considered, not all of the above listed factors affect considerably the infection dynamics. Additional analysis is need to better understand the control factors operating during distinct stages of HIV-1 infection. Finally, in order to eradicate most efficiently the infection the powerful mathematical methods of optimal control of the system dynamics and the clinical data-based parameter estimation are required to be a part of the modelling-guided studies. Some progress made in this direction can be found in [1], [4], [3], [11], [10].