Stability of a viral infection model with state-dependent delay, CTL and antibody immune responses

A virus dynamics model with intracellular state-dependent delay and nonlinear infection rate of Beddington-DeAngelis functional response is studied. The technique of Lyapunov functionals is used to analyze stability of an interior infection equilibrium which describes the case of both CTL and antibody immune responses activated. We consider first a particular biologically motivated class of discrete state-dependent delays. Next, the general case is investigated.


Introduction
We are interested in mathematical models of infectious diseases. The diseases are caused by pathogenic microorganisms, such as bacteria, viruses, parasites or fungi. According to World Health Organization, many viruses (as Ebola virus, Zika virus, HIV, HBV, HCV and others) continue to be a major global public health issues.
In our research we concentrate on models of viral infections. There are variety of models with and without delays which describe dynamics between different viral infections and immune responses. Delays could be concentrated or distributed. We will not describe the historical evolution of such models, just mention that early models [13,12] contain three variables: susceptible host cells, infected cells and free virus. Next step was to take into account, as written in [25], that "one particular part of the immune system that is very important in the fight against viral infections are the killer T cells or cytotoxic T lymphocytes (CTL)." See also [29] and references therein. There is another adaptive immune response by antibodies. The relative balance of both types of immune response "can be a decisive factor that determines whether patients are asymptomatic or whether pathology is observe" [24]. These lead to introduction of two additional variables of both immune responses [24,25] (see also [27] and references therein).
The model under consideration contains five variables: susceptible (noninfected) host cells T , infected cells T * , free virus V , a CTL response Y , and an antibody response A. In case of bilinear nonlinearities and one constant concentrated delay (see, for example, [26]) it has the following form Here the dot over a function denotes the time derivative i.g,Ṫ (t) = dT (t) dt , all the constants λ, d, k, δ, p, N, c, q, β, γ, g, b, ω are positive. As for the immune responses, the fourth equation describes the regulation of CTL response with pY (t)T * (t) (in the second equation) being the rate of In the above model, the standard bilinear incidence rate is used according to the principle of mass action. For more details and references on the models of infectious diseases with more general types of nonlinear incidence rates see e.g. [9,5]. In paper [28], following [8,22,26], authors assume that the infection rate of the virus dynamics models is given by the Beddington-DeAngelis functional response [1,2], kT V 1+k1T +k2V , where k, k 1 , k 2 ≥ 0 are constants. The Lyapunov asymptotic stability [11] of points of equilibrium is studied for the following model with constant The functional response was introduced by Beddington [1] and DeAngelis et al. [2].
It is evident that the constancy of the delay is an extra assumption which essentially simplifies the analysis, but is not motivated by the biological background of the model. It was a reason (see e.g. [23]) to discuss distributed delay models as an alternative to discrete constant delay ones. We propose another approach.
Our goal is to remove the restriction of the constancy of the delay and investigate the wellposedness and Lyapunov stability of the following virus infection model with Beddington-DeAngelis functional response and state dependent delay. It appears that the analysis essentially differs from the constant delay case. To the best of our knowledge, such models have not been considered before. It is well known that differential equations with state dependent delay are always nonlinear by its nature (see the review [7] for more details and discussion).
As usual in a delay system with (maximal) delay with the functional response f (T, V ) given by (3). We denote by F the right-hand side of (4) to write the system shortly asu(t) = F (u t ).
The paper is organized as follows. In Section 2, we discuss and choose different sets of initial data and prove the existence and uniqueness of solutions. Next we prove that the sets are invariant.
Section 3 is devoted to the stability properties of a stationary solution. We study the stability of an interior equilibrium which describes the case when both CTL and antibody immune responses are activated. We believe this infection equilibrium is only biologically meaningful in the study of the disease. We use the technique of Lyapunov functionals [11] and consider first a particular (biologically motivated) case of state-dependent delay. We complete the paper by investigation of the general case.

