Viral infection model with diffusion and state-dependent delay: stability of classical solutions

A class of reaction-diffusion virus dynamics models with intracellular state-dependent delay and a general non-linear infection rate functional response is investigated. We are interested in classical solutions with Lipschitz in-time initial functions which are adequate to the discontinuous change of parameters due to, for example, drug administration. The Lyapunov functions technique is used to analyse stability of interior infection equilibria which describe the cases of a chronic disease.


Introduction
In our research we are interested in mathematical models of viral diseases. 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. Particularly, in the recent The Global hepatitis report (WHO, April 2017) we find [40] "a large number of people -about 325 million worldwide in 2015 -are carriers of hepatitis B or C virus infections, which can remain asymptomatic for decades." and "Viral hepatitis caused 1.34 million deaths in 2015, a number comparable to deaths caused by tuberculosis and higher than those caused by HIV." In such a situation any steps toward understanding viral diseases are important.
There are variety of models with and without delays which describe dynamics of different viral infections. Delays could be concentrated or distributed, constant, time-dependent or statedependent.
We notice that classical models [19,22] contain ordinary differential equations (without delay) for three variables: susceptible host cells T , infected host cells T * and free virus particles V . The intracellular delay is an important property of the biological problem, so we formulate the delay problem       Ṫ (t) = λ − dT (t) − f (T (t), V (t)), V (t) = N δT * (t) − cV (t).
( 1) In (1), susceptible cells T are produced at a rate λ, die at rate dT , and become infected at rate f (T, V ). Properties and examples of incidence function f are discussed below. Infected cells T * die at rate δT * , free virions V are produced by infected cells at rate N δT * and are removed at rate cV (t). In (1) h denotes the delay between the time a virus particle contacts a target cell and the time the cell becomes actively infected (start producing new virions). It is clear that the constancy of the delay is an extra assumption which essentially simplifies the analysis, but has no biological background.
To the best of our knowledge, viral infection models with state-dependent delay (SDD) have been considered for the first time in [29] (see also [30]). It is well known that differential equations with state dependent delay are always non-linear by its nature (see the review [9] for more details and discussion).
As usual in a delay system with (maximal) delay h > 0 [8,13,5], for a function v(t), t ∈ The ODEs delay system (1) is extended to the state-dependent one T * (t) = e −ωh f (T (t − η(u t )), V (t − η(u t ))) − δT * (t), Here u(t) = (T (t), T * (t), V (t)). System (2) is a particular case of the system with state-dependent delay studied in [29,30]. The ODE system is formulated assuming host cells do not move and the diffusion of free virus particles is very quick, so they are mixed enough to consider homogeneous distribution over the spatial domain in a host organ. Similar situation is in case of all cells and free virions are well mixed (e.g., in case of HIV and other infections targeting blood cells). To consider more realistic nonhomogeneous situation one introduces spatial coordinate x ∈ Ω and allow the unknowns to depend on it, i.e. Consider a connected bounded domain Ω ⊂ R n with a smooth boundary ∂Ω. Now we are ready to present the PDEs system under consideration Here the dot over a function denotes the partial time derivative i.g,Ṫ (t, x) = ∂T (t,x) ∂t , all the constants λ, d, δ, N, c, ω are positive while d i , i = 1, 2, 3 (diffusion coefficients) are non negative. We consider a general functional response f (T, V ) satisfying natural assumptions presented below. In earlier models (with constant or without delay) the study was started in case of bilinear f (T, V ) = const · T V and then extended to more general classes of non-linearities, see Remark 2 below.
Boundary conditions are of Neumann type for the corresponding unknown if d i = 0 i.e.

