Modelling immune system dynamics : the interaction of HIV and recombinant virus

We investigate the dynamical behavior of a mathematical model of HIV and recombinant rabies virus (RV), designed to infect only the lymphocytes previously infected by HIV. This model is described by five ordinary differential equations with two discrete delays. The effect of two time delays on stability of the equilibria of the system has been studied. Stability switches and Hopf bifurcations when time delays cross through some critical values are found. Numerical simulations are performed to illustrate the theoretical results.


Problem formulation
According to Word Health Organization, "there were approximately 35 million people living with HIV at the end of 2013" and it "continues to be a major global public health issue, having claimed more than 39 million lives so far" [10].While antiretroviral therapy is progressing every year, there is still a need to improve life expectancy and quality of HIV infected patients.Lately the concept of a recombinant rabies virus (RV), which would be genetically modified to infect only the lymphocytes previously infected by HIV and thus could prevent further HIV production, was introduced [4].
Here we examine a delay differential equations (DDE) model with two discrete delays, representing HIV and RV dynamics in vivo.
The basic HIV model contains three ordinary differential equations (ODE) representing the dynamics of healthy and infected CD4 + T lymphocytes and HIV load.The model and motivation behind its analysis are briefly yet informatively reviewed in [5].The paper also covers the process of introducing the antiretroviral drugs into the system.Revilla and García-Ramos [6] suggested another way of expanding the basic model by adding the recombinant virus and double (HIV and RV) infected CD4 + T lymphocytes.In the article the basic reproduction numbers for HIV and RV were given as well as numerical simulations of various cases.The local and global asymptotic stability of the expanded system's equilibria was analyzed in [2].Tian et al. [9] incorporate in the system a single delay τ , which represents the time needed for virus to infect a healthy lymphocyte, and the exponential probability that the lymphocyte will still be infected after the delay period: here x(t), y(t) and z(t) are healthy, infected and double (by both HIV and RV) infected CD4 + T lymphocytes, respectively, v(t) and w(t) are the HIV and RV load.Parameters λ and d are healthy lymphocytes reproduction and clearance rates, k and p are HIV production and clearance rates, c and q are RV production and clearance rates.β and α are the HIV and RV infection rates, while a and b are the clearance rates of single and double infected lymphocytes, respectively.The nonnegativeness of the solution and the local and global stability of equilibria are proven in [9].Also the Hopf bifurcation conditions and the critical delay value at which it occurs are given.
As mentioned, here we investigate five DDE model with two discrete delays: τ 1 correspond to delay τ in (1) and τ 2 is the time period in which virus is produced and released from the infected lymphocyte.The model is biologically meaningful: the nonnegativeness and boundedness of the solution can be easily proven using basic DDE theory and Theorem 2.1 in [7].The main purpose of this paper is to study the influence of the time delay τ 2 on the dynamics of (2).

Stability investigation
Calculating the equilibria points we get the same expressions as given in [2,6]: E 0 -the situation when no virus is present in the system, E 1 -the case when only HIV survives and E 2 -when both HIV and RV survive.Following [1], we calculated the single HIV infection reproduction number R 1 := βkλ adp , which also agrees with the expression obtained in [2,6].It is easy to see that E 1 is biologically meaningful (nonnegative) ⇔ R 1 > 1. Similarly it can be shown that E 2 is biologically meaningful ⇔ double HIV and RV infection reproduction number , where R * := 1 + bqβk αcdp .Procedure of finding local and global asymptotic stability intervals of equilibria E 0 , E 1 and local asymptotic stability interval of E 2 is given in [2,9].However, because of the second delay, linearization about the equilibria points and corresponding characteristic equations have additional members and the stability proofs have to be expanded accordingly.We are going to analyze only the periodic nature of changes between local asymptotic stability of the equilibrium E 2 of system (2) and occurring Hopf bifurcation when variating the two delays.Linearizing the system (2) about the point Then in order to find a purely imaginary conjugate eigenvalue ξ = ±i̟, ̟ ∈ R, we define the real (R(̟, θ)) and imaginary (̟I(̟, θ)) parts of characteristic equation (F (i̟, θ)), corresponding to (3), where and τ := τ 1 +τ 2 , θ := ̟τ , constants A i , B j , i = 0, 4, j = 1, 2, 3, depend on system's (2) parameters, except the two delays τ 1 , τ 2 .The presence of periodic trigonometric functions in (4), ( 5) equations is common for DDE.Due to them, solving R(̟, θ) = 0, I(̟, θ) = 0, we get infinitely many solutions at which the stability changes between asymptotically stable equilibrium point and Hopf bifurcation resulting in periodic oscillations.According to our notation, every such solution gives aggregated delay and, as we can see from Fig. 1, these values of τ can again give infinitely many combinations of two delays (τ 1 , τ 2 ).It means that in every stability region we can have a case that is mathematically invariant, but, due to different delay values, may be biologically distinct.

Numerical results
In this chapter numerical simulations are given to illustrate the theoretical results.In order to compare obtained data, we fix the parameters to the same values as in [9]: , and choose delays τ 1 , τ 2 as bifurcation parameters.Then equations (4), (5) are and graphically solving R(̟, θ) = 0, I(̟, θ) = 0 we get Fig. 2. The values of the first few stability changing points are given in Table 1.In [3] the mean intracellular delay is estimated to be 0.92 ± 0.43 days, so only the first stability changing point could be biologically meaningful, while others may be purely mathematically interesting.We will check Hopf bifurcation conditions (see i.e. [8, 90 p.]) at the first point (̟ 1 , τ 0 1 ) corresponding with (θ 0 1 , ̟ 1 , ): • the partial derivative of characteristic equation is not equal to zero: Fig. 3. Hopf bifurcations resulting in periodic oscillations (blue, solid) and locally asymptotically stable equilibrium E 2 (green, dashed).The two bottom graphs are corresponding phase portraits: a limit cycle and a stable focus.Dotted red color marks the results when simulating system (2) without the RV.
So, system (2) undergoes Hopf bifurcation at the critical value of aggregated delay (6) τ 0 1 = 0.5942700703 (see Fig. 3), which is smaller than the critical Hopf bifurcation value obtained in [9] for a system with single delay (1).Apart from Hopf bifurcation, two other cases, obtained with the same values of delays τ 1 , τ 2 , are shown in Fig. 3: locally asymptotically stable equilibrium E 2 and simulation of system (2) without the recombinant virus.When comparing them, it can be seen that a system without recombinant virus shows smaller numbers of healthy CD4 + T lymphocytes and larger numbers of both infected CD4 + T lymphocytes and HIV.

Discussion and conclusions
Even a simplified model with two delays indicates interesting dynamics from both mathematical and biological standpoint.Periodic stability switches could be further investigated, comparing the amplitude and period of oscillations, which rise from Hopf bifurcation.The direction of such bifurcation could also be proved.
Our results suggest that the behavior of the system (2) is defined by the aggregated delay rather than by delays τ 1 and τ 2 separately.And a simple numeric example comparing the count of healthy and infected lymphocytes and HIV agrees with the theoretical results given in [2] for a system without delays.

Table 1 .
The first four stability changing points (̟, θ) and the corresponding value of τ .