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

TL;DR

This paper analyzes the stability of classical solutions in a reaction-diffusion virus model incorporating intracellular state-dependent delay and nonlinear infection rates, using Lyapunov functions to study interior infection equilibria relevant to chronic diseases.

Contribution

It introduces a reaction-diffusion virus model with state-dependent delay and nonlinear infection response, providing stability analysis for solutions with realistic initial conditions.

Findings

Lyapunov functions establish stability of infection equilibria.

Model accommodates discontinuous parameter changes like drug administration.

Provides insights into chronic disease dynamics.

Abstract

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.

Equations

$$\left\{\begin{array}[]{l}\dot{T}(t)=\lambda-dT(t)-f(T(t),V(t)),\\ \dot{T}^{*}(t)=e^{-\omega h}f(T(t-h),V(t-h))-\delta T^{*}(t),\\ \dot{V}(t)=N\delta T^{*}(t)-cV(t).\\ \end{array}\right.$$

$$\left\{\begin{array}[]{l}\dot{T}(t)=\lambda-dT(t)-f(T(t),V(t)),\\ \dot{T}^{*}(t)=e^{-\omega h}f(T(t-\eta(u_{t})),V(t-\eta(u_{t})))-\delta T^{*}(t),\\ \dot{V}(t)=N\delta T^{*}(t)-cV(t).\\ \end{array}\right.$$

$$\left\{\begin{array}[]{l}\dot{T}(t,x)=\lambda-dT(t,x)-f(T(t,x),V(t,x))+d^{1}\Delta T(t,x),\\ \dot{T}^{*}(t,x)=e^{-\omega h}f(T(t-\eta(u_{t}),x),V(t-\eta(u_{t}),x))-\delta T^{*}(t,x)+d^{2}\Delta T^{*}(t,x),\\ \dot{V}(t,x)=N\delta T^{*}(t,x)-cV(t,x)+d^{3}\Delta V(t,x).\\ \end{array}\right.$$

$$\frac{d}{dt}u(t)+\mathcal{A}u(t)=F(u_t),t>0.$$

$$F(\varphi)=F(\varphi)(x)=\left(\begin{array}[]{l}\lambda-d\varphi^{1}(t,x)-f(\varphi^{1}(t,x),\varphi^{3}(t,x))\\ e^{-\omega h}f(\varphi^{1}(-\eta(\varphi),x),\varphi^{3}(-\eta(\varphi),x))-\delta\varphi^{2}(t,x)\\ N\delta\varphi^{2}(t,x)-c\varphi^{3}(t,x)\\ \end{array}\right).$$

$$\phi \in Lip([-h,0];X)\equiv \left\{\psi \in C:\underset{s\ne t}{\sup }\frac{||\psi (s)-\psi (t)|{|}_X}{|s-t|}<\infty \right\},\phi (0)\in D(\mathcal{A}).$$

$$u(t)=e^{-\mathcal{A}t}\phi (0)+{\int }_0^t{e}^{-\mathcal{A}(t-s)}F(u_s)ds,t\in [0,T).$$

$$(\mathbf{H}{\mathbf{f}}_1)\text{there exists }\mu >0\text{ such that }|f(T,V)|\le \mu |T|\text{for all }T,V\in \mathbb{R},$$

$$\frac{dv(t)}{dt}+\mathcal{A}v(t)=g(t),v(0)=x\in X,$$

$${\Omega }_{Lip}\equiv \{\phi =({\phi }^1,{\phi }^2,{\phi }^3)\in Lip([-h,0];X))\subset C,\phi (0)\in D(\mathcal{A}):0\le {\phi }^1(\theta )\le \frac{\lambda }{d},$$

$$0\le {\phi }^2(\theta )\le \frac{\lambda \mu }{d\delta }{e}^{-\omega h},0\le {\phi }^3(\theta )\le \frac{N\lambda \mu }{dc}{e}^{-\omega h},\theta \in [-h,0]\},$$

$${\bf(Hf_{1}+)}\quad\left\{\begin{array}[]{l}f(T,0)=f(0,V)=0,\quad\text{ and }\quad f(T,V)>0\text{ for all }T>0,V>0;\\ f\text{ is strictly increasing in both coordinates for all }T>0,V>0;\\ \text{there exists }\mu>0\,\text{ such that }|f(T,V)|\leq\mu|T|\,\,\text{for all }T,V\in\mathbb{R}.\end{array}\right.$$

$$\left\{\begin{array}[]{l}0=\lambda-d\widehat{T}-f(\widehat{T},\widehat{V}),\qquad 0=e^{-\omega h}f(\widehat{T},\widehat{V})-\delta\widehat{T^{*}},\\ 0=N\delta\widehat{T^{*}}-c\widehat{V}.\\ \end{array}\right.$$

$$h_f(s)\equiv f\left(\frac{\lambda }{d}-\frac{\delta }{d}{e}^{\omega h}\cdot s,\frac{N\delta }{c}\cdot s\right)-\delta e^{\omega h}\cdot s.$$

$$\forall \mu \in (0,1)\forall s\in (1-\mu ,1+\mu )\text{mbox}onehas\frac{(s-1{)}^2}{2(1+\mu )}\le v(s)\le \frac{(s-1{)}^2}{2(1-\mu )}.$$

$$(\mathbf{H}{\mathbf{f}}_3)\left(\frac{V}{\hat{V}}-\frac{f(T,V)}{f(T,\hat{V})}\right)\cdot \left(\frac{f(T,V)}{f(T,\hat{V})}-1\right)>0.$$

$$[f(T,\hat{V}){]}^{-1}\ge C_f^1+C_f^2\frac{1}{T},T>0,C_f^i=C_f^i(\hat{V})\ge 0,i=1,2.$$

$$M_F\equiv \left\{\phi \in C^1([-h,0];X),\phi (0)\in D(\mathcal{A}),\dot{\phi }(0)+\mathcal{A}\phi (0)=F(\phi )\right\}.$$

$$\frac{d}{dt}\eta (u_t)=[(D\eta )(u_t)]({\dot{u}}_t),$$

$$\left|\frac{d}{dt}\eta (u_t)\right|\le {\alpha }_{\epsilon },\text{ with }{\alpha }_{\epsilon }\to 0,\text{ as }\epsilon \to 0.$$

$$U^{\mathrm{sdd}-\mathrm{x}}(t,x)\equiv \left(T(t,x)-\hat{T}-{\int }_{\hat{T}}^{T(t,x)}\frac{f(\hat{T},\hat{V})}{f(\theta ,\hat{V})}d\theta \right)e^{-\omega h}+\widehat{T^{\ast }}\cdot v\left(\frac{T^{\ast }(t,x)}{\widehat{T^{\ast }}}\right)$$

$$+\frac{\hat{V}}{N}\cdot v\left(\frac{V(t,x)}{\hat{V}}\right)+\delta \widehat{T^{\ast }}{\int }_{t-\eta (u_t)}^{t}v\left(\frac{f(T(\theta ,x),V(\theta ,x))}{f(\hat{T},\hat{V})}\right)d\theta .$$

$$U^{\mathrm{sdd}}(t)\equiv {\int }_{\Omega }{U}^{\mathrm{sdd}-\mathrm{x}}(t,x)dx.$$

$$\frac{\partial }{\partial t}\left[{\int }_{t-\eta (u_t)}^{t}v\left(\frac{f(T(\theta ,x),V(\theta ,x))}{f(\hat{T},\hat{V})}\right)d\theta \right]$$

$$=v\left(\frac{f(T(t,x),V(t,x))}{f(\hat{T},\hat{V})}\right)-v\left(\frac{f(T(t-\eta (u_t),x),V(t-\eta (u_t),x))}{f(\hat{T},\hat{V})}\right)\left(1-\frac{d}{dt}\eta (u_t)\right)$$

$$=v\left(\frac{f(T(t,x),V(t,x))}{f(\hat{T},\hat{V})}\right)-v\left(\frac{f(T(t-\eta (u_t),x),V(t-\eta (u_t),x))}{f(\hat{T},\hat{V})}\right)+S^{\mathrm{sdd}}(t,x),$$

$$S^{\mathrm{sdd}}(t,x)\equiv v\left(\frac{f(T(t-\eta (u_t),x),V(t-\eta (u_t),x))}{f(\hat{T},\hat{V})}\right)\cdot \frac{d}{dt}\eta (u_t).$$

$$\frac{\partial U^{\mathrm{sdd}-\mathrm{x}}(t,x)}{\partial t}=\left(1-\frac{f(\hat{T},\hat{V})}{f(T(t,x),\hat{V})}\right)e^{-\omega h}\cdot \frac{\partial T(t,x)}{\partial t}+\left(1-\frac{\widehat{T^{\ast }}}{T^{\ast }(t,x)}\right)\cdot \frac{\partial \widehat{T^{\ast }}(t,x)}{\partial t}$$

$$+\frac{1}{N}\cdot \left(1-\frac{\hat{V}}{V(t,x)}\right)\cdot \frac{\partial V(t,x)}{\partial t}+\delta \widehat{T^{\ast }}v\left(\frac{f(T(t,x),V(t,x))}{f(\hat{T},\hat{V})}\right)$$

$$-\delta \widehat{T^{\ast }}v\left(\frac{f(T(t-\eta (u_t),x),V(t-\eta (u_t),x))}{f(\hat{T},\hat{V})}\right)+\delta \widehat{T^{\ast }}{S}^{\mathrm{sdd}}(t,x).$$

Full text

2010 Mathematics Subject Classification: 93C23, 34K20,35K57, 97M60.

1 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 state-dependent.

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

[TABLE]

In (1), susceptible cells $T$ are produced at a rate $\lambda$, 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 $\delta T^{*}$, free virions $V$ are produced by infected cells at rate $N\delta 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\in[a-h,b]\subset\mathbb{R},b>a$, we denote the history segment $v_{t}=v_{t}(\theta)\equiv v(t+\theta),\theta\in[-h,0],t\in[a,b].$

The ODEs delay system (1) is extended to the state-dependent one

[TABLE]

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\in\Omega$ and allow the unknowns to depend on it, i.e. $T(t,x),T^{*}(t,x),V(t,x)$. Now $T(t,x),T^{*}(t,x),V(t,x)$ represent the densities of uninfected cells, infected cells and free virions at position $x$ at time $t$.

Consider a connected bounded domain $\Omega\subset\mathbb{R}^{n}$ with a smooth boundary $\partial\Omega$. Now we are ready to present the PDEs system under consideration

[TABLE]

Here the dot over a function denotes the partial time derivative i.g, $\dot{T}(t,x)={\partial T(t,x)\over\partial t}$, all the constants $\lambda,d,\delta,N,c,\omega$ 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)=\mbox{const}\cdot TV$ 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}\neq 0$ i.e. ${\partial T(t,x)\over\partial n}|_{\partial\Omega}=0$ if $d^{1}\neq 0$ and similarly for $T^{*}(t,x)$ and $V(t,x)$. Here ${\partial\over\partial n}$ is the outward normal derivative on $\partial\Omega$. 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 $\Omega\subset\mathbb{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 $\eta$ and discuss the structure of $\eta$ 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 state-dependent 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.

2 Basic properties of the model

Define the following linear operator  $-\mathcal{A}^{0}$ $=diag\,(d^{1}\Delta,d^{2}\Delta,d^{3}\Delta)$ in $C(\overline{\Omega};\mathbb{R}^{3})$ with $D(\mathcal{A}^{0})\equiv D(d^{1}\Delta)\times D(d^{2}\Delta)\times D(d^{3}\Delta)$. Here, for $d^{i}\neq 0$ we set $D(d^{i}\Delta)\equiv\{v\in C^{2}(\overline{\Omega}):{\partial v(x)\over\partial n}|_{\partial\Omega}=0\}$ and $D(d^{j}\Delta)\equiv C(\overline{\Omega})$ for $d^{j}=0$. We omit the space coordinate $x$, for short, for unknown $u(t)=(T(t),T^{*}(t),V(t))\in X\equiv[C(\overline{\Omega})]^{3}\equiv C(\overline{\Omega};\mathbb{R}^{3})$. It is well-known that the closure $-\mathcal{A}$ (in $X$) of the operator $-\mathcal{A}^{0}$ generates a $C_{0}$-semigroup $e^{-\mathcal{A}t}$ on $X$ which is analytic and nonexpansive [16, p.5]. We denote the space of continuous functions by $C\equiv C([-h,0];X)$ equipped with the sup-norm $||\psi||_{C}\equiv\max_{\theta\in[-h,0]}||\psi(\theta)||_{X}$.

We write, the system (3) in abstract form

[TABLE]

The non-linear continuous mapping $F:C\to X$ is defined by

[TABLE]

Here $\varphi=(\varphi^{1},\varphi^{2},\varphi^{3})\in C$. 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(\theta,x)=\varphi(\theta,x)=(T(\theta,x),T^{*}(\theta,x),V(\theta,x)),\theta\in[-h,0]$ for the delay problem (4)

[TABLE]

In our study we use the standard (c.f. [21, Def. 2.3, p.106] and [21, Def. 2.1, p.105])

Definition 2.1

A function $u\in C([-h,T];X)$ is called a mild solution on $[-h,T)$ of the initial value problem (4), (6) if it satisfies (6) and

[TABLE]

A function $u\in C([-h,T);X)\bigcap C^{1}((0,T);X)$ is called a classical solution on $[-h,T)$ of the initial value problem (4), (6) if it satisfies (6), $u(t)\in D(A)$ for $0<t<T$ and (4) is satisfied on $(0,T)$.

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:\mathbb{R}^{2}\to\mathbb{R}$ is Lipschitz continuous and satisfies

[TABLE]

We have the following result

Proposition 1

Let nonlinear function $f$ be Lipschitz and satisfy ${\bf(Hf_{1})}$ (see (8)), state-dependent delay $\eta:C\to[0,h]$ is locally Lipschitz. Then the initial value problem (4), (6) has a unique classical solution which is global in time i.e. defined for all $t\geq 0$.

Proof of Proposition 1. We start with discussion of mild solutions. Since the semigroup generated by the linear part $-\mathcal{A}$ is not necessarily compact (in cases when at least one constant $d^{i}\neq 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 $\varphi$ 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 extension of $\bar{\varphi}(t)=\varphi(t)$ for $t\in[-h,0]$ and $\bar{\varphi}(t)=e^{-\mathcal{A}t}\varphi(0)$ for $t\geq 0$. Next we change variable $u(t)=\bar{\varphi}(t)+y(t)$ and consider complete metric space $A(\alpha,\beta,\gamma)\equiv\{y\in C([-h,\alpha];X),y_{0}\equiv 0,\max_{t\in[0,\alpha]}||y(t)||_{X}\leq\beta,\sup_{s\neq t}||y(s)-y(t)||_{X}\cdot|s-t|^{-1}\leq\gamma\}$ endowed by the metrics of the space of continuous functions. The operator $\mathcal{F}:A(\alpha,\beta,\gamma)\to C([-h,\alpha];X)$ is defined as $\mathcal{F}(y)(t)\equiv 0$ for $t\in[-h,0]$ and $\mathcal{F}(y)(t)\equiv\int^{t}_{0}e^{\mathcal{A}(t-\tau)}F(\bar{\varphi}_{\tau}+y_{\tau})\,d\tau$ for $t\in(0,\alpha]$. It is not difficult to check that our non-linear mapping F, defined by (5), satisfies (see the estimate $|f(T,V)|\leq\mu|T|$ in ${\bf(Hf_{1})}$) $||F(\psi)||_{X}\leq n_{1}+n_{2}||\psi||_{C}$ and is locally almost Lipschitz on $A(\alpha,\beta,\gamma)$ by the terminology of [15]. The last means $||F(\psi^{1})-F(\psi^{2})||_{X}\leq L_{F}(\gamma)||\psi^{1}-\psi^{2}||_{C}$. Standard computations show that operator $\mathcal{F}$ maps $A(\alpha,\beta,\gamma)$ into itself provided $\alpha,\beta,\gamma$ satisfy $\alpha(n_{1}+n_{2}(||\varphi||+\beta))\leq\beta,n_{1}+n_{2}\beta\leq\gamma$. Additional condition $\alpha L_{F}(\gamma)<1$ guarantees the contraction of $\mathcal{F}$. The classical Banach Fixed Point Theorem gives the unique fixed point $\hat{y}$ and hence the unique mild solution $u=\bar{\varphi}+\hat{y}$. The linear growth bound of $F$ implies the global continuation of the mild solution.

Our next step is to show that any mild solution is classical. Let us fix any mild solution $u$ to (4), (6) and define $g(t)\equiv F(u_{t}),t\geq 0.$ For any $t^{0}>0$, mapping $g$ is continuous on $[0,t^{0}]$ since $F$ and $u$ are continuous. We notice that by construction, the solution is Lipschitz in time on $[0,t^{0}]$ (see also restrictions in (6)). Hence, $||g(t)-g(s)||=||F(u_{t})-F(u_{s})||\leq L_{F}\max_{\theta\in[-h,0]}||u(t+\theta)-u(s+\theta)||\leq L_{F}L^{[0,t^{0}]}_{u}\cdot|t-s|$. Here we use the almost Lipschitz property of $F$. Now we consider the following (non-delayed) initial value problem

[TABLE]

which has a unique solution. The solution of (9) is $v=u$ in case $x=u(0)$.

We remind that $C_{0}$-semigroup $e^{-\mathcal{A}t}$ is analytic on $X$ [16, p.5]. Hence theorem 3.5 [21, p.114] implies that the mild solution (of (9) and hence of (4), (6)) is classical for $t\geq 0.$ The proof of Proposition 1 is complete.

Define the set (c.f. (6))

[TABLE]

[TABLE]

where $\mu$ is defined in $(Hf_{1})$ and all the inequalities hold pointwise w.r.t. $x\in\overline{\Omega}.$

We need further assumptions (which include $(Hf_{1})$) on Lipschitz function $f$ :

[TABLE]

We have the following result

Proposition 2

Let non-linear function $f$ satisfy $(Hf_{1}+)$ (see (11)), state-depen-dent delay $\eta:C\to[0,h]$ is locally Lipschitz. Then $\Omega_{Lip}$ is invariant i.e. for any $\varphi\in\Omega_{Lip}$ the unique solution to problem (4), (6) satisfies $u_{t}\in\Omega_{Lip}$ for all $t\geq 0$.

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 Lipschitz property of nonlinearity $F$. The estimates (for the subtangential condition) are the same as for the constant delay case, see e.g. [17, Theorem 2.2]. We do not repeat it here. It is important to notice that the solutions are classic for all $t\geq 0$ (but not for $t\geq h$ as could be in the case of merely continuous initial functions $\varphi\in C$). The proof of Proposition 2 is complete.

2.1 Stationary solutions

Let us discuss stationary solutions of (3). By such solutions we mean time independent $\widehat{u}$ which, in general, may depend on $x\in\overline{\Omega}$. Consider the system (3) with $u(t)=u(t-\eta(u_{t}))=\widehat{u}$ and denote the coordinates of a stationary solution by $(\widehat{T},\widehat{T^{*}},\widehat{V})=\widehat{u}\equiv\widehat{\varphi}(\theta),\,\theta\in[-h,0]$. Since stationary solutions of (3) do not depend on the type of delay (state-dependent or constant) we have (see e.g. [17])

[TABLE]

Equations hold pointwise w.r.t. $x\in\overline{\Omega}.$

It is easy to see that the trivial stationary solution $(\lambda d^{-1},0,0)$ always exists. We are interested in nontrivial disease stationary solutions of (3). Using (12), we have $\widehat{T}=(\lambda-\delta\widehat{T^{*}}e^{\omega h})d^{-1}$ and $\widehat{V}={N\delta\over c}\widehat{T^{*}}$. It gives the condition on the coordinate $\widehat{T^{*}}$ which should belong to $(0,\lambda e^{\omega h}\delta^{-1}]$. Denote (c.f. [17])

[TABLE]

Assume $f$ satisfies

${\bf(Hf_{2})}\hskip 28.45274pth_{f}(s)=0$ has at least one and at most finite roots on $(0,\lambda e^{\omega h}\delta^{-1}]$.

We denote an arbitrary root of $h_{f}(s)=0$ by $\widehat{T^{*}}$ and define the corresponding $\widehat{T}=(\lambda-\delta\widehat{T^{*}}e^{\omega h})d^{-1}$ and $\widehat{V}={N\delta\over c}\widehat{T^{*}}$. The point $(\widehat{T},\widehat{T^{*}},\widehat{V})$ satisfies (12), so it is a disease stationary solutions of (3).

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\in\Omega$. We remind that $\Omega$ is a connected set, so a function $v\in C(\overline{\Omega})$ may take either one or continuum values. Assumption ${\bf(Hf_{2})}$ implies $\widehat{T^{*}}(x)\equiv\widehat{T^{*}}\in\mathbb{R}$, so $(\widehat{T},\widehat{T^{*}},\widehat{V})$ is independent of $x\in\overline{\Omega}.$

Remark 2

Below we mention some well-known examples of non-linear functions $f$ when we have exactly one root of $h_{f}(s)=0$. The first one is the DeAngelis-Bendington [1, 4] functional response $f(T,V)={kTV\over 1+k_{1}T+k_{2}V}$, with $k,k_{1}\geq 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)={kTV\over 1+k_{2}V}$. Another example of the nonlinearity is the Crowley-Martin incidence rate $f(T,V)={kTV\over(1+k_{1}T)(1+k_{2}V)}$, with $k\geq 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 ${\bf(Hf_{2})}$.

3 Stability of disease stationary solutions

The following Volterra function $v(s)=s-1-\ln s:(0,+\infty)\to\mathbb{R}_{+}$ plays an important role in construction of Lyapunov functionals [12, 17]. One can see that $v(s)\geq 0$ and $v(s)=0$ if and only if $s=1$. The derivative equals $\dot{v}(s)=1-{1\over s}$, which is obviously negative for $x\in(0,1)$ and positive for $x>1$. The graph of $v$ explains the use of the composition $v\left({s\over s^{0}}\right)$ in the study of the stability properties of an equilibrium $s^{0}$. Another important property is the following [29] estimate

[TABLE]

To check it, one simply observes that all three functions vanish at $s=1$ and $\left|{d\over ds}\left({(s-1)^{2}\over 2(1+\mu)}\right)\right|\leq|{d\over ds}v(s)|\leq\left|{d\over ds}\left({(s-1)^{2}\over 2(1-\mu)}\right)\right|$ in the $\mu$-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 ${\bf(Hf_{2})}$).

[TABLE]

This property simply means that the value ${f(T,V)\over f(T,\widehat{V})}$ is always strictly between $1$ and ${V\over\widehat{V}}$ for any $T\geq 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

${\bf(Hf_{4})}$ Function $f$ is either differentiable with respect to its first coordinate or satisfies

[TABLE]

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)

[TABLE]

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 ($\varphi\in\Omega_{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.

Theorem 3.1

Let the nonlinear function $f$ satisfy $(Hf_{1}+),(Hf_{2}),(Hf_{3}),(Hf_{4})$ (see (11), (16), (15)) and state-dependent delay $\eta:C\to[0,h]$ be locally Lipschitz in $C$ and continuously differentiable in a neighbourhood of equilibrium $\widehat{\varphi}\equiv(\widehat{T},\widehat{T^{*}},\widehat{V})$. Then the stationary solution $\widehat{\varphi}$ is locally asymptotically stable (in $M_{F}$).

Remark 4

Similar to ODE case, described in [29, Remark 13], we have the following property. For any $u\in C^{1}([-h,b);X)$ one has for $t\in[0,b)$

[TABLE]

where $[(D\eta)(u_{t})](\cdot)$ is the Fréchet derivative of $\eta$ at point $u_{t}$. Hence, (for a solution in $\varepsilon$-neighborhood of the stationary solution $\widehat{\varphi}$) the estimate $|{d\over dt}\eta(u_{t})|\leq||(D\eta)(u_{t})||_{L(C;R)}\cdot||\dot{u}_{t}||_{C}\leq\varepsilon\,||(D\eta)(u_{t})||_{L(C;R)}$ guarantees the property

[TABLE]

due to the boundedness of $||(D\eta)(\psi)||_{L(C;R)}$ as $\varepsilon\to 0$ (here $||\psi-\widehat{\varphi}||_{C}<\varepsilon$).

Proof of Theorem 3.1. Let us consider (point-wise) the following auxiliary functional

[TABLE]

[TABLE]

Now we can introduce the following Lyapunov functional with state-dependent delay along a solution of (3)

[TABLE]

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^{\mathrm{sdd-x}}(t,x)$ defined in (19). This time derivative is considered along classical solutions of (3). It gives the possibility to consider ${\partial T(t,x)\over\partial t},{\partial\widehat{T^{*}}(t,x)\over\partial t},{\partial V(t,x)\over\partial 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 integral term

[TABLE]

[TABLE]

[TABLE]

where we denoted for short

[TABLE]

Remark 5

The term $S^{\mathrm{sdd}}$ appears due to the presence of the state-dependent delay. It makes the technical calculations more challenging. The sign of $S^{\mathrm{sdd}}$ is undefined, so we propose below (see also [29, 30]) a way to compensate/bound $S^{\mathrm{sdd}}$ by other positive defined terms in ${\partial U^{\mathrm{sdd-x}}\over\partial t}$ to have the time derivative of the Lyapunov functional (along a solution) negative defined relative to the equilibrium.

Now we differentiate

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Calculations, using (12), particularly, $\lambda=d\widehat{T}+f(\widehat{T},\widehat{V})$ give

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where, for short, we collected some terms as $C^{1}$. It is written as follows

[TABLE]

[TABLE]

[TABLE]

Calculations show that

[TABLE]

[TABLE]

[TABLE]

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$):

[TABLE]

[TABLE]

which is split on the sum of seven logarithms to write shortly, using the Volterra function $v$

[TABLE]

[TABLE]

[TABLE]

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 $\Omega$. Denote them, for short, as

[TABLE]

[TABLE]

In case of differentiable $f$ $(Hf_{4})$ (see (16)) need the following simple

Proposition 3

Let $p:\mathbb{R}\to\mathbb{R}$ be differentiable. Then $\int_{\Omega}p(u(x))\Delta u(x)\,dx=-\int_{\Omega}p^{\prime}(u)||\nabla u||^{2}\,dx$ for any $u\in C^{2}(\overline{\Omega})$ satisfying ${\partial u(x)\over\partial n}|_{\partial\Omega}=0$.

Proof of Proposition 3. We use the classical Gauss-Ostrogradsky theorem. Consider the vector field $E\equiv p(u)\nabla u$. Hence $\mathrm{div}\,E=p^{\prime}(u)||\nabla u||^{2}+p(u)\Delta u$. One has $\int_{\Omega}\mathrm{div}\,E\,dx=\int_{\Omega}p^{\prime}(u)||\nabla u||^{2}\,dx+\int_{\Omega}p(u)\Delta u\,dx=\int_{\partial\Omega}p(u)(\nabla u,n)\,dS=0.$ The last equality due to the Neumann boundary conditions. Finally, $\int_{\Omega}p(u)\Delta u\,dx=-\int_{\Omega}p^{\prime}(u)||\nabla u||^{2}\,dx.$ It completes the proof of Proposition 3.

Now we apply Proposition 3 to show that $D^{\mathrm{diff-3}}(t)\leq 0.$ Let us start with the first term in $D^{\mathrm{diff-3}}(t,x)$, see (24), and show that $\int_{\Omega}\left(1-{f(\widehat{T},\widehat{V})\over f(T(t,x),\widehat{V})}\right)\,\Delta T(t,x)\,dx\leq 0.$ For this we set $p(T)=\left(1-{f(\widehat{T},\widehat{V})\over f(T(t,x),\widehat{V})}\right)$ and check that $p^{\prime}(T)=f^{\prime}_{1}(T,\widehat{V})f(\widehat{T},\widehat{V})\times[f(T,\widehat{V})]^{-2}\geq 0$ due to $f^{\prime}_{1}(T,\cdot)\geq 0$ by the assumption on $f$. Similar considerations with the second and third terms in (24) show that

[TABLE]

[TABLE]

[TABLE]

Remark 6

In case of nondifferentiable $f$ we prove $D^{\mathrm{diff-3}}(t)\leq 0$, using alternative (geometrical) conditions on $f$ given in $(Hf_{4})$ (see (16)).

Then $D^{\mathrm{diff-3}}(t)\leq 0$ along any classical solution.

Now we combine the arguments above to study the Lyapunov functional $U^{\mathrm{sdd}}(t)$, see (20). We have the following equality (c.f. (22))

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Here $D^{\mathrm{diff-3}}(t)$ is defined in (24) and transformed in (25), $C^{1}$ is presented as in (23) and $S^{\mathrm{sdd}}$ is defined in (21). We remind (see (12)) that $\delta\widehat{T^{*}}=e^{-\omega h}f(\widehat{T},\widehat{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

[TABLE]

[TABLE]

[TABLE]

[TABLE]

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^{\mathrm{diff-3}}(t)\leq 0$ is given in (25). To show that

[TABLE]

we use the property $(Hf_{3})$ of function $f$ (see (15)).

Now we plan to prove that ${d\over dt}U^{\mathrm{sdd}}(t)\leq 0$ in a small neighbourhood of the stationary solution with the equality only in case of $(T,T^{*},V)=(\widehat{T},\widehat{T^{*}},\widehat{V})$. In the particular case of constant delay, one has $S^{\mathrm{sdd}}(t,x)=0$ which may lead to the global stability of $(\widehat{T},\widehat{T^{*}},\widehat{V})$.

We rewrite, for short, (26) as

[TABLE]

where $D^{\mathrm{sdd}}(t,x)$ contains all the terms except the last one in (26). As proved above $\int_{\Omega}D^{\mathrm{sdd}}(t,x)\,dx\geq 0$. Let us start with an analysis of the zero-sets $D^{\mathrm{sdd}}(t,x)=0,S^{\mathrm{sdd}}(t,x)=0$ and ${d\over dt}U^{\mathrm{sdd}}(t)=0$.

We start with $D^{\mathrm{sdd}}(t,x)=0$. One sees from (26) that $T=\widehat{T}$. Since $v(s)=0$ iff $s=1$, we see from (15) that $V=\widehat{V}$. Hence $T^{*}=\widehat{T^{*}}$. One also sees $f(T(t-\eta(u_{t}),x),V(t-\eta(u_{t}),x))=f(\widehat{T},\widehat{V})$. Moreover, $D^{\mathrm{diff-3}}(t)=0$, means (see Proposition 3 and (25)) that $T,T^{*}$ and $V$ are independent of $x$. The zero set $S^{\mathrm{sdd}}(t,x)=0$ is described (see (21) by $f(T(t-\eta(u_{t}),x),V(t-\eta(u_{t}),x))=f(\widehat{T},\widehat{V})$ or ${d\over dt}\eta(u_{t})=0$ along a solution. It is important for us that the zero-set $D^{\mathrm{sdd}}(t,x)=0$ is a singleton $(T,T^{*},V)=(\widehat{T},\widehat{T^{*}},\widehat{V})$ and is a subset of $S^{\mathrm{sdd}}(t,x)=0$. The rest of the proof that in a small neighbourhood of $(\widehat{T},\widehat{T^{*}},\widehat{V})$ one has $|S^{\mathrm{sdd}}(t,x)|<D^{\mathrm{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\over dt}U^{\mathrm{sdd}}(t)\leq 0$ in a small neighbourhood of the stationary solution with the equality only in case of $(T,T^{*},V)=(\widehat{T},\widehat{T^{*}},\widehat{V})$ completes the proof of Theorem 3.1.

It is interesting to notice that $\varphi\in M_{F}$ (see (17)) is not a necessary condition for our approach. Now we consider a wider set $\Omega_{Lip}$ (see (10)). Let us discuss a particular simple form of the delay (c.f. examples in [29])

[TABLE]

with a locally Lipschitz $\xi$. To check the property (18) we calculate

[TABLE]

Hence, in the $\varepsilon$-neighborhood of the stationary solution $\hat{u}$, one has

[TABLE]

Here $L_{\xi,\varepsilon}$ is the Lipschitz constant of $\xi$. More general delay terms could be used

[TABLE]

with a differentiable $\rho:\mathbb{R}\to[0,h]$. The example (28) is a particular case of (29) with $\rho(s)\equiv s$ and $\kappa(s)\equiv 1.$

The discussion above shows that property (18) of the state-dependent delay (29) allows to use the proof of Theorem 3.1 to get the following result in $\Omega_{Lip}$

Theorem 3.2

Let non-linear function $f$ satisfy $(Hf_{1}+),(Hf_{2}),(Hf_{3}),(Hf_{4})$ (see (11), (16), (15)) and state-dependent delay $\eta:C\to[0,h]$ be of the form (29). Then the stationary solution $\widehat{\varphi}$ is locally asymptotically stable.

Acknowledgments. The author is thankful to anonymous referees for useful comments and suggestions. This work was supported in part by GA CR under project 16-06678S.