Delayed Modeling Approach to Forecast the Periodic Behavior of SARS-2

1. Introduction

In the field of biological sciences, the delayed processes takes place, not only at macro scale but also at micro scale. The computational framework for such problems can help in understanding the dynamics in a more cost effective manner (Yan, 2007; Fang et al., 2020; Sohail and Nutini, 2020).

Recently, the World Health Organization (WHO) has declared the novel corona virus (2019-nCoV) outbreak a Public Health Emergency of International Concern (PHEIC) It is named as severe acute respiratory syndrome corona virus 2 (SARS-CoV-2). The SARS-CoV-2 has been determined as the seventh member of the corona viruses infected humans (Zhu et al., 2020).

The antiviral immune response, also in the case of SARS-CoV2, behaves according to two characteristics: “lytic” and “non-lytic.” These characteristics emerge from the action of two fundamental elements of the immune system in the case of an anti-viral response: antibodies (non-lytic response) and cytotoxic T cells also called CTL—Cytotoxic Lymphocites (lytic response).

At the moment the epitopes that can be identified to analyze these responses, are not clear. There are strong homologies with the SARS-CoV virus (75.5%). Furthermore, there is a strong alteration of antigenicity compared to SARS-CoV2 even if important epitopes in the Spike protein have been identified through the analysis of the localization of the RDB sequences, which can therefore become specific targets for drugs and vaccines (Zheng and Song, 2020).

The non-lytic response occurs against adaptive humeral immunity through the production of antibodies by B cells whose action is to neutralize the virion through direct connection with the same; specific antibodies are important in the defense against viruses during the early stages of infection when the virus is still extra-cellular; neutralizing antibodies facing the virus bind to the capsid or viral pericapside proteins, preventing their adhesion to the cell surface and therefore entry into the cells. Opsonizing antibodies can potentate the elimination of viral particles by phagocytosis.

In the case of SARS-CoV2, virus-specific IgG reached 100% approximately 17–19 days after symptom onset, while the proportion of patients with positive virus-specific IgM reached a peak of 94.1% approximately 20–22 days after symptom onset while during the first 3 weeks after symptom onset, there was increase in virus-specific IgG and IgM antibody titters.

Three types of seroconversion are possible: synchronous seroconversion of IgG and IgM, IgM seroconversion earlier than that of IgG and IgM seroconversion later than that of IgG (Long et al., 2020); specific data on the production of IgG and IgM is crucial to allow the rapid identification of the infection (di Mauro Gabriella et al., 2020). The lytic response depends on the action of the CTL cells which, in order to take place efficiently, must result from a cooperation between CD4⁺ and CD8⁺ lymphocytes. CD8⁺ cells recognize endogenously synthesized viral antigens in association with MHC class I molecules on all cell types; CD4⁺ cells recognize viral antigens in association with MHC class II antigens and the complete differentiation of CD8⁺ cells requires the intervention of CD4⁺ T helper cells. The antiviral effect of CTL is due to the lysis of infected cells, the activation of endocellular enzymes (endonucleases) which cause the degradation of the viral genome and the secretion of cytokines with interferon activity.

CTL epitopes of SARS-CoV-2 have been predicted by several studies, which can be used effectively to understand the pathogenesis (Kumar et al., 2020).

The fundamental mechanism by which the advancement of the virus in an organism is contrasted is given by a balance between the two mechanisms (lytic and non lytic). During multiple stages of the infection, the initial coordinated action of the antibodies join the lytic activity of the CTL cells to eradicate the virus.

The current research is motivated by the possible lack of the “control” and “coordination” between the two types of immune response, such that, a strong CTL response induces a limited and insufficient antibody action (or the reverse); is interlinked with the strong decrease in viral load induced by one of the two types of response.

During this imbalanced process, the coordinative response of the CD4⁺ helper T cells also varies considerably: in SARS-CoV-2 infected patients, it has been reported that the analysis showed activation and reduction in CD4⁺ and CD8⁺ T cell counts (Li et al., 2020).

The purpose of this paper, therefore, is to quantify a possible mathematical model that investigates the possible non-coordination between lytic and non-lytic response and indicate a mathematical quantification of the best antiviral immune response in the case of SARS-CoV2.

2. Materials and the Methods

2.1. Basic Model of Virus Transmission

Consider the basic transmission model, with the coefficients, as listed in Table 1.

TABLE 1

Table 1. Parameters description.

Model

$\begin{array}{ll}\begin{array}{l}\frac{dU}{dt}=\rho -dU-\mu VU,\\ \frac{dW}{dt}=\mu VU-aW-pWC,\\ \frac{dV}{dt}=\kappa W-qVA-\varphi V,\\ \frac{dA}{dt}=gVA-hA,\\ \frac{dC}{dt}=cWC-bC.\end{array}& (1)\end{array}$

The description of compartments U(t), W(t), V(t), A(t), and C(t) are presented in Table 2. The basic concept of modeling was adapted from the work of Wodarz (2005).

TABLE 2

Table 2. Compartments and their description.

2.1.1. Positivity of Solution

For the model (1) transmission to be epidemiologically feasible, it is necessary to show that for all time every state variable is non-negative. Thus, the solutions of the model with non-negative initial value for all time t > 0 will remain non-negative.

Theorem 2.1. Suppose that the model (1) consists of all feasible solutions with non-negative initial result then it remains non-negative for all time t.

Proof: Let the model (1) satisfy the initial non-negative solution, i.e.,

U(0) ≥ 0, W(0) ≥ 0, V(0) ≥ 0, A(0) ≥ 0 and C(0) ≥ 0.