Preliminaries
We first study the basic questions of the existence and uniqueness of solutions to the problem (4).
Since two functions T and V are used in (4) at different time moments (current time t and delayed time t − η(u t )), we should consider initial values T (θ), V (θ) for θ ∈ [−h, 0]. As usual for such a biological system, one should check the non-negativeness and boundedness of all the coordinates provided initial values are non-negative.
We will study the system (4) with an initial function where R + ≡ [0, +∞), We consider the following assumption on the state-dependent delay Remark 1. We notice that even more restrictive assumption η(ψ) > 0 for all ψ ∈ C + is biologically well motivated. On the other hand, even this restriction (the so-called "non-vanishing delay") does not guarantee the uniqueness of solutions with merely continuous data (see [4]).
The first result is the following (ii) If additionally, η satisfies (H1 η ), then for any initial function ϕ ≡ (T 0 , T * 0 , V 0 , A 0 , Y 0 ) ∈ Ω C such that T 0 , V 0 are Lipschitz functions, the problem (4), (5) has a unique solution. The solution is globally Lipschitz in time and satisfies Proof of Theorem 2. (i) The existence of continuous solutions is guaranteed by the continuity of the right-hand side of (4) and classical results on delay equations [6,3].
(ii) Since T 0 , V 0 are Lipschitz functions, the uniqueness of continuous solutions follows from the general results on differential equations with state-dependent delay (see the review on ordinary equations [7] for details and references and also [14,15,16,18] for PDEs)). Let us show that the set Ω C is invariant i.e. any solution starting from ϕ ∈ Ω C remains in Ω C for all t ≥ 0.
We notice that in the case of constant delay the non-negativeness of all coordinates of a solution follows from the quasi-positivity property of the right-hand side of (4) (see e.g. [ (4). We could use the corresponding extension to the state-dependent delay case [17], but we propose another way here.
To prove the non-negativeness of all coordinates of a solution we use the direct analysis of each coordinate. It is easy to see that T (t) → 0+ impliesṪ (t) → λ > 0 which makes impossible for T to become negative. The direct integration shows that coordinates satisfy Equations (10) show that For the constant delay case, equations (8), (9) would imply the similar result for T * (t), V (t), but for the state-dependent delay we need more care. First, (8) shows the property (for some t 1 ≥ 0) Similarly, (9) gives Now, let us assume that the non-negativity of T * or V falls. Properties (11), (12) show that T * and V should change the sign simultaneously i.e. there exist a (smallest possible) time moment By this property and (8), one has T * (t) ≥ 0 for t ∈ (t 1 , t 1 + δ 2 ). It contradicts the choice of t 1 and completes the proof of the non-negativity of all coordinates.
Let us prove the upper bounds on the coordinates, given in (6). To save the space we formulate an easy variant of the Gronwall's lemma.
Proof of Lemma 3 is easy multiplication of the inequality by e c2t and integration over [a, t]. It It completes the proof of Lemma 3. Since f is non-negative for non-negative arguments (see (3)), we get from the first equation of and Lemma 3 implies the needed bound for T * in (6). The bound for T * and the third equation in (4) giveV (t) ≤ N δT * (t) − cV (t) ≤ N kλ dk2 e −ωh − cV (t). Lemma 3 proves the estimate for V in (6). Next, we use the second and the fourth equations in (4) to geṫ Lemma 3 proves the bound for T * (t) + p β Y (t) in (6). In the similar way, using the third and fifth equations in (4), one getṡ Lemma 3 implies the last estimate in (6). All solutions are global (defined for all t ≥ −h). It completes the proof of Theorem 2.
Remark 4. We notice the our invariant set Ω C differs from the absorbing set Γ, used in [28] for the constant delay system. Let us denote by Ω ε C the set where all the upper bounds in (6) are increased by ε. Then the second part of the Lemma 3 implies that for any ε > 0 the set Ω ε C is absorbing for any solution (not necessary starting in Ω C ). Another difference is that all the five coordinates of ϕ ∈ Ω C are continuous functions in contrast to the constant delay case [28], where the second, fourth and fifth coordinates belong to R + .
Remark 5. It is well known that continuous solutions to differential equations with statedependent delay may be non-unique (see examples in [4]). There are two ways to insure the uniqueness of solutions. The first one is to restrict the set of initial functions [7]. The second one is to restrict the class of state-dependent delays [14,16] and work with continuous initial functions.
If one is interested in continuously differentiable solutions, we could also apply the solution manifold approach [21,7] to the initial value problem (4), (5). We remind the short notation for the system (4) asu(t) = F (u t ). Let us introduce the following subset of Ω C (c.f. (6)) The following result is a corollary of Theorem 2.
One can see that any continuous solution u started at ϕ ∈ Ω C , satisfies u t ∈ Ω F for t > h.

Stationary solutions
We look for stationary solutions of (4). Consider the system with u(t) = u(t − η(u t )) = u and denote the coordinates of a stationary solution by ( T , Since the stationary solutions of (4) are the same as for a constant delay system, we proceed as follows (see e.g. [28]). The last two equations in (14) give T * = γ β , V = b g . Using this, the third equation implies A = N δγg−βcb βqb . To insure the positivity of A we assume the constants in the system (4) satisfy Now we substitute the value V = b g into the first equation of (14) and use (3) to get the quadratic equation for T One can easily see that the discriminant of the above equation is positive and there are two roots of different signs. We denote below by T the unique positive root. The first two equations in (14) give (remind that T * is already known) Y = λ−d T −e ωh δ T * e ωh p T * . The positivity of Y follows from the following assumption where T is the unique positive root of the quadratic equation (16).
We notice that, from biological point of view, (H2), (H3) are standard assumptions on reproduction numbers, which are given here in a short form. We could summarize the above calculations in the following Lemma 7. Let assumptions (H2) and (H3) be satisfied (see (15), (17)). Then the system (14) has a unique solution ( T , T * , V , Y , A) (the unique stationary solution of (4)). All the coordinates are positive, T is the unique positive root of the quadratic equation (16) and coordinates satisfy   These equations connecting the coordinates of the stationary solution will be used below in the study of the stability properties (c.f. (14)).

Stability properties
The equalsv(x) = 1− 1 x , which is evidently negative for x ∈ (0, 1) and positive for x > 1. The graph of v explains the use of the composition v x x 0 in the study of the stability properties of an equilibrium x 0 . Another important property is the following estimate To check it, one simply observes that all three functions vanish at x = 1 and d dx in the δ-neighborhood of x = 1.

Particular case of a state-dependent delay
As before, we denote u(t) = (T (t), T * (t), V (t), Y (t), A(t)). Consider arbitrary ϕ = (ϕ 1 , ϕ 2 , ϕ 3 , ϕ 4 , ϕ 5 ) ∈ C. We are interested in the following particular form of the state-dependent It means that η(u t ) = F (T (t), V (t)) which looks natural since the delay appears in the nonlinearity f which depends on T and V only (see the second equation in (4)).
Proof of Theorem 8. Consider the Lyapunov functional We use the same notations as in [28] to simplify for the reader the comparison of the computations.
In spite of the same Lyapunov functional as in [28], the time derivative of U 1 (t) along a solution u of (4) is different due to the state-dependence of the delay in the system. It reads as follows Opening parenthesis, grouping similar terms and canceling some of them, we obtain To save the space we omit long computations where we intensively used equations (18), for example, For the sum {...} above, the particular form of the function f (see (3)) and computations give Now we add into [...] above and substitute (23) to obtain The first four terms in [...] above suggest to split the logarithm as follows Substitution of (25) into (24) implies Here we used the function v(x) = x − 1 − ln x to save the space. Next, we can rewrite the first term in (26), using (3), We substitute the last equality into (26) to get where One can see, using v(x) ≥ 0, that D 1 (t) ≥ 0. The sign of S 1 (t) is undefined. Our goal is to show that (−D 1 (t) + S 1 (t)) ≤ 0, i.e. d dt U 1 (t) ≤ 0 (see (27)) and d dt U 1 (t) = 0 at the stationary point only. To estimate the term S 1 (t) we notice that the functional response f (T, V ), given by (3), is Lipschitz Remark 10. Both coordinates T (t) and V (t) of a solution u(t) of (4) are Lipschitz in time.
We denote the corresponding Lipschitz constants for arbitrary solution as L T u , L V u . It is easy to see that for any δ-neighborhood of the stationary point ϕ the Lipschitz constants of arbitrary solution We continue, using the assumptions on the state-dependent delay η (see (20) and (21)), This and (29) give the estimate Now we can choose small enough δ (see remark) to make the coefficient e −ωh 1 − T * T * (t) · L f 1 L T,δ + L f 2 L V,δ · c η in (30) arbitrary small, which implies (see the form of D 1 (t) (28)) the desired property d dt U 1 (t) = −D 1 (t) + S 1 (t) < 0. Remark 11. It is easy to see from the calculations above that the small value of |T * (t) − T * | alone gives d dt U 1 (t) < 0. No need to ask the values of |T (t) − T |, |V (t) − V | to be small. On the other hand, the small values of The previous remark suggests that, alternatively, the smaller value of constant c η (see (21) and (30)) the bigger the set where d dt U 1 (t) < 0 holds. The latter implies that in case c η is sufficiently small, the stationary solution is globally stable. Here the LaSalle invariance principle (see e.g. [20, Theorem 5.17, p.80]) is applied. To show that the maximal invariant subset of d dt U 1 (t) = 0 is the stationary solution only, we use the form of d dt U 1 (t) (we have shown that d dt U 1 (t) = 0 iff D 1 (t) = 0). We notice that T (t) = T , T * (t) = T * , V (t) = V follows immediately from (28). Hence the second equation of (4) implies 0 =Ṫ * (t) = e −ωh f ( T , V ) − δ T * − pY (t) T * . Similarly, the third equation of (4) gives 0 =V (t) = N δ T * − c V − qA(t) V . Clearly, (14) implies Y (t) ≡ Y , A(t) ≡ A.
The proof of Theorem 8 is complete.
T * · f (x (4) , x (5) ) The sum of all terms in (35) is needed to bound |S (5) (x)|. To see it, one should compare the sets where each functional vanishes. Denote the zero-sets as Z S (5) and Z rhs (for the right-hand side of (38)). Then one sees that Z S (5) Z rhs . Moreover, in any neighborhood of the point (x (2) , x (4) , x (5) ) = ( T * , T , V ) ∈ R 3 one can find points where the right-hand side of (38) is zero, while the the left-hand side is positive. Clearly, the coordinates of such points should satisfy f (x (4) , x (5) ) = f ( T , V ), T * · f (x (4) , x (5) ) = x (2) · f ( T , V ).