Global Properties of a Diffusive SARS-CoV-2 Infection Model with Antibody and Cytotoxic T-Lymphocyte Immune Responses

Abstract

A severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) infection can lead to morbidity and mortality. SARS-CoV-2 infects the epithelial cells of the respiratory tract and causes coronavirus disease 2019 (COVID-19). The immune system’s response plays a significant role in viral progression. This article develops and analyzes a system of partial differential equations (PDEs), which describe the in-host dynamics of SARS-CoV-2 under the effect of cytotoxic T-lymphocyte (CTL) and antibody immune responses. The model characterizes the interplay between six compartments, healthy epithelial cells (ECs), latent infected ECs, active infected ECs, free SARS-CoV-2 particles, CTLs, and antibodies. We consider the logistic growth of healthy ECs. We first investigate the properties of the model’s solutions, then, we calculate all steady states and determine the conditions of their existence and global stability. The global asymptotic stability is examined by constructing Lyapunov functions. The analytical findings are supported via numerical simulations.

1. Introduction

Recently, many viruses have spread that infect the human body, which can lead to illness and death and thus have a great impact on health and the global economy. These viruses include human viral infections such as human T-cell lymphotropic virus (HTLV), human immunodeficiency virus (HIV), Ebola virus, hepatitis B virus, hepatitis C virus, influenza virus, chikungunya virus, Middle East Respiratory Syndrome coronavirus (MERS-CoV), Zika virus, and dengue virus. At the end of 2019, the world witnessed the emergence of a new respiratory virus in Wuhan, China, called severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), which causes coronavirus disease 2019 (COVID-19). Within a few months, this viral infection had spread to most countries in the world and infected many people of different ages. Since the beginning of the spread of the virus, scientists from many disciplines have focused on finding ways to confront this pandemic such as by applying suitable control methods to reduce viral transmission, synthesizing antiviral drugs, and synthesizing vaccinations [1,2]. The WHO approved eleven vaccines for COVID-19 for emergency use [3]. Disease progression and outcome in COVID-19 are highly dependent on the host immune response, particularly in the elderly in whom immunosenescence may predispose them to increased risk of infection [4]. Immunosenescence enhances the susceptibility to viral infections and renders vaccinations less effective [5].

SARS-CoV-2 is a single-stranded positive-sense RNA virus and is a member of the Coronaviridae family. It is an airborne-transmitted virus and can infect the upper and lower respiratory tracts [6]. The epithelial cells (ECs) of the host respiratory tract are the target for SARS-CoV-2 particles [7]. The in-host dynamics of SARS-CoV-2 were mathematically modeled in [8] as