It can be deduced from the model (1)

$\begin{array}{ll}\frac{dU}{dt}=\rho -dU-\mu VU,& (2)\end{array}$

as the solution of variable U can be computed by following result,

$\begin{array}{ll}U=U(0)e^G+{\displaystyle {\int }_0^t}\pi e^{H(k)}d(k)\ge 0,& (3)\end{array}$

where $G=-{\int }_0^{t}G(k)d(k)$ and $H(k)=-{\int }_0^{t}G(l)d(l)$. This gives non-negativity of U(0) ≥ 0 for all t ≥ 0. The non-negativity of the rest of the variables in the model (1) is given below.

$\begin{array}{ll}\begin{array}{l}\frac{dW}{dt}=\mu VU-aW-pWC,\\ \frac{dV}{dt}=\kappa W-qVA-\varphi V,\\ \frac{dA}{dt}=gVA-hA,\\ \frac{dC}{dt}=cWC-bC.\end{array}& (4)\end{array}$

In term of matrix the above can be expressed as,

$\begin{array}{ll}\frac{dF(t)}{dt}=H(t)+MF(t)& (5)\end{array}$

where

$\begin{array}{ll}F(t)=\left(\begin{array}{l}W(t)\\ V(t)\\ A(t)\\ C(t)\end{array}\right),H(t)=\left(\begin{array}{l}a\\ 0\\ 0\\ 0\end{array}\right)& (6)\end{array}$

and M matrix.

$\begin{array}{ll}M=\left(\begin{array}{cccc}-a& \frac{\mu \rho }{d}& 0& 0\\ \kappa & -\varphi & 0& 0\\ 0& 0& -h& 0\\ 0& 0& 0& -b\end{array}\right).& (7)\end{array}$

The matrix M is a Matzler matrix so by the result presented in Smith (1996) the model is monotone. The fact that $R_{+}^4$ is invariant with respect to stream of model (1).

2.1.2. Qualitative Analysis

Stability analysis is an essential key to validate the models in the field of science and technology. The remarkable contributions by Tunc (2002); Tunç (2004, 2008, 2013) and other researchers in the cross references can not be denied, since their contribution help the researchers to deal with the highly nonlinear models.

Here we present the stability analysis of given mathematical model (1). The model (1) is locally asymptotically stable at uninfected and infected equilibrium points. For uninfected equilibrium, the model is locally stable, if the value of reproduction number ℝ₀ < 1, whereas for infected equilibrium the model is stable if the value of the basic reproduction number ℝ₀ > 1. Furthermore, we will investigate the model (1) is locally stable at uninfected and infected. The uninfected equilibrium point E⁰ = (U⁰, 0, 0, 0, 0) and infected equilibrium point E^* = (U^*, W^*, V^*, A^*, C^*) of model (1) are constructed by following theorems.

Theorem 2.2. The uninfected equilibrium point E⁰ of the model (1) is given by

$\begin{array}{l}{E}^0=(U^0,0,0,0,0),\end{array}$

where:

$\begin{array}{l}{U}^0=\frac{\rho }{d}.\end{array}$

Proof: By putting the equations of model (1) equal to zero, the dynamics is determined as follows:

$\begin{array}{r}\begin{array}{r}\rho -dU-\mu VU=0,\\ \mu VU-aW-pWC=0,\\ \kappa W-qVA-\varphi V=0,\\ gVA-hA=0,\\ cWC-bC=0,\end{array}\end{array}$

and after further algebraic manipulation, we got

$\begin{array}{ll}{U}^0=\frac{\rho }{d},& (8)\end{array}$

which completes the proof.

Theorem 2.3. The model (1) admits a unique infected equilibrium E^* = (U^*, W^*, V^*, A^*, C^*) if and only if ℝ₀ > 1.

Proof: The infected equilibrium point is given as:

$\begin{array}{l}{U}^{*}=\frac{g\rho }{dg+\mu h},\\ W^{*}=\frac{b}{c},\\ V^{*}=\frac{h}{g},\\ A^{*}=\frac{bg\kappa -ch\varphi }{chq},\\ C^{*}=\frac{-abdg-ab\mu h+\mu ch\rho }{bp(dg+\mu h)}.& (9)\end{array}$

2.1.3. Reproduction Number ℝ₀ and Stability Analysis

The basic reproduction number ℝ₀ is formulated by evaluating infection matrix F and transmission matrix V constructed from Jacobian matrix J of model (1). The reproduction number ℝ₀ is the spectral radius of matrix FV⁻¹. The Jacobian matrix for J₀ for uninfected equilibrium point is

$\begin{array}{ll}{J}_0=\left(\begin{array}{ccccc}-d& 0& -\frac{\mu \rho }{d}& 0& 0\\ 0& -a& \frac{\mu \rho }{d}& 0& 0\\ 0& \kappa & -\varphi & 0& 0\\ 0& 0& 0& -h& 0\\ 0& 0& 0& 0& -b\end{array}\right)& (10)\end{array}$

The infection matrix F and rest of transmission matrix V are give as follow:

$\begin{array}{ll}\text{V}=\left(\begin{array}{ccccc}d& 0& 0& 0& 0\\ 0& a& 0& 0& 0\\ 0& -\kappa & \varphi & 0& 0\\ 0& 0& 0& h& 0\\ 0& 0& 0& 0& b\end{array}\right),& (11)\end{array}$

$\begin{array}{ll}F=\left(\begin{array}{ccccc}0& 0& -\frac{\mu \rho }{d}& 0& 0\\ 0& 0& \frac{\mu \rho }{d}& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\end{array}\right).& (12)\end{array}$

