<figure><img src="image_1.png" /></figure>

HAL
open science

# Nonlinear Spatiotemporal Viral Infection Model with CTL Immunity: Mathematical Analysis

Jaouad Danane, Karam Allali, Leon Tine, Vitaly Volpert

► To cite this version:

Jaouad Danane, Karam Allali, Leon Tine, Vitaly Volpert. Nonlinear Spatiotemporal Viral Infec-
tion Model with CTL Immunity: Mathematical Analysis. Mathematics, 2020, 8 (1), pp.1-13.
10.3390/math8010052. hal-03029312

HAL Id: hal-03029312

https://hal.science/hal-03029312v1

Submitted on 11 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

Article

# Nonlinear Spatiotemporal Viral Infection Model with CTL Immunity: Mathematical Analysis †

Jaouad Danane <fnref n="1" />, Karam Allali <fnref n="1" />, Léon Matar Tine <fnref n="2" />,<fnref n="3" />,* and Vitaly Volpert <fnref n="2" />,<fnref n="3" />,<fnref n="4" />

Received: 27 November 2019; Accepted: 13 December 2019; Published: 1 January 2020

<figure><img src="image_4.png" /></figure>

check for
updates

**Abstract:** A mathematical model describing viral dynamics in the presence of the latently infected cells and the cytotoxic T-lymphocytes cells (CTL), taking into consideration the spatial mobility of free viruses, is presented and studied. The model includes five nonlinear differential equations describing the interaction among the uninfected cells, the latently infected cells, the actively infected cells, the free viruses, and the cellular immune response. First, we establish the existence, positivity, and boundedness for the suggested diffusion model. Moreover, we prove the global stability of each steady state by constructing some suitable Lyapunov functionals. Finally, we validated our theoretical results by numerical simulations for each case.

**Keywords:** viral infection; diffusion; Lyapunov functional; convergence

## 1. Introduction

Viral infections represent a major cause of morbidity with important consequences for patient's health and for the society. Among the most dangerous, let us cite the human immunodeficiency virus (HIV) that attacks immune cells leading to the deficiency of the immune system [1,2], the human papillomavirus (HPV) that infects basal cells of the cervix [3,4], and the hepatitis B virus (HBV) and the hepatitis C virus (HCV) that attack liver cells [5-8]. Mathematical modeling becomes an important tool for the understanding and predicting the spread of viral infection, and for the development of efficient strategies to control its dynamics [9-12]. One of the basic models of viral infection suggested by Nowak in 1996 describes the interactions among uninfected cells, infected cells, and free viruses. Nowadays, modeling of viral infection actively develops with a variety of new models and methods [9-15] (see also the monograph [16] and the references therein). The action of immune system was introduced into the basic model with cytotoxic T-lymphocytes cells (CTL) killing infected cells [12-15,17,18]. The impact of CTL cells with a saturated incidence function was considered in [12]:

<footnote n="1">¹ Laboratory of Mathematics and Applications, Faculty of Sciences and Technologies, University Hassan II of Casablanca, P.O.Box 146, Mohammed, Morocco; jaouaddanane@gmail.com (J.D.); allali@hotmail.com (K.A.)</footnote>

<footnote n="2">² CNRS UMR 5208 Institut Camille Jordan, University Claude Bernard Lyon 1, University de Lyon, 69622 Villeurbanne CEDEX, France; volpert@math.univ-lyon1.fr</footnote>

<footnote n="3">³ INRIA Team Dracula, INRIA Lyon La Doua, 69603 Villeurbanne, France</footnote>

<footnote n="4">⁴ S.M. Nikol'skii Mathematical Institute, Peoples' Friendship University of Russia (RUDN University), 6 Miklukho-Maklaya St, 117198 Moscow, Russia</footnote>

<footnote>* Correspondence: leon-matar.tine@univ-lyon1.fr</footnote>