∂T (t,x)
∂n | ∂Ω = 0 if d 1 = 0 and similarly for T * (t, x) and V (t, x). Here ∂ ∂n is the outward normal derivative on ∂Ω. In case d i = 0, no boundary conditions are needed for the corresponding unknown(s).
Our main goals are to present the existence and uniqueness results for the model (3) in the sense of classical solutions, and to study the local asymptotic stability of non-trivial diseased equilibria.
We apply the Lyapunov approach [14] to the state-dependent delay PDE model and allow, but not require, diffusion terms in each state equation.
There is a number of works studying the case d 1 = d 2 = 0, d 3 > 0 (see e.g. [38,37,39] for models without delay and [17,10] with constant delay; see also references therein). In the mentioned works authors assume that the host cells (healthy and infected) do not move or are well mixed, while viral particles diffuse freely. Let us discuss the cases when an infection affects one particular organ as, for example, liver in case of HBV, HCV. In such cases the spatial domain Ω ⊂ R 3 represents the organ. The Neumann boundary conditions say that viral particles do not leave the organ. It is not relevant from the biological point of view since viral particles circulate together with the blood stream in and out the organ (e.g. liver). For the mathematical system to cover such cases one could assume d 3 = 0 and no boundary conditions for V . Taking into account the high speed of the blood stream, this means the viral particles are well mixed. Even more interesting case is d i > 0, i = 1, 2, d 3 = 0. To the best of our knowledge, this case has not been considered before. The case d 2 > 0 may reflect the cell-to-cell transmission of the infection when viral particles cross the membranes of the nearest cells (see [2] for more discussion and references; c.f. [39]). The infection spreads similar to diffusion to cells in a neighbourhood of an infected cell.
The case d 1 > 0 may reflect natural division of healthy cells in order to fill the space previously occupied by infected cells (after the death of the last ones). In cases d 1 > 0, d 2 > 0, the host cells (both healthy and infected) do not leave the organ, so Neumann boundary conditions are quite relevant.
In study of state-dependent delay equations the choice of the set of initial functions is particularly important and non-trivial (see review [9] for ODE case and works [23,24,25,3] for PDEs).
We are interested in classical solutions with Lipschitz in-time initial functions which are adequate to the discontinuous change of parameters due to, for example, drug administration (for more discussion and references see [30]). The main motivation here is the situation (see e.g. [31,20]) when the drug effectiveness is decreased in a stepwise manner. In terms of system (3), the parameter N could change its value in a discontinuous way (see equation (2) in [31, p.920]). It is clear that at any time moment of discontinuity of (any) parameter, the solution is continuous, but not differentiable (c.f. figure 2-B in [31, p.921] and also fig.1 in [20, p.23]).
Since delay is a central part of the paper, it would be interesting to present examples of SDD η and discuss the structure of η from biological point of view. Unfortunately, up to now, the biological side of virus dynamics is not fully understood. Even current in vitro study does not provide enough information. In vivo study is essentially more complicated, and up to now, there are no technical (biological/medical) tools for the real time monitoring of disease dynamics available. In such a situation we present a rather general class of SDD (see (28), (29) below). Delays of the form (28), (29) take into account all the prehistory u t by integrating a solution over [t − h, t].
For general facts on PDEs with constant delay see e.g. [34,16,41] and PDEs with statedependent delay [23,24,25,26,27,28,3]. We also mention that the case of all d i > 0 is, in a sense, easier from mathematical point of view since the linear part generates a compact semi-group.
We use the Lyapunov functions technique [14] to analyse stability of interior infection equilibria which describe the cases of chronic disease. To the best of our knowledge, viral infection models with diffusion and state-dependent delay have not been considered before.