The matrix K = FV⁻¹ is constructed as follows:

$\begin{array}{ll}K=\left(\begin{array}{ccccc}0& -\frac{\mu \kappa \rho }{ad\varphi }& -\frac{\mu \rho }{d\varphi }& 0& 0\\ 0& \frac{\mu \kappa \rho }{ad\varphi }& \frac{\mu \rho }{d\varphi }& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\end{array}\right).& (13)\end{array}$

The spectral radius of matrix K = FV⁻¹ give reproductive number of model (1) is

$\begin{array}{ll}{ℝ}_0=\frac{\mu \kappa \rho }{ad\varphi }.& (14)\end{array}$

The stability analysis of system (1) is presented in following Theorem 2.4.

Theorem 2.4. Uninfected equilibrium point E₀, will be locally asymptotically stable if ℝ₀ < 1 and otherwise it will be unstable.

Proof: The Jacobian matrix J(E₀) of model (1) for uninfected equilibrium point is formulated as follows:

$\begin{array}{ll}J(E_0)=\left(\begin{array}{ccccc}-d& 0& -\frac{\mu \rho }{d}& 0& 0\\ 0& -a& \frac{\mu \rho }{d}& 0& 0\\ 0& \kappa & -\varphi & 0& 0\\ 0& 0& 0& -h& 0\\ 0& 0& 0& 0& -b\end{array}\right).& (15)\end{array}$

The eigenvalues of Jacobian matrix J(E₀) is evaluated by characteristic equation that is det(J(E₀) − Iλ_i) = 0, for i = 1 : 5. The eigenvalues are given as follows:

λ₁ = −b, λ₂ = −d, λ₃ = −h, ${\lambda }_4=-\frac{\left(\sqrt{d}\right)(a+\varphi )+\sqrt{d{(a-\varphi )}^2+4\mu \kappa \rho }}{2\sqrt{d}}$ and ${\lambda }_5=-\frac{\sqrt{d}(a+\varphi )-\sqrt{d{(a-\varphi )}^2+4\mu \kappa \rho }}{2\sqrt{d}}$.

Thus, all eigenvalues λ_i are strictly negative for i = 1 : 5. Hence, model (1) is locally asymptotically stable.

Lemma 2.5. If R₀ > 1, the infected equilibrium point of system is locally asymptotically stable, otherwise it is unstable.

Proof: The results can by obtained by the same procedure of (2.4).

2.1.4. Sensitivity Analysis

The sensitivity of basis reproductive

$\begin{array}{ll}{ℝ}_0=\frac{\mu \kappa \rho }{ad\varphi }& (16)\end{array}$

is analyze with respect to each parameters is follows:

$\begin{array}{ll}\begin{array}{l}\frac{\partial {ℝ}_0}{\partial \mu }=\frac{\kappa \rho }{ad\varphi }>0,\\ \frac{\partial {ℝ}_0}{\partial \kappa }=\frac{\mu \rho }{ad\varphi }>0,\\ \frac{\partial {ℝ}_0}{\partial \rho }=\frac{\mu \kappa }{ad\varphi }>0,\\ \frac{\partial {ℝ}_0}{\partial a}=-\frac{\mu \kappa \rho }{a^{2}d\varphi }<0,\\ \frac{\partial {ℝ}_0}{\partial d}=-\frac{\mu \kappa \rho }{a{d}^2\varphi }<0,\\ \frac{\partial {ℝ}_0}{\partial \varphi }=-\frac{\mu \kappa \rho }{ad{\varphi }^2}<0.\end{array}& (17)\end{array}$

It is clear that reproductive number is directly proportional to the replication rate of the virus, rate of infected cells produced in free virus and rate of generation of susceptible host cells. Decrease with the death rate of infected cells, decays rate in free virus and the death rate of suspected host cells.

The basic reproductive number is ${ℝ}_0=\frac{\mu \kappa \rho }{ad\varphi }>1$ in case of infection. In the absence of an immune responses the model (1) converges to the following equilibrium points E⁽⁰⁾ = (U⁽⁰⁾, W⁽⁰⁾, V⁽⁰⁾, A⁽⁰⁾, C⁽⁰⁾) which is defined as

$\begin{array}{ll}{U}^{(0)}=(\frac{\rho }{d},0,0,0,0).& (18)\end{array}$

We assumed that the immune responses is present. this desires following conditions: cW⁽⁰⁾ > b and gV⁽⁰⁾ > h. In this case following can be observed.

1. The anti body response can't established and the CTL response develops. Because the CTL response is strong and decrease virus load to levels which are vary low to stimulate antibody response. It has following equilibrium points,

E⁽¹⁾ = (U⁽¹⁾, W⁽¹⁾, V⁽¹⁾, A⁽¹⁾, C⁽¹⁾) where,

$\begin{array}{ll}\begin{array}{l}{U}^{(1)}=\frac{\text{c}\varphi \rho }{b\kappa \mu +cd\varphi },\\ W^{(1)}=\frac{b}{c},\\ V^{(1)}=\frac{b\kappa }{c\varphi },\\ A^{(1)}=0,\\ C^{(1)}=\frac{-ab\kappa \mu -acd\varphi +c\kappa \mu \rho }{p(b\kappa \mu +cd\varphi )}.\end{array}& (19)\end{array}$

This is obtained if $\frac{bg\kappa }{c\varphi }<h$ and $\frac{ch\rho \mu }{a(dg+h\mu )}>b$.

2. The sustained CTL response zero and antibody response develop, because the antibody response is strong relative to CTL response and decrease virus load to levels which are vary low to stimulate CTL. It has following equilibrium points.