$$\left\{\begin{array}{c}\frac{dE\left(t\right)}{dt}=-\stackrel{\mathrm{SARS-CoV-}2\mathrm{infectious}\mathrm{transmission}}{\overbrace{\psi E\left(t\right)S\left(t\right)}},\hfill \\ \frac{dW\left(t\right)}{dt}=\stackrel{\mathrm{SARS-CoV-}2\mathrm{infectious}\mathrm{transmission}}{\overbrace{\psi E\left(t\right)S\left(t\right)}}-\stackrel{\mathrm{latent}\mathrm{transition}}{\overbrace{\vartheta W\left(t\right)}},\hfill \\ \frac{dP\left(t\right)}{dt}=\stackrel{\mathrm{latent}\mathrm{transition}}{\overbrace{\vartheta W\left(t\right)}}-\stackrel{\mathrm{natural}\mathrm{death}}{\overbrace{{\beta }_{P}P\left(t\right)}},\hfill \\ \frac{dS\left(t\right)}{dt}=\stackrel{\mathrm{SARS-CoV-}2\mathrm{production}}{\overbrace{\varrho P\left(t\right)}}-\stackrel{\mathrm{natural}\mathrm{death}}{\overbrace{{\beta }_{S}S\left(t\right)}},\hfill \end{array}\right.$$

(1)

where E, W, P, and S represent the concentrations of healthy ECs, latent infected ECs, active infected ECs, and SARS-CoV-2 particles, respectively. Li et al. [9] considered the regeneration and death of uninfected ECs

$$\frac{dE\left(t\right)}{dt}={\beta }_E(E\left(0\right)-E\left(t\right))-\psi E\left(t\right)S\left(t\right),$$

where $E\left(0\right)$is the initial concentration of healthy ECs.

Several extensions and modifications have been made on the SARS-CoV-2 infection models presented in [8,9] by taking into consideration the effect of different drug therapies [10,11,12] and the influence of time delay [13]. The innate immune response represents the first line of defense that recognizes the antigens and activates the adaptive immune response. B cells and T cells are important components of the adaptive immune response. It has been reported in [14] that both B- and T-cell responses against SARS-CoV-2 are detected in the blood around 1 week after the onset of COVID-19 symptoms. B cells have the ability to produce antibodies specific to SARS-CoV-2, such as IgA, IgG, and IgM, in order to neutralize the virus at the infection site [15,16,17,18]. CD8${}^{+}$ T cells (also known as cytotoxic T lymphocytes (CTLs)) are important for directly attacking and killing SARS-CoV-2-infected epithelial cells, whereas CD4${}^{+}$ T cells are crucial to prime both B cells and CTLs. CD4${}^{+}$ T cells are also responsible for cytokine production to drive immune cell recruitment [14]. The impact of the immune response was introduced in the SARS-CoV-2 infection model in [19,20,21,22,23,24,25]. Elaiw et al. [13] added a logistic term for the proliferation of healthy ECs as

$$\frac{dE\left(t\right)}{dt}=\varphi -{\beta }_{E}E\left(t\right)+\upsilon E\left(t\right)\left(1-\frac{E\left(t\right)}{E_{max}}\right)-\psi E\left(t\right)S\left(t\right).$$

It was assumed that healthy ECs are regenerated with a constant $\varphi$ and are proliferated with a logistic growth rate $\upsilon E\left(1-\frac{E}{E_{max}}\right)$, where $\upsilon$ is the rate of growth and $E_{max}$ is the maximum capacity of healthy ECs in the human body.

Stability analysis for models describing the in-host dynamics of SARS-CoV-2 infection was conducted in [23,24,25,26,27,28]. Hattaf and Yousfi [23] studied the global stability of an in-host SARS-CoV-2 infection model with cell-to-cell transmission and the cytotoxic T-lymphocyte (CTL) immune response. The model included both lytic (destruction of infected cells [29]) and nonlytic (inhibition of viral replication by soluble mediators secreted by immune cells [29]) immune responses. A SARS-CoV-2 infection model with both CTL and antibody immunities was developed and analyzed in [24]. Mathematical analysis of the model presented in [9] was studied in [26]. Both local and global stability analyses of the model’s equilibria were established. Almocera et al. [25] studied the stability of the two-dimensional SARS-CoV-2 dynamics model with the immune response presented in [8]. Elaiw et al. [13] studied the global stability of a delayed SARS-CoV-2 dynamics model with the logistic growth of uninfected ECs and antibody immunity. In very recent works, the Lyapunov method was used to establish the global stability of co-infection models including SARS-CoV-2/HIV-1 [27], SARS-CoV-2/Influenza A virus [28], and SARS-CoV-2/malaria [30].

In all the above mentioned works, it was assumed that the viruses and cells are homogeneously distributed in the human body. However, such an assumption is unrealistic because the diffusion of viruses and cells causes spatial variations within the body. Spatial structure plays a major role in describing the dynamical behaviors of SARS-CoV-2 infection within a host. The effect of spatial structure on the SARS-CoV-2 dynamics was addressed in [31,32]. In these papers, it was assumed that SARS-CoV-2 infection is resisted by antibodies, whereas the effect of CTL immunity was neglected. Moreover, the proliferation of healthy ECs was not included. In the present article, we construct and analyze a diffusive SARS-CoV-2 infection model with both antibody and CTL immune responses. We prove the non-negativity and boundedness of the solutions. We consider the logistic growth of healthy ECs. We derive the threshold parameters that determine the existence and stability of the model’s steady states. The global stability of all steady states is established by constructing Lyapunov functions and applying LaSalle’s invariance principle (LIP). We illustrate our theoretical results with some numerical simulations.

2. Model Formulation

We propose a diffusive SARS-CoV-2 model with antibody and CTL immune responses. Let $E(x,t)$, $W(x,t)$, $P(x,t)$, $S(x,t)$, $H(x,t)$, and $M(x,t)$ be the concentrations of healthy ECs, latent infected cells, active infected ECs, SARS-CoV-2 particles, antibodies, and CTLs, respectively, at spatial location x and time t.

$$\begin{array}{cc}\hfill \frac{\partial E(x,t)}{\partial t}& =D_E\Delta E(x,t)+\varphi -{\beta }_{E}E(x,t)+\upsilon E(x,t)\left(1-\frac{E(x,t)}{E_{max}}\right)-\psi E(x,t)S(x,t),\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill \frac{\partial W(x,t)}{\partial t}& =D_W\Delta W(x,t)+\gamma \psi S(x,t)E(x,t)-\vartheta W(x,t)-{\beta }_{W}W(x,t),\hfill \end{array}$$

(3)

$$\begin{array}{cc}\hfill \frac{\partial P(x,t)}{\partial t}& =D_P\Delta P(x,t)+(1-\gamma )\psi S(x,t)E(x,t)+\vartheta W(x,t)-{\beta }_{P}P(x,t)-\sigma P(x,t)M(x,t),\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill \frac{\partial S(x,t)}{\partial t}& =D_S\Delta S(x,t)+\varrho P(x,t)-{\beta }_{S}S(x,t)-\lambda S(x,t)H(x,t),\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill \frac{\partial H(x,t)}{\partial t}& =D_H\Delta H(x,t)+\mu S(x,t)H(x,t)-{\beta }_{H}H(x,t),\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill \frac{\partial M(x,t)}{\partial t}& =D_M\Delta M(x,t)+\eta P(x,t)M(x,t)-{\beta }_{M}M(x,t),\hfill \end{array}$$

(7)

where $x=\left(x_1,x_2,\cdots ,x_m\right)\in \varpi$ and $t>0$. The spatial domain $\varpi \subset {\mathbb{R}}^m,$ $m\ge 1$ is connected and bounded and the boundary $\partial \varpi$ is smooth. $\Delta ={\displaystyle \frac{{\partial }^2}{\partial x^2}}$ is the Laplacian operator and $D_{\mathcal{Q}}$ is the diffusion coefficient of $\mathcal{Q}\in \{E,W,P,S,H,M\}$. Parameter $\gamma \in (0,1)$ is the part of the healthy ECs that enters the latent state and ${\beta }_{W}W$ is the death rate of the latent infected ECs. The active infected ECs are killed by CTLs at rate $\sigma MP$. The SARS-CoV-2 particles are neutralized by antibodies at the rate $\lambda HS$. The terms $\mu HS$ and $\eta MP$ represent the proliferation rates of antibodies and CTLs, respectively. The terms ${\beta }_{H}H$ and ${\beta }_{M}M$ denote the death rate of antibodies and CTLs, respectively. The initial conditions are given by

$$\begin{array}{cc}\hfill E(x,0)& ={\delta }_1\left(x\right),W(x,0)={\delta }_2\left(x\right),P(x,0)={\delta }_3\left(x\right),S(x,0)={\delta }_4\left(x\right),\hfill \\ \hfill H(x,0)& ={\delta }_5\left(x\right),M(x,0)={\delta }_6\left(x\right),x\in \overline{\varpi },\hfill \end{array}$$

(8)

where ${\delta }_i\left(x\right)$, $i=1,\cdots ,6,$ are non-negative and continuous functions. In addition, we consider the following homogeneous Neumann boundary conditions:

$$\frac{\partial E}{\partial \overrightarrow{n}}=\frac{\partial W}{\partial \overrightarrow{n}}=\frac{\partial P}{\partial \overrightarrow{n}}=\frac{\partial S}{\partial \overrightarrow{n}}=\frac{\partial H}{\partial \overrightarrow{n}}=\frac{\partial M}{\partial \overrightarrow{n}}=0,t>0,x\in \partial \varpi ,$$

(9)

where ${\displaystyle \frac{\partial }{\partial \overrightarrow{n}}}$ is the outward normal derivative on the boundary $\partial \varpi$.

3. Properties of Solutions

This section examines the existence, non-negativity, and boundedness of the solutions of Systems (2)–(7) with the initial conditions in (8) and the boundary conditions in (9). Moreover, it determines the existence conditions of the spatially homogeneous steady states of the model.

Theorem 1.

Assume that $D_E=D_W=D_P=D_S=D_H=D_M=D$. Then, Models (2)–(7) have unique, non-negative, and bounded solutions defined on $\overline{\varpi }\times [0,+\infty )$ for any given initial data satisfying (8).

Proof. 

Let $\mathbb{X}=BUM\left(\overline{\varpi },{\mathbb{R}}^6\right)$ be the set of all bounded and uniformly continuous functions from $\overline{\varpi }$ to ${\mathbb{R}}^6$, and ${\mathbb{X}}_{+}=BUM\left(\overline{\varpi },{\mathbb{R}}_{+}^6\right)\subset \mathbb{X}$. Hence, the positive cone ${\mathbb{X}}_{+}$ induces a partial order on $\mathbb{X}$. Define ${\parallel \theta \parallel }_{\mathbb{X}}=\underset{x\in \overline{\varpi }}{sup}\left|\theta \left(x\right)\right|$, where $|·|$ is the Euclidean norm on ${\mathbb{R}}^6$. It follows that the space $\left({\mathbb{X},\parallel ·\parallel }_{\mathbb{X}}\right)$ is a Banach lattice [33,34]. □

For any initial data $\delta ={({\delta }_1,{\delta }_2,{\delta }_3,{\delta }_4,{\delta }_5,{\delta }_6)}^{\prime }\in {\mathbb{X}}_{+}$, we define $\Gamma ={({\Gamma }_1,{\Gamma }_2,{\Gamma }_3,{\Gamma }_4,{\Gamma }_5,{\Gamma }_6)}^{\prime }:{\mathbb{X}}_{+}\to \mathbb{X}$ by

$$\begin{array}{cc}\hfill {\Gamma }_1\left(\delta \right)\left(x\right)& =\varphi -{\beta }_E{\delta }_1\left(x\right)+\upsilon {\delta }_1\left(x\right)\left(1-\frac{{\delta }_1\left(x\right)}{E_{max}}\right)-\psi {\delta }_1\left(x\right){\delta }_4\left(x\right),\hfill \\ \hfill {\Gamma }_2\left(\delta \right)\left(x\right)& =\gamma \psi {\delta }_1\left(x\right){\delta }_4\left(x\right)-\vartheta {\delta }_2\left(x\right)-{\beta }_W{\delta }_2\left(x\right),\hfill \\ \hfill {\Gamma }_3\left(\delta \right)\left(x\right)& =(1-\gamma )\psi {\delta }_1\left(x\right){\delta }_4\left(x\right)+\vartheta {\delta }_2\left(x\right)-{\beta }_P{\delta }_3\left(x\right)-\sigma {\delta }_3\left(x\right){\delta }_6\left(x\right),\hfill \\ \hfill {\Gamma }_4\left(\delta \right)\left(x\right)& =\varrho {\delta }_3\left(x\right)-{\beta }_S{\delta }_4\left(x\right)-\lambda {\delta }_4\left(x\right){\delta }_5\left(x\right),\hfill \\ \hfill {\Gamma }_5\left(\delta \right)\left(x\right)& =\mu {\delta }_4\left(x\right){\delta }_5\left(x\right)-{\beta }_H{\delta }_5\left(x\right),\hfill \\ \hfill {\Gamma }_6\left(\delta \right)\left(x\right)& =\eta {\delta }_3\left(x\right){\delta }_6\left(x\right)-{\beta }_M{\delta }_6\left(x\right).\hfill \end{array}$$

We observe that $\Gamma$ is locally Lipschitz on ${\mathbb{X}}_{+}$. Models (2)–(7) with Conditions (8) and (9) can be written as the following abstract functional DE

$$\begin{array}{c}{\displaystyle \frac{d\mathcal{Q}}{dt}}=D\mathcal{Q}+\Gamma \left(\mathcal{Q}\right),t>0,\hfill \\ \mathcal{Q}\left(0\right)=\delta \in {\mathbb{X}}_{+},\hfill \end{array}$$

where $\mathcal{Q}={(E,W,P,S,H,M)}^{\prime }$ and $D\mathcal{Q}={\left(D_E\Delta E,D_W\Delta W,D_P\Delta P,D_S\Delta S,D_H\Delta H,D_M\Delta M\right)}^{\prime }.$ can be shown that

$$\underset{h\to 0^{+}}{lim}\frac{1}{h}\mathrm{dist}\left(\delta \left(0\right)+h\Gamma \left(\delta \right),{\mathbb{X}}_{+}\right)=0,\mathrm{for}\mathrm{all}\delta \in {\mathbb{X}}_{+}.$$

It follows from [33,34,35] that for any $\delta \in {\mathbb{X}}_{+}$, Systems (2)–(7) with (8) and (9) have unique non-negative mild solutions $(E\left(x,t\right),W\left(x,t\right),P\left(x,t\right),S\left(x,t\right),H\left(x,t\right),M\left(x,t\right))$ defined on $\overline{\varpi }\times [0,{\mathcal{T}}_{max})$, where $[0,{\mathcal{T}}_{max})$ is the maximal existence time interval on which the solution exists. Further, this solution is also a classical solution for the given problem.

Now, we define

$$\Phi \left(x,t\right)=E\left(x,t\right)+W\left(x,t\right)+P\left(x,t\right)+\frac{{\beta }_P}{2\varrho }S\left(x,t\right)+\frac{{\beta }_P\lambda }{2\varrho \mu }H\left(x,t\right)+\frac{\sigma }{\eta }M\left(x,t\right).$$

Since $D_E=D_W=D_P=D_S=D_H=D_M=D$, then, using Systems (2)–(7) we obtain

$$\begin{array}{cc}\hfill \frac{\partial \Phi (x,t)}{\partial t}-D\Delta \Phi (x,t)& =\varphi -{\beta }_{E}E(x,t)+\upsilon E(x,t)\left(1-\frac{E(x,t)}{E_{max}}\right)-{\beta }_{W}W(x,t)-\frac{{\beta }_P}{2}P(x,t)\hfill \\ \hfill & -\frac{{\beta }_P{\beta }_S}{2\varrho }S(x,t)-\frac{{\beta }_P{\beta }_H\lambda }{2\varrho \mu }H(x,t)-\frac{{\beta }_M\sigma }{\eta }M(x,t)\hfill \\ \hfill & =-\frac{\upsilon }{E_{max}}{E}^2(x,t)+\upsilon E(x,t)+\varphi -{\beta }_{E}E(x,t)-{\beta }_{W}W(x,t)\hfill \\ \hfill & -\frac{{\beta }_P}{2}P(x,t)-\frac{{\beta }_P{\beta }_S}{2\varrho }S(x,t)-\frac{{\beta }_P{\beta }_H\lambda }{2\varrho \mu }H(x,t)-\frac{{\beta }_M\sigma }{\eta }M(x,t).\hfill \end{array}$$

Let us define $\Pi \left(E\right)=-{\displaystyle \frac{\upsilon }{E_{max}}}{E}^2+\upsilon E+\varphi$ and calculate ${\Pi }^{\prime }\left(E\right)$ as

$${\Pi }^{\prime }\left(E\right)=-\frac{2\upsilon }{E_{max}}E+\upsilon =0\Rightarrow E=\frac{E_{max}}{2}$$

and

$${\Pi }^{\prime \prime }\left(E\right)=-\frac{2\upsilon }{E_{max}}<0.$$

Then,

$$\Pi \left(\frac{E_{max}}{2}\right)=-\frac{\upsilon }{E_{max}}{\left(\frac{E_{max}}{2}\right)}^2+\upsilon \left(\frac{E_{max}}{2}\right)+\varphi =\frac{\upsilon E_{max}}{4}+\varphi .$$

Let $N_1={\displaystyle \frac{\upsilon E_{max}+4\varphi }{4}}>0$ and $q_1=min\{{\beta }_E,{\beta }_W,{\displaystyle \frac{{\beta }_P}{2}},{\beta }_S,{\beta }_H,{\beta }_M\}$, then, $\Phi (x,t)$ satisfies the following system:

$$\left\{\begin{array}{c}{\displaystyle \frac{\partial \Phi (x,t)}{\partial t}}-D\Delta \Phi (x,t)\le N_1-q_1\Phi (x,t),\hfill \\ {\displaystyle \frac{\partial \Phi }{\partial \overrightarrow{n}}}=0,\hfill \\ \Phi (x,0)={\delta }_1\left(x\right)+{\delta }_2\left(x\right)+{\delta }_3\left(x\right)+{\displaystyle \frac{{\beta }_P}{2\varrho }}{\delta }_4\left(x\right)+{\displaystyle \frac{{\beta }_P\lambda }{2\varrho \mu }}{\delta }_5\left(x\right)+{\displaystyle \frac{\sigma }{\eta }}{\delta }_6\left(x\right)\ge 0.\hfill \end{array}\right.$$

Let $\tilde{\Phi }\left(t\right)$ be a solution to the following ODE system

$$\left\{\begin{array}{c}{\displaystyle \frac{d\tilde{\Phi }\left(t\right)}{dt}}=N_1-q_1\tilde{\Phi }\left(t\right),\hfill \\ \tilde{\Phi }\left(0\right)=\underset{x\in \overline{\varpi }}{max}\Phi (x,0).\hfill \end{array}\right.$$

This gives that $\tilde{\Phi }\left(t\right)\le max\left\{{\displaystyle \frac{N_1}{q_1}},\underset{x\in \overline{\varpi }}{max}\Phi (x,0)\right\}$. According to the comparison principle [36], we have $\Phi (x,t)\le \tilde{\Phi }\left(t\right)$. Then, we obtain

$$\Phi (x,t)\le max\left\{\frac{N_1}{q_1},\underset{x\in \overline{\varpi }}{max}\Phi (x,0)\right\},$$

which implies that $E\left(x,t\right),$ $W\left(x,t\right),$ $P\left(x,t\right),$ $S\left(x,t\right),$ $H\left(x,t\right),$ and $M\left(x,t\right)$ are bounded on $\overline{\varpi }\times [0,{\mathcal{T}}_{max})$. From the standard theory for semi-linear parabolic systems, we deduce that ${\mathcal{T}}_{max}=+\infty$ [37]. This shows that solution $\left(E\left(x,t\right),W\left(x,t\right),P\left(x,t\right),S\left(x,t\right),H\left(x,t\right),M\left(x,t\right)\right)$ is defined for all $x\in \varpi ,$ $t>0$ and is also non-negative and unique.

4. Steady States

This section computes all steady states of Systems (2)–(7) and the threshold parameters that guarantee the existence of these steady states. Let $SS=(E,W,P,S,H,M)$ be any steady state of Systems (2)–(7) satisfying the following system of nonlinear equations:

$$\begin{array}{cc}\hfill 0& =\varphi -{\beta }_{E}E+\upsilon E\left(1-\frac{E}{E_{max}}\right)-\psi ES,\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill 0& =\gamma \psi ES-(\vartheta +{\beta }_W)W,\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill 0& =(1-\gamma )\psi ES+\vartheta W-{\beta }_{P}P-\sigma MP,\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill 0& =\varrho P-{\beta }_{S}S-\lambda HS,\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill 0& =\mu HS-{\beta }_{H}H,\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill 0& =\eta MP-{\beta }_{M}M.\hfill \end{array}$$

(15)

By solving Systems (10)–(15), we find that Systems (2)–(7) have the following steady states:

(i) Healthy steady state $S{S}_0=(E_0,0,0,0,0,0)$, where $E_0$ is the positive root of $\varphi -{\beta }_{E}E+\upsilon E\left(1-{\displaystyle \frac{E}{E_{max}}}\right)=0$ and is given by

$$E_0=\frac{E_{max}}{2\upsilon }\left[\upsilon -{\beta }_E+\sqrt{{(\upsilon -{\beta }_E)}^2+\frac{4\upsilon \varphi }{E_{max}}}\right].$$

(16)

The basic reproduction number ${\mathcal{R}}_0$ for Systems (2)–(7) is given by

$${\mathcal{R}}_0=\frac{\varrho \psi E_0}{{\beta }_P{\beta }_S}\rho .$$

where, $\rho ={\displaystyle \frac{\vartheta \gamma }{\vartheta +{\beta }_W}}+1-\gamma$.

(ii) Infected steady state with inactive immune responses $S{S}_1=(E_1,W_1,P_1,S_1,0,0)$, where

$$\begin{array}{cc}\hfill E_1& =\frac{{\beta }_P{\beta }_S}{\varrho \psi \rho }=\frac{E_0}{{\mathcal{R}}_0},\hfill \\ \hfill W_1& =\frac{\gamma }{\vartheta +{\beta }_W}\psi E_1{S}_1,\hfill \\ \hfill P_1& =\frac{{\beta }_S}{\varrho }{S}_1,\hfill \\ \hfill S_1& =\frac{\varphi \varrho \rho }{{\beta }_P{\beta }_S}+\frac{\upsilon }{\psi }-\left(\frac{{\beta }_E}{\psi }+\frac{\upsilon {\beta }_P{\beta }_S}{\varrho {\psi }^2{E}_{max}\rho }\right).\hfill \end{array}$$

Assume that ${\beta }_E-\upsilon +{\displaystyle \frac{\upsilon E_1}{E_{max}}}>0$, then, we obtain

$${\beta }_E-\upsilon +\frac{\upsilon }{E_{max}}\frac{{\beta }_P{\beta }_S}{\varrho \psi \rho }>0.$$

(17)

We note that

$$\begin{array}{cc}\hfill {\mathcal{R}}_0& >1⟺\frac{E_{max}}{2\upsilon }\left[\upsilon -{\beta }_E+\sqrt{{(\upsilon -{\beta }_E)}^2+\frac{4\upsilon \varphi }{E_{max}}}\right]>\frac{{\beta }_P{\beta }_S}{\varrho \psi \rho }\hfill \\ \hfill & ⟺\sqrt{{(\upsilon -{\beta }_E)}^2+\frac{4\upsilon \varphi }{E_{max}}}>\frac{2\upsilon {\beta }_P{\beta }_S}{\varrho \psi E_{max}\rho }-(\upsilon -{\beta }_E).\hfill \end{array}$$

From Inequality (17), we have ${\displaystyle \frac{2\upsilon {\beta }_P{\beta }_S}{\varrho \psi E_{max}\rho }}-(\upsilon -{\beta }_E)>0$. Then,

$$\begin{array}{cc}\hfill {\mathcal{R}}_0>1& ⟺\frac{4\upsilon \varphi }{E_{max}}>\frac{4{\upsilon }^2{\beta }_P{\beta }_S}{{\varrho }^2{\psi }^2{E}_{max}^2{\rho }^2}-\frac{4\upsilon {\beta }_P{\beta }_S}{\varrho \psi E_{max}\rho }(\upsilon -{\beta }_E)\hfill \\ \hfill & ⟺\upsilon \varphi >\frac{{\upsilon }^2{\beta }_P{\beta }_S}{{\varrho }^2{\psi }^2{E}_{max}{\rho }^2}-\frac{{\upsilon }^2{\beta }_P{\beta }_S}{\varrho \psi \rho }+\frac{\upsilon {\beta }_E{\beta }_P{\beta }_S}{\varrho \psi \rho }\hfill \\ \hfill & ⟺\frac{\varphi \varrho \rho }{{\beta }_P{\beta }_S}+\frac{\upsilon }{\psi }-\left(\frac{{\beta }_E}{\psi }+\frac{\upsilon {\beta }_P{\beta }_S}{\varrho {\psi }^2{E}_{max}\rho }\right)>0\hfill \\ \hfill & ⟺S_1>0.\hfill \end{array}$$

Thus, $S{S}_1$ exists when ${\mathcal{R}}_0>1$ and ${\beta }_E-\upsilon +{\displaystyle \frac{\upsilon E_1}{E_{max}}}>0$. At this steady state, the virus exists while the immune response is inhibited.

(iii) Infected steady state with only active antibody immunity $S{S}_2=(E_2,W_2,P_2,S_2,H_2,0)$, where

$$\begin{array}{cc}\hfill E_2& =\frac{E_{max}}{2\upsilon }\left[\upsilon -{\beta }_E-\frac{{\beta }_H\psi }{\mu }+\sqrt{{\left(\upsilon -{\beta }_E-\frac{{\beta }_H\psi }{\mu }\right)}^2+\frac{4\upsilon \varphi }{E_{max}}}\right],\hfill \\ \hfill W_2& =\frac{{\beta }_H\gamma \psi }{\mu (\vartheta +{\beta }_W)}\frac{E_{max}}{2\upsilon }\left[\upsilon -{\beta }_E-\frac{{\beta }_H\psi }{\mu }+\sqrt{{\left(\upsilon -{\beta }_E-\frac{{\beta }_H\psi }{\mu }\right)}^2+\frac{4\upsilon \varphi }{E_{max}}}\right]=\frac{{\beta }_H\gamma \psi E_2}{\mu (\vartheta +{\beta }_W)},\hfill \\ \hfill P_2& =\frac{{\beta }_H\psi }{\mu {\beta }_P}\frac{E_{max}}{2\upsilon }\left[\upsilon -{\beta }_E-\frac{{\beta }_H\psi }{\mu }+\sqrt{{\left(\upsilon -{\beta }_E-\frac{{\beta }_H\psi }{\mu }\right)}^2+\frac{4\upsilon \varphi }{E_{max}}}\right]\rho =\frac{{\beta }_H\psi E_2}{\mu {\beta }_P}\rho ,\hfill \\ \hfill S_2& =\frac{{\beta }_H}{\mu },\hfill \\ \hfill H_2& =\frac{{\beta }_S}{\lambda }\left[\frac{\varrho \psi }{{\beta }_P{\beta }_S}\frac{E_{max}}{2\upsilon }\left(\upsilon -{\beta }_E-\frac{{\beta }_H\psi }{\mu }+\sqrt{{\left(\upsilon -{\beta }_E-\frac{{\beta }_H\psi }{\mu }\right)}^2+\frac{4\upsilon \varphi }{E_{max}}}\right)\rho -1\right]\hfill \\ \hfill & =\frac{{\beta }_S}{\lambda }\left(\frac{\varrho \psi E_2}{{\beta }_P{\beta }_S}\rho -1\right).\hfill \end{array}$$

We note that $S{S}_2$ exists when ${\displaystyle \frac{\varrho \psi E_2}{{\beta }_P{\beta }_S}}\rho >1$. Let ${\mathcal{R}}_1^H$ be the antibody immunity activation number defined as

$${\mathcal{R}}_1^H=\frac{\varrho \psi E_2}{{\beta }_P{\beta }_S}\rho .$$

which determines when the antibody immunity is activated. Thus, $H_2={\displaystyle \frac{{\beta }_S}{\lambda }}\left({\mathcal{R}}_1^H-1\right)$. We note that $H_2>0$ when ${\mathcal{R}}_1^H>1$. Thus, $S{S}_2$ exists when ${\mathcal{R}}_1^H>1$.

(iv) Infected steady state with only active CTL immunity $S{S}_3=(E_3,W_3,P_3,S_3,0,M_3)$, where

$$\begin{array}{cc}\hfill E_3& =\frac{E_{max}}{2\upsilon }\left[\upsilon -{\beta }_E-\frac{{\beta }_M\varrho \psi }{{\beta }_S\eta }+\sqrt{{\left(\upsilon -{\beta }_E-\frac{{\beta }_M\varrho \psi }{{\beta }_S\eta }\right)}^2+\frac{4\upsilon \varphi }{E_{max}}}\right],\hfill \\ \hfill W_3& =\frac{{\beta }_M\gamma \varrho \psi }{{\beta }_S\eta (\vartheta +{\beta }_W)}\frac{E_{max}}{2\upsilon }\left[\upsilon -{\beta }_E-\frac{{\beta }_M\varrho \psi }{{\beta }_S\eta }+\sqrt{{\left(\upsilon -{\beta }_E-\frac{{\beta }_M\varrho \psi }{{\beta }_S\eta }\right)}^2+\frac{4\upsilon \varphi }{E_{max}}}\right]\hfill \\ \hfill & =\frac{{\beta }_M\gamma \varrho \eta \psi E_3}{{\beta }_S\eta (\vartheta +{\beta }_W)},\hfill \\ \hfill P_3& =\frac{{\beta }_M}{\eta },S_3=\frac{{\beta }_M\varrho }{{\beta }_S\eta },\hfill \\ \hfill M_3& =\frac{{\beta }_P}{\sigma }\left[\frac{\varrho \psi E_3}{{\beta }_P{\beta }_S}\left(\frac{\vartheta \gamma }{\vartheta +{\beta }_W}+1-\gamma \right)-1\right]=\frac{{\beta }_P}{\sigma }\left(\frac{\varrho \psi E_3}{{\beta }_P{\beta }_S}\rho -1\right).\hfill \end{array}$$

We note that $M_3>0$ when ${\displaystyle \frac{\varrho \psi E_3}{{\beta }_P{\beta }_S}}\rho >1$. Then, we define the CTL immunity activation number ${\mathcal{R}}_1^M$ as

$${\mathcal{R}}_1^M=\frac{\varrho \psi E_3}{{\beta }_P{\beta }_S}\rho .$$

Consequently, $S{S}_3$ exists when ${\mathcal{R}}_1^M>1$.

(v) Infected steady state with both active antibody and CTL immunities $S{S}_4=(E_4,W_4,P_4,S_4,H_4,M_4)$, where

$$\begin{array}{cc}\hfill E_4& =\frac{E_{max}}{2\upsilon }\left[\upsilon -{\beta }_E-\frac{{\beta }_H\psi }{\mu }+\sqrt{{\left(\upsilon -{\beta }_E-\frac{{\beta }_H\psi }{\mu }\right)}^2+\frac{4\upsilon \varphi }{E_{max}}}\right]=E_2,\hfill \\ \hfill W_4& =\frac{{\beta }_H\gamma \psi }{\mu (\vartheta +{\beta }_W)}\frac{E_{max}}{2\upsilon }\left[\upsilon -{\beta }_E-\frac{{\beta }_H\psi }{\mu }+\sqrt{{\left(\upsilon -{\beta }_E-\frac{{\beta }_H\psi }{\mu }\right)}^2+\frac{4\upsilon \varphi }{E_{max}}}\right]=\frac{{\beta }_H\gamma \psi E_4}{\mu (\vartheta +{\beta }_W)},\hfill \\ \hfill P_4& =\frac{{\beta }_M}{\eta }=P_3,S_4=\frac{{\beta }_H}{\mu }=S_2,\hfill \\ \hfill H_4& =\frac{{\beta }_M\varrho \mu }{{\beta }_H\eta \lambda }-\frac{{\beta }_S}{\lambda }=\frac{{\beta }_S}{\lambda }\left[\frac{{\beta }_M\varrho \mu }{{\beta }_S{\beta }_H\eta }-1\right],\hfill \\ \hfill M_4& =\frac{\eta }{{\beta }_M\sigma }\left[\frac{{\beta }_H\psi E_4}{\mu }\left(\frac{\vartheta \gamma }{\vartheta +{\beta }_W}+1-\gamma \right)-\frac{{\beta }_P{\beta }_M}{\eta }\right]=\frac{{\beta }_P}{\sigma }\left[\frac{{\beta }_H\eta \psi E_4}{{\beta }_P{\beta }_M\mu }\rho -1\right].\hfill \end{array}$$

We see that, $H_4$ and $M_4$ exist when ${\displaystyle \frac{{\beta }_M\varrho \mu }{{\beta }_S{\beta }_H\eta }}>1$ and ${\displaystyle \frac{{\beta }_H\eta \psi E_4}{{\beta }_P{\beta }_M\mu }}\rho >1$. Now, we define

$${\mathcal{R}}_2^M=\frac{{\beta }_H\eta \psi E_4}{{\beta }_P{\beta }_M\mu }\rho .$$

Hence, $H_4$ and $M_4$ can be rewritten as

$$H_4=\frac{{\beta }_S}{\lambda }\left[\frac{{\mathcal{R}}_1^H}{{\mathcal{R}}_2^M}-1\right],M_4=\frac{{\beta }_P}{\sigma }\left[{\mathcal{R}}_2^M-1\right].$$

Therefore, $S{S}_4$ exists when ${\mathcal{R}}_1^H>{\mathcal{R}}_2^M$ and ${\mathcal{R}}_2^M>1$. Here, ${\mathcal{R}}_2^M$ refers to the competing CTL immunity number.

Lemma 1.

Systems (2)–(7) have four threshold parameters ${\mathcal{R}}_0>0$, ${\mathcal{R}}_1^H>0$, ${\mathcal{R}}_1^M>0$, and ${\mathcal{R}}_2^M>0$, such that:

(i)

if ${\mathcal{R}}_0\le 1$, then there exists only one steady state $S{S}_0$;

(ii)

if ${\mathcal{R}}_1^H\le 1<{\mathcal{R}}_0$, ${\mathcal{R}}_1^M\le 1$ and ${\beta }_E-\upsilon +{\displaystyle \frac{\upsilon E_1}{E_{max}}}>0$, then there exist two steady states $S{S}_0$ and $S{S}_1$;

(iii)

if ${\mathcal{R}}_1^H>1$ and ${\mathcal{R}}_2^M\le 1$, then there exist only three steady states $S{S}_0$, $S{S}_1$, and $S{S}_2$;

(iv)

if ${\mathcal{R}}_1^M>1$ and ${\mathcal{R}}_1^H\le {\mathcal{R}}_2^M$, then there exist only three steady states $S{S}_0$, $S{S}_1$, and $S{S}_3$;

(v)

if ${\mathcal{R}}_1^H>{\mathcal{R}}_2^M>1$, then there exist five steady states $S{S}_0$, $S{S}_1$, $S{S}_2$, $S{S}_3$, and $S{S}_4$.

5. Global Properties

In this section, we prove the global asymptotic stability of all five steady states $S{S}_i$, $i=0,1,2,3,4$ by constructing Lyapunov functions [38]. We use the arithmetic–geometric mean inequality

$$\frac{1}{n}\sum _{i=1}^n{{\rm Y}}_i\ge \sqrt[n]{\prod _{i=1}^n{{\rm Y}}_i},{{\rm Y}}_i\ge 0,i=1,2,\cdots$$

(18)

Using Inequality (18), we obtain the following relations:

$$\begin{array}{cc}\hfill \frac{E_j}{E}-\frac{W_{j}ES}{W{E}_j{S}_j}-\frac{P_{j}W}{P{W}_j}-\frac{S_{j}P}{S{P}_j}& \ge 4,\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill \frac{E_j}{E}-\frac{P_{j}ES}{P{E}_j{S}_j}-\frac{S_{j}P}{S{P}_j}& \ge 3,j=1,\cdots ,4,\hfill \end{array}$$

(20)

Let a function ${\mathcal{G}}_j(E,W,P,S,H,M)$ and define

$${\tilde{\mathcal{G}}}_j\left(t\right)={\int }_{\varpi }{\mathcal{G}}_j(x,t)dx,j=0,1,\cdots ,4.$$

Let ${\tilde{\Xi }}_j$ be the largest invariant subset of

$${\Xi }_j=\left\{(E,W,P,S,H,M):\frac{d{\tilde{\mathcal{G}}}_j}{dt}=0\right\},j=0,1,\cdots ,4.$$

From the Neumann boundary conditions in (9) and the Divergence Theorem we obtain

$$\begin{array}{cc}\hfill 0& ={\int }_{\partial \varpi }\nabla \mathcal{Q}·\overrightarrow{n}dx={\int }_{\varpi }\mathrm{div}(\nabla \mathcal{Q})dx={\int }_{\varpi }\Delta \mathcal{Q}dx,\hfill \\ \hfill 0& ={\int }_{\partial \varpi }\frac{1}{\mathcal{Q}}\nabla \mathcal{Q}·\overrightarrow{n}dx={\int }_{\varpi }\mathrm{div}(\frac{1}{\mathcal{Q}}\nabla \mathcal{Q})dx={\int }_{\varpi }\left(\frac{\Delta \mathcal{Q}}{\mathcal{Q}}-\frac{{\parallel ▿\mathcal{Q}\parallel }^2}{{\mathcal{Q}}^2}\right)dx,\hfill \end{array}$$

for $\mathcal{Q}\in \{E,W,P,S,H,M\}.$ Thus, we obtain

$$\begin{array}{cc}\hfill {\int }_{\varpi }\Delta \mathcal{Q}dx& =0,\hfill \\ \hfill {\int }_{\varpi }\frac{\Delta \mathcal{Q}}{\mathcal{Q}}dx& ={\int }_{\varpi }\frac{{\parallel ▿\mathcal{Q}\parallel }^2}{{\mathcal{Q}}^2}dx,\mathrm{for}\mathcal{Q}\in \{E,W,P,S,H,M\}.\hfill \end{array}$$

(21)

For convenience, we drop the input notation, i.e., $(E,W,P,S,H,M)=\left(E\right(x,t),W(x,t),$

$P(x,t),S(x,t),H(x,t),M(x,t))$. Define a function $\mathcal{H}\left(\theta \right)=\theta -1-ln\theta \ge 0$ for all $\theta >0$.

The following result suggests that when ${\mathcal{R}}_0\le 1$, the SARS-CoV-2 infection is predicted to be removed regardless of the initial conditions.

Theorem 2.

The steady state $S{S}_0$ is globally asymptotically stable (GAS) when ${\mathcal{R}}_0\le 1$.

Proof. 

Define a function ${\mathcal{G}}_0(x,t)$ as

$${\mathcal{G}}_0=\rho E_0\mathcal{H}\left(\frac{E}{E_0}\right)+\frac{\vartheta }{\vartheta +{\beta }_W}W+P+\frac{{\beta }_P}{\varrho }S+\frac{{\beta }_P\lambda }{\varrho \mu }H+\frac{\sigma }{\eta }M.$$

Clearly, ${\mathcal{G}}_0(E,W,P,S,H,M)>0$ for all $E,W,P,S,H,M>0$, and ${\mathcal{G}}_0(E_0,0,0,0,0,0)=0$. We calculate ${\displaystyle \frac{\partial {\mathcal{G}}_0}{\partial t}}$ along the solutions of Systems (2)–(7) as

$$\begin{array}{cc}\hfill \frac{\partial {\mathcal{G}}_0}{\partial t}& =\rho \left(1-\frac{E_0}{E}\right)\frac{\partial E}{\partial t}+\frac{\vartheta }{\vartheta +{\beta }_W}\frac{\partial W}{\partial t}+\frac{\partial P}{\partial t}+\frac{{\beta }_P}{\varrho }\frac{\partial S}{\partial t}+\frac{{\beta }_P\lambda }{\varrho \mu }\frac{\partial H}{\partial t}+\frac{\sigma }{\eta }\frac{\partial M}{\partial t}\hfill \\ \hfill & =\rho \left(1-\frac{E_0}{E}\right)\left[D_E\Delta E+\varphi -{\beta }_{E}E+\upsilon E\left(1-\frac{E}{E_{max}}\right)-\psi ES\right]\hfill \\ \hfill & +\frac{\vartheta }{\vartheta +{\beta }_W}\left[D_W\Delta W+\gamma \psi ES-(\vartheta +{\beta }_W)W\right]+D_P\Delta P+(1-\gamma )\psi ES+\vartheta W-{\beta }_{P}P-\sigma MP\hfill \\ \hfill & +\frac{{\beta }_P}{\varrho }\left[D_S\Delta S+\varrho P-{\beta }_{S}S-\lambda HS\right]+\frac{{\beta }_P\lambda }{\varrho \mu }\left[D_H\Delta H+\mu HS-{\beta }_{H}H\right]\hfill \\ \hfill & +\frac{\sigma }{\eta }\left[D_M\Delta M+\eta MP-{\beta }_{M}M\right].\hfill \end{array}$$

(22)

Collecting the terms of Equation (22), we obtain

$$\begin{array}{cc}\hfill \frac{\partial {\mathcal{G}}_0}{\partial t}& =\rho \left(1-\frac{E_0}{E}\right)\left[\varphi -{\beta }_{E}E+\upsilon E\left(1-\frac{E}{E_{max}}\right)\right]\hfill \\ \hfill & +\rho \psi E_{0}S-\frac{{\beta }_P{\beta }_S}{\varrho }S-\frac{{\beta }_P{\beta }_H\lambda }{\varrho \mu }H-\frac{{\beta }_M\sigma }{\eta }M+\rho D_E\left(1-\frac{E_0}{E}\right)\Delta E+\frac{\vartheta D_W}{\vartheta +{\beta }_W}\Delta W\hfill \\ \hfill & +D_P\Delta P+\frac{{\beta }_P{D}_S}{\varrho }\Delta S+\frac{{\beta }_P\lambda D_H}{\varrho \mu }\Delta H+\frac{\sigma D_M}{\eta }\Delta M.\hfill \end{array}$$

Form the steady-state conditions of $S{S}_0$, we obtain $\varphi ={\beta }_E{E}_0-\upsilon E_0\left(1-{\displaystyle \frac{E_0}{E_{max}}}\right)$ and thus,

$$\varphi -{\beta }_{E}E+\upsilon E\left(1-\frac{E}{E_{max}}\right)=(E_0-E)\left({\beta }_E-\upsilon +\frac{\upsilon E_0}{E_{max}}+\frac{\upsilon E}{E_{max}}\right).$$

Therefore, we deduce that

$$\begin{array}{cc}\hfill \frac{\partial {\mathcal{G}}_0}{\partial t}& =-\rho \left({\beta }_E-\upsilon +\frac{\upsilon E_0}{E_{max}}+\frac{\upsilon E}{E_{max}}\right)\frac{{(E-E_0)}^2}{E}+\rho \psi E_{0}S-\frac{{\beta }_P{\beta }_S}{\varrho }S-\frac{{\beta }_P{\beta }_H\lambda }{\varrho \mu }H-\frac{{\beta }_M\sigma }{\eta }M\hfill \\ \hfill & +\rho D_E\left(1-\frac{E_0}{E}\right)\Delta E+\frac{\vartheta D_W}{\vartheta +{\beta }_W}\Delta W+D_P\Delta P+\frac{{\beta }_P{D}_S}{\varrho }\Delta S+\frac{{\beta }_P\lambda D_H}{\varrho \mu }\Delta H+\frac{\sigma D_M}{\eta }\Delta M.\hfill \end{array}$$

The following inequality will be used for $i=0,1,2,3,4$

$$-\rho \left({\beta }_E-\upsilon +\frac{\upsilon E_i}{E_{max}}+\frac{\upsilon E}{E_{max}}\right)\frac{{(E-E_i)}^2}{E}\le -\rho \left({\beta }_E-\upsilon +\frac{\upsilon E_i}{E_{max}}\right)\frac{{(E-E_i)}^2}{E}.$$

(23)

Using Inequality (23) in case of $i=0$, we obtain

$$\begin{array}{cc}\hfill \frac{\partial {\mathcal{G}}_0}{\partial t}& \le -\rho \left({\beta }_E-\upsilon +\frac{\upsilon E_0}{E_{max}}\right)\frac{{(E-E_0)}^2}{E}+\frac{{\beta }_P{\beta }_S}{\varrho }\left(\frac{\varrho \psi E_0}{{\beta }_P{\beta }_S}\rho -1\right)S-\frac{{\beta }_P{\beta }_H\lambda }{\varrho \mu }H-\frac{{\beta }_M\sigma }{\eta }M\hfill \\ \hfill & +\rho D_E\left(1-\frac{E_0}{E}\right)\Delta E+\frac{\vartheta D_W}{\vartheta +{\beta }_W}\Delta W+D_P\Delta P+\frac{{\beta }_P{D}_S}{\varrho }\Delta S+\frac{{\beta }_P\lambda D_H}{\varrho \mu }\Delta H+\frac{\sigma D_M}{\eta }\Delta M.\hfill \end{array}$$

Using the definition of ${\mathcal{R}}_0$ we obtain

$$\begin{array}{cc}\hfill \frac{\partial {\mathcal{G}}_0}{\partial t}& \le -\rho \left({\beta }_E-\upsilon +\frac{\upsilon E_0}{E_{max}}\right)\frac{{(E-E_0)}^2}{E}+\frac{{\beta }_P{\beta }_S}{\varrho }\left({\mathcal{R}}_0-1\right)S-\frac{{\beta }_P{\beta }_H\lambda }{\varrho \mu }H-\frac{{\beta }_M\sigma }{\eta }M\hfill \\ \hfill & +\rho D_E\left(1-\frac{E_0}{E}\right)\Delta E+\frac{\vartheta D_W}{\vartheta +{\beta }_W}\Delta W+D_P\Delta P+\frac{{\beta }_P{D}_S}{\varrho }\Delta S+\frac{{\beta }_P\lambda D_H}{\varrho \mu }\Delta H+\frac{\sigma D_M}{\eta }\Delta M.\hfill \end{array}$$

Consequently, we calculate ${\displaystyle \frac{d{\tilde{\mathcal{G}}}_0}{dt}}$ as follows:

$$\begin{array}{cc}\hfill \frac{d{\tilde{\mathcal{G}}}_0}{dt}& ={\int }_{\varpi }\frac{\partial {\mathcal{G}}_0}{\partial t}dx\hfill \\ \hfill & \le -\rho \left({\beta }_E-\upsilon +\frac{\upsilon E_0}{E_{max}}\right){\int }_{\varpi }\frac{{(E-E_0)}^2}{E}dx+\frac{{\beta }_P{\beta }_S}{\varrho }\left({\mathcal{R}}_0-1\right){\int }_{\varpi }Sdx-\frac{{\beta }_P{\beta }_H\lambda }{\varrho \mu }{\int }_{\varpi }Hdx\hfill \\ \hfill & -\frac{{\beta }_M\sigma }{\eta }{\int }_{\varpi }Mdx+\rho D_E{\int }_{\varpi }\left(1-\frac{E_0}{E}\right)\Delta Edx+\frac{\vartheta D_W}{\vartheta +{\beta }_W}{\int }_{\varpi }\Delta Wdx+D_P{\int }_{\varpi }\Delta Pdx\hfill \\ \hfill & +\frac{{\beta }_P{D}_S}{\varrho }{\int }_{\varpi }\Delta Sdx+\frac{{\beta }_P{D}_H\lambda }{\varrho \mu }{\int }_{\varpi }\Delta Hdx+\frac{\sigma D_M}{\eta }{\int }_{\varpi }\Delta Mdx.\hfill \end{array}$$

Using Equality (21) we obtain

$$\begin{array}{cc}\hfill \frac{d{\tilde{\mathcal{G}}}_0}{dt}& \le -\rho \left({\beta }_E-\upsilon +\frac{\upsilon E_0}{E_{max}}\right){\int }_{\varpi }\frac{{(E-E_0)}^2}{E}\mathrm{d}x+\frac{{\beta }_P{\beta }_S}{\varrho }\left({\mathcal{R}}_0-1\right){\int }_{\varpi }Sdx-\frac{{\beta }_P{\beta }_H\lambda }{\varrho \mu }{\int }_{\varpi }Hdx\hfill \\ \hfill & -\frac{{\beta }_M\sigma }{\eta }{\int }_{\varpi }Mdx-\rho D_E{E}_0{\int }_{\varpi }\frac{{\parallel ▿E\parallel }^2}{E^2}dx.\hfill \end{array}$$

At the steady state $S{S}_0$, we have $\varphi ={\beta }_E{E}_0-\upsilon E_0\left(1-{\displaystyle \frac{E_0}{E_{max}}}\right)$, which implies that ${\beta }_E-\upsilon +{\displaystyle \frac{\upsilon E_0}{E_{max}}}>0$. It follows that ${\displaystyle \frac{d{\tilde{\mathcal{G}}}_0}{dt}}\le 0$ when ${\mathcal{R}}_0\le 1$. In addition, ${\displaystyle \frac{d{\tilde{\mathcal{G}}}_0}{dt}}=0$ when $E=E_0$, $S=0$, $H=0$ and $M=0$. The solutions of Systems (2)–(7) converge to ${\tilde{\Xi }}_0$. For any elements in ${\tilde{\Xi }}_0$ we have $E=E_0$ and $S=H=M=0$, and then ${\displaystyle \frac{\partial S}{\partial t}}=\Delta S=0$. From Equation (5), we obtain $0={\displaystyle \frac{\partial S}{\partial t}}=\varrho P$, which gives $P=0$ and thus, ${\displaystyle \frac{\partial P}{\partial t}}=\Delta P=0$. From Equation (4), we obtain $0={\displaystyle \frac{\partial P}{\partial t}}=\vartheta W$, which gives $W=0$. It follows that ${\tilde{\Xi }}_0=\left\{S{S}_0\right\}$. By LIP [39], we conclude that $S{S}_0$ is GAS when ${\mathcal{R}}_0\le 1$.

The next result establishes that when ${\mathcal{R}}_1^H\le 1<{\mathcal{R}}_0$, ${\mathcal{R}}_1^M\le 1$ and ${\beta }_E-\upsilon +{\displaystyle \frac{\upsilon E_1}{E_{max}}}>0$, the SARS-CoV-2 infection with inactive immune responses is always established regardless of the initial conditions.

Theorem 3.

The steady state $S{S}_1$ is GAS when ${\mathcal{R}}_1^H\le 1<{\mathcal{R}}_0$, ${\mathcal{R}}_1^M\le 1$ and ${\beta }_E-\upsilon +{\displaystyle \frac{\upsilon E_1}{E_{max}}}>0$.

Proof. 

Define a function ${\mathcal{G}}_1(x,t)$ as

$${\mathcal{G}}_1=\rho E_1\mathcal{H}\left(\frac{E}{E_1}\right)+\frac{\vartheta }{\vartheta +{\beta }_W}{W}_1\mathcal{H}\left(\frac{W}{W_1}\right)+P_1\mathcal{H}\left(\frac{P}{P_1}\right)+\frac{{\beta }_P}{\varrho }{S}_1\mathcal{H}\left(\frac{S}{S_1}\right)+\frac{{\beta }_P\lambda }{\varrho \mu }H+\frac{\sigma }{\eta }M.$$

Clearly, ${\mathcal{G}}_1(E,W,P,S,H,M)>0$ for all $E,W,P,S,H,M>0$, and ${\mathcal{G}}_1(E_1,W_1,P_1,S_1,0,0)=0$. Hence, ${\displaystyle \frac{\partial {\mathcal{G}}_1}{\partial t}}$ is given by

$$\begin{array}{cc}\hfill \frac{\partial {\mathcal{G}}_1}{\partial t}& =\rho \left(1-\frac{E_1}{E}\right)\left[D_E\Delta E+\varphi -{\beta }_{E}E+\upsilon E\left(1-\frac{E}{E_{max}}\right)-\psi ES\right]\hfill \\ \hfill & +\frac{\vartheta }{\vartheta +{\beta }_W}\left(1-\frac{W_1}{W}\right)\left[D_W\Delta W+\gamma \psi ES-(\vartheta +{\beta }_W)W\right]\hfill \\ \hfill & +\left(1-\frac{P_1}{P}\right)\left[D_P\Delta P+(1-\gamma )\psi ES+\vartheta W-{\beta }_{P}P-\sigma MP\right]\hfill \\ \hfill & +\frac{{\beta }_P}{\varrho }\left(1-\frac{S_1}{S}\right)\left[D_S\Delta S+\varrho P-{\beta }_{S}S-\lambda HS\right]+\frac{{\beta }_P\lambda }{\varrho \mu }\left[D_H\Delta H+\mu HS-{\beta }_{H}H\right]\hfill \\ \hfill & +\frac{\sigma }{\eta }\left[D_M\Delta M+\eta MP-{\beta }_{M}M\right].\hfill \end{array}$$

(24)

Collecting the terms of Equation (24), we obtain

$$\begin{array}{cc}\hfill \frac{\partial {\mathcal{G}}_1}{\partial t}& =\rho \left(1-\frac{E_1}{E}\right)\left[\varphi -{\beta }_{E}E+\upsilon E\left(1-\frac{E}{E_{max}}\right)\right]+\rho \psi E_{1}S-\frac{\vartheta \gamma \psi }{\vartheta +{\beta }_W}\frac{W_{1}ES}{W}+\vartheta W_1\hfill \\ \hfill & -(1-\gamma )\psi \frac{P_{1}ES}{P}-\vartheta \frac{P_{1}W}{P}+{\beta }_P{P}_1+\sigma P_{1}M-\frac{{\beta }_P{\beta }_S}{\varrho }S-{\beta }_P\frac{S_{1}P}{S}+\frac{{\beta }_P{\beta }_S}{\varrho }{S}_1+\frac{{\beta }_P\lambda }{\varrho }{S}_{1}H\hfill \\ \hfill & -\frac{{\beta }_P{\beta }_H\lambda }{\varrho \mu }H-\frac{{\beta }_M\sigma }{\eta }M+\rho D_E\left(1-\frac{E_1}{E}\right)\Delta E+\frac{\vartheta D_W}{\vartheta +{\beta }_W}\left(1-\frac{W_1}{W}\right)\Delta W\hfill \\ \hfill & +D_P\left(1-\frac{P_1}{P}\right)\Delta P+\frac{{\beta }_P{D}_S}{\varrho }\left(1-\frac{S_1}{S}\right)\Delta S+\frac{{\beta }_P\lambda D_H}{\varrho \mu }\Delta H+\frac{\sigma D_M}{\eta }\Delta M.\hfill \end{array}$$

Using the steady-state conditions at $S{S}_1$

$$\begin{array}{cc}\hfill & \varphi ={\beta }_E{E}_1-\upsilon E_1\left(1-\frac{E_1}{E_{max}}\right)+\psi E_1{S}_1,\hfill \\ \hfill & \vartheta W_1=\frac{\vartheta \gamma }{\vartheta +{\beta }_W}\psi E_1{S}_1,\hfill \\ \hfill & {\beta }_P{P}_1=\rho \psi E_1{S}_1=\frac{{\beta }_P{\beta }_S}{\varrho }{S}_1,\hfill \end{array}$$

we obtain

$$\varphi -{\beta }_{E}E+\upsilon E\left(1-\frac{E}{E_{max}}\right)=(E_1-E)\left({\beta }_E-\upsilon +\frac{\upsilon E_1}{E_{max}}+\frac{\upsilon E}{E_{max}}\right)+\psi E_1{S}_1.$$

Further, we obtain

$$\begin{array}{cc}\hfill \frac{\partial {\mathcal{G}}_1}{\partial t}& =-\rho \left({\beta }_E-\upsilon +\frac{\upsilon E_1}{E_{max}}+\frac{\upsilon E}{E_{max}}\right)\frac{{(E-E_1)}^2}{E}+\rho \left(1-\frac{E_1}{E}\right)\psi E_1{S}_1+2\rho \psi E_1{S}_1\hfill \\ \hfill & +\rho \psi E_{1}S-\frac{{\beta }_P{\beta }_S}{\varrho }S-\frac{\vartheta \gamma }{\vartheta +{\beta }_W}\psi E_1{S}_1\frac{W_{1}ES}{W{E}_1{S}_1}+\frac{\vartheta \gamma }{\vartheta +{\beta }_W}\psi E_1{S}_1-(1-\gamma )\psi E_1{S}_1\frac{P_{1}ES}{P{E}_1{S}_1}\hfill \\ \hfill & -\frac{\vartheta \gamma }{\vartheta +{\beta }_W}\psi E_1{S}_1\frac{P_{1}W}{P{W}_1}-\rho \psi E_1{S}_1\frac{S_{1}P}{S{P}_1}+\frac{{\beta }_P\lambda }{\varrho }\left(S_1-\frac{{\beta }_H}{\mu }\right)H+\sigma \left(P_1-\frac{{\beta }_M}{\eta }\right)M\hfill \\ \hfill & +\rho D_E\left(1-\frac{E_1}{E}\right)\Delta E+\frac{\vartheta D_W}{\vartheta +{\beta }_W}\left(1-\frac{W_1}{W}\right)\Delta W+D_P\left(1-\frac{P_1}{P}\right)\Delta P\hfill \\ \hfill & +\frac{{\beta }_P{D}_S}{\varrho }\left(1-\frac{S_1}{S}\right)\Delta S+\frac{{\beta }_P\lambda D_H}{\varrho \mu }\Delta H+\frac{\sigma D_M}{\eta }\Delta M.\hfill \end{array}$$

Using Inequality (23) in case of $i=1$, we obtain

$$\begin{array}{cc}\hfill \frac{\partial {\mathcal{G}}_1}{\partial t}& \le -\rho \left({\beta }_E-\upsilon +\frac{\upsilon E_1}{E_{max}}\right)\frac{{(E-E_1)}^2}{E}+3\rho \psi E_1{S}_1-\rho \psi E_1{S}_1\frac{E_1}{E}+\left(\rho \psi E_1-\frac{{\beta }_P{\beta }_S}{\varrho }\right)S\hfill \\ \hfill & -\frac{\vartheta \gamma }{\vartheta +{\beta }_W}\psi E_1{S}_1\frac{W_{1}ES}{W{E}_1{S}_1}+\frac{\vartheta \gamma }{\vartheta +{\beta }_W}\psi E_1{S}_1-(1-\gamma )\psi E_1{S}_1\frac{P_{1}ES}{P{E}_1{S}_1}\hfill \\ \hfill & -\frac{\vartheta \gamma }{\vartheta +{\beta }_W}\psi E_1{S}_1\frac{P_{1}W}{P{W}_1}-\rho \psi E_1{S}_1\frac{S_{1}P}{S{P}_1}+\frac{{\beta }_P\lambda }{\varrho }\left(S_1-\frac{{\beta }_H}{\mu }\right)H+\sigma \left(P_1-\frac{{\beta }_M}{\eta }\right)M\hfill \\ \hfill & +\rho D_E\left(1-\frac{E_1}{E}\right)\Delta E+\frac{\vartheta D_W}{\vartheta +{\beta }_W}\left(1-\frac{W_1}{W}\right)\Delta W+D_P\left(1-\frac{P_1}{P}\right)\Delta P\hfill \\ \hfill & +\frac{{\beta }_P{D}_S}{\varrho }\left(1-\frac{S_1}{S}\right)\Delta S+\frac{{\beta }_P\lambda D_H}{\varrho \mu }\Delta H+\frac{\sigma D_M}{\eta }\Delta M.\hfill \end{array}$$

We have $\rho \psi E_1-{\displaystyle \frac{{\beta }_P{\beta }_S}{\varrho }}=0$, then,

$$\begin{array}{cc}\hfill \frac{\partial {\mathcal{G}}_1}{\partial t}& \le -\rho \left({\beta }_E-\upsilon +\frac{\upsilon E_1}{E_{max}}\right)\frac{{(E-E_1)}^2}{E}+\frac{\vartheta \gamma }{\vartheta +{\beta }_W}\psi E_1{S}_1\left(4-\frac{E_1}{E}-\frac{W_{1}ES}{W{E}_1{S}_1}-\frac{P_{1}W}{P{W}_1}-\frac{S_{1}P}{S{P}_1}\right)\hfill \\ \hfill & +(1-\gamma )\psi E_1{S}_1\left(3-\frac{E_1}{E}-\frac{P_{1}ES}{P{E}_1{S}_1}-\frac{S_{1}P}{S{P}_1}\right)+\frac{{\beta }_P\lambda }{\varrho }\left(S_1-\frac{{\beta }_H}{\mu }\right)H+\sigma \left(P_1-\frac{{\beta }_M}{\eta }\right)M\hfill \\ \hfill & +\rho D_E\left(1-\frac{E_1}{E}\right)\Delta E+\frac{\vartheta D_W}{\vartheta +{\beta }_W}\left(1-\frac{W_1}{W}\right)\Delta W+D_P\left(1-\frac{P_1}{P}\right)\Delta P\hfill \\ \hfill & +\frac{{\beta }_P{D}_S}{\varrho }\left(1-\frac{S_1}{S}\right)\Delta S+\frac{{\beta }_P\lambda D_H}{\varrho \mu }\Delta H+\frac{\sigma D_M}{\eta }\Delta M.\hfill \end{array}$$

Hence, $\frac{d{\tilde{\mathcal{G}}}_1}{dt}$ is calculated as

$$\begin{array}{cc}\hfill \frac{d{\tilde{\mathcal{G}}}_1}{dt}& ={\int }_{\varpi }\frac{\partial {\mathcal{G}}_1}{\partial t}dx\hfill \\ \hfill & \le -\rho \left({\beta }_E-\upsilon +\frac{\upsilon E_1}{E_{max}}\right){\int }_{\varpi }\frac{{(E-E_1)}^2}{E}dx\hfill \\ \hfill & +\frac{\vartheta \gamma }{\vartheta +{\beta }_W}\psi E_1{S}_1{\int }_{\varpi }\left(4-\frac{E_1}{E}-\frac{W_{1}ES}{W{E}_1{S}_1}-\frac{P_{1}W}{P{W}_1}-\frac{S_{1}P}{S{P}_1}\right)dx\hfill \\ \hfill & +(1-\gamma )\psi E_1{S}_1{\int }_{\varpi }\left(3-\frac{E_1}{E}-\frac{P_{1}ES}{P{E}_1{S}_1}-\frac{S_{1}P}{S{P}_1}\right)dx+\frac{{\beta }_P\lambda }{\varrho }\left(S_1-\frac{{\beta }_H}{\mu }\right){\int }_{\varpi }Hdx\hfill \\ \hfill & +\sigma \left(P_1-\frac{{\beta }_M}{\eta }\right){\int }_{\varpi }Mdx-\rho D_E{E}_1{\int }_{\varpi }\frac{{\parallel ▿E\parallel }^2}{E^2}dx-\frac{\vartheta D_W{W}_1}{\vartheta +{\beta }_W}{\int }_{\varpi }\frac{{\parallel ▿W\parallel }^2}{W^2}dx\hfill \\ \hfill & -D_P{P}_1{\int }_{\varpi }\frac{{\parallel ▿P\parallel }^2}{P^2}dx-\frac{{\beta }_P{D}_S{S}_1}{\varrho }{\int }_{\varpi }\frac{{\parallel ▿S\parallel }^2}{S^2}dx.\hfill \end{array}$$

Since ${\beta }_E-\upsilon +{\displaystyle \frac{\upsilon E_1}{E_{max}}}>0$, then, we obtain

$$\begin{array}{cc}\hfill {\beta }_E-\upsilon +\frac{\upsilon E_1}{E_{max}}& ={\beta }_E-\upsilon +\frac{\upsilon {\beta }_P{\beta }_S}{\varrho \psi E_{max}\rho }>0\hfill \\ \hfill & ⟹{\beta }_E-\upsilon +\frac{{\beta }_H\psi }{\mu }+\frac{2\upsilon {\beta }_P{\beta }_S}{\varrho \psi E_{max}\rho }>0\hfill \\ \hfill & ⟹\frac{2\upsilon {\beta }_P{\beta }_S}{\varrho \psi E_{max}\rho }-\left(\upsilon -{\beta }_E-\frac{{\beta }_H\psi }{\mu }\right)>0.\hfill \end{array}$$

We obtain

$$\begin{array}{cc}\hfill {\mathcal{R}}_1^H\le 1& ⟺\frac{E_{max}}{2\upsilon }\left[\upsilon -{\beta }_E-\frac{{\beta }_H\psi }{\mu }+\sqrt{{\left(\upsilon -{\beta }_E-\frac{{\beta }_H\psi }{\mu }\right)}^2+\frac{4\upsilon \varphi }{E_{max}}}\right]\le \frac{{\beta }_P{\beta }_S}{\varrho \psi \rho }\hfill \\ \hfill & ⟺\sqrt{{\left(\upsilon -{\beta }_E-\frac{{\beta }_H\psi }{\mu }\right)}^2+\frac{4\upsilon \varphi }{E_{max}}}<\frac{2\upsilon {\beta }_P{\beta }_S}{\varrho \psi E_{max}\rho }-\left(\upsilon -{\beta }_E-\frac{{\beta }_H\psi }{\mu }\right)\hfill \\ \hfill & ⟺\frac{4\upsilon \varphi }{E_{max}}<\frac{4{\upsilon }^2{\beta }_P^2{\beta }_S^2}{{\varrho }^2{\psi }^2{E}_{max}^2{\rho }^2}-\frac{4\upsilon {\beta }_P{\beta }_S}{\varrho \psi E_{max}\rho }\left(\upsilon -{\beta }_E-\frac{{\beta }_H\psi }{\mu }\right)\hfill \\ \hfill & ⟺\upsilon \varphi <\frac{{\upsilon }^2{\beta }_P^2{\beta }_S^2}{{\varrho }^2{\psi }^2{E}_{max}{\rho }^2}-\frac{{\upsilon }^2{\beta }_P{\beta }_S}{\varrho \psi \rho }+\frac{\upsilon {\beta }_E{\beta }_P{\beta }_S}{\varrho \psi \rho }+\frac{\upsilon {\beta }_P{\beta }_S{\beta }_H}{\varrho \mu \rho }\hfill \\ \hfill & ⟺\frac{\varphi \varrho \rho }{{\beta }_P{\beta }_S}+\frac{\upsilon }{\psi }-\left(\frac{{\beta }_E}{\psi }+\frac{\upsilon {\beta }_P{\beta }_S}{\varrho {\psi }^2{E}_{max}\rho }\right)<\frac{{\beta }_H}{\mu }\hfill \\ \hfill & ⟺S_1<\frac{{\beta }_H}{\mu }.\hfill \end{array}$$

in addition, since ${\beta }_E-\upsilon +{\displaystyle \frac{\upsilon E_1}{E_{max}}}>0$, then we obtain

$$\begin{array}{cc}\hfill {\beta }_E-\upsilon +\frac{\upsilon E_1}{E_{max}}& ={\beta }_E-\upsilon +\frac{\upsilon {\beta }_P{\beta }_S}{\varrho \psi E_{max}\rho }>0\hfill \\ \hfill & ⟹{\beta }_E-\upsilon +\frac{{\beta }_M\varrho \psi }{{\beta }_S\eta }+\frac{2\upsilon {\beta }_P{\beta }_S}{\varrho \psi E_{max}\rho }>0\hfill \\ \hfill & ⟹\frac{2\upsilon {\beta }_P{\beta }_S}{\varrho \psi E_{max}\rho }-\left(\upsilon -{\beta }_E-\frac{{\beta }_M\varrho \psi }{{\beta }_S\eta }\right)>0.\hfill \end{array}$$

We obtain

$$\begin{array}{cc}\hfill {\mathcal{R}}_1^M\le 1& ⟺\frac{E_{max}}{2\upsilon }\left[\upsilon -{\beta }_E-\frac{{\beta }_M\varrho \psi }{{\beta }_S\eta }+\sqrt{{\left(\upsilon -{\beta }_E-\frac{{\beta }_M\varrho \psi }{{\beta }_S\eta }\right)}^2+\frac{4\upsilon \varphi }{E_{max}}}\right]\le \frac{{\beta }_P{\beta }_S}{\varrho \psi \rho }\hfill \\ \hfill & ⟺\sqrt{{\left(\upsilon -{\beta }_E-\frac{{\beta }_M\varrho \psi }{{\beta }_S\eta }\right)}^2+\frac{4\upsilon \varphi }{E_{max}}}<\frac{2\upsilon {\beta }_P{\beta }_S}{\varrho \psi E_{max}\rho }-\left(\upsilon -{\beta }_E-\frac{{\beta }_M\varrho \psi }{{\beta }_S\eta }\right)\hfill \\ \hfill & ⟺\frac{4\upsilon \varphi }{E_{max}}<\frac{4{\upsilon }^2{\beta }_P^2{\beta }_S^2}{{\varrho }^2{\psi }^2{E}_{max}^2{\rho }^2}-\frac{4\upsilon {\beta }_P{\beta }_S}{\varrho \psi E_{max}\rho }\left(\upsilon -{\beta }_E-\frac{{\beta }_M\varrho \psi }{{\beta }_S\eta }\right)\hfill \\ \hfill & ⟺\upsilon \varphi <\frac{{\upsilon }^2{\beta }_P^2{\beta }_S^2}{{\varrho }^2{\psi }^2{E}_{max}{\rho }^2}-\frac{{\upsilon }^2{\beta }_P{\beta }_S}{\varrho \psi \rho }+\frac{\upsilon {\beta }_E{\beta }_P{\beta }_S}{\varrho \psi \rho }+\frac{\upsilon {\beta }_P{\beta }_M}{\eta \rho }\hfill \\ \hfill & ⟺\frac{{\beta }_S}{\varrho }\left[\frac{\varphi \varrho \upsilon }{{\beta }_S}+\frac{{\upsilon }^2{\beta }_P}{\psi \rho }-\left(\frac{\upsilon {\beta }_E{\beta }_P}{\psi \rho }+\frac{{\upsilon }^2{\beta }_P^2{\beta }_S}{\varrho {\psi }^2{E}_{max}{\rho }^2}\right)\right]<\frac{\upsilon {\beta }_P{\beta }_M}{\eta \rho }\hfill \\ \hfill & ⟺\frac{{\beta }_S}{\varrho }\left[\frac{\varphi \varrho \rho }{{\beta }_P{\beta }_S}+\frac{\upsilon }{\psi }-\left(\frac{{\beta }_E}{\psi }+\frac{\upsilon {\beta }_P{\beta }_S}{\varrho {\psi }^2{E}_{max}\rho }\right)\right]<\frac{{\beta }_M}{\eta }\hfill \\ \hfill & ⟺P_1<\frac{{\beta }_M}{\eta }.\hfill \end{array}$$

If ${\mathcal{R}}_1^H\le 1$, ${\mathcal{R}}_1^M\le 1$ and ${\beta }_E-\upsilon +{\displaystyle \frac{\upsilon E_1}{E_{max}}}>0$, then, using Inequalities (19) and (20) we get ${\displaystyle \frac{d{\tilde{\mathcal{G}}}_1}{dt}}\le 0$. Moreover, ${\displaystyle \frac{d{\tilde{\mathcal{G}}}_1}{dt}}=0$ when $E=E_1$, $W=W_1$, $P=P_1$, $S=S_1$, $H=0$ and $M=0$. Thus, ${\tilde{\Xi }}_1=\left\{S{S}_1\right\}$ and by LIP, $S{S}_1$ is GAS when ${\mathcal{R}}_1^H\le 1<{\mathcal{R}}_0$, ${\mathcal{R}}_1^M\le 1$ and ${\beta }_E-\upsilon +{\displaystyle \frac{\upsilon E_1}{E_{max}}}>0$. □

The next result illustrates that when ${\mathcal{R}}_1^H>1$, ${\mathcal{R}}_2^M\le 1$ and ${\beta }_E-\upsilon +{\displaystyle \frac{\upsilon E_2}{E_{max}}}>0$, the SARS-CoV-2 infection with only active antibody immunity is always established regardless of the initial conditions.

Theorem 4.

The steady state $S{S}_2$ is GAS when ${\mathcal{R}}_1^H>1$, ${\mathcal{R}}_2^M\le 1$ and ${\beta }_E-\upsilon +{\displaystyle \frac{\upsilon E_2}{E_{max}}}>0$.

Proof. 

Define a function ${\mathcal{G}}_2(x,t)$ as

$$\begin{array}{cc}\hfill {\mathcal{G}}_2& =\rho E_2\mathcal{H}\left(\frac{E}{E_2}\right)+\frac{\vartheta }{\vartheta +{\beta }_W}{W}_2\mathcal{H}\left(\frac{W}{W_2}\right)+P_2\mathcal{H}\left(\frac{P}{P_2}\right)+\frac{{\beta }_P}{\varrho }{S}_2\mathcal{H}\left(\frac{S}{S_2}\right)+\frac{{\beta }_P\lambda }{\varrho \mu }{H}_2\mathcal{H}\left(\frac{H}{H_2}\right)\hfill \\ \hfill & +\frac{\sigma }{\eta }M.\hfill \end{array}$$

Clearly, ${\mathcal{G}}_2(E,W,P,S,H,M)>0$ for all $E,W,P,S,H,M>0$, and ${\mathcal{G}}_2(E_2,W_2,P_2,S_2,H_2,0)=0$. We calculate ${\displaystyle \frac{\partial {\mathcal{G}}_2}{\partial t}}$ as

$$\begin{array}{cc}\hfill \frac{\partial {\mathcal{G}}_2}{\partial t}& =\rho \left(1-\frac{E_2}{E}\right)\left[D_E\Delta E+\varphi -{\beta }_{E}E+\upsilon E\left(1-\frac{E}{E_{max}}\right)-\psi ES\right]\hfill \\ \hfill & +\frac{\vartheta }{\vartheta +{\beta }_W}\left(1-\frac{W_2}{W}\right)\left[D_W\Delta W+\gamma \psi ES-(\vartheta +{\beta }_W)W\right]\hfill \\ \hfill & +\left(1-\frac{P_2}{P}\right)\left[D_P\Delta P+(1-\gamma )\psi ES+\vartheta W-{\beta }_{P}P-\sigma MP\right]\hfill \\ \hfill & +\frac{{\beta }_P}{\varrho }\left(1-\frac{S_2}{S}\right)\left[D_S\Delta S+\varrho P-{\beta }_{S}S-\lambda HS\right]+\frac{{\beta }_P\lambda }{\varrho \mu }\left(1-\frac{H_2}{H}\right)\left[D_H\Delta H+\mu HS-{\beta }_{H}H\right]\hfill \\ \hfill & +\frac{\sigma }{\eta }\left[D_M\Delta M+\eta MP-{\beta }_{M}M\right].\hfill \end{array}$$

(25)

Collecting the terms of (25), we obtain

$$\begin{array}{cc}\hfill \frac{\partial {\mathcal{G}}_2}{\partial t}& =\rho \left(1-\frac{E_2}{E}\right)\left[\varphi -{\beta }_{E}E+\upsilon E\left(1-\frac{E}{E_{max}}\right)\right]+\rho \psi E_{2}S\hfill \\ \hfill & -\frac{\vartheta \gamma \psi }{\vartheta +{\beta }_W}\frac{W_{2}ES}{W}+\vartheta W_2-(1-\gamma )\psi \frac{P_{2}ES}{P}-\vartheta \frac{P_{2}W}{P}+{\beta }_P{P}_2+\sigma P_{2}M\hfill \\ \hfill & -\frac{{\beta }_P{\beta }_S}{\varrho }S-{\beta }_P\frac{S_{2}P}{S}+\frac{{\beta }_P{\beta }_S}{\varrho }{S}_2+\frac{{\beta }_P\lambda }{\varrho }{S}_{2}H-\frac{{\beta }_P{\beta }_H\lambda }{\varrho \mu }H-\frac{{\beta }_P\lambda }{\varrho }{H}_{2}S\hfill \\ \hfill & +\frac{{\beta }_P{\beta }_H\lambda }{\varrho \mu }{H}_2-\frac{{\beta }_M\sigma }{\eta }M+\rho D_E\left(1-\frac{E_2}{E}\right)\Delta E+\frac{\vartheta D_W}{\vartheta +{\beta }_W}\left(1-\frac{W_2}{W}\right)\Delta W\hfill \\ \hfill & +D_P\left(1-\frac{P_2}{P}\right)\Delta P+\frac{{\beta }_P{D}_S}{\varrho }\left(1-\frac{S_2}{S}\right)\Delta S+\frac{{\beta }_P\lambda D_H}{\varrho \mu }\left(1-\frac{H_2}{H}\right)\Delta H+\frac{\sigma D_M}{\eta }\Delta M.\hfill \end{array}$$

At the steady state $S{S}_2$ we have

$$\begin{array}{cc}\hfill & \varphi ={\beta }_E{E}_2-\upsilon E_2\left(1-\frac{E_2}{E_{max}}\right)+\psi E_2{S}_2,\hfill \\ \hfill & \vartheta W_2=\frac{\vartheta \gamma }{\vartheta +{\beta }_W}\psi E_2{S}_2,\hfill \\ \hfill & {\beta }_P{P}_2=\rho \psi E_2{S}_2=\frac{{\beta }_P{\beta }_S}{\varrho }{S}_2+\frac{{\beta }_P\lambda }{\varrho }{H}_2{S}_2,\hfill \\ \hfill & S_2=\frac{{\beta }_H}{\mu },\hfill \end{array}$$

and we obtain

$$\varphi -{\beta }_{E}E+\upsilon E\left(1-\frac{E}{E_{max}}\right)=(E_2-E)\left({\beta }_E-\upsilon +\frac{\upsilon E_2}{E_{max}}+\frac{\upsilon E}{E_{max}}\right)+\psi E_2{S}_2.$$

By using the above conditions, we obtain

$$\begin{array}{cc}\hfill \frac{\partial {\mathcal{G}}_2}{\partial t}& =-\rho \left({\beta }_E-\upsilon +\frac{\upsilon E_2}{E_{max}}+\frac{\upsilon E}{E_{max}}\right)\frac{{(E-E_2)}^2}{E}+\rho \left(1-\frac{E_2}{E}\right)\psi E_2{S}_2+2\rho \psi E_2{S}_2+\rho \psi E_{2}S\hfill \\ \hfill & -\frac{{\beta }_P{\beta }_S}{\varrho }S-\frac{{\beta }_P\lambda }{\varrho }{H}_{2}S-\frac{\vartheta \gamma }{\vartheta +{\beta }_W}\psi E_2{S}_2\frac{W_{2}ES}{W{E}_2{S}_2}+\frac{\vartheta \gamma }{\vartheta +{\beta }_W}\psi E_2{S}_2-(1-\gamma )\psi E_2{S}_2\frac{P_{2}ES}{P{E}_2{S}_2}\hfill \\ \hfill & -\frac{\vartheta \gamma }{\vartheta +{\beta }_W}\psi E_2{S}_2\frac{P_{2}W}{P{W}_2}-\rho \psi E_2{S}_2\frac{S_{2}P}{S{P}_2}+\sigma \left(P_2-\frac{{\beta }_M}{\eta }\right)M+\rho D_E\left(1-\frac{E_2}{E}\right)\Delta E\hfill \\ \hfill & +\frac{\vartheta D_W}{\vartheta +{\beta }_W}\left(1-\frac{W_2}{W}\right)\Delta W+D_P\left(1-\frac{P_2}{P}\right)\Delta P+\frac{{\beta }_P{D}_S}{\varrho }\left(1-\frac{S_2}{S}\right)\Delta S\hfill \\ \hfill & +\frac{{\beta }_P\lambda D_H}{\varrho \mu }\left(1-\frac{H_2}{H}\right)\Delta H+\frac{\sigma D_M}{\eta }\Delta M.\hfill \end{array}$$

Using Inequality (23) in the case of $i=2$, we obtain

$$\begin{array}{cc}\hfill \frac{\partial {\mathcal{G}}_2}{\partial t}& \le -\rho \left({\beta }_E-\upsilon +\frac{\upsilon E_2}{E_{max}}\right)\frac{{(E-E_2)}^2}{E}+3\rho \psi E_2{S}_2-\rho \psi E_2{S}_2\frac{E_2}{E}\hfill \\ \hfill & +\left(\rho \psi E_2-\frac{{\beta }_P{\beta }_S}{\varrho }-\frac{{\beta }_P\lambda }{\varrho }{H}_2\right)S-\frac{\vartheta \gamma }{\vartheta +{\beta }_W}\psi E_2{S}_2\frac{W_{2}ES}{W{E}_2{S}_2}+\frac{\vartheta \gamma }{\vartheta +{\beta }_W}\psi E_2{S}_2\hfill \\ \hfill & -(1-\gamma )\psi E_2{S}_2\frac{P_{2}ES}{P{E}_2{S}_2}-\frac{\vartheta \gamma }{\vartheta +{\beta }_W}\psi E_2{S}_2\frac{P_{2}W}{P{W}_2}-\rho \psi E_2{S}_2\frac{S_{2}P}{S{P}_2}+\sigma \left(P_2-\frac{{\beta }_M}{\eta }\right)M\hfill \\ \hfill & +\rho D_E\left(1-\frac{E_2}{E}\right)\Delta E+\frac{\vartheta D_W}{\vartheta +{\beta }_W}\left(1-\frac{W_2}{W}\right)\Delta W+D_P\left(1-\frac{P_2}{P}\right)\Delta P\hfill \\ \hfill & +\frac{{\beta }_P{D}_S}{\varrho }\left(1-\frac{S_2}{S}\right)\Delta S+\frac{{\beta }_P\lambda D_H}{\varrho \mu }\left(1-\frac{H_2}{H}\right)\Delta H+\frac{\sigma D_M}{\eta }\Delta M.\hfill \end{array}$$

The conditions of $S{S}_2$ imply that

$$\rho \psi E_2-\frac{{\beta }_P{\beta }_S}{\varrho }-\frac{{\beta }_P\lambda }{\varrho }{H}_2=0.$$

We obtain

$$\begin{array}{cc}\hfill \frac{\partial {\mathcal{G}}_2}{\partial t}& \le -\rho \left({\beta }_E-\upsilon +\frac{\upsilon E_2}{E_{max}}\right)\frac{{(E-E_2)}^2}{E}+\frac{\vartheta \gamma }{\vartheta +{\beta }_W}\psi E_2{S}_2\left(4-\frac{E_2}{E}-\frac{W_{2}ES}{W{E}_2{S}_2}-\frac{P_{2}W}{P{W}_2}-\frac{S_{2}P}{S{P}_2}\right)\hfill \\ \hfill & +(1-\gamma )\psi E_2{S}_2\left(3-\frac{E_2}{E}-\frac{P_{2}ES}{P{E}_2{S}_2}-\frac{S_{2}P}{S{P}_2}\right)+\sigma \left(P_2-\frac{{\beta }_M}{\eta }\right)M+\rho D_E\left(1-\frac{E_2}{E}\right)\Delta E\hfill \\ \hfill & +\frac{\vartheta D_W}{\vartheta +{\beta }_W}\left(1-\frac{W_2}{W}\right)\Delta W+D_P\left(1-\frac{P_2}{P}\right)\Delta P+\frac{{\beta }_P{D}_S}{\varrho }\left(1-\frac{S_2}{S}\right)\Delta S\hfill \\ \hfill & +\frac{{\beta }_P\lambda D_H}{\varrho \mu }\left(1-\frac{H_2}{H}\right)\Delta H+\frac{\sigma D_M}{\eta }\Delta M.\hfill \end{array}$$

Using the definition of ${\mathcal{R}}_2^M$ we obtain

$$\begin{array}{cc}\hfill \frac{\partial {\mathcal{G}}_2}{\partial t}& \le -\rho \left({\beta }_E-\upsilon +\frac{\upsilon E_2}{E_{max}}\right)\frac{{(E-E_2)}^2}{E}+\frac{\vartheta \gamma }{\vartheta +{\beta }_W}\psi E_2{S}_2\left(4-\frac{E_2}{E}-\frac{W_{2}ES}{W{E}_2{S}_2}-\frac{P_{2}W}{P{W}_2}-\frac{S_{2}P}{S{P}_2}\right)\hfill \\ \hfill & +(1-\gamma )\psi E_2{S}_2\left(3-\frac{E_2}{E}-\frac{P_{2}ES}{P{E}_2{S}_2}-\frac{S_{2}P}{S{P}_2}\right)+\frac{{\beta }_M\sigma }{\eta }\left({\mathcal{R}}_2^M-1\right)M+\rho D_E\left(1-\frac{E_2}{E}\right)\Delta E\hfill \\ \hfill & +\frac{\vartheta D_W}{\vartheta +{\beta }_W}\left(1-\frac{W_2}{W}\right)\Delta W+D_P\left(1-\frac{P_2}{P}\right)\Delta P+\frac{{\beta }_P{D}_S}{\varrho }\left(1-\frac{S_2}{S}\right)\Delta S\hfill \\ \hfill & +\frac{{\beta }_P\lambda D_H}{\varrho \mu }\left(1-\frac{H_2}{H}\right)\Delta H+\frac{\sigma D_M}{\eta }\Delta M.\hfill \end{array}$$

Then, $\frac{d{\tilde{\mathcal{G}}}_2}{dt}$ is calculated as

$$\begin{array}{cc}\hfill \frac{d{\tilde{\mathcal{G}}}_2}{dt}& ={\int }_{\varpi }\frac{\partial {\mathcal{G}}_2}{\partial t}dx\hfill \\ \hfill & \le -\rho \left({\beta }_E-\upsilon +\frac{\upsilon E_2}{E_{max}}\right){\int }_{\varpi }\frac{{(E-E_2)}^2}{E}dx\hfill \\ \hfill & +\frac{\vartheta \gamma }{\vartheta +{\beta }_W}\psi E_2{S}_2{\int }_{\varpi }\left(4-\frac{E_2}{E}-\frac{W_{2}ES}{W{E}_2{S}_2}-\frac{P_{2}W}{P{W}_2}-\frac{S_{2}P}{S{P}_2}\right)dx\hfill \\ \hfill & +(1-\gamma )\psi E_2{S}_2{\int }_{\varpi }\left(3-\frac{E_2}{E}-\frac{P_{2}ES}{P{E}_2{S}_2}-\frac{S_{2}P}{S{P}_2}\right)dx+\frac{{\beta }_M\sigma }{\eta }\left({\mathcal{R}}_2^M-1\right){\int }_{\varpi }Mdx\hfill \\ \hfill & -\rho D_E{E}_2{\int }_{\varpi }\frac{{\parallel ▿E\parallel }^2}{E^2}dx-\frac{\vartheta D_W{W}_2}{\vartheta +{\beta }_W}{\int }_{\varpi }\frac{{\parallel ▿W\parallel }^2}{W^2}dx-D_P{P}_2{\int }_{\varpi }\frac{{\parallel ▿P\parallel }^2}{P^2}dx\hfill \\ \hfill & -\frac{{\beta }_P{D}_S{S}_2}{\varrho }{\int }_{\varpi }\frac{{\parallel ▿S\parallel }^2}{S^2}dx-\frac{{\beta }_P\lambda D_H{H}_2}{\varrho \mu }{\int }_{\varpi }\frac{{\parallel ▿H\parallel }^2}{H^2}dx.\hfill \end{array}$$

If ${\mathcal{R}}_1^H>1$, ${\mathcal{R}}_2^M\le 1$ and ${\beta }_E-\upsilon +{\displaystyle \frac{\upsilon E_2}{E_{max}}}>0$, then, using Inequalities (19) and (20) we obtain ${\displaystyle \frac{d{\tilde{\mathcal{G}}}_2}{dt}}\le 0$. Moreover, ${\displaystyle \frac{d{\tilde{\mathcal{G}}}_2}{dt}}=0$ when, $E=E_2$, $W=W_2$, $P=P_2$, $S=S_2$ and $M=0$. The solutions of Systems (2)–(7) tend toward ${\tilde{\Xi }}_2$. For each element in ${\tilde{\Xi }}_2$ we have $S=S_2$, then, ${\displaystyle \frac{\partial S}{\partial t}}=\Delta S=0$ and from Equation (5), we have

$$0={\displaystyle \frac{\partial S}{\partial t}}=\varrho P_2-{\beta }_S{S}_2-\lambda H{S}_2⟹H=H_2.$$

It follows that ${\tilde{\Xi }}_2=\left\{S{S}_2\right\}$. By LIP, $S{S}_2$ is GAS when ${\mathcal{R}}_1^H>1$, ${\mathcal{R}}_2^M\le 1$ and ${\beta }_E-\upsilon +{\displaystyle \frac{\upsilon E_2}{E_{max}}}>0$. □

The next theorem shows that when ${\mathcal{R}}_1^M>1$, ${\mathcal{R}}_1^H\le {\mathcal{R}}_2^M$ and ${\beta }_E-\upsilon +{\displaystyle \frac{\upsilon E_3}{E_{max}}}>0$, the SARS-CoV-2 infection with only active CTL immunity is always established regardless of the initial conditions.

Theorem 5.

The steady state $S{S}_3$ is GAS when ${\mathcal{R}}_1^M>1$, ${\mathcal{R}}_1^H\le {\mathcal{R}}_2^M$ and ${\beta }_E-\upsilon +{\displaystyle \frac{\upsilon E_3}{E_{max}}}>0$.

Proof. 

Define ${\mathcal{G}}_3(x,t)$ as

$$\begin{array}{cc}\hfill {\mathcal{G}}_3& =\rho E_3\mathcal{H}\left(\frac{E}{E_3}\right)+\frac{\vartheta }{\vartheta +{\beta }_W}{W}_3\mathcal{H}\left(\frac{W}{W_3}\right)+P_3\mathcal{H}\left(\frac{P}{P_3}\right)+\frac{({\beta }_P+\sigma M_3)}{\varrho }{S}_3\mathcal{H}\left(\frac{S}{S_3}\right)\hfill \\ \hfill & +\frac{({\beta }_P+\sigma M_3)\lambda }{\varrho \mu }H+\frac{\sigma }{\eta }{M}_3\mathcal{H}\left(\frac{M}{M_3}\right).\hfill \end{array}$$

We have ${\mathcal{G}}_3(E,W,P,S,H,M)>0$ for all $E,W,P,S,H,M>0$, and ${\mathcal{G}}_3(E_3,W_3,P_3,S_3,0,M_3)=0$. Then, we calculate ${\displaystyle \frac{\partial {\mathcal{G}}_3}{\partial t}}$ as

$$\begin{array}{cc}\hfill \frac{\partial {\mathcal{G}}_3}{\partial t}& =\rho \left(1-\frac{E_3}{E}\right)\left[D_E\Delta E+\varphi -{\beta }_{E}E+\upsilon E\left(1-\frac{E}{E_{max}}\right)-\psi ES\right]\hfill \\ \hfill & +\frac{\vartheta }{\vartheta +{\beta }_W}\left(1-\frac{W_3}{W}\right)\left[D_W\Delta W+\gamma \psi ES-(\vartheta +{\beta }_W)W\right]\hfill \\ \hfill & +\left(1-\frac{P_3}{P}\right)\left[D_P\Delta P+(1-\gamma )\psi ES+\vartheta W-{\beta }_{P}P-\sigma MP\right]\hfill \\ \hfill & +\frac{({\beta }_P+\sigma M_3)}{\varrho }\left(1-\frac{S_3}{S}\right)\left[D_S\Delta S+\varrho P-{\beta }_{S}S-\lambda HS\right]\hfill \\ \hfill & +\frac{({\beta }_P+\sigma M_3)\lambda }{\varrho \mu }\left[D_H\Delta H+\mu HS-{\beta }_{H}H\right]+\frac{\sigma }{\eta }\left(1-\frac{M_3}{M}\right)\left[D_M\Delta M+\eta MP-{\beta }_{M}M\right].\hfill \end{array}$$

(26)

Collecting the terms of Equation (26), we obtain

$$\begin{array}{cc}\hfill \frac{\partial {\mathcal{G}}_3}{\partial t}& =\rho \left(1-\frac{E_3}{E}\right)\left[\varphi -{\beta }_{E}E+\upsilon E\left(1-\frac{E}{E_{max}}\right)\right]+\rho \psi E_{3}S-\frac{\vartheta \gamma \psi }{\vartheta +{\beta }_W}\frac{W_{3}ES}{W}\hfill \\ \hfill & +\vartheta W_3-(1-\gamma )\psi \frac{P_{3}ES}{P}-\vartheta \frac{P_{3}W}{P}+{\beta }_P{P}_3+\sigma P_{3}M-\frac{({\beta }_P+\sigma M_3){\beta }_S}{\varrho }S\hfill \\ \hfill & -({\beta }_P+\sigma M_3)\frac{S_{3}P}{S}+\frac{({\beta }_P+\sigma M_3){\beta }_S}{\varrho }{S}_3+\frac{({\beta }_P+\sigma M_3)\lambda }{\varrho }{S}_{3}H-\frac{({\beta }_P+\sigma M_3){\beta }_H\lambda }{\varrho \mu }H\hfill \\ \hfill & -\frac{{\beta }_M\sigma }{\eta }M+\frac{{\beta }_M\sigma }{\eta }{M}_3+\rho D_E\left(1-\frac{E_3}{E}\right)\Delta E+\frac{\vartheta D_W}{\vartheta +{\beta }_W}\left(1-\frac{W_3}{W}\right)\Delta W+D_P\left(1-\frac{P_3}{P}\right)\Delta P\hfill \\ \hfill & +\frac{({\beta }_P+\sigma M_3)D_S}{\varrho }\left(1-\frac{S_3}{S}\right)\Delta S+\frac{({\beta }_P+\sigma M_3)\lambda D_H}{\varrho \mu }\Delta H+\frac{\sigma D_M}{\eta }\left(1-\frac{M_3}{M}\right)\Delta M.\hfill \end{array}$$

The components of $S{S}_3$ satisfy

$$\begin{array}{cc}\hfill & \varphi ={\beta }_E{E}_3-\upsilon E_3\left(1-\frac{E_3}{E_{max}}\right)+\psi E_3{S}_3,\hfill \\ \hfill & \vartheta W_3=\frac{\vartheta \gamma }{\vartheta +{\beta }_W}\psi E_3{S}_3,\hfill \\ \hfill & ({\beta }_P+\sigma M_3)P_3=\rho \psi E_3{S}_3=\frac{({\beta }_P+\sigma M_3){\beta }_S}{\varrho }{S}_3,\hfill \\ \hfill & P_3=\frac{{\beta }_M}{\sigma },\varrho P_3={\beta }_S{S}_3,\hfill \end{array}$$

and we obtain

$$\varphi -{\beta }_{E}E+\upsilon E\left(1-\frac{E}{E_{max}}\right)=(E_3-E)\left({\beta }_E-\upsilon +\frac{\upsilon E_3}{E_{max}}+\frac{\upsilon E}{E_{max}}\right)+\psi E_3{S}_3.$$

By using the above conditions we obtain

$$\begin{array}{cc}\hfill \frac{\partial {\mathcal{G}}_3}{\partial t}& =-\rho \left({\beta }_E-\upsilon +\frac{\upsilon E_3}{E_{max}}+\frac{\upsilon E}{E_{max}}\right)\frac{{(E-E_3)}^2}{E}+\rho \left(1-\frac{E_3}{E}\right)\psi E_3{S}_3+2\rho \psi E_3{S}_3\hfill \\ \hfill & \rho \psi E_{3}S-\frac{({\beta }_P+\sigma M_3){\beta }_S}{\varrho }S-\frac{\vartheta \gamma }{\vartheta +{\beta }_W}\psi E_3{S}_3\frac{W_{3}ES}{W{E}_3{S}_3}+\frac{\vartheta \gamma }{\vartheta +{\beta }_W}\psi E_3{S}_3\hfill \\ \hfill & -(1-\gamma )\psi E_3{S}_3\frac{P_{3}ES}{P{E}_3{S}_3}-\frac{\vartheta \gamma }{\vartheta +{\beta }_W}\psi E_3{S}_3\frac{P_{3}W}{P{W}_3}-\rho \psi E_3{S}_3\frac{S_{3}P}{S{P}_3}\hfill \\ \hfill & +\frac{({\beta }_P+\sigma M_3)\lambda }{\varrho }\left(S_3-\frac{{\beta }_H}{\mu }\right)H+\rho D_E\left(1-\frac{E_3}{E}\right)\Delta E+\frac{\vartheta D_W}{\vartheta +{\beta }_W}\left(1-\frac{W_3}{W}\right)\Delta W\hfill \\ \hfill & +D_P\left(1-\frac{P_3}{P}\right)\Delta P+\frac{({\beta }_P+\sigma M_3)D_S}{\varrho }\left(1-\frac{S_3}{S}\right)\Delta S+\frac{({\beta }_P+\sigma M_3)\lambda D_H}{\varrho \mu }\Delta H\hfill \\ \hfill & +\frac{\sigma D_M}{\eta }\left(1-\frac{M_3}{M}\right)\Delta M.\hfill \end{array}$$

Using Inequality (23) in the case of $i=3$, we obtain

$$\begin{array}{cc}\hfill \frac{\partial {\mathcal{G}}_3}{\partial t}& \le -\rho \left({\beta }_E-\upsilon +\frac{\upsilon E_3}{E_{max}}\right)\frac{{(E-E_3)}^2}{E}+3\rho \psi E_3{S}_3-\rho \psi E_3{S}_3\frac{E_3}{E}\hfill \\ \hfill & +\left(\rho \psi E_3-\frac{({\beta }_P+\sigma M_3){\beta }_S}{\varrho }\right)S-\frac{\vartheta \gamma }{\vartheta +{\beta }_W}\psi E_3{S}_3\frac{W_{3}ES}{W{E}_3{S}_3}+\frac{\vartheta \gamma }{\vartheta +{\beta }_W}\psi E_3{S}_3\hfill \\ \hfill & -(1-\gamma )\psi E_3{S}_3\frac{P_{3}ES}{P{E}_3{S}_3}-\frac{\vartheta \gamma }{\vartheta +{\beta }_W}\psi E_3{S}_3\frac{P_{3}W}{P{W}_3}-\rho \psi E_3{S}_3\frac{S_{3}P}{S{P}_3}\hfill \\ \hfill & +\frac{({\beta }_P+\sigma M_3)\lambda }{\varrho }\left(S_3-\frac{{\beta }_H}{\mu }\right)H+\rho D_E\left(1-\frac{E_3}{E}\right)\Delta E+\frac{\vartheta D_W}{\vartheta +{\beta }_W}\left(1-\frac{W_3}{W}\right)\Delta W\hfill \\ \hfill & +D_P\left(1-\frac{P_3}{P}\right)\Delta P+\frac{({\beta }_P+\sigma M_3)D_S}{\varrho }\left(1-\frac{S_3}{S}\right)\Delta S+\frac{({\beta }_P+\sigma M_3)\lambda D_H}{\varrho \mu }\Delta H\hfill \\ \hfill & +\frac{\sigma D_M}{\eta }\left(1-\frac{M_3}{M}\right)\Delta M.\hfill \end{array}$$

We have $\rho \psi E_3-{\displaystyle \frac{({\beta }_P+\sigma M_3){\beta }_S}{\varrho }}=0$. Then, we obtain

$$\begin{array}{cc}\hfill \frac{\partial {\mathcal{G}}_3}{\partial t}& \le -\rho \left({\beta }_E-\upsilon +\frac{\upsilon E_3}{E_{max}}\right)\frac{{(E-E_3)}^2}{E}+\frac{\vartheta \gamma }{\vartheta +{\beta }_W}\psi E_3{S}_3\left(4-\frac{E_3}{E}-\frac{W_{3}ES}{W{E}_3{S}_3}-\frac{P_{3}W}{P{W}_3}-\frac{S_{3}P}{S{P}_3}\right)\hfill \\ \hfill & +(1-\gamma )\psi E_3{S}_3\left(3-\frac{E_3}{E}-\frac{P_{3}ES}{P{E}_3{S}_3}-\frac{S_{3}P}{S{P}_3}\right)+\frac{({\beta }_P+\sigma M_3)\lambda }{\varrho }\left(S_3-\frac{{\beta }_H}{\mu }\right)H\hfill \\ \hfill & +\rho D_E\left(1-\frac{E_3}{E}\right)\Delta E+\frac{\vartheta D_W}{\vartheta +{\beta }_W}\left(1-\frac{W_3}{W}\right)\Delta W+D_P\left(1-\frac{P_3}{P}\right)\Delta P\hfill \\ \hfill & +\frac{({\beta }_P+\sigma M_3)D_S}{\varrho }\left(1-\frac{S_3}{S}\right)\Delta S+\frac{({\beta }_P+\sigma M_3)\lambda D_H}{\varrho \mu }\Delta H+\frac{\sigma D_M}{\eta }\left(1-\frac{M_3}{M}\right)\Delta M.\hfill \end{array}$$

Using the definition of ${\mathcal{R}}_1^H$ and ${\mathcal{R}}_2^M$ we obtain

$$\begin{array}{cc}\hfill \frac{\partial {\mathcal{G}}_3}{\partial t}& \le -\rho \left({\beta }_E-\upsilon +\frac{\upsilon E_3}{E_{max}}\right)\frac{{(E-E_3)}^2}{E}+\frac{\vartheta \gamma }{\vartheta +{\beta }_W}\psi E_3{S}_3\left(4-\frac{E_3}{E}-\frac{W_{3}ES}{W{E}_3{S}_3}-\frac{P_{3}W}{P{W}_3}-\frac{S_{3}P}{S{P}_3}\right)\hfill \\ \hfill & +(1-\gamma )\psi E_3{S}_3\left(3-\frac{E_3}{E}-\frac{P_{3}ES}{P{E}_3{S}_3}-\frac{S_{3}P}{S{P}_3}\right)+\frac{({\beta }_P+\sigma M_3){\beta }_H\lambda }{\varrho \mu }\left(\frac{{\mathcal{R}}_1^H}{{\mathcal{R}}_2^M}-1\right)H\hfill \\ \hfill & +\rho D_E\left(1-\frac{E_3}{E}\right)\Delta E+\frac{\vartheta D_W}{\vartheta +{\beta }_W}\left(1-\frac{W_3}{W}\right)\Delta W+D_P\left(1-\frac{P_3}{P}\right)\Delta P\hfill \\ \hfill & +\frac{({\beta }_P+\sigma M_3)D_S}{\varrho }\left(1-\frac{S_3}{S}\right)\Delta S+\frac{({\beta }_P+\sigma M_3)\lambda D_H}{\varrho \mu }\Delta H+\frac{\sigma D_M}{\eta }\left(1-\frac{M_3}{M}\right)\Delta M.\hfill \end{array}$$

(27)

Therefore, $\frac{d{\tilde{\mathcal{G}}}_3}{dt}$ is given by

$$\begin{array}{cc}\hfill \frac{d{\tilde{\mathcal{G}}}_3}{dt}& ={\int }_{\varpi }\frac{\partial {\mathcal{G}}_3}{\partial t}dx\hfill \\ \hfill & \le -\rho \left({\beta }_E-\upsilon +\frac{\upsilon E_3}{E_{max}}\right){\int }_{\varpi }\frac{{(E-E_3)}^2}{E}dx\hfill \\ \hfill & +\frac{\vartheta \gamma }{\vartheta +{\beta }_W}\psi E_3{S}_3{\int }_{\varpi }\left(4-\frac{E_3}{E}-\frac{W_{3}ES}{W{E}_3{S}_3}-\frac{P_{3}W}{P{W}_3}-\frac{S_{3}P}{S{P}_3}\right)dx\hfill \\ \hfill & +(1-\gamma )\psi E_3{S}_3{\int }_{\varpi }\left(3-\frac{E_3}{E}-\frac{P_{3}ES}{P{E}_3{S}_3}-\frac{S_{3}P}{S{P}_3}\right)dx\hfill \\ \hfill & +\frac{({\beta }_P+\sigma M_3){\beta }_H\lambda }{\varrho \mu }\left(\frac{{\mathcal{R}}_1^H}{{\mathcal{R}}_2^M}-1\right){\int }_{\varpi }Hdx-\rho D_E{E}_3{\int }_{\varpi }\frac{{\parallel ▿E\parallel }^2}{E^2}dx\hfill \\ \hfill & -\frac{\vartheta D_W{W}_3}{\vartheta +{\beta }_W}{\int }_{\varpi }\frac{{\parallel ▿W\parallel }^2}{W^2}dx-D_P{P}_3{\int }_{\varpi }\frac{{\parallel ▿P\parallel }^2}{P^2}dx-\frac{({\beta }_P+\sigma M_3)D_S{S}_3}{\varrho }{\int }_{\varpi }\frac{{\parallel ▿S\parallel }^2}{S^2}dx\hfill \\ \hfill & -\frac{\sigma D_M{M}_3}{\eta }{\int }_{\varpi }\frac{{\parallel ▿M\parallel }^2}{M^2}dx.\hfill \end{array}$$

If ${\mathcal{R}}_1^M>1$, ${\displaystyle \frac{{\mathcal{R}}_1^H}{{\mathcal{R}}_2^M}}\le 1$ and ${\beta }_E-\upsilon +{\displaystyle \frac{\upsilon E_3}{E_{max}}}>0$ then, using Inequalities (19) and (20) we obtain ${\displaystyle \frac{d{\tilde{\mathcal{G}}}_3}{dt}}\le 0$. In addition, ${\displaystyle \frac{d{\tilde{\mathcal{G}}}_3}{dt}}=0$ when, $E=E_3$, $W=W_3$, $P=P_3$, $S=S_3$, $H=0$ and $M=M_3$. It follows that ${\tilde{\Xi }}_3=\left\{S{S}_3\right\}$. By LIP, $S{S}_3$ is GAS when ${\mathcal{R}}_1^M>1$, ${\displaystyle \frac{{\mathcal{R}}_1^H}{{\mathcal{R}}_2^M}}\le 1$ and ${\beta }_E-\upsilon +{\displaystyle \frac{\upsilon E_3}{E_{max}}}>0$. □

The next theorem demonstrates that, when ${\mathcal{R}}_1^H>{\mathcal{R}}_2^M>1$ and ${\beta }_E-\upsilon +{\displaystyle \frac{\upsilon E_4}{E_{max}}}>0$, the SARS-CoV-2 infection with both active antibody and CTL immune responses is always established regardless of the initial conditions.

Theorem 6.

The steady state $S{S}_4$ is GAS when ${\mathcal{R}}_1^H>{\mathcal{R}}_2^M>1$ and ${\beta }_E-\upsilon +{\displaystyle \frac{\upsilon E_4}{E_{max}}}>0$.

Proof. 

Define a function ${\mathcal{G}}_4(x,t)$ as

$$\begin{array}{cc}\hfill {\mathcal{G}}_4& =\rho E_4\mathcal{H}\left(\frac{E}{E_4}\right)+\frac{\vartheta }{\vartheta +{\beta }_W}{W}_4\mathcal{H}\left(\frac{W}{W_4}\right)+P_4\mathcal{H}\left(\frac{P}{P_4}\right)+\frac{({\beta }_P+\sigma M_4)}{\varrho }{S}_4\mathcal{H}\left(\frac{S}{S_4}\right)\hfill \\ \hfill & +\frac{({\beta }_P+\sigma M_4)\lambda }{\varrho \mu }{H}_4\mathcal{H}\left(\frac{H}{H_4}\right)+\frac{\sigma }{\eta }{M}_4\mathcal{H}\left(\frac{M}{M_4}\right).\hfill \end{array}$$

We have ${\mathcal{G}}_4(E,W,P,S,H,M)>0$ for all $E,W,P,S,H,M>0$, and ${\mathcal{G}}_4(E_4,W_4,P_4,S_4,H_4,M_4)=0$. Then, we calculate ${\displaystyle \frac{\partial {\mathcal{G}}_4}{\partial t}}$ as

$$\begin{array}{cc}\hfill \frac{\partial {\mathcal{G}}_4}{\partial t}& =\rho \left(1-\frac{E_4}{E}\right)\left[D_E\Delta E+\varphi -{\beta }_{E}E+\upsilon E\left(1-\frac{E}{E_{max}}\right)-\psi ES\right]\hfill \\ \hfill & +\frac{\vartheta }{\vartheta +{\beta }_W}\left(1-\frac{W_4}{W}\right)\left[{\beta }_W\Delta W+\gamma \psi ES-(\vartheta +{\beta }_W)W\right]\hfill \\ \hfill & +\left(1-\frac{P_4}{P}\right)\left[D_P\Delta P+(1-\gamma )\psi ES+\vartheta W-{\beta }_{P}P-\sigma MP\right]\hfill \\ \hfill & +\frac{({\beta }_P+\sigma M_4)}{\varrho }\left(1-\frac{S_4}{S}\right)\left[D_S\Delta S+\varrho P-{\beta }_{S}S-\lambda HS\right]\hfill \\ \hfill & +\frac{({\beta }_P+\sigma M_4)\lambda }{\varrho \mu }\left(1-\frac{H_4}{H}\right)\left[D_H\Delta H+\mu HS-{\beta }_{H}H\right]\hfill \\ \hfill & +\frac{\sigma }{\eta }\left(1-\frac{M_4}{M}\right)\left[D_M\Delta M+\eta MP-{\beta }_{M}M\right].\hfill \end{array}$$

(28)

Collecting the terms of Equation (28), we obtain

$$\begin{array}{cc}\hfill \frac{\partial {\mathcal{G}}_4}{\partial t}& =\rho \left(1-\frac{E_4}{E}\right)\left[\varphi -{\beta }_{E}E+\upsilon E\left(1-\frac{E}{E_{max}}\right)\right]+\rho \psi E_{4}S-\frac{\vartheta \gamma \psi }{\vartheta +{\beta }_W}\frac{W_{4}ES}{W}\hfill \\ \hfill & +\vartheta W_4-(1-\gamma )\psi \frac{P_{4}ES}{P}-\vartheta \frac{P_{4}W}{P}+{\beta }_P{P}_4+\sigma P_{4}M-\frac{({\beta }_P+\sigma M_4){\beta }_S}{\varrho }S-({\beta }_P+\sigma M_4)\frac{S_{4}P}{S}\hfill \\ \hfill & +\frac{({\beta }_P+\sigma M_4){\beta }_S}{\varrho }{S}_4+\frac{({\beta }_P+\sigma M_4)\lambda }{\varrho }{S}_{4}H-\frac{({\beta }_P+\sigma M_4){\beta }_H\lambda }{\varrho \mu }H-\frac{({\beta }_P+\sigma M_4)\lambda }{\varrho }{H}_{4}S\hfill \\ \hfill & +\frac{({\beta }_P+\sigma M_4){\beta }_H\lambda }{\varrho \mu }{H}_4-\frac{{\beta }_M\sigma }{\eta }M+\frac{{\beta }_M\sigma }{\eta }{M}_4+\rho D_E\left(1-\frac{E_4}{E}\right)\Delta E\hfill \\ \hfill & +\frac{\vartheta D_W}{\vartheta +{\beta }_W}\left(1-\frac{W_4}{W}\right)\Delta W+D_P\left(1-\frac{P_4}{P}\right)\Delta P+\frac{({\beta }_P+\sigma M_4)D_S}{\varrho }\left(1-\frac{S_4}{S}\right)\Delta S\hfill \\ \hfill & +\frac{({\beta }_P+\sigma M_4)\lambda D_H}{\varrho \mu }\left(1-\frac{H_4}{H}\right)\Delta H+\frac{\sigma D_M}{\eta }\left(1-\frac{M_4}{M}\right)\Delta M.\hfill \end{array}$$

The steady state $S{S}_4$ satisfies the following:

$$\begin{array}{cc}\hfill & \varphi ={\beta }_E{E}_4-\upsilon E_4\left(1-\frac{E_4}{E_{max}}\right)+\psi E_4{S}_4,\hfill \\ \hfill & \vartheta W_4=\frac{\vartheta \gamma }{\vartheta +{\beta }_W}\psi E_4{S}_4,\hfill \\ \hfill & ({\beta }_P+\sigma M_4)P_4=\rho \psi E_4{S}_4=\frac{({\beta }_P+\sigma M_4){\beta }_S}{\varrho }{S}_4+\frac{({\beta }_P+\sigma M_4)\lambda }{\varrho }{S}_4{H}_4,\hfill \\ \hfill & S_4=\frac{{\beta }_H}{\mu },P_4=\frac{{\beta }_M}{\eta },\hfill \end{array}$$

and we obtain

$$\varphi -{\beta }_{E}E+\upsilon E\left(1-\frac{E}{E_{max}}\right)=(E_4-E)\left({\beta }_E-\upsilon +\frac{\upsilon E_4}{E_{max}}+\frac{\upsilon E}{E_{max}}\right)+\psi E_4{S}_4.$$

By using the above conditions, we obtain

$$\begin{array}{cc}\hfill \frac{\partial {\mathcal{G}}_4}{\partial t}& =-\rho \left({\beta }_E-\upsilon +\frac{\upsilon E_4}{E_{max}}+\frac{\upsilon E}{E_{max}}\right)\frac{{(E-E_4)}^2}{E}+\rho \left(1-\frac{E_4}{E}\right)\psi E_4{S}_4+2\rho \psi E_4{S}_4\hfill \\ \hfill & +\rho \psi E_{4}S-\frac{({\beta }_P+\sigma M_4){\beta }_S}{\varrho }S-\frac{({\beta }_P+\sigma M_4)\lambda }{\varrho }{H}_{4}S\hfill \\ \hfill & -\frac{\vartheta \gamma }{\vartheta +{\beta }_W}\psi E_4{S}_4\frac{W_{4}ES}{W{E}_4{S}_4}+\frac{\vartheta \gamma }{\vartheta +{\beta }_W}\psi E_4{S}_4-(1-\gamma )\psi E_4{S}_4\frac{P_{4}ES}{P{E}_4{S}_4}-\frac{\vartheta \gamma }{\vartheta +{\beta }_W}\psi E_4{S}_4\frac{P_{4}W}{P{W}_4}\hfill \\ \hfill & -\rho \psi E_4{S}_4\frac{S_{4}P}{S{P}_4}+\rho D_E\left(1-\frac{E_4}{E}\right)\Delta E+\frac{\vartheta D_W}{\vartheta +{\beta }_W}\left(1-\frac{W_4}{W}\right)\Delta W+D_P\left(1-\frac{P_4}{P}\right)\Delta P\hfill \\ \hfill & +\frac{({\beta }_P+\sigma M_4)D_S}{\varrho }\left(1-\frac{S_4}{S}\right)\Delta S+\frac{({\beta }_P+\sigma M_4)\lambda D_H}{\varrho \mu }\left(1-\frac{H_4}{H}\right)\Delta H\hfill \\ \hfill & +\frac{\sigma D_M}{\eta }\left(1-\frac{M_4}{M}\right)\Delta M.\hfill \end{array}$$

Using Inequality (23) in the case of $i=4$, we obtain

$$\begin{array}{cc}\hfill \frac{\partial {\mathcal{G}}_4}{\partial t}& \le -\rho \left({\beta }_E-\upsilon +\frac{\upsilon E_4}{E_{max}}\right)\frac{{(E-E_4)}^2}{E}+3\rho \psi E_4{S}_4-\rho \psi E_4{S}_4\frac{E_4}{E}\hfill \\ \hfill & +\left(\rho \psi E_4-\frac{({\beta }_P+\sigma M_4){\beta }_S}{\varrho }-\frac{({\beta }_P+\sigma M_4)\lambda }{\varrho }{H}_4\right)S\hfill \\ \hfill & -\frac{\vartheta \gamma }{\vartheta +{\beta }_W}\psi E_4{S}_4\frac{W_{4}ES}{W{E}_4{S}_4}+\frac{\vartheta \gamma }{\vartheta +{\beta }_W}\psi E_4{S}_4-(1-\gamma )\psi E_4{S}_4\frac{P_{4}ES}{P{E}_4{S}_4}-\frac{\vartheta \gamma }{\vartheta +{\beta }_W}\psi E_4{S}_4\frac{P_{4}W}{P{W}_4}\hfill \\ \hfill & -\rho \psi E_4{S}_4\frac{S_{4}P}{S{P}_4}+\rho D_E\left(1-\frac{E_4}{E}\right)\Delta E+\frac{\vartheta D_W}{\vartheta +{\beta }_W}\left(1-\frac{W_4}{W}\right)\Delta W+D_P\left(1-\frac{P_4}{P}\right)\Delta P\hfill \\ \hfill & +\frac{({\beta }_P+\sigma M_4)D_S}{\varrho }\left(1-\frac{S_4}{S}\right)\Delta S+\frac{({\beta }_P+\sigma M_4)\lambda D_H}{\varrho \mu }\left(1-\frac{H_4}{H}\right)\Delta H\hfill \\ \hfill & +\frac{\sigma D_M}{\eta }\left(1-\frac{M_4}{M}\right)\Delta M.\hfill \end{array}$$

The steady-state conditions of $S{S}_4$ imply that

$$\rho \psi E_4-{\displaystyle \frac{({\beta }_P+\sigma M_4){\beta }_S}{\varrho }}-{\displaystyle \frac{({\beta }_P+\sigma M_4)\lambda }{\varrho }}{H}_4=0.$$

We obtain

$$\begin{array}{cc}\hfill \frac{\partial {\mathcal{G}}_4}{\partial t}& \le -\rho \left({\beta }_E-\upsilon +\frac{\upsilon E_4}{E_{max}}\right)\frac{{(E-E_4)}^2}{E}+\frac{\vartheta \gamma }{\vartheta +{\beta }_W}\psi E_4{S}_4\left(4-\frac{E_4}{E}-\frac{W_{4}ES}{W{E}_4{S}_4}-\frac{P_{4}W}{P{W}_4}-\frac{S_{4}P}{S{P}_4}\right)\hfill \\ \hfill & +(1-\gamma )\psi E_4{S}_4\left(3-\frac{E_4}{E}-\frac{P_{4}ES}{P{E}_4{S}_4}-\frac{S_{4}P}{S{P}_4}\right)+\rho D_E\left(1-\frac{E_4}{E}\right)\Delta E\hfill \\ \hfill & +\frac{\vartheta D_W}{\vartheta +{\beta }_W}\left(1-\frac{W_4}{W}\right)\Delta W+D_P\left(1-\frac{P_4}{P}\right)\Delta P+\frac{({\beta }_P+\sigma M_4)D_S}{\varrho }\left(1-\frac{S_4}{S}\right)\Delta S\hfill \\ \hfill & +\frac{({\beta }_P+\sigma M_4)\lambda D_H}{\varrho \mu }\left(1-\frac{H_4}{H}\right)\Delta H+\frac{\sigma D_M}{\eta }\left(1-\frac{M_4}{M}\right)\Delta M.\hfill \end{array}$$

Calculating $\frac{d{\tilde{\mathcal{G}}}_4}{dt}$ as

$$\begin{array}{cc}\hfill \frac{d{\tilde{\mathcal{G}}}_4}{dt}& ={\int }_{\varpi }\frac{\partial {\mathcal{G}}_4}{\partial t}dx\hfill \\ \hfill & \le -\rho \left({\beta }_E-\upsilon +\frac{\upsilon E_4}{E_{max}}\right){\int }_{\varpi }\frac{{(E-E_4)}^2}{E}dx\hfill \\ \hfill & +\frac{\vartheta \gamma }{\vartheta +{\beta }_W}\psi E_4{S}_4{\int }_{\varpi }\left(4-\frac{E_4}{E}-\frac{W_{4}ES}{W{E}_4{S}_4}-\frac{P_{4}W}{P{W}_4}-\frac{S_{4}P}{S{P}_4}\right)dx\hfill \\ \hfill & +(1-\gamma )\psi E_4{S}_4{\int }_{\varpi }\left(3-\frac{E_4}{E}-\frac{P_{4}ES}{P{E}_4{S}_4}-\frac{S_{4}P}{S{P}_4}\right)dx-\rho D_E{E}_4{\int }_{\varpi }\frac{{\parallel ▿E\parallel }^2}{E^2}dx\hfill \\ \hfill & -\frac{\vartheta D_W{W}_4}{\vartheta +{\beta }_W}{\int }_{\varpi }\frac{{\parallel ▿W\parallel }^2}{W^2}dx-D_P{P}_4{\int }_{\varpi }\frac{{\parallel ▿P\parallel }^2}{P^2}dx-\frac{({\beta }_P+\sigma M_4)D_S{S}_4}{\varrho }{\int }_{\varpi }\frac{{\parallel ▿S\parallel }^2}{S^2}dx\hfill \\ \hfill & -\frac{({\beta }_P+\sigma M_4)\lambda D_H{H}_4}{\varrho \mu }{\int }_{\varpi }\frac{{\parallel ▿H\parallel }^2}{H^2}dx-\frac{\sigma D_M{M}_4}{\eta }{\int }_{\varpi }\frac{{\parallel ▿M\parallel }^2}{M^2}dx.\hfill \end{array}$$

We see that ${\displaystyle \frac{d{\tilde{\mathcal{G}}}_4}{dt}}\le 0$ when ${\mathcal{R}}_1^H>{\mathcal{R}}_2^M>1$, ${\beta }_E-\upsilon +{\displaystyle \frac{\upsilon E_4}{E_{max}}}>0$ and using Inequalities (19) and (20). Moreover, ${\displaystyle \frac{d{\tilde{\mathcal{G}}}_4}{dt}}=0$ when, $E=E_4$, $W=W_4$, $P=P_4$, $S=S_4$, $H=H_4$ and $M=M_4$. It follows that ${\tilde{\Xi }}_4=\left\{S{S}_4\right\}$. By LIP, $S{S}_4$ is GAS when ${\mathcal{R}}_1^H>{\mathcal{R}}_2^M>1$ and ${\beta }_E-\upsilon +{\displaystyle \frac{\upsilon E_4}{E_{max}}}>0$. □

Based on the above findings, we summarize the existence and global stability conditions for all steady state points in

Table 1. Steady states and their global stability conditions for Models (2)–(7).

Steady State	Global Stability Conditions
$S{S}_0=(E_0,0,0,0,0,0)$	${\mathcal{R}}_0\le 1$
$S{S}_1=(E_1,W_1,P_1,S_1,0,0)$	${\mathcal{R}}_1^H\le 1<{\mathcal{R}}_0$, ${\mathcal{R}}_1^M\le 1$ and ${\beta }_E-\upsilon +\frac{\upsilon E_1}{E_{max}}>0$
$S{S}_2=(E_2,W_2,P_2,S_2,H_2,0)$	${\mathcal{R}}_1^H>1$, ${\mathcal{R}}_2^M\le 1$ and ${\beta }_E-\upsilon +\frac{\upsilon E_2}{E_{max}}>0$
$S{S}_3=(E_3,W_3,P_3,S_3,0,M_3)$	${\mathcal{R}}_1^M>1$, ${\mathcal{R}}_1^H\le {\mathcal{R}}_2^M$ and ${\beta }_E-\upsilon +\frac{\upsilon E_3}{E_{max}}>0$
$S{S}_4=(E_4,W_4,P_4,S_4,H_4,M_4)$	${\mathcal{R}}_1^H>{\mathcal{R}}_2^M>1$ and ${\beta }_E-\upsilon +\frac{\upsilon E_4}{E_{max}}>0$

.

6. Numerical Simulations

In this section, we execute numerical simulations to support the results of Theorems 2–6. The MATLB PDE solver (pdepe) is used to solve Systems (2)–(7). We consider the spatial domain as $\varpi =[0,2]$, with a step size of $0.02$. The time step size is chosen as $0.1$. Moreover, we consider the following initial conditions:

$$\begin{array}{cc}\hfill & E(x,0)=5\left[1+0.2{cos}^2\left(\pi x\right)\right],W(x,0)=0.005\left[1+0.8{cos}^2\left(\pi x\right)\right],\hfill \\ \hfill & P(x,0)=0.5\left[1+0.2{cos}^2\left(\pi x\right)\right],S(x,0)=0.05\left[1+0.3{cos}^2\left(\pi x\right)\right],\hfill \\ \hfill & H(x,0)=3\left[1+0.3{cos}^2\left(\pi x\right)\right],M(x,0)=0.01\left[1+0.1{cos}^2\left(\pi x\right)\right],x\in [0,2].\hfill \end{array}$$

(29)

In addition, we consider the homogeneous Neumann boundary conditions:

$$\frac{\partial E}{\partial \overrightarrow{n}}=\frac{\partial W}{\partial \overrightarrow{n}}=\frac{\partial P}{\partial \overrightarrow{n}}=\frac{\partial S}{\partial \overrightarrow{n}}=\frac{\partial H}{\partial \overrightarrow{n}}=\frac{\partial M}{\partial \overrightarrow{n}}=0,t>0,x=0,2.$$

(30)

The parameters $(\psi ,\mu ,\eta )$ of Models (2)–(7) are taken as free parameters. The other parameters are fixed, as shown in

Table 2. Model parameters.

Parameter	Value	Parameter	Value	Parameter	Value	Parameter	Value
$\varphi$	$0.1$	$\sigma$	$0.8$	${\beta }_W$	$0.001$	$D_W$	$0.01$
$\upsilon$	$0.01$	$\varrho$	$0.24$	${\beta }_P$	$0.1$	$D_P$	$0.01$
$E_{max}$	12	$\lambda$	$0.5$	${\beta }_S$	$4.36$	$D_S$	$0.01$
$\psi$	varied	$\mu$	varied	${\beta }_H$	$0.05$	$D_H$	$0.01$
$\gamma$	$0.5$	$\eta$	varied	${\beta }_M$	$0.1$	$D_M$	$0.01$
$\vartheta$	$4.08$	${\beta }_E$	$0.01$	$D_E$	$0.01$

. To illustrate the global stability of the five steady states of Models (2)–(7), we have the following cases:

Case 1 (Stability of $S{S}_0$): $(\psi ,\mu ,\eta )=(0.05,0.5,0.09)$. For these values, we obtain ${\mathcal{R}}_0=0.3015<1$. According to Theorem 2, $S{S}_0$ is GAS and the SARS-CoV-2 infection is predicted to be completely cleared from the body. From Figure 1, we can see that the numerical results agree with the results of Theorem 2. We observe that the concentration of healthy ECs is increased and converged to its normal value $E_0=10.9495$, whereas the concentrations of other compartments are reducing and tending to zero.

Case 2 (Stability of $S{S}_1$): $(\psi ,\mu ,\eta )=(0.5,0.5,0.09)$. Using these values, we compute ${\mathcal{R}}_1^H=0.5332<1$, ${\mathcal{R}}_1^M=0.8314<1<{\mathcal{R}}_0=3.0146$ and ${\beta }_E-\upsilon +{\displaystyle \frac{\upsilon E_1}{E_{max}}}=0.003>0$. It means that the conditions of Theorem 3 are valid and thus $S{S}_1$ is GAS. From Figure 2, we see that there is an agreement between the numerical simulations and the results of Theorem 3. Further, the states of the system converge to the steady state $S{S}_1=(3.6335,0.0109,0.8899,0.049,0,0)$. In this case, the SARS-CoV-2 exists in the body but without any response from the immune system.

Case 3 (Stability of $S{S}_2$): $(\psi ,\mu ,\eta )=(0.5,2.3,0.09)$. Consequently, ${\mathcal{R}}_1^H=1.7137>1,{\mathcal{R}}_2^M=0.6091<1$ and ${\beta }_E-\upsilon +{\displaystyle \frac{\upsilon E_2}{E_{max}}}=0.0052>0$. This shows that the conditions of Theorem 4 are fulfilled and thus $S{S}_2$ is GAS. The numerical solutions displayed in Figure 3 are consistent with the results of Theorem 4. Further, the states of the system converge to the steady state $S{S}_2=(6.2268,0.0083,0.6776,0.0218,6.2213,0)$. In this case, only the antibody immunity is activated.

Case 4 (Stability of $S{S}_3$): $(\psi ,\mu ,\eta )=(0.5,1,0.13)$. We compute ${\mathcal{R}}_1^M=1.1203>1,{\displaystyle \frac{{\mathcal{R}}_1^H}{{\mathcal{R}}_2^M}}=0.8469<1$ and ${\beta }_E-\upsilon +{\displaystyle \frac{\upsilon E_3}{E_{max}}}=0.0034>0$. According to Theorem 5, $S{S}_3$ is GAS and this is shown in Figure 4. We can see that the states of the system converge to the steady state $S{S}_3=(4.0730,0.0106,0.7713,0.0425,0,0.0150)$. For this case, the CTLs are activated, whereas the antibody immune response is unstimulated.

Case 5 (Stability of $S{S}_4$): $(\psi ,\mu ,\eta )=(0.5,1.6,0.13)$. Hence, we compute ${\displaystyle \frac{{\mathcal{R}}_1^H}{{\mathcal{R}}_2^M}}=1.355>1,{\mathcal{R}}_2^M=1.0243>1$ and ${\beta }_E-\upsilon +{\displaystyle \frac{\upsilon E_4}{E_{max}}}=0.0042>0$. According to Theorem 6, $S{S}_4$ is GAS and this is clarified numerically in Figure 5. The states of the system converge to the steady state $S{S}_4=(5.0433,0.0097,0.7691,0.0313,3.0927,0.0031)$. In this situation, both antibodies and CTLs are activated against the viral infection.

7. Discussion

SARS-CoV-2 infection represents a real concern worldwide. Therefore, the modeling and analysis of SARS-CoV-2 are needed to understand the dynamics of this virus within a host. In this paper, we develop a diffusive SARS-CoV-2 infection model with antibody and CTL immune responses. We study the dynamical behavior of the model. We establish that the model has five steady states and we prove the following:

• The healthy steady state $S{S}_0$ always exists and it is GAS when ${\mathcal{R}}_0\le 1$. In this case, the patient will be recoveredfrom COVID-19. From a control viewpoint, making ${\mathcal{R}}_0\le 1$ will be a good strategy. This can be obtained by reducing the parameters $\psi$ and $\varrho$. Let us consider the effect of two types of antiviral drugs, one for blocking the infection with drug efficacy ${ϵ}_1\in [0,1]$ and the other for blocking the production of SARS-CoV-2 with drug efficacy ${ϵ}_2\in [0,1]$ [8]. Modeling the two antiviral drugs will change parameters $\psi$ and $\varrho$ to $(1-{ϵ}_1)\psi$ and $(1-{ϵ}_2)\varrho$ [21]. Let us consider ${ϵ}_1={ϵ}_2=ϵ$ and the other parameters are fixed, then, ${\mathcal{R}}_0$ can be given as functions of $ϵ$ as follows:

  $$\begin{array}{cc}\hfill {\mathcal{R}}_0\left(ϵ\right)& =\frac{{(1-ϵ)}^2\varrho \psi E_{max}}{2\upsilon {\beta }_P{\beta }_S}\left({\displaystyle \frac{\vartheta \gamma }{\vartheta +{\beta }_W}}+1-\gamma \right)\left(\upsilon -{\beta }_E+\sqrt{{(\upsilon -{\beta }_E)}^2+\frac{4\upsilon \varphi }{E_{max}}}\right)\hfill \\ \hfill & ={(1-ϵ)}^2{\mathcal{R}}_0\left(0\right).\hfill \end{array}$$
  To make ${\mathcal{R}}_0\le 1$, the effectiveness $ϵ$ has to satisfy

  $${ϵ}^{min}\le ϵ\le 1,{ϵ}^{min}=max\left\{0,1-\frac{1}{\sqrt{{\mathcal{R}}_0\left(0\right)}}\right\},$$

  where ${ϵ}^{min}$ is the minimum drug efficacy required to eradicate SARS-CoV-2 from the body.
  We note that ${\mathcal{R}}_0$ does not depend on the immune response parameters $\sigma$, $\lambda$, $\mu$, and $\eta$. Therefore, both CTL and antibody immune responses can control the viral infection but they do not play the role of clearing the viruses.
• The infected steady state with inactive immune responses $S{S}_1$ exists when ${\mathcal{R}}_0>1$ and ${\beta }_E-\upsilon +{\displaystyle \frac{\upsilon E_1}{E_{max}}}>0$. Further, $S{S}_1$ is GAS when ${\mathcal{R}}_1^H\le 1<{\mathcal{R}}_0$, ${\mathcal{R}}_1^M\le 1$ and ${\beta }_E-\upsilon +{\displaystyle \frac{\upsilon E_1}{E_{max}}}>0$. This result suggests that starting from any disease stage, the COVID-19 patient will still have a SARS-CoV-2 infection but without immune responses.
• The infected steady state with only active antibody immunity $S{S}_2$ exists when ${\mathcal{R}}_1^H>1$. Moreover, $S{S}_2$ is GAS, when ${\mathcal{R}}_1^H>1$, ${\mathcal{R}}_2^M\le 1$ and ${\beta }_E-\upsilon +{\displaystyle \frac{\upsilon E_2}{E_{max}}}>0$. This result suggests that starting from any disease stage, the COVID-19 patient will still have a SARS-CoV-2 infection but with only an active antibody immune response.
• The infected steady state with only active CTL immunity $S{S}_3$ exists when ${\mathcal{R}}_1^M>1$, whereas it is GAS when ${\mathcal{R}}_1^M>1$, ${\mathcal{R}}_1^H\le {\mathcal{R}}_2^M$ and ${\beta }_E-\upsilon +{\displaystyle \frac{\upsilon E_3}{E_{max}}}>0$. This result suggests that starting from any disease stage, the COVID-19 patient will still have a SARS-CoV-2 infection but with only an active CTL immune response.
• The infected steady state with both active antibody and CTL immunities $S{S}_4$ exists when ${\mathcal{R}}_1^H>{\mathcal{R}}_2^M$ and ${\mathcal{R}}_2^M>1$. Further, $S{S}_4$ is GAS when ${\mathcal{R}}_1^H>{\mathcal{R}}_2^M>1$ and ${\beta }_E-\upsilon +{\displaystyle \frac{\upsilon E_4}{E_{max}}}>0$. This result suggests that starting from any disease stage, the COVID-19 patient will still have a SARS-CoV-2 infection despite both antibody and CTL immune responses being active.
  We performed the numerical simulations for the model and showed that both the numerical and theoretical results are consistent.

We presented some numerical results and showed that they agreed with the theoretical results. The main limitation of the present work is that we did not validate the model using real data from COVID-19 patients (such as the concentrations of SARS-CoV-2, CTLs, antibodies, etc.). In fact, the real data on SARS-CoV-2 infections are still very limited. Collecting such data from SARS-CoV-2-infected patients is not an easy task and needs further experimental studies.

8. Conclusions

Mathematical models can be helpful for understanding the complex behavior of viral infections and the reaction of the immune system. We noted that the great majority of works on within-host SARS-CoV-2 infection models are based on the assumption that the viruses and cells are homogeneously distributed in the human body. This is a poor approach because the diffusion of viruses and cells causes spatial variations within the body. In this paper, we constructed a diffusive SARS-CoV-2 infection model with antibody and CTL immune responses. The model was given by a system of PDEs that describes the interaction of six compartments: healthy ECs, latent infected ECs, active infected ECs, free SARS-CoV-2 particles, CTLs, and antibodies. We considered a logistic term for healthy ECs. The non-negativity and boundedness of the solutions of the model were proven. Further, we derived four threshold parameters that determine the existence and stability of the five steady states of the model. The global stability of all steady states of the model was investigated by constructing Lyapunov functions and applying LIP. We performed numerical simulations for the model and showed that both the numerical and theoretical results are consistent.

Our model suggest that both CTL and antibody immune responses can play a role in controlling SARS-CoV-2 infection but not in clearing the virus from the body.

We discussed the effect of the antiviral treatment of the SARS-CoV-2 dynamics. Our model can help pharmaceutical companies and biologists to develop effective antiviral drugs for COVID-19 patients that make the basic reproduction number ${\mathcal{R}}_0$ for patients less than or equal to one. This will lead to the clearance of SARS-CoV-2 from the patient’s body.