Basic properties of the model
We denote the space of continuous functions by We write, the system (3) in abstract form The non-linear continuous mapping F : C → X is defined by Here Mapping F is not Lipschitz on the space C which is typical for a mapping which includes discrete state-dependent delays (see review [9] for ODE case and works [23,24,25,3] for PDEs).
We need initial conditions u(θ,  (4), (6) if it satisfies (6) and (4) is satisfied on In the study below we are mainly interested in classical solutions which preserve the regularity of the Lipschitzian initial data (see (6)).
Assume the non-linear function f : R 2 → R is Lipschitz continuous and satisfies We have the following result Proposition 1 Let nonlinear function f be Lipschitz and satisfy (Hf 1 ) (see (8)), state-dependent Then the initial value problem (4), (6) has a unique classical solution which is global in time i.e. defined for all t ≥ 0.
Proof of Proposition 1. We start with discussion of mild solutions. Since the semigroup generated by the linear part −A is not necessarily compact (in cases when at least one constant d i = 0), see e.g. [16], we cannot directly use results of [23,24,25,3]. On the other hand, as mentioned above, non-linearity F is not Lipschitz on C, so we cannot directly apply the existence result of [16]. Moreover, the extension provided in [27] cannot be directly applied to our case since we do not assume here the ignoring condition on the state-dependent delay (see more details in [24,26,27]). Nevertheless, the restrictions on initial function ϕ posed by (6) give the possibility to prove the existence of a (unique) mild solution to initial-value problem (3), (6) using the standard line based on Banach Fixed Point Theorem (in a complete metric space) as in the ODE case. We outline only main steps of the proof. First we consider the following . It is not difficult to check that our non-linear mapping F, defined by (5), satisfies (see the estimate |f (T, V )| ≤ µ|T | in (Hf 1 )) ||F (ψ)|| X ≤ n 1 + n 2 ||ψ|| C and is locally almost Lipschitz on A(α, β, γ) by the terminology of [15].
Proof of Proposition 2. The existence and uniqueness of solution is proven in Proposition 1.
The proof of the invariance part follows the invariance result of [16] with the use of the almost to notice that the solutions are classic for all t ≥ 0 (but not for t ≥ h as could be in the case of merely continuous initial functions ϕ ∈ C). The proof of Proposition 2 is complete.

Stationary solutions
Let us discuss stationary solutions of (3). By such solutions we mean time independent u which, in general, may depend on x ∈ Ω. Consider the system (3) with u(t) = u(t − η(u t )) = u and denote the coordinates of a stationary solution by ( T , T * , V ) = u ≡ ϕ(θ), θ ∈ [−h, 0]. Since stationary solutions of (3) do not depend on the type of delay (state-dependent or constant) we have (see e.g.

Remark 1
We notice that the finiteness of roots (which are obviously isolated) does not allow the existence of equilibria which depend on spatial coordinate x ∈ Ω. We remind that Ω is a connected set, so a function v ∈ C(Ω) may take either one or continuum values. Assumption (Hf 2 ) implies 1+k1T +k2V , with k, k 1 ≥ 0, k 2 > 0. We also mention that the functional response includes as a special case (k 1 = 0) the saturated incidence rate f (T, V ) = kT V 1+k2V . Another example of the nonlinearity is the Crowley-Martin incidence rate f (T, V ) = kT V (1+k1T )(1+k2V ) , with k ≥ 0, k 1 , k 2 > 0 (see e.g. [42]). For more general class of functions f see, e.g. [17,10,30], where under additional conditions, one has exactly one root of h f (s) = 0. We notice that, in contrast to [17,10], we do not assume here the differentiability of f .

Remark 3
It is important to mention that usually in study of stability properties of stationary solutions (for viral dynamics problems) one uses conditions on the so-called reproduction numbers.
These conditions are used to separate the case of a unique stationary solution. Then the global stability of the equilibrium is investigated. In our study, taking into account the state-dependence of the delay, we discuss the local stability. As a consequence, it allows the co-existence of multiple equilibria. We believe this framework provides a way to model more complicated situations with rich dynamics (in contrast to a globally stable equilibrium). The conditions on the reproduction numbers do not appear explicitly here, but could be seen as particular sufficient conditions for (Hf 2 ).

Stability of disease stationary solutions
The following Volterra function v(s) = s − 1 − ln s : (0, +∞) → R + plays an important role in construction of Lyapunov functionals [12,17]. One can see that v(s) ≥ 0 and v(s) = 0 if and only if s = 1. The derivative equalsv(s) = 1 − 1 s , which is obviously negative for x ∈ (0, 1) and positive for x > 1. The graph of v explains the use of the composition v s s 0 in the study of the stability properties of an equilibrium s 0 . Another important property is the following [29] estimate To check it, one simply observes that all three functions vanish at s = 1 and d ds in the µ-neighbourhood of s = 1.
In this section we use the following local assumptions on f in a small neighbourhood of a disease equilibrium (given by (Hf 2 )).
This property simply means that the value f (T,V ) f (T, V ) is always strictly between 1 and V V for any T ≥ 0 (c.f. with the non-strict property [17, p.74]). The strict inequality in (15) will be needed to handle the state-dependence of the delay. In the particular case of constant delay, the non-strict property is enough.
We will also use the following assumption (Hf 4 ) Function f is either differentiable with respect to its first coordinate or satisfies For simplicity of presentation we start with stability analysis for smooth initial data belonging to the so-called solution manifold (see e.g. [35,9] for ODE case and [28] for PDEs) The equation in (17), called the compatibility conditions, is an equality in X. Below (see Theorem 3.2) we return to more general case of Lipschitz initial functions (ϕ ∈ Ω Lip , not necessarily continuously differentiable) which are important to cover the cases of drug administration when the time derivative may be discontinuous, see [30] for more discussion. Remark 4 Similar to ODE case, described in [29,Remark 13], we have the following property.
Proof of Theorem 3.1. Let us consider (point-wise) the following auxiliary functional Now we can introduce the following Lyapunov functional with state-dependent delay along a solution of (3) The form of the functional is standard except the low limit of the last integral in (19) which is state-dependent. This state-dependence was first considered in [29] (see also [30]). For the constant delay case, see e.g. [17].
Now, for the simplicity of presentation, we consider the point-wise time derivative of the functional U sdd−x (t, x) defined in (19). This time derivative is considered along classical solutions of (3). It gives the possibility to consider ∂T (t,x) ∂t , for any t > 0. The computations below are in a sense close to the ones in [17], but here we have two additional diffusion terms and the state-dependence in both the system (3) and the Lyapunov functional. First we consider the where we denoted for short Remark 5 The term S sdd appears due to the presence of the state-dependent delay. It makes the technical calculations more challenging. The sign of S sdd is undefined, so we propose below (see also [29,30]) a way to compensate/bound S sdd by other positive defined terms in ∂U sdd−x ∂t to have the time derivative of the Lyapunov functional (along a solution) negative defined relative to the equilibrium.

Now we differentiate
Calculations, using (12), particularly, where, for short, we collected some terms as C 1 . It is written as follows Calculations show that In the above expression we see two positive and five negative fraction terms, so we write 3 = −2 + 5 and add the following zero term (0 = ln 1): which is split on the sum of seven logarithms to write shortly, using the Volterra function v As before, we use v(s) = s − 1 − ln s. Now we discuss the diffusion terms (the ones with coefficients d i ) in (22). More precisely, we are interested in the sign of these terms after integration by x in Ω. Denote them, for short, as In case of differentiable f (Hf 4 ) (see (16)) need the following simple for any u ∈ C 2 (Ω) satisfying ∂u(x) ∂n | ∂Ω = 0.
Proof of Proposition 3. We use the classical Gauss-Ostrogradsky theorem. Consider the vector field E ≡ p(u)∇u. Hence div E = p ′ (u)||∇u|| 2 + p(u)∆u. One has Ω div E dx = Ω p ′ (u)||∇u|| 2 dx + Ω p(u)∆u dx = ∂Ω p(u)(∇u, n) dS = 0. The last equality due to the Neumann boundary conditions. Finally, Ω p(u)∆u dx = − Ω p ′ (u)||∇u|| 2 dx. It completes the proof of Proposition 3. Now we apply Proposition 3 to show that D diff−3 (t) ≤ 0. Let us start with the first term in x), see (24), and show that ≥ 0 by the assumption on f . Similar considerations with the second and third terms in (24) show that
Now we combine the arguments above to study the Lyapunov functional U sdd (t), see (20). We have the following equality (c.f. (22)) Here D diff−3 (t) is defined in (24) and transformed in (25), C 1 is presented as in (23) and S sdd is defined in (21). We remind (see (12)) that δ T * = e −ωh f ( T , V ) which leads to cancellation of the first and second terms in the last integral with the corresponding terms in C 1 (see (23)). We continue calculations We will show that all the terms in (26) are non-positive except for the last one which, in general, may change sign. The first term in (26) is non-positive due to monotonicity of f with respect to the first coordinate. The property D diff−3 (t) ≤ 0 is given in (25). To show that we use the property (Hf 3 ) of function f (see (15)).
Now we plan to prove that d dt U sdd (t) ≤ 0 in a small neighbourhood of the stationary solution with the equality only in case of (T, T * , V ) = ( T , T * , V ). In the particular case of constant delay, one has S sdd (t, x) = 0 which may lead to the global stability of ( T , T * , V ).
Moreover, D diff−3 (t) = 0, means (see Proposition 3 and (25)) that T, T * and V are independent of x. The zero set S sdd (t, x) = 0 is described (see (21) by f (T (t − η(u t ), x), V (t − η(u t ), x)) = f ( T , V ) or d dt η(u t ) = 0 along a solution. It is important for us that the zero-set D sdd (t, x) = 0 is a singleton (T, T * , V ) = ( T , T * , V ) and is a subset of S sdd (t, x) = 0. The rest of the proof that in a small neighbourhood of ( T , T * , V ) one has |S sdd (t, x)| < D sdd (t, x) follows the streamline of the proof [29,Theorem 12] (see also [30,Theorem 3.3]). It relies on property (18), auxiliary quadratic functionals due to property (14) of Volterra function v and the change of variables to the polar ones (see [29, (33)-(35)]). We do not repeat the calculations here. The property d dt U sdd (t) ≤ 0 in a small neighbourhood of the stationary solution with the equality only in case of (T, T * , V ) = ( T , T * , V ) completes the proof of Theorem 3.1.
It is interesting to notice that ϕ ∈ M F (see (17)) is not a necessary condition for our approach. Now we consider a wider set Ω Lip (see (10)). Let us discuss a particular simple form of the delay Hence, in the ε-neighborhood of the stationary solutionû, one has d dt η(u t ) ≤ |ξ(u(t)) − ξ(u(t − h))| ≤ 2ε L ξ,ε ≡ α ε → 0 as ε → 0.
The discussion above shows that property (18) ) and state-dependent delay η : C → [0, h] be of the form (29). Then the stationary solution ϕ is locally asymptotically stable.