$\begin{array}{l}{E}^{(2)}=(U^{(2)},W^{(2)},V^{(2)},A^{(2)},C^{(2)}),\end{array}$

which is defined as follows:

$\begin{array}{ll}\begin{array}{l}{U}^{(2)}=\frac{g\rho }{dg+h\mu },\\ W^{(2)}=\frac{h\mu \rho }{a(dg+h\mu )},\\ V^{(2)}=\frac{h}{g},\\ A^{(2)}=\frac{-adg\varphi -ah\mu \varphi +g\kappa \mu \rho }{aq(dg+h\mu )},\\ C^{(2)}=0.\end{array}& (20)\end{array}$

This is obtained if $\frac{bgk}{c\varphi }>h$ and $\frac{ch\rho \mu }{a(dg+h\mu )}<b$.

3. CTL and anti body response develops. This equilibrium points are as follows:

$\begin{array}{ll}{E}^{(3)}=(U^{(3)},W^{(3)},V^{(3)},A^{(3)},C^{(3)})& (21)\end{array}$

where,

$\begin{array}{ll}\begin{array}{l}{U}^{(3)}=\frac{g\rho }{dg+h\mu },\\ W^{(3)}=\frac{b}{c},\\ V^{(3)}=\frac{h}{g},\\ A^{(3)}=\frac{bg\kappa -ch\varphi }{chq},\\ C^{(3)}=\frac{ch\mu \rho -abdg-abh\mu }{bp(dg+h\mu )}.\end{array}& (22)\end{array}$

This is obtained if $\frac{bgk}{c\varphi }>h$ and $\frac{ch\rho \mu }{a(dg+h\mu )}>b$.

2.2. Modified Modeling Approach

$\begin{array}{ll}\begin{array}{l}\frac{dU}{dt}=\rho -dU-\mu \frac{VU}{1+\eta U},\\ \frac{dW}{dt}=\mu \frac{VU}{1+\eta U}-aW-pWC,\\ \frac{dV}{dt}=\kappa W-qVA-\varphi V,\\ \frac{dA}{dt}=gVA-hA,\\ \frac{dC}{dt}=cWC-bC,\end{array}& (23)\end{array}$

where the description of compartments U(t), W(t), V(t), A(t), and C(t) are presented in Table 2.

2.2.1. Positivity of Solution

For the model (23) transmission to be epidemiologically feasible, it is necessary to show that for all time every state variable is non-negative. Thus, the solutions of the model with non-negative initial value for all time t > 0 will remain non-negative.

Theorem 2.6. Suppose that the model(23) consists of all feasible solutions with non-negative initial result then it remains non-negative for all time t.

Proof: Let the model (23) satisfy that the initial non-negative solution, i.e.,

U(0) ≥ 0, W(0) ≥ 0, V(0) ≥ 0, A(0) ≥ 0, and C(0) ≥ 0.

It can be deduced from model (23) that is,

$\begin{array}{ll}\frac{dU}{dt}=\rho -dU-\mu \frac{VU}{1+\eta U}.& (24)\end{array}$

As the solution of variable U can be computed by following result,

$\begin{array}{ll}U=U(0)e^G+{\displaystyle {\int }_0^t}\pi e^{H(k)}d(k)\ge 0& (25)\end{array}$

where $G=-{\int }_0^{t}G(k)d(k)$ and $H(k)=-{\int }_0^{t}G(l)d(l)$. This gives non-negativity of U(0) ≥ 0 for all t ≥ 0. The non-negativity of rest variables in the model (23) is given as follows,

$\begin{array}{ll}\begin{array}{l}\frac{dW}{dt}=\mu \frac{VU}{1+\eta U}-aW-pWC,\\ \frac{dV}{dt}=\kappa W-qVA-\varphi V,\\ \frac{dA}{dt}=gVA-hA,\\ \frac{dC}{dt}=cWC-bC.\end{array}& (26)\end{array}$

In term of matrix the above can be expressed as,

$\begin{array}{ll}\frac{dF(t)}{dt}=H(t)+MF(t)& (27)\end{array}$

where F(t) and H(t) are as follows:

$\begin{array}{ll}F(t)=\left(\begin{array}{c}W(t)\\ V(t)\\ A(t)\\ C(t)\end{array}\right),H(t)=\left(\begin{array}{c}a\\ 0\\ 0\\ 0\end{array}\right)& (28)\end{array}$

and we have matrix M

$\begin{array}{ll}M=\left(\begin{array}{cccc}-a& \frac{\mu \rho }{d\left(\frac{\eta \rho }{d}+1\right)}& 0& 0\\ \kappa & -\varphi & 0& 0\\ 0& 0& -h& 0\\ 0& 0& 0& -b\end{array}\right).& (29)\end{array}$

The matrix M is a Matzler matrix so by the result presented in Smith (1996) the model is monotone. the fact that $R_{+}^4$ is invariant with respect to stream of model (23).

2.2.2. Qualitative Analysis

Here we present the stability analysis of (23). The model (23) is locally asymptotically stable at uninfected and infected equilibrium points. We investigate the model (23) is locally stable at uninfected and infected. The uninfected equilibrium point E⁰ = (U⁰, 0, 0, 0, 0) and infected equilibrium point E^* = (U^*, W^*, V^*, A^*, C^*) of model (23) are constructed by following theorems.

Theorem 2.7. The uninfected equilibrium point E⁰ of the model (23) is given by

$\begin{array}{l}{E}^0=(U^0,0,0,0,0),\end{array}$

where:

$\begin{array}{l}{U}^0=\frac{\rho }{d}.\end{array}$