<footnote>† Fully documented templates are available in the elsarticle package on CTAN.</footnote>
$$
\left\{
\begin{aligned}
\dot{H} &= \lambda - d_1 H - \frac{k_1 H V}{H + V'}, \\
\dot{S} &= \frac{k_1 H V}{H + V} - d_2 S - k_2 S, \\
\dot{Y} &= k_2 S - d_3 Y - p Y Z, \\
\dot{V} &= a Y - d_4 V, \\
\dot{Z} &= c Y Z - b Z.
\end{aligned}
\right.
\qquad (1)
$$

Here $H, S, Y, V$, and $Z$ represent the densities of uninfected cells, exposed cells, infected cells, free virus, and CTL cells, respectively. Our model uses a more realistic saturated incidence function $\frac{k_1HV}{H+V}$ [10,11]. This saturated incidence functional describes the infection rate taking into consideration the effect of free viruses crowd near the healthy cells. The parameters of the system in Equation (1) are described in Table 1.

<table><caption>**Table 1.** The parameters of the mathematical model and their descriptions.</caption><thead><tr><th>Coefficient</th><th>Description</th></tr></thead><tbody><tr><td>λ</td><td>The birth rate of the uninfected cells</td></tr><tr><td>k<sub>1</sub></td><td>The rate of infection</td></tr><tr><td>d<sub>1</sub></td><td>The natural mortality of the susceptible cells</td></tr><tr><td>d<sub>2</sub></td><td>The death rate of exposed cells</td></tr><tr><td>k<sub>2</sub></td><td>The average that exposed cells become infected</td></tr><tr><td>d<sub>3</sub></td><td>The death rate of infected cells, not by CTL killing</td></tr><tr><td>a</td><td>The rate of production the virus by infected cells</td></tr><tr><td>d<sub>4</sub></td><td>The rate of viral clearance</td></tr><tr><td>p</td><td>Clearance rate of infection</td></tr><tr><td>c</td><td>Activation rate CTL cells</td></tr><tr><td>b</td><td>Death rate of CTL cells</td></tr><tr><td>d</td><td>Diffusion coefficient</td></tr></tbody></table>

The majority of mathematical models of viral infection ignores the spatial movement of viruses and cells, assuming that the virus and cell populations are well mixed [19]. However, their mobility and a nonuniform spatial distribution can play an important role for the infection development [20]. Thus far, few studies have been devoted to the influence of spatial structure on the dynamics of the virus [21,22]. Reaction-diffusion waves of infection spreading were studied in [23-25].

In this work, we consider the previous model in Equation (1) taking into account virus diffusion:

$$
\left\{
\begin{aligned}
\frac{\partial H}{\partial t} &= \lambda - d_1 H - \frac{k_1 H V}{H + V'}, \\
\frac{\partial S}{\partial t} &= \frac{k_1 H V}{H + V} - d_2 S - k_2 S, \\
\frac{\partial Y}{\partial t} &= k_2 S - d_3 Y - p Y Z, \\
\frac{\partial V}{\partial t} &= d \Delta V + a Y - d_4 V, \\
\frac{\partial Z}{\partial t} &= c Y Z - b Z.
\end{aligned}
\right.
\qquad (2)
$$

As in (1), $H, S, Y, V$, and $Z$ represent the densities of uninfected cells, latently infected cells, infected cells, free virus, and CTLs cells, depending now on the space coordinate location $x$ and on time $t$, $x \in \Omega$, where $\Omega$ is an open bounded set in $\mathbb{R}^n$. The positive constant $d$ is the virus diffusion coefficient. All parameters of the system in Equation (2) have the same biological meanings as in the
model in Equation (1). In this paper, we consider the system in Equation (2) with the homogeneous Neumann boundary conditions:

$$ \frac{\partial V}{\partial \nu} = 0, \quad \text{on } \partial\Omega \times [0, +\infty) \qquad (3) $$

and the initial conditions:

$$ H(x, 0) = \varphi_1(x) \geq 0; S(x, 0) = \varphi_2(x) \geq 0; Y(x, 0) = \varphi_3(x) \geq 0; \qquad (4) $$

$$ V(x, 0) = \varphi_4(x) \geq 0; Z(x, 0) = \varphi_5(x) \geq 0 \quad \forall x \in \bar{\Omega}. \qquad (5) $$

The paper is organized as follows. The next section is devoted to the well-posedness of the model, followed in Section 3 by the global stability analysis. In Section 4, we illustrate the results with the numerical simulations.

## 2. Well-Posedness of Model

In this section, we investigate the well-posedness of the model in Equation (2) proving the global existence, the positivity and the boundedness of solutions.

**Proposition 1.** For any initial condition satisfying Equations (4)–(5), there exists a unique solution to the problem in Equations (2)–(3) defined on $t \in (0, +\infty)$. Moreover, this solution stays non-negative and bounded for all $t > 0$.

**Proof.** Let consider the set

$$ X_T = \{V \in L^2([0, T] \times \Omega) \cap C([0, T]; L^2(\Omega)) ; 0 \leq V(x, t) \leq C_T\}, $$

where $T$ is a fixed positive constant and $C_T = \| V_0 \|_{L^\infty(\Omega)} + a \frac{\lambda}{m} T (V_0 = V(x, 0))$, where $m = \min(d_1, d_2, d_3, pb/c)$.

To prove the existence of the solution, we define the following map

$$ \Psi : X_T \to L^2([0, T] \times \Omega) \\ \bar{V} \mapsto \Psi(\bar{V}) = V $$

such that

$$ \left\{ \begin{array}{ll} \displaystyle \frac{\partial V}{\partial t} - d\Delta V = aY - d_4V, & \text{in } \Omega \\ \displaystyle \frac{\partial V}{\partial \nu} = 0, & \text{in } [0, T] \times \partial\Omega \\ V(x, 0) = \varphi_4(x), & \forall x \in \bar{\Omega}, \end{array} \right. \qquad (6) $$

where $Y$ is the third component of the solution vector of the following subsystem

$$ \left\{ \begin{array}{l} \displaystyle \frac{\partial H(x,t)}{\partial t} = \lambda - d_1 H(x,t) - \frac{k_1 H(x,t) \bar{V}(x,t)}{H(x,t) + \bar{V}(x,t)}, \\ \displaystyle \frac{\partial S(x,t)}{\partial t} = \frac{k_1 H(x,t) \bar{V}(x,t)}{H(x,t) + \bar{V}(x,t)} - d_2 S(x,t) - k_2 S(x,t), \\ \displaystyle \frac{\partial Y(x,t)}{\partial t} = k_2 S(x,t) - d_3 Y(x,t) - pY(x,t)Z(x,t), \\ \displaystyle \frac{\partial Z(x,t)}{\partial t} = cY(x,t)Z(x,t) - bZ(x,t), \end{array} \right. \qquad (7) $$
with the initial data

$$ 
\begin{aligned}
H(x,0) &= \varphi_1(x) \ge 0; S(x,0) = \varphi_2(x) \ge 0; Y(x,0) = \varphi_3(x) \ge 0; \\
Z(x,0) &= \varphi_5(x) \ge 0 \quad \forall x \in \bar{\Omega}.
\end{aligned}
$$

Then, the system in Equation (7) can be written abstractly in $X = \Omega^4$ by the following form

$$ 
\begin{aligned}
U'(t) &= AU(t) + F(U(t)), \quad \forall t > 0, \\
U(0) &= U_0 \in X,
\end{aligned}
$$

with $U = (H, S, Y, Z)^T$, $U_0 = (\varphi_1, \varphi_2, \varphi_3, \varphi_5)^T$, and

$$ F(U(t)) = \begin{pmatrix} \lambda - d_1 H(x,t) - \frac{k_1 H(x,t) \bar{V}(x,t)}{H(x,t) + \bar{V}(x,t)} \\ \frac{k_1 H(x,t) \bar{V}(x,t)}{H(x,t) + \bar{V}(x,t)} - d_2 S(x,t) - k_2 S(x,t) \\ k_2 S(x,t) - d_3 Y(x,t) - p Y(x,t) Z(x,t) \\ c Y(x,t) Z(x,t) - b Z(x,t) \end{pmatrix}. $$

It is clear that $F$ is locally Lipschitz in $U$. Using the theorem of Cauchy-Lipschitz, we deduce that the system in Equation (7) admits a unique local solution on $[0, \tau]$, where $\tau \le T$. In addition, the system in Equation (7) can be written of the form

$$ 
\left\{ 
\begin{aligned}
\frac{\partial H}{\partial t} &= F_1(H, S, Y, Z), \\
\frac{\partial S}{\partial t} &= F_2(H, S, Y, Z), \\
\frac{\partial Y}{\partial t} &= F_3(H, S, Y, Z), \\
\frac{\partial Z}{\partial t} &= F_4(H, S, Y, Z).
\end{aligned}
\right. 
$$

It is easy to see that the functions $F_i(H, S, Y, Z)$, $1 \le i \le 4$ are continuously differentiable, verifying $F_1(0, S, Y, Z) = \lambda \ge 0$, $F_2(H, 0, Y, Z) = \frac{k_1 H \bar{V}}{H + \bar{V}} \ge 0$, $F_3(H, S, 0, Z) = k_2 S \ge 0$ and $F_4(H, S, Y, 0) = 0$ for all $H, S, Y, \bar{V}, Z \ge 0$. Since the initial data of the system in Equation (7) are nonnegative, we obtain the positivity of $H, S, Y$, and $Z$ thanks to the quasi-reversibility principle.

Now, we show the boundedness of solution. Let

$$ 
\begin{aligned}
\mathcal{T}(x,t) &= H(x,t) + S(x,t) + Y(x,t) + \frac{p}{c} Z(x,t) \\
\frac{\partial \mathcal{T}}{\partial t} &= \lambda - d_1 H(x,t) - d_2 S(x,t) - d_3 Y(x,t) - p \frac{b}{c} Z(x,t) \\
&\le \lambda - m \mathcal{T}(x,t),
\end{aligned}
$$

Then,

$$ \mathcal{T}(x,t) \le \mathcal{T}(x,0)e^{-mt} + \frac{\lambda}{m}(1 - e^{-mt}), $$

For biological reasons, we assume that the problem initial data are upper-bounded by the carrying capacity. This means that $\mathcal{T}(x,0) \le \frac{\lambda}{m}$. We deduce that
$$ 
\begin{aligned}
H(x,t) &\le \frac{2\lambda}{m}, \quad S(x,t) \le \frac{2\lambda}{m}, \\
Y(x,t) &\le \frac{2\lambda}{m}, \quad Z(x,t) \le \frac{2\lambda}{m}.
\end{aligned}
$$

Thus, H, S, Y, and Z are bounded.

Let us recast the system in Equation (6) as follows

$$ 
\left\{ 
\begin{array}{ll}
\displaystyle \frac{\partial V}{\partial t} - d\Delta V + d_4 V = aY, & \text{on } \Omega, \\
\\
\displaystyle \frac{\partial V}{\partial v} = 0, & \text{on } \partial\Omega.
\end{array}
\right.
$$

We know that $0 \le Y(x,t) \le \frac{\lambda}{m}$, then from the proposition 2.1 in [26], we deduce for all $V_0 \in L^2(\Omega)$ the existence and the uniqueness of the solution $V \in L^2([0, T]; H^1(\Omega)) \cap C([0, T]; L^2(\Omega))$ such that $\partial_t V \in L^2([0, T]; H^1(\Omega))'$. Furthermore, if $V_0 \in L^\infty(\Omega)$, by using the maximum principle relation, we have

$$ 0 \le V(t,x) \le \| V_0 \|_{L^\infty(\Omega)} + a \frac{\lambda}{m} T. $$

We note that $\Psi$ is well defined and continuous and $\Psi(X_T)$ is compact, thus $\Psi$ admits a fixed point. Then, we conclude the existence of the solution $V$ of (6) and it is positive and bounded. $\square$

## Equilibria and Basic Reproduction Number

The system in Equation (2) has an infection-free equilibrium $E_f = (\frac{\lambda}{d_1}, 0, 0, 0, 0)$, corresponding to the total absence of viral infection. The basic reproduction number of the system in Equation (2) is given by

$$ R_0 = k_1 \frac{k_2}{d_2 + k_2} \frac{a}{d_3} \frac{1}{d_4}, \qquad (8) $$

with $\frac{k_2}{d_2 + k_2}$ the ratio of exposed cells that will become infected, $\frac{a}{d_3}$ the average of free virus produced by an infected cell, and $\frac{1}{d_4}$ the lifespan of the virus. The biological interpretation of $R_0$ represents the rate of secondary infections generated by an infected cell when it is introduced into a population of uninfected cells.

In addition to the disease free equilibrium, our system (Equation (2)) admits three endemic equilibria. The first of them is $E_1 = (H_1, S_1, Y_1, V_1, Z_1)$, where

$$ 
\begin{aligned}
H_1 &= \frac{\lambda}{d_1 + k_1(1 - \frac{1}{R_0})}, & Z_1 &= 0, \\
S_1 &= \frac{k_1 \lambda R_0 (1 - \frac{1}{R_0})}{(d_2 + k_2)(d_1 + k_1(1 - \frac{1}{R_0}))((1 - \frac{1}{R_0})R_0 + 1)}, & \\
Y_1 &= \frac{d_4 \lambda R_0 (1 - \frac{1}{R_0})}{ad_1 + ak_1(1 - \frac{1}{R_0})}, & V_1 &= \frac{\lambda R_0 (1 - \frac{1}{R_0})}{d_1 + k_1(1 - \frac{1}{R_0})}.
\end{aligned}
$$

This endemic steady state is specified as endemic equilibrium without cellular immunity. The second endemic steady state is $E_2 = (H_2, S_2, Y_2, V_2, Z_2)$, where
$$H_2 = \frac{-abd_1 - abk_1 + \lambda cd_4 + \sqrt{A}}{2cd_1d_4},$$

$$S_2 = \frac{d_3 R_0 b(-abd_1 - abk_1 + \lambda cd_4 + \sqrt{A})}{k_2 c(abd_1 - abk_1 + \lambda cd_4 + \sqrt{A})},$$

$$Y_2 = \frac{b}{c}, \quad V_2 = \frac{ba}{cd_4},$$

$$Z_2 = \frac{d_3((R_0 - 1)(-ak_1b + \lambda cd_4 + \sqrt{A}) - abd_1(R_0 + 1))}{p(abd_1 - ak_1b + \lambda cd_4 + \sqrt{A})}.$$

This endemic steady state is specified as endemic equilibrium with cellular immunity. It is also called interior equilibrium. The third endemic steady state is $E_3 = (H_3, S_3, Y_3, V_3, Z_3)$, where

$$H_3 = -\frac{abd_1 + abk_1 - \lambda cd_4 + \sqrt{A}}{2cd_1d_4},$$

$$S_3 = \frac{k_1 H_3 V_3}{(k_2 + d_2)(H_3 + V_3)},$$

$$Y_3 = \frac{b}{c}, \quad V_3 = \frac{ba}{cd_4},$$

$$Z_3 = \frac{k_2 S_3 - d_3 Y_3}{p Y_3},$$

with $A = (abk_1 - \lambda cd_4)^2 + a^2 b^2 d_1^2 + 2a^2 b^2 d_1 k_1 + 2\lambda abcd_1 d_4$.

Noting that $H_3 < 0$, this is not biologically relevant, thus the steady state $E_3$ is not considered. When $R_0 > 1$, the equilibrium $E_1$ exists. We define the reproduction rate of the CTL $R_{CTL}$ immune response by

$$R_{CTL} = \frac{cY_1}{b} = \frac{cd_4\lambda R_0(1 - \frac{1}{R_0})}{abd_1 + abk_1(1 - \frac{1}{R_0})}.$$

Note that the endemic state $E_2$ exists when $R_{CTL} > 1$. Indeed, if one considers $R_0 > 1$ then, in total absence of CTL immune response, the infected cell loaded per unit time is $\frac{d_4\lambda R_0(1 - \frac{1}{R_0})}{ad_1 + ak_1(1 - \frac{1}{R_0})}$. From the system in Equation (2), we have the CTL cells reproduced due to infected cells stimulating per unit time is $\frac{cd_4\lambda R_0(1 - \frac{1}{R_0})}{ad_1 + ak_1(1 - \frac{1}{R_0})} = cY_1$. The CTL charge during the lifespan of a CTL cell is $\frac{cd_4\lambda R_0(1 - \frac{1}{R_0})}{abd_1 + abk_1(1 - \frac{1}{R_0})} = R_{CTL}$. Thus, if $\frac{cd_4\lambda R_0(1 - \frac{1}{R_0})}{abd_1 + abk_1(1 - \frac{1}{R_0})} > 1$, we deduce the existence of $E_2$.

### 3. Global Stability

To prove the global stability of the uninfected and the infected steady states, we use the method of construction of Lyapunov functions developed in [12] and can claim the following result

**Theorem 1.** The disease-free equilibrium $E_f$ of the model in Equation (2) is globally asymptotically stable when $R_0 < 1$.

**Proof.** We define the function $G_f$ by

$$G_f(x, t) = S + \frac{d_2 + k_2}{k_2} Y + \frac{d_3(d_2 + k_2)}{ak_2} V + \frac{p}{c} \frac{d_2 + k_2}{k_2} Z.$$
Then, by using the equations of the system in Equation (2), the time derivative of $G_f$ verifies

$$ \frac{\partial G_f}{\partial t} \leq \frac{d_3 d_4 (d_2 + k_2)}{a k_2} (R_0 - 1) V. $$

Now, we define a Lyapunov function as follows

$$ L_f = \int_{\Omega} G_f dx. $$

Calculating the time derivative of $L_f$ along the positive solutions of the model in Equation (2), we obtain

$$ \begin{aligned} \frac{dL_f}{dt} &= \int_{\Omega} \frac{\partial G_f}{\partial t} dx \\ &\leq \int_{\Omega} \left( \frac{d_3 d_4 (d_2 + k_2)}{a k_2} (R_0 - 1) V \right) dx. \end{aligned} $$

Thus, if $R_0 < 1$ implies that $\frac{dL_f}{dt} \leq 0$. The largest compact invariant is

$$ E = \{(H, S, Y, V, Z) | V = 0\}, $$

according to LaSalle's invariance principle, $\lim_{t \to +\infty} V(x, t) = 0$, the limit system of equations is

$$ \left\{ \begin{aligned} \frac{\partial H(x,t)}{\partial t} &= \lambda - d_1 H(x,t), \\ \frac{\partial S(x,t)}{\partial t} &= -d_2 S(x,t) - k_2 S(x,t), \\ \frac{\partial Y(x,t)}{\partial t} &= k_2 S(x,t) - d_3 Y(x,t) - pY(x,t)Z(x,t), \\ \frac{\partial Z(x,t)}{\partial t} &= cY(x,t)Z(x,t) - bZ(x,t). \end{aligned} \right. $$

For simplicity, we use the same notation,

$$ G_f(H, S, Y, Z) = \frac{1}{H_0} \left( H - H_0 - H_0 \ln \frac{H}{H_0} \right) + S + \frac{d_2 + k_2}{k_2} Y + \frac{p}{c} \frac{d_2 + k_2}{k_2} Z. $$

Since $H_0 = \frac{\lambda}{d_1}$,

$$ \frac{\partial G_f}{\partial t} (H, S, Y, Z) = d_1 \left( 2 - \frac{H}{H_0} - \frac{H_0}{H} \right) - \frac{d_3(d_2 + k_2)}{k_2} Y - \frac{pb}{c} \frac{d_2 + k_2}{k_2} Z. $$

We define another Lyapunov function

$$ L_f = \int_{\Omega} G_f dx, $$

then, the time derivative of $L_f$ satisfies

$$ \begin{aligned} \frac{dL_f}{dt} &= \int_{\Omega} \frac{\partial G_f}{\partial t} dx \\ &\leq \int_{\Omega} \left( d_1 \left( 2 - \frac{H}{H_0} - \frac{H_0}{H} \right) - \frac{d_3(d_2 + k_2)}{k_2} Y - \frac{pb}{c} \frac{d_2 + k_2}{k_2} Z \right) dx. \end{aligned} $$
Since the arithmetic mean is greater than or equal to the geometric mean, it follows

$$2 - \frac{H}{H_0} - \frac{H_0}{H} \le 0,$$

therefore $\frac{dL_f}{dt} \le 0$ and the equality holds if $H = H_0$ and $S = Y = Z = 0$, which completes the proof. □

Now, we are interested in the stability of the infected steady state $E_1$. Let us state the following theorem.

**Theorem 2.** The infected steady state $E_1$ of the model in Equation (2) is globally asymptotically stable when $R_{CTL} \le 1 < R_0$. In this case, the other infected steady state $E_2$ does not exist.

**Proof.** Firstly, we define the function

$$G_1(x,t) = H - H_1 - \int_{H_1}^{H} \frac{(d_2 + k_2)S_1}{k_1 u V_1} du + S - S_1 - S_1 \ln \frac{S}{S_1} \\ + \frac{d_2 + k_2}{k_2} \left( Y - Y_1 - Y_1 \ln \left( \frac{Y}{Y_1} \right) \right) + \frac{d_3(d_2 + k_2)}{ak_2} \left( V - V_1 - V_1 \ln \left( \frac{V}{V_1} \right) \right) \\ + \frac{p}{c} \frac{d_2 + k_2}{k_2} Z.$$

Using the same technique proposed in [12], we obtain

$$\frac{\partial G_1}{\partial t} = - \frac{d_1 H_1}{H(H_1 + V_1)} (H - H_1)^2 \\ - (d_2 + k_2) S_1 \left( \frac{H(V - V_1)^2}{V_1(H + V_1)(H + V)} \right) \\ + (d_2 + k_2) S_1 \left( 5 - \frac{H_1}{H} \frac{H + V_1}{H_1 + V_1} - \frac{S_1}{S} \frac{HV}{H_1 V_1} \frac{H_1 + V_1}{H + V} - \frac{SY_1}{S_1 Y} - \frac{YV_1}{Y_1 V} - \frac{H + V}{H + V_1} \right) \\ + pZ \frac{d_2 + k_2 b}{k_2 c} (R_{CTL} - 1).$$

Now, let us consider the following Lyapunov function

$$L_1 = \int_{\Omega} G_1 dx,$$

then we deduce

$$\frac{dL_1}{dt} = \int_{\Omega} \frac{\partial G_1}{\partial t} dx \\ = \int_{\Omega} \left( - \frac{d_1 H_1}{H(H_1 + V_1)} (H - H_1)^2 - (d_2 + k_2) S_1 \left( \frac{H(V - V_1)^2}{V_1(H + V_1)(H + V)} \right) \right. \\ + (d_2 + k_2) S_1 \left( 5 - \frac{H_1}{H} \frac{H + V_1}{H_1 + V_1} - \frac{S_1}{S} \frac{HV}{H_1 V_1} \frac{H_1 + V_1}{H + V} - \frac{SY_1}{S_1 Y} - \frac{YV_1}{Y_1 V} - \frac{H + V}{H + V_1} \right) \\ \left. + pZ \frac{d_2 + k_2 b}{k_2 c} (R_{CTL} - 1) \right) dx.$$
Again, since the arithmetic mean is greater than or equal to the geometric mean, it follows

$$5 - \frac{H_1}{H} \frac{H + V_1}{H_1 + V_1} - \frac{S_1}{S} \frac{HV}{H_1 V_1} \frac{H_1 + V_1}{H + V} - \frac{SY_1}{S_1 Y} - \frac{YV_1}{Y_1 V} - \frac{H + V}{H + V_1} \le 0.$$

In addition, when $R_{CTL} < 1$, which means that $\frac{dL_1}{dt} \le 0$.

Therefore, by Lyapunov-LaSalle invariance theorem, $E_1$ is globally asymptotically stable when $R_0 > 1$ and $R_{CTL} \le 1$. □

To prove the stability of $E_2$ equilibrium, let us state the following theorem

**Theorem 3.** The infected steady state $E_2$ of the model in Equation (2) is globally asymptotically stable when $R_0 > 1$ and $R_{CTL} > 1$. In this case, the other infected steady state $E_1$ is unstable.

**Proof.** We consider the following function

$$G_2(x,t) = H - H_2 - \int_{H_2}^{H} \frac{(d_2 + k_2)S_2}{\frac{k_1 u V_2}{u + V_2}} du + S - S_2 - S_2 \ln \frac{S}{S_2} \\ + \frac{d_2 + k_2}{k_2} \left( Y - Y_2 - Y_2 \ln \left( \frac{Y}{Y_2} \right) \right) + \frac{d_3(d_2 + k_2) + (d_2 + k_2)pZ_2}{ak_2} \\ \times \left( V - V_2 - V_2 \ln \left( \frac{V}{V_2} \right) \right) + \frac{p}{c} \frac{d_2 + k_2}{k_2} \left( Z - Z_2 - Z_2 \ln \left( \frac{Z}{Z_2} \right) \right).$$

Then, we have

$$\frac{\partial G_2}{\partial t} = - \frac{d_1 V_2}{H(H_2 + V_2)} (H - H_2)^2 \\ - (d_2 + k_2) S_2 \left( \frac{H(V - V_2)^2}{V_2(H + V_2)(H + V)} \right) \\ + (d_2 + k_2) S_2 \left( 5 - \frac{H_2}{H} \frac{H + V_2}{H_2 + V_2} - \frac{S_2}{S} \frac{HV}{H_2 V_2} \frac{H_2 + V_2}{H + V} - \frac{SY_2}{S_2 Y} - \frac{YV_2}{Y_2 V} - \frac{H + V}{H + V_2} \right).$$

As a result, we define a Lyapunov function as follows

$$L_2 = \int_{\Omega} G_2 dx,$$

then,

$$\frac{dL_2}{dt} = \int_{\Omega} \frac{\partial G_2}{\partial t} dx \\ = \int_{\Omega} \left( - \frac{d_1 V_2}{H(H_2 + V_2)} (H - H_2)^2 - (d_2 + k_2) S_2 \left( \frac{H(V - V_2)^2}{V_2(H + V_2)(H + V)} \right) \right. \\ \left. + (d_2 + k_2) S_2 \left( 5 - \frac{H_2}{H} \frac{H + V_2}{H_2 + V_2} - \frac{S_2}{S} \frac{HV}{H_2 V_2} \frac{H_2 + V_2}{H + V} - \frac{SY_2}{S_2 Y} - \frac{YV_2}{Y_2 V} - \frac{H + V}{H + V_2} \right) \right) dx,$$

since the arithmetic mean is greater than or equal to the geometric mean, it follows

$$5 - \frac{H_2}{H} \frac{H + V_2}{H_2 + V_2} - \frac{S_2}{S} \frac{HV}{H_2 V_2} \frac{H_2 + V_2}{H + V} - \frac{SY_2}{S_2 Y} - \frac{YV_2}{Y_2 V} - \frac{H + V}{H + V_2} \le 0,$$

which means that $\frac{dL_2}{dt} \le 0$, and the equality holds when $H = H_2$, $S = S_2$, $Y = Y_2$, $V = V_2$, and $Z = Z_2$. By the LaSalle invariance principle, the endemic point $E_2$ is globally stable. □
4. Numerical Simulations

In this section, we present the results of numerical simulations to validate the theoretical results of the previous section. We used the finite difference numerical method with Euler explicit scheme. The convergence of our numerical method was tested by successively decreasing the time and space steps. The values of parameters are given in Appendix A.

We considered the one-dimensional interval $0 \le x \le L$ and time $0 < t \le T$, where $L = 20$ (dimensionless space units) and $T = 200$ days. The initial conditions were chosen space dependent to illustrate behavior of spatially inhomogeneous solutions. We used an explicit numerical method with the space step $h_x = 0.01$ and time step $h_t = 0.1$. The program was implemented with Matlab (2014a, MathWorks, Natick, MA, USA).

Figure 1 shows spatiotemporal dynamics of uninfected cells (left) and the maximal and the minimal values of the virus concentration in space as a function of time (right). For the values of parameters considered in this example, the basic reproduction number $R_0 = 0.22 < 1$, which implies that the virus-free equilibrium is stable. Therefore, as expected, solution converges toward the equilibrium $E_f = (8.27 \times 10^2, 0, 0, 0, 0)$.

<figure><img src="image_5.png" /><figcaption>**Figure 1.** Dynamics of solution for $\lambda = 10, d_1 = 0.0139, k_1 = 0.04, d_2 = 0.0495, k_2 = 1.1, d_3 = 0.5776, a = 2, d_4 = 0.6, p = 0.0024, c = 0.15$, and $b = 0.5$. The concentration of uninfected cells is shown as a function of $x$ and $t$ (left). The maximal and the minimal value of the virus concentration with respect to $x$ are shown as functions of time (right).</figcaption></figure>

To illustrate convergence to the endemic equilibrium point $E_1$, we considered the values of parameters presented in Figure 2. In this case, the basic reproduction number is greater than 1, $R_0 = 11.05 > 1$, and the immune reproduction number is less than 1, $R_{CTL} = 3.596 \times 10^{-1} < 1$. Therefore, the free-immune endemic equilibrium $E_1 = (19.96, 5.98 \times 10^{-1}, 1.14, 199.78, 0)$ is globally stable.
<figure><img src="image_6.png" /><figcaption>**Figure 2.** Dynamics of solution for $\lambda = 1$, $d_1 = 0.0139$, $k_1 = 0.04$, $d_2 = 0.0495$, $k_2 = 1.1$, $d_3 = 0.5776$, $a = 100$, $d_4 = 0.6$, $p = 0.0024$, $c = 0.15$, and $b = 0.5$. The concentration of uninfected cells is shown as a function of $x$ and $t$ (**left**). The maximal and the minimal value of the virus concentration with respect to $x$ are shown as functions of time (**right**).</figcaption></figure>

In the case considered in Figure 3, we obtain $R_0 = 11.05 > 1$ and $R_{CTL} = 4.13 > 1$. Consequently, we observe that the endemic equilibrium $E_2 = (285.12, 6.55, 3.33, 555.55, 660.86)$ is stable. This, numerical results support the theoretical findings.

<figure><img src="image_7.png" /><figcaption>**Figure 3.** Dynamics of solution for $\lambda = 10$, $d_1 = 0.0139$, $k_1 = 0.04$, $d_2 = 0.0495$, $k_2 = 1.1$, $d_3 = 0.5776$, $a = 100$, $d_4 = 0.6$, $p = 0.0024$, $c = 0.15$, and $b = 0.5$. The concentration of uninfected cells is shown as a function of $x$ and $t$ (**left**). The maximal and the minimal value of the virus concentration with respect to $x$ are shown as functions of time (**right**).</figcaption></figure>

**5. Discussion and Conclusions**

In this work, we study a model of viral infection in the presence of CTL cells and latently infected cells. We take into consideration not only the variation in time, but also spatial variation of virus distribution and its diffusion, where uninfected cells, latently infected cells, infected cells, and CTL cells do not exhibit any such spatial mobility. As a result, we show the stability of the free-disease equilibrium using the construction of a Lyapunov function when the reproduction number is less than one ($R_0 < 1$). If $R_0 > 1$, then two scenarios are established. If $R_{CTL} < 1$, then the equilibrium point $E_1$ is globally asymptotically stable, while, for $R_{CTL} > 1$, the equilibrium point $E_2$ is globally asymptotically stable. We also performed numerical simulations of our model to illustrate behavior of solutions and to confirm the theoretical results. It was established that spatial diffusion of free viruses
has no effect on the stability of the steady states. However, the effect appears only for the first days of infection observation.

**Author Contributions:** All authors contributed equally to this work. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by Centre National de la Recherche Scientifique, grant number PICS 244832.

**Acknowledgments:** The work was partially supported by the France–Morocco program PICS 244832, “RUDN University Program 5-100” and the French–Russian project PRC2307.

**Conflicts of Interest:** The authors declare no conflict of interest.

## Appendix A. The Values of Parameters Used in the Numerical Simulations

The values of parameters used in the numerical simulations are given in Table A1:

<table><caption>Table A1. The values, units and meaning of all used parameters.</caption><thead><tr><th>Parameters</th><th>Units</th><th>Meaning</th><th>Value</th><th>References</th></tr></thead><tbody><tr><td>λ</td><td>cells µL<sup>-1</sup> day<sup>-1</sup></td><td>Source rate of CD4+ T cells</td><td>[0, 10]</td><td>[27]</td></tr><tr><td>k<sub>1</sub></td><td>µL virion<sup>-1</sup> day<sup>-1</sup></td><td>Average of infection</td><td>[2.5 × 10<sup>-4</sup>, 0.5]</td><td>[12]</td></tr><tr><td>d<sub>1</sub></td><td>day<sup>-1</sup></td><td>Decay rate of healthy cells</td><td>0.0139</td><td>[12]</td></tr><tr><td>d<sub>2</sub></td><td>day<sup>-1</sup></td><td>Death rate of exposed CD4+ T cells</td><td>0.0495</td><td>[12]</td></tr><tr><td>k<sub>2</sub></td><td>day<sup>-1</sup></td><td>The rate that exposed cells become infected CD4+ T cells</td><td>1.1</td><td>[12]</td></tr><tr><td>d<sub>3</sub></td><td>day<sup>-1</sup></td><td>Death rate of infected CD4+ T cells, not by CTL killing</td><td>0.5776</td><td>[12]</td></tr><tr><td>a</td><td>day<sup>-1</sup></td><td>The rate of production the virus by infected CD4+ T cells</td><td>[2, 1250]</td><td>[12]</td></tr><tr><td>d<sub>4</sub></td><td>day<sup>-1</sup></td><td>Clearance rate of virus</td><td>[0.3466, 2.4]</td><td>[12]</td></tr><tr><td>p</td><td>µL cell<sup>-1</sup> day<sup>-1</sup></td><td>Clearance rate of infection</td><td>0.0024</td><td>[28]</td></tr><tr><td>c</td><td>cells cell<sup>-1</sup> day<sup>-1</sup></td><td>Activation rate CTL cells</td><td>0.15</td><td>[28]</td></tr><tr><td>b</td><td>day<sup>-1</sup></td><td>Death rate of CTL cells</td><td>0.5</td><td>[28]</td></tr><tr><td>d</td><td>mm<sup>2</sup>day<sup>-1</sup></td><td>Diffusion coefficient</td><td>0.01</td><td>–</td></tr></tbody></table>

## References

1. Perelson, A.S.; Neumann, A.U.; Markowitz, M.; Leonard, J.M.; Ho, D.D. HIV-1 dynamics in vivo: Virion clearance rate, infected cell life-span, and viral generation time. *Science* **1996**, *271*, 1582–1586. [CrossRef] [PubMed]

2. Adams, B.M.; Banks, H.T.; Davidian, M.; Kwon, H.D.; Tran, H.T.; Wynne, S.N.; Rosenberg, E.S. HIV dynamics: Modeling, data analysis, and optimal treatment protocols. *J. Comput. Appl. Math.* **2005**, *184*, 10–49. [CrossRef]

3. Burchell, A.N.; Winer, R.L.; de Sanjosé, S.; Franco, E.L. Epidemiology and transmission dynamics of genital HPV infection. *Vaccine* **2006**, *24*, S52–S61. [CrossRef] [PubMed]

4. Elbasha, E.H.; Dasbach, E.J.; Insinga, R.P. A multi-type HPV transmission model. *Bull. Math. Biol.* **2008**, *70*, 2126–2176. [CrossRef] [PubMed]

5. Li, M.; Zu, J. The review of differential equation models of HBV infection dynamics. *J. Virol. Methods* **2019**, *266*, 103–113. [CrossRef] [PubMed]

6. Meskaf, A.; Allali, K.; Tabit, Y. Optimal control of a delayed hepatitis B viral infection model with cytotoxic T-lymphocyte and antibody responses. *Int. J. Dyn. Control.* **2017**, *5*, 893–902. [CrossRef]

7. Wodarz, D. Hepatitis C virus dynamics and pathology: The role of CTL and antibody responses. *J. Gen. Virol.* **2003**, *84*, 1743–1750. [CrossRef]

8. Layden, J.E.; Layden, T.J. How can mathematics help us understand HCV? *Gastroenterology* **2001**, *120*, 1546–1549. [CrossRef]

9. Nowak, M.A.; Bangham, C.R.M. Population dynamics of immune responses to persistent viruses. *Science* **1996**, *272*, 74–79. [CrossRef]

10. Sun, Q.; Min, L.; Kuang, Y. Global stability of infection-free state and endemic infection state of a modified human immunodeficiency virus infection model. *IET Syst. Biol.* **2015**, *9*, 95–103. [CrossRef]

11. Sun, Q.; Min, L. Dynamics Analysis and Simulation of a Modified HIV Infection Model with a Saturated Infection Rate. *Comput. Math. Methods Med.* **2014**. [CrossRef] [PubMed]

12. Allali, K.; Danane, J.; Kuang, Y. Global Analysis for an HIV Infection Model with CTL Immune Response and Infected Cells in Eclipse Phase. *Appl. Sci.* **2017**, *7*, 861. [CrossRef]
13. Smith, H.L.; De Leenheer, P. Virus dynamics: A global analysis. *SIAM J. Appl. Math.* **2003**, *63*, 1313–1327.
[CrossRef]

14. Daar, E.S.; Moudgil, T.; Meyer, R.D.; Ho, D.D. Transient highlevels of viremia in patients with primary human immunodeficiency virus type 1. *N. Engl. J. Med.* **1991**, *324*, 961–964. [CrossRef]

15. Kahn, J.O.; Walker, B.D. Acute human immunodeficiency virus type 1 infection. *N. Engl. J. Med.* **1998**, *339*, 33–39. [CrossRef]

16. Bocharov, G.; Volpert, V.; Ludewig, B.; Meyerhans, A. *Mathematical Immunology of Virus Infections*; Springer International Publishing: Cham, Switzerland, 2018; doi:10.1007/978-3-319-72317-4; [CrossRef]

17. Kaufmann, G.R.; Cunningham, P.; Kelleher, A.D.; Zauders, J.; Carr, A.; Vizzard, J.; Law, M.; Cooper, D.A. The Sydney Primary HIV Infection Study Group. Patterns of viral dynamics during primary human immunodeficiency virus type 1 infection. *J. Infect. Dis.* **1998**, *178*, 1812–1815. [CrossRef]

18. Schacker, T.; Collier, A.; Hughes, J.; Shea, T.; Corey, L. Clinical and epidemiologic features of primary HIV infection. *Ann. Int. Med.* **1996**, *125*, 257–264. [CrossRef]

19. Britton, N.F. *Essential Mathematical Biology*; Springer: London, UK, 2003.

20. Grebennikov, D.; Bouchnita, A.; Volpert, V.; Bessonov, N.; Meyerhans, A.; Bocharov, G. Spatial lymphocyte dynamics in lymph nodes predicts the CTL frequency needed for HIV infection control. *Front. Immunol.* **2019**, *10*, 1213. [CrossRef]

21. Maziane, M.; Hattaf, K.; Yousfi, N. Dynamics of a class of HIV infection models with cure of infected cells in eclipse stage. *Acta Biotheor.* **2015**, *63*, 363–380. [CrossRef]

22. Hattaf, K.; Yousfi, N. A numerical method for delayed partial differential equations describing infectious diseases. *Comput. Math. Appl.* **2016**, *72*, 2741–2750. [CrossRef]

23. Trofimchuk, S.; Volpert, V. Traveling waves for a bistable reaction-diffusion equation with delay. *SIAM J. Math. Anal.* **2018**, *50*, 1175–1199. [CrossRef]

24. Bocharov, G.; Meyerhans, A.; Bessonov, N.; Trofimchuk, S.; Volpert, V. Interplay between reaction and diffusion processes in governing the dynamics of virus infections. *J. Theor. Biol.* **2018**, *457*, 221–236. [CrossRef] [PubMed]

25. Bessonov, N.; Bocharov, G.; Touaoula, T.M.; Trofimchuk, S.; Volpert, V. Delay reaction-diffusion equation for infection dynamics. *Discret. Contin. Dyn. Syst.-B* **2019**, *24*, 2073–2091. [CrossRef]

26. Goudon, T.; Lagoutiere, F.; Tine, L.M. The Lifschitz-Slyozov equation with space-diffusion of monomers. *Kinet. Relat. Model.* **2012**, *5*, 325–355. [CrossRef]

27. Wang, Y.; Zhou, Y.; Wu, J.; Heffernan, J. Oscillatory viral dynamics in a delayed HIV pathogenesis model. *Math. Biosci.* **2009**, *219*, 104–112. [CrossRef]

28. Zhu, H.; Luo, Y.; Chen, M. Stability and Hopf bifurcation of a HIV infection model with CTL-response delay. *Comput. Math. Appl.* **2011**, *62*, 3091–3102. [CrossRef]

<figure><img src="image_8.png" /></figure>