Proof: By putting the equations in (23) equal to zero, the dynamics is determined as follows,

$\begin{array}{r}\begin{array}{r}\rho -dU-\mu \frac{VU}{1+\eta U}=0,\\ \mu \frac{VU}{1+\eta U}-aW-pWC=0,\\ \kappa W-qVA-\varphi V=0\\ gVA-hA=0,\\ cWC-bC=0.\end{array}\end{array}$

By solving these equations, we got

$\begin{array}{ll}{U}^0=\frac{\rho }{d}.& (30)\end{array}$

This completes the proof.

Theorem 2.8. The model (23) admits a unique infected equilibrium E^* = (U^*, W^*, V^*, A^*, C^*) if and only if ℝ₀ > 1.

Proof: Calculating the infected equilibrium point, we obtain

$\begin{array}{ll}\begin{array}{l}{U}^{*}=\frac{-\sqrt{4d\eta g^2\rho +{(dg-\eta g\rho +h\mu )}^2}}{2dg\eta }\\ +\frac{\eta g\rho -dg-h\mu }{2dg\eta },\\ W^{*}=\frac{b}{c},\\ V^{*}=\frac{h}{g},\\ A^{*}=\frac{bg\kappa -ch\varphi }{chq},\\ C^{*}=\frac{-ab-cd{U}^{*}+c\rho }{bp}.\end{array}& (31)\end{array}$

2.2.3. Reproduction Number ℝ₀ and Stability Analysis

The basic reproduction number ℝ₀ is formulated by evaluating infection matrix F and transmission matrix V constructed from Jacobian matrix J of system (23). The reproduction number ℝ₀ is the spectral radius of matrix FV⁻¹. The Jacobian matrix for J₀ for uninfected equilibrium point is,

$\begin{array}{ll}{J}_0=\left(\begin{array}{ccccc}0& 0& -\frac{\mu \rho }{d\left(\frac{\eta \rho }{d}+1\right)}& 0& 0\\ 0& -a& \frac{\mu \rho }{d\left(\frac{\eta \rho }{d}+1\right)}& 0& 0\\ 0& \kappa & -\varphi & 0& 0\\ 0& 0& 0& -h& 0\\ 0& 0& 0& 0& -b\end{array}\right).& (32)\end{array}$

The infection and transmission matrices, F & V are give as follow:

$\begin{array}{ll}\text{V}=\left(\begin{array}{ccccc}0& 0& \frac{\mu \rho }{d\left(\frac{\eta \rho }{d}+1\right)}& 0& 0\\ 0& a& -\frac{\mu \rho }{d\left(\frac{\eta \rho }{d}+1\right)}& 0& 0\\ 0& -\kappa & \varphi & 0& 0\\ 0& 0& 0& h& 0\\ 0& 0& 0& 0& 0\end{array}\right),& (33)\end{array}$

$\begin{array}{ll}F=\left(\begin{array}{cccc}0& -\frac{\mu \rho }{d\left(\frac{\eta \rho }{d}+1\right)}& 0& 0\\ 0& \frac{\mu \rho }{d\left(\frac{\eta \rho }{d}+1\right)}& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right).& (34)\end{array}$

The matrix K = FV⁻¹ is constructed as follows:

$\begin{array}{ll}K=\left(\begin{array}{cccc}-\frac{\kappa \mu \rho }{ad\left(\frac{\eta \rho }{d}+1\right)\varphi }& -\frac{\mu \rho }{d\left(\frac{\eta \rho }{d}+1\right)\varphi }& 0& 0\\ \frac{\kappa \mu \rho }{ad\left(\frac{\eta \rho }{d}+1\right)\varphi }& \frac{\mu \rho }{d\left(\frac{\eta \rho }{d}+1\right)\varphi }& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right).& (35)\end{array}$

The spectral radius of matrix K = FV⁻¹ finally provides the reproductive number for given model (23):

$\begin{array}{ll}{ℝ}_0=\frac{\mu \rho (a-\kappa )}{a\varphi (d+\eta \rho )}.& (36)\end{array}$

The stability analysis of system (23) is presented by follows theorem.

Theorem 2.9. For uninfected equilibrium point E₀ if the real part of eigenvalues of Jacobian matrix J of the system (23) are strictly negative, then the system (23) is locally stable otherwise it is unstable.

Proof: The Jacobian matrix of model (23) for uninfected equilibrium point J(E₀) is formulated as follows:

$\begin{array}{ll}J(E_0)=\left(\begin{array}{ccccc}0& 0& -\frac{\mu \rho }{d\left(\frac{\eta \rho }{d}+1\right)}& 0& 0\\ 0& -a& \frac{\mu \rho }{d\left(\frac{\eta \rho }{d}+1\right)}& 0& 0\\ 0& \kappa & -\varphi & 0& 0\\ 0& 0& 0& -h& 0\\ 0& 0& 0& 0& -b\end{array}\right).& (37)\end{array}$

The eigenvalues of Jacobian matrix J(E₀) is evaluated by characteristic equation det(J(E₀) − Iλ_i) = 0, for i = 1 : 5. The eigenvalues are given as follows:

λ₁ = 0, λ₂ = −b, λ₃ = −h,

${\lambda }_4=\frac{1}{2}\left(\sqrt{{(a-\varphi )}^2+\frac{4\kappa \mu \rho }{d+\eta \rho }}+a+\varphi \right)$

and

${\lambda }_5=\frac{1}{2}\left(-\sqrt{{(a-\varphi )}^2+\frac{4\kappa \mu \rho }{d+\eta \rho }}+a+\varphi \right)$.

Since the all eigenvalues λ_i, satisfy the criteria (where as i = 1 : 5). Hence, the model (23) is locally stable.

2.2.4. Sensitivity Analysis

The sensitivity of basis reproductive

$\begin{array}{ll}{ℝ}_0=\frac{\mu \rho (a-\kappa )}{a\varphi (d+\eta \rho )}& (38)\end{array}$

is analyze with respect to each parameters is as follows:

$\begin{array}{l}\begin{array}{l}\frac{\partial {ℝ}_0}{\partial a}=\frac{\kappa \mu \rho }{a^2\varphi (d+\eta \rho )}>0,\\ \frac{\partial {ℝ}_0}{\partial \mu }=\frac{a\rho -\kappa \rho }{ad\varphi +a\eta \rho \varphi }>0,\end{array}\\ \begin{array}{l}\frac{\partial {ℝ}_0}{\partial \rho }=\frac{d\mu (a-\kappa )}{a\varphi {(d+\eta \rho )}^2}>0,\\ \frac{\partial {ℝ}_0}{\partial \kappa }=-\frac{\mu \rho }{ad\varphi +a\eta \rho \varphi }<0\\ \frac{\partial {ℝ}_0}{\partial d}=-\frac{\mu \rho (a-\kappa )}{a\varphi {(d+\eta \rho )}^2}<0,\\ \frac{\partial {ℝ}_0}{\partial \varphi }=-\frac{\mu \rho (a-\kappa )}{a{\varphi }^2(d+\eta \rho )}<0,\\ \frac{\partial {ℝ}_0}{\partial \eta }=-\frac{\mu {\rho }^2(a-\kappa )}{a\varphi {(d+\eta \rho )}^2}<0.\end{array}& (39)\end{array}$

Therefore, the reproductive number ℝ₀ is directly proportional to a, μ, and ρ and inversely proportional to κ, d, ϕ and η.

The basic reproductive number is ${ℝ}_0=\frac{\mu \kappa \rho }{ad\varphi }>1$ in case of infection.

In the absence of an immune responses the model (23) converges to the following equilibrium points: E⁽⁰⁾ = (U⁽⁰⁾, W⁽⁰⁾, V⁽⁰⁾, A⁽⁰⁾, C⁽⁰⁾) which is defined as $U^{(0)}=\frac{\rho }{d},W^{(0)}=0,V^{(0)}=0,A^{(0)}=0,C^{(0)}=0$. We assumed that the immune responses is present. This desires following conditions: cW⁽⁰⁾ > b and gV⁽⁰⁾ > h.

In this case, the following can be observed:

1. The anti body response can not be established and the CTL response develops.

Because of strong CTL response and low virus load to levels, it has the following equilibrium points:

$\begin{array}{ll}\begin{array}{l}{E}^{(1)}=(U^{(1)},W^{(1)},V^{(1)}{A}^{(1)},C^{(1)}),\\ E^{(1)}=(X,\frac{b}{c},\frac{b\kappa }{c\varphi },0,Y),\end{array}& (40)\end{array}$

where:

$\begin{array}{l}{U}^{(1)}=\frac{-b\kappa \mu -c\varphi (d-\eta \rho )}{2cd\eta \varphi },\\ -\frac{\sqrt{{(b\kappa \mu +c\varphi (d-\eta \rho ))}^2+4c^{2}d\eta \rho {\varphi }^2}}{2cd\eta \varphi },\\ W^{(1)}=\frac{b}{c},\\ V^{(1)}=\frac{b\kappa }{c\varphi },\\ A^{(1)}=0,\\ C^{(1)}=\frac{-ab-cd{U}^{(1)}+c\rho }{bp}.& (41)\end{array}$

This equilibrium points is obtain if it follows the following condition:

$\frac{bg\kappa }{c\varphi }<h$ and $\frac{c\sqrt{4d\eta g^2\rho +{(-dg+\eta g\rho -h\mu )}^2}+cdg+c\eta g\rho +ch\mu }{2ag\eta }>b$.

2. For null sustained CTL response, the antibody response develops. Because of strong antibody response to relative CTL response and low virus load to levels which is too small for stimulate CTL. It has the following equilibrium points:

$\begin{array}{ll}{E}^{(2)}=(U^{(2)},W^{(2)},V^{(2)},A^{(2)},C^{(2)}),& (42)\end{array}$

which is defined as follows:

$\begin{array}{l}\begin{array}{l}{U}^{(2)}=\frac{-dg+\eta g\rho -h\mu }{2dg\eta },\\ -\frac{\sqrt{4d\eta g^2\rho +{(dg-\eta g\rho +h\mu )}^2}}{2dg\eta },\\ W^{(2)}=\frac{g\rho }{a}-U^{(2)},\end{array}\end{array}$

$\begin{array}{ll}\begin{array}{l}{V}^{(2)}=\frac{h}{g},\\ A^{(2)}=\frac{g\kappa (\rho -d{U}^{(2)})-ah\varphi }{ahq},\\ C^{(2)}=0.\end{array}& (43)\end{array}$

The equilibrium point satisfies the following criteria: $\frac{bg\kappa }{c\varphi }>h$ and $\frac{c\sqrt{4d\eta g^2\rho +{(-dg+\eta g\rho -h\mu )}^2}+cdg+c\eta g\rho +ch\mu }{2ag\eta }<b$.

3. CTL and anti body response develops. The equilibrium point is given as:

$\begin{array}{ll}{E}^{*}=(U^{*},W^{*},V^{*},A^{*},C^{*})& (44)\end{array}$

where:

$\begin{array}{l}{U}^{*}=\frac{-\sqrt{4d\eta g^2\rho +{(dg-\eta g\rho +h\mu )}^2}}{2dg\eta },\\ +\frac{\eta g\rho -dg-h\mu }{2dg\eta },\\ W^{*}=\frac{b}{c},\\ V^{*}=\frac{h}{g},\\ A^{*}=\frac{bg\kappa -ch\varphi }{chq},\\ C^{*}=\frac{-ab-cd{U}^{*}+c\rho }{bp},& (45)\end{array}$

under the criteria: $\frac{bg\kappa }{c\varphi }>h$ and $\frac{c\sqrt{4d\eta g^2\rho +{(-dg+\eta g\rho -h\mu )}^2}+cdg+c\eta g\rho +ch\mu }{2ag\eta }>b$.

2.3. Delayed Model

Different biological models with delayed transmission are available in the literature (Maiti et al., 2008; Lv and Yuan, 2009; Kuniya and Nakata, 2012). The schematic 2 presents the transmission of the virus from one compartment to the other. The delay takes place after the interaction of virus with the target cells at τ₁, it is further assumed that at delay τ₂, the infected cells and the CTLs interact and the antibodies and the virus interact.

$\begin{array}{ll}\begin{array}{l}\frac{dU(t)}{dt}=\rho -dU(t)-\mu \frac{V(t)U(t)}{1+\eta U(t)},\\ \frac{dW(t)}{dt}=\frac{\mu V(t-{\tau }_1)U(t-{\tau }_1)}{1+\eta U(t-{\tau }_1)}-aW(t),\\ -pW(t)C(t),\\ \frac{dV(t)}{dt}=\kappa W(t)-qV(t)A(t)-\varphi V(t),\\ \frac{dA(t)}{dt}=gV(t-{\tau }_2)A(t-{\tau }_2)-hA(t),\\ \frac{dC(t)}{dt}=cW(t-{\tau }_2)C(t-{\tau }_2)-bC(t).\end{array}& (46)\end{array}$

2.3.1. Local Stability for Delay Differential Equations

In ODE's, the local stability of consistent condition relies upon the area of underlying foundations of characteristic function, that are polynomial in shape. The unfaltering condition is steady if the majority of the roots having −ve real part. The outstanding Routh-Hurwitz criteria provide exact situation for arbitrary polynomials. For DDE's, nearby stability is likewise controlled by the area of trademark work, yet for this situation, this capacity appears as an alleged quasi polynomial, that is supernatural. Hence, there are vastly numerous roots. Moreover, the Routh-Hurwitz criteria are not pertinent. Numerous methodologies decide the stability of steady states delay equations.

2.3.2. Domain Subdivision

D-subdivision or Domain subdivision, utilizes basic details about the actions of the roots of characteristic functions, as a parameter modifies to split parameter space into parts where the no. of roots with +ve real parts is constant. The roots' position depends consistently on the models' parameter and when the parameters are altered, another root rises if imaginary root exists for a set of parameters.

Now we subdivide the parametric domain by hypersurfaces comprising of parameter routines for which at least one simply imaginary roots exist. When the districts are limited by these hypersurfaces, the quantity of roots with +ve real part is constant. Obviously, the locales where the number is zero and their supplements are more interesting. This technique is especially simple to picture when the framework in question relies upon two parameters, so the area is ℝ², and, the hypersurfaces are curve.

3. Results and Discussion

During this research, we have used the Matlab^TM delay differential equations bifurcation analysis tools. For the validity of the computational model, extensive numerical experiments were conducted using simulink toolbox.

3.1. Analysis of SARS-CoV2 and Antibodies Interaction

With the aid of mathematical model we have concluded the following:

1. The analysis of lytic vs. non lytic immune response plays an important role in infection control.

2. The Hill function is important in kinetic modeling and the Hill coefficient is important parameter to forecast a complete cycle of infection.

3. The analytic approach and numerical Matcont bifurcation analysis proved to be efficient in parametric approximation for such complex dynamics.

Figure 1 presents the schematic to understand the interaction of key players. Figure 2 provides the phase space portraits to explore the interaction between the CTL's, the anti bodies and the virus, for three different values of q. We can see that for increasing values of q, there is reduction in the concentration of CTL, as well as the length of the cycle increases over time (top panel). On the other hand, when delay was considered, the dynamics were different. Figure 3 provides analysis relative to the replication rate of virus. Figure 4 presents the dynamics relative to the parameter g. We can see a twist in the phase space portraits when the delay was taken into account. In the supplementary figures, (Figures S1 and S2, we can see the dynamics more clearly relative to change in parametric values). We can see that the dynamics are more visible to witness the rapid action of SARS-2. Thus, a model without delay, and with η = 0 will not be able to demonstrate the dynamics well. We thus conclude that the disease transmission and the immune response depends on time delay as well as nonlinear Hill formalism.

FIGURE 1

Figure 1. Schematic of the model.

FIGURE 2

Figure 2. For three different values of antibody-SARS-2 interaction rates (q). (Top) without, (Bottom) with delay τ₁ = 0.5, τ₂ = 0.08.

FIGURE 3

Figure 3. For three different values of the replication rate of virus. (Top) without, (Bottom) with delay τ₁ = 0.5, τ₂ = 0.08.

FIGURE 4

Figure 4. For three different values of g. (Top) without, (Bottom) (with) delay τ₁ = 0.5, τ₂ = 0.08.

4. Conclusions and Future Work

The manuscript presents a state of the art model, with delay, from one compartment to the next due to transition. In nature, there is always a delay in the onset of infections. The proposed mathematical model quantifies and analyzes this imbalance and describes the temporal trend of the phenomenon, leaving its application open to possible direct therapies in that sense.

• Mathematical analysis of the non-lytic and lytic action of the immune reaction to SARS-CoV2.

• Construction of a model that describes the balance between antibody reaction and cellular reaction mediated by CTL cells.

• Analysis of the imbalance between non-lytic and lytic action of the immune response.

• Description and quantification of the model related to the infection and functionality of CTL cells over time.

• Evidence of delay in disease transmission.

Data Availability Statement

The original contributions generated for the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

Author Contributions

ZY and AS did the modeling. AN and RA did the literature review and simulations. All of the authors equally contributed to the manuscript and participated in results and discussion.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmolb.2020.585245/full#supplementary-material

References

di Mauro Gabriella, S. C., Concetta, R., Francesco, R., and Annalisa, C. (2020). Sars-cov-2 infection: response of human immune system and possible implications for the rapid test and treatment. Int. Immunopharmacol. 84:106519. doi: 10.1016/j.intimp.2020.106519

PubMed Abstract | CrossRef Full Text | Google Scholar

Fang, J., Liu, C., Simos, T., and Famelis, I. T. (2020). Neural network solution of single-delay differential equations. Mediterr. J. Math. 17:30. doi: 10.1007/s00009-019-1452-5

CrossRef Full Text | Google Scholar

Kumar, S., Maurya, V. K., Prasad, A. K., Bhatt, M. L., and Saxena, S. K. (2020). Structural, glycosylation and antigenic variation between 2019 novel coronavirus (2019-nCov) and SARS coronavirus (SARS-CoV). VirusDisease 31, 13–21. doi: 10.1007/s13337-020-00571-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Kuniya, T., and Nakata, Y. (2012). Permanence and extinction for a nonautonomous seirs epidemic model. Appl. Math. Comput. 218, 9321–9331. doi: 10.1016/j.amc.2012.03.011

CrossRef Full Text | Google Scholar

Li, X., Geng, M., Peng, Y., Meng, L., and Lu, S. (2020). Molecular immune pathogenesis and diagnosis of COVID-19. J. Pharm. Anal. 10, 102–108. doi: 10.1016/j.jpha.2020.03.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Long, Q.-X., Liu, B.-Z., Deng, H.-J., Wu, G.-C., Deng, K., Chen, Y.-K., et al. (2020). Antibody responses to SARS-CoV-2 in patients with COVID-19. Nat. Med. 26, 845–848. doi: 10.1038/s41591-020-0897-1

CrossRef Full Text | Google Scholar

Lv, C., and Yuan, Z. (2009). Stability analysis of delay differential equation models of hiv-1 therapy for fighting a virus with another virus. J. Math. Anal. Appl. 352, 672–683. doi: 10.1016/j.jmaa.2008.11.026

CrossRef Full Text | Google Scholar

Maiti, A., Pal, A., and Samanta, G. (2008). Effect of time-delay on a food chain model. Appl. Math. Comput. 200, 189–203. doi: 10.1016/j.amc.2007.11.011

CrossRef Full Text | Google Scholar

Smith, H. L. (1996). Monotone dynamical systems: an introduction to the theory of competitive and cooperative systems. Bull. Am. Math. Soc. 33, 203–209. doi: 10.1090/S0273-0979-96-00642-8

CrossRef Full Text | Google Scholar

Sohail, A., and Nutini, A. (2020). Forecasting the timeframe of coronavirus and human cells interaction with reverse engineering. Prog. Biophys. Mol. Biol. 155, 29–35. doi: 10.1016/j.pbiomolbio.2020.04.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Tunc, C. (2002). A study of the stability and boundedness of the solutions of nonlinear differential equations of the fifth order. Indian J. Pure Appl. Math. 33, 519–529.

Google Scholar

Tunç, C. (2004). A note on the stability and boundedness results of solutions of certain fourth order differential equations. Appl. Math. Comput. 155, 837–843. doi: 10.1016/S0096-3003(03)00819-1

CrossRef Full Text | Google Scholar

Tunç, C. (2008). On the stability of solutions to a certain fourth-order delay differential equation. Nonlin. Dyn. 51, 71–81. doi: 10.1007/s11071-006-9192-z

CrossRef Full Text | Google Scholar

Tunç, C. (2013). New results on the stability and boundedness of nonlinear differential equations of fifth order with multiple deviating arguments. Bull. Malays. Math. Sci. Soc. 36, 671-682.

Google Scholar

Wodarz, D. (2005). Mathematical models of immune effector responses to viral infections: virus control versus the development of pathology. J. Comput. Appl. Math. 184, 301–319. doi: 10.1016/j.cam.2004.08.016

CrossRef Full Text | Google Scholar

Yan, X.-P. (2007). Stability and hopf bifurcation for a delayed prey-predator system with diffusion effects. Appl. Math. Comput. 192, 552–566. doi: 10.1016/j.amc.2007.03.033

CrossRef Full Text | Google Scholar

Zheng, M., and Song, L. (2020). Novel antibody epitopes dominate the antigenicity of spike glycoprotein in SARS-CoV-2 compared to SARS-CoV. Cell. Mol. Immunol. 17, 536–538. doi: 10.1038/s41423-020-0385-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhu, N., Zhang, D., Wang, W., Li, X., Yang, B., Song, J., et al. (2020). A novel coronavirus from patients with pneumonia in china, 2019. N. Engl. J. Med. 382, 727–733. doi: 10.1056/NEJMoa2001017

CrossRef Full Text | Google Scholar