Within-host delay differential model for SARS-CoV-2 kinetics with saturated antiviral responses

1.   Introduction

Since January 2020, a deadly pandemic known as COVID-19 has spread from person to person, wreaking havoc across the globe. The virus responsible for the illness is SARS-CoV-2, also known as the 2019-ncov ^[1]. Multiple cases of pneumonia, dry cough, fever, fatigue, breathing difficulty and bilateral lung infiltration have been reported due to the virus's spread within the human body ^[2,3]. According to the World Health Organization (WHO), as of March 18, 2023 there have been 760 million cases reported worldwide, with 6.8 million fatalities. Vaccines like Pfizer-BioNTech, Moderna (mRNA-1273), and COVAXIN are available for hospitality and to control the pandemic.

Mathematical models give more insight and knowledge into the dynamics of infectious diseases, which may help control the pandemic. Several researchers ^{[4,5,6,7,8,9]} have investigated the human-to-human transmission dynamics of the COVID-19 pandemic using mathematical modeling. Mathematical models incorporating time delays are also used for understanding infectious disease dynamics because they are more realistic and accurate due to the memory effect. Yang et al. ^[10] proposed a COVID-19 model with a time delay that characterized the viral infection cycle and treatment time. They fitted their model with the data from Beijing and Wuhan and suggested that early detection and isolation are the most important ways to prevent the spread of the epidemic. In ^[11], the effect of the time delay of the immune response on COVID-19 transmission dynamics has been investigated. A mathematical model has been designed in ^[12] to explore the importance of time lag in precautionary measures to control the COVID-19 pandemic. Babasola et al. ^[13] explored the role of time delay with a convex incidence rate in the COVID-19 epidemic. They observed that time delay destabilizes the system and generates periodic solutions. The effects of time delay in using vaccinations for COVID-19 spread dynamics has been investigated in ^[14].

Moreover, controlling the COVID-19 pandemic also requires a thorough understanding of the within-host dynamics of the SARS-CoV-2 virus. For this purpose, numerous modeling studies have been performed to predict the behavior of the SARS-CoV-2 virus at the cellular level. Both viral and host factors play a significant role in SARS-CoV-2 infections. Host factors during infection trigger an immune response that combats the virus. The host's innate immune system can identify viral infections by utilizing pattern recognition receptors to identify pathogen-associated molecular patterns. As the innate immune system detects the viruses, it typically triggers the pulmonary and systemic inflammatory reactions linked to the SARS-CoV-2 virus. According to a modeling study ^[15], the innate immune response may be able to clear the virus more effectively if the adaptive immune response is temporarily suppressed. The kinetics of the human SARS-CoV-2 infection have been mathematically modeled in ^[16]. Prakash et al. ^[17] proposed a multi-scale model that incorporates both the intra-host and inter-host dynamics of COVID-19. Their findings suggest that treatment with antiviral drugs, immunotherapies, and improved sanitation and social isolation was proposed as the most effective method for lowering the virus's basic reproduction number and environmental spread in a population. Chhetri et al. ^[18] studied the role of various drugs at various stages of COVID-19 pathogenesis using a mathematical model of the complex interaction of virus replication and the host immune response. The effector T cell response to SARS-CoV-2 and the dynamics of the virus were studied in ^[19]. The dynamics of SARS-CoV-2 in infected patients were simulated in ^[20]. Chowdhury et al. ^[21] analyzed a mathematical model to investigate the impact of natural killer cells and T cells on SARS-CoV-2 kinetics. Li et al. ^[22] proposed a within-host SARS-CoV-2 model to observe drug effects on virus growth and the immunity effect of patients. In ^[23], dynamical analysis of the model ^[22] was studied. Ghosh et al. ^[24] observed the within-host dynamics of SARS-CoV-2 in the presence of innate and adaptive immune responses and then observed the effect of vaccination and antiviral treatments. Stochastic dynamics of the within-host COVID-19 epidemic with discrete time delay and noise have been studied in ^[25] and showed the impact of delay tactics and noise on the extinction of the disease. Staroverov et al. ^[26] proposed a novel mathematical model for SARS-CoV-2 dynamics that explicitly modeled intracellular events like the exhaustion of cellular resources needed for virus production and innate immune responses. Pillis et al. ^[27] showed a mathematical model of how SARS-CoV-2 spreads within a host and how neutralizing antibodies react to vaccination.

It is found in the above literature that during the infection by the virus, the epithelial cells resist the viral infection by activating the immune system. Thus, considering this fact, we modified the model discussed in ^[22,23]. So, in this study, we consider a model incorporating the antiviral responses of uninfected epithelial cells in terms of Michaelis-Menten, which is suitably described as the "saturated responses". In addition, it takes some time for uninfected epithelial cells to develop a suitable antiviral response after contracting a virus. Therefore, we introduce a time delay into the saturated response term. The rest of the manuscript is assembled as follows: In Section 2, we have formulated our model. The basic characteristics of the solution of the model have been discussed in Section 3. The stability of the model and estimation of the time lag for preserving a stable limit cycle have been investigated in Section 4. The numerical simulation has been performed in Section 5. A concluding remark has been made in Section 6.

2.   Model formulation

This section will formulate a model that describes the within-host SARS-CoV-2 viral kinetics. The considered model was originally proposed by Li et al. ^[22], and the model is described by the following equations:

The numbers of pulmonary epithelial cells that are uninfected, infected and the virus concentration are represented here as $E(t)$, $I(t)$, and $v(t)$, respectively. The rate at which virus-free epithelial cells are infected by the virus is indicated by the symbol $\beta$. The virus-free pulmonary epithelial cells are produced at a rate of $d_{1}E(0)$ and die at a rate of $d_{1}$. The term $d_{2}$ represents the death rate of virus-infected epithelial cells. Also, $\mu$ and $d_{3}$ represent the production rate and death rate of viruses, respectively. The epithelial cells are infected once the virus comes into contact with a human being. In addition, epithelial cells can fight off the virus during infection by triggering the immune system to control antiviral responses. Since uninfected epithelial cells resist viral infection by secreting immune cytokines and chemokines, we modified the model (2.1) to describe this behavior by the nonlinear term $\frac{\beta E(t)v(t)}{1+v(t)}$ ^[28,29,30]. In addition, uninfected epithelial cells need time to develop a proper antiviral response after infection by the virus. Therefore, there is a time lag to regulate the antiviral response and this time delay can be introduced as a discrete-time delay, $\tau$ into the term, $\frac{\beta E(t)v(t)}{1+v(t)}$. Thus, the considered model becomes,

subject to the initial conditions:

where

where the continuous functions $\phi_{i}$ are defined on the interval $\theta \in [-\tau, 0]$.

3.   Nonnegativity and boundedness

In this section, we will discuss the nonnegativity and boundedness of solutions for the considered model.

Proposition 3.1. Corresponding to initial conditions (2.3), all solutions of the system (2.2) are in $\mathcal{R}_{+}^{3}$ and bounded for all $t > 0$.

Proof. The epithelial cells are subdivided into uninfected and infected epithelial cells. Thus the first two equations of the system (2.2) give us

where $\delta = min\{d_{1}, d_{2}\}$. Thus, an upper bound always exists for uninfected and infected epithelial cells. It is easy to see from the model's third equation that the free virus population $v$ is also bounded above.

Now, the system (2.2) can be re-written in vector notation as

with $X = [E(t), I(t), v(t)]^{T}\in \mathcal{R}_{+}^{3}$ and

where, $\Gamma \in \mathcal{C}^{\infty}(\mathcal{R}_{+}^{3})$ and $\Gamma:\mathcal{R}_{+}^{3} \rightarrow \mathcal{R}_{+}^{3}$. The righthand side of system (2.2) is locally Lipchitz, thus the derivatives are bounded and satisfy

Also, using the second lemma of ^[31], the solution of system (2.2) remains positive throughout the domain $\mathcal{R}_{+}^{3}$, $\forall t > 0$ corresponding to initial conditions (2.3). Thus, the solutions of model (2.2) are nonnegative and bounded for time $t > 0$. This completes the proof. □

4.   Stability analysis

4.1.   Basic reproduction number

The basic reproduction number in virus dynamics determines whether the disease will persist or die. It also controls disease spread in the population. We will use the method described in Driessche and Watmough ^[32] to calculate the basic reproduction number of model (2.2).

The two infected compartments in the model (2.2) are $I$ and $v$. If $\mathcal{F}_{i}$ and $\mathcal{V}_{i}$ denote the appearance rate of new infections in the $i^{th}$ compartment and the transfer rate of individuals from the $i^{th}$ compartment, respectively, then we have two matrices $\mathcal{F}$ and $\mathcal{V}$ as follows:

The Jacobian matrices of $\mathcal{F}$ and $\mathcal{V}$ at the virus-free equilibrium point, denoted by $F$ and $V$ respectively, are evaluated as

Now, we compute the next generation matrix $FV^{-1}$ as follows

Following Driessche and Watmough ^[32], the basic reproduction number $\chi_{0}$ of the model (2.2) is the spectral radius of the matrix $FV^{-1}$ which is given by $\chi_0 = \frac{\mu\beta E(0)}{d_{2}d_{3}}$.

4.2.   Equilibrium points and local stability

In this section, we will investigate the existence of the equilibrium points for the model (2.2) and study their local stability. As the time delay $\tau$ neither affects the number nor the type of equilibrium points; the existing equilibrium points of the model (2.2) are

● $P^{\prime}(E(0), 0, 0)$, uninfected or virus-free equilibrium point, which always exists.

● $P^{*}(E^{*}, I^{*}, v^{*})$, coexisting equilibrium point; where,

Clearly, the coexisting equilibrium $P^{*}(E^{*}, I^{*}, v^{*})$ exists if $\chi_{0} > 1$.

For local stability analysis, we consider two cases:

Case 1: $\tau = 0$

$({\rm{i}})$ The Jacobian matrix of the system (2.2) at uninfected or virus-free equilibrium point $P^{\prime}(E(0), 0, 0)$ is

The characteristic equation of the matrix $J_{P^{\prime}}^{0}$ is

where

The Routh-Hurwitz stability criterion states the equilibrium point $P^{\prime}(E(0), 0, 0)$ is locally asymptotically stable if $A_{11} > 0$, $A_{12} > 0$, $A_{13} > 0$ and $A_{11}A_{12}-A_{13} > 0$. These inequalities imply that $P^{\prime}(E(0), 0, 0)$ is locally asymptotically stable if $\chi_{0} < 1 \implies \mu < \frac{d_{2}d_{3}}{\beta E(0)}$; otherwise, it is unstable. Thus, if the production rate of the virus is less than the ratio between the product of the decay rate of infected epithelial cells and virus concentration and the product of the infection rate of the virus and density of uninfected epithelial cells, then the virus-free equilibrium point $P^{\prime}(E(0), 0, 0)$ is locally asymptotically stable.

$({\rm{ii}})$ Corresponding to the co-existing equilibrium point $P^{*}(E^{*}, I^{*}, v^{*})$, the Jacobian matrix of system (2.2) is

The characteristic equation of the matrix $J_{P^{*}}^{0}$ is

where

According to the Routh-Hurwitz stability criteria, whenever $\chi_{0} > 1$, the equilibrium point $P^{*}(E^{*}, I^{*}, v^{*})$ is locally asymptotically stable if $B_{11} > 0$, $B_{12} > 0$, $B_{13} > 0$ and $B_{11}B_{12}-B_{13} > 0$; otherwise it is unstable. However, clearly $B_{11} > 0$, $B_{12} > 0$, $B_{13} > 0$ and $B_{11}B_{12}-B_{13} > 0$ if $\chi_{0} > 1$. Thus, for the stability of the coexisting equilibrium point $P^{*}(E^{*}, I^{*}, v^{*})$, the system (2.2) must have the production rate of the virus greater than the ratio between the product of the decay rate of infected epithelial cells and virus concentration and the product of the infection rate of the virus and density of uninfected epithelial cells.

Case 2: $\tau\neq 0$

$({\rm{i}})$ The Jacobian matrix of system (2.2) at uninfected or virus-free equilibrium point $P^{\prime}(E(0), 0, 0)$ is

The one eigenvalue of the matrix $J_{P^{\prime}}$ is $\lambda_{P^{\prime}, 1} = -d_{1} < 0$, and the other two eigenvalues will be calculated from the following quadratic equation:

Let $\tau > 0$. Since Eq (4.5) is transcendental with an infinite number of solutions, we assume that $\lambda = \iota\omega$, and without loss of generality we may assume that $\omega = 0$ is a root of (4.5). Therefore, putting $\lambda = \iota\omega$ in (4.5), we get

Separating the real and imaginary parts of Eq (4.6), we have the following system:

Squaring both the equations of (4.7) and then adding them, we obtain

Since $\omega\neq 0$, Eq (4.8) is equivalent to

If $\chi_{0} < 1$, then

and

so that

Therefore, we have $\omega^{2} < 0$, which is a contradiction. Hence, whenever $\chi_{0} < 1$, according to Rouché's theorem ^[33,34], Eq (4.5) cannot have pure imaginary roots and the virus-free equilibrium $P^{\prime}(E(0), 0, 0)$ is locally asymptotically stable, for any strictly positive time-delay $\tau$.

Suppose $\chi_{0} > 1$. We already know that the one eigenvalue of the matrix $J_{P}$ is $\lambda_{P^{\prime}, 1} = -d_{1} < 0$, and the other two eigenvalues are the zeroes of the following expression:

As $\lambda_{P, 1} = -d_{1} < 0$, for stability, we need to check whether the other two eigenvalues have negative real parts or not. However, $\Gamma(0) = d_{2}d_{3}-\mu\beta E(0) < 0$ for $\chi_{0} > 1$ and also ${\lim_{\lambda \to +\infty}} \Gamma(\lambda) = +\infty$. Therefore, by continuity of $\Gamma(\lambda)$, there is at least one positive zero for the expression (4.10). Hence, we conclude that the virus-free equilibrium $P^{\prime}(E(0), 0, 0)$ is unstable when $\chi_{0} > 1$ i.e., also for the case of $\tau\neq0$, the virus-free equilibrium $P^{\prime}(E(0), 0, 0)$ is stable if $\chi_{0} < 1 \implies \mu < \frac{d_{2}d_{3}}{\beta E(0)}$.

$({\rm{ii}})$ At coexisting equilibrium point $P^{*}(E^{*}, I^{*}, v^{*})$, the Jacobian matrix of the system (2.2) is

The characteristics equation of the matrix $J_{P^{*}}$ is

where

Since Eq (4.11) is a transcendental equation with an infinite number of solutions, we cannot apply the classical Routh-Hurwitz criterion to Eq (4.11). To determine the stability of $P^{*}$, we consider $\lambda = \pm \iota\Omega$ (where, $\Omega > 0$), then Eq (4.11) can be written as

Separating real and imaginary parts, we get

Squaring both the equations of (4.12) and then adding them, we obtain,

For simplification, we can rewrite it as

Here,

It can be verified that $P_{11}^{2}-2P_{12}-Q_{11}^{2} = 28.6852 > 0$, and $P_{13}^{2}-Q_{13}^{2} = -0.010076 < 0$ for the parameter values prescribed in Table (1). Hence, there exists at least one nonnegative real root $\Omega_{0}$ for Eq (4.13). Therefore, we found a purely imaginary root $\pm \Omega_{0}\iota$ for Eq (4.11). So, the stability of system (2.2) may switch at $P^{*}$ as $\tau$ changes.

We eliminate $\text{sin}(\Omega\tau)$ from both the equations of (4.12) and obtain

Then, $\tau_{n}^{*}$ corresponding to $\Omega_{0}$ is given by

For $\tau = 0$, $P^{*}(E^{*}, I^{*}, v^{*})$ is stable. Hence, $P^{*}$ will remain stable for $\tau < \tau_{0}$, where $\tau_{0} = \tau_{0}^{*}$ and $n = 0$ ^[35].

Now, we will verify the transversality condition

for the occurrence Hopf bifurcation.

As Eq (4.11) has purely imaginary roots $\lambda_{0} = \iota\Omega_{0}$ and $\bar{\lambda}_{0} = -\iota\Omega_{0}$ at $\tau_{0}$

So, the implicit function theorem states that in a neighborhood of $\tau_{0}$, there must exist a complex function $\lambda = \lambda(\tau)$ such that $\lambda_{0}(\tau_{0}) = \lambda_{0}$ and $\Gamma(\lambda(\tau), \tau) = 0$, and $\frac{d\lambda}{d\tau} = -\frac{\partial \Gamma (\lambda(\tau), \tau)/\partial \tau}{\partial\Gamma(\lambda(\tau), \tau)/\partial\Phi}$ for $\tau$ in a neighborhood of $\tau_{0}$. Thus,

By using Eqs (4.11) and (4.12) in (4.14), we obtain the following expression,

Let $\lambda(\tau) = \text{Re}(\lambda(\tau)) + \iota\text{Im}(\lambda(\tau))$.

Then, from Eq (4.15), we have

By using the prescribed parameter values in Table (1), we get $\bigg[\frac{d\text{Re}(\lambda(\tau))}{d\tau}\bigg]_{\tau = \tau_{0}}^{-1} = 27.255 > 0$, which verifies the transversality condition for the Hopf bifurcation.

Therefore, at $\tau = \tau_{0}$ system (2.2) undergoes a Hopf bifurcation around the point $P^{*}$. It implies that an isolated periodic orbit emerges in the neighborhood $P^{*}$. Thus, $P^{*}$ is locally asymptotically stable for $\tau = 0$. Moreover, if there exists a threshold value of $\tau$ (say $\tau_{0}$) then $P^{*}$ is locally asymptotically stable for $0 < \tau < \tau_{0}$ and unstable for $\tau > \tau_{0}$. Also, a Hopf bifurcation occurs around the point $P^{*}$ when $\tau = \tau_{0}$.

Now, we will investigate the dynamics of appearing periodic solutions and compute the time lag to preserve the limit cycles. In order to do this, we consider model (2.2) that satisfies the initial conditions (2.3) on the interval $[-\tau, 0)$ and also the space of all continuous real-valued functions defined on $[-\tau, +\infty)$. Linearizing model (2.2) around $P^{*} (E^{*}, I^{*}, v^{*})$, we get

By using Laplace transformation of (4.16), leading to

with

where $L_{E} (\eta)$, $L_{I} (\eta)$, and $L_{v} (\eta)$ are the respective Laplace transformations of $E(t)$, $I(t)$, and $v(t)$. According to the well-known theory described in ^[35,36] and using classical Nyquist criteria stated in ^[37], the coexistence equilibria $P^{*}$ is asymptotically stable, for

with

and $\xi_{0} > 0$ is the smallest root of the expressions (4.18) and (4.19).

We can rewrite the expressions (4.18) and (4.19) as

which gives sufficient conditions for the stability of coexisting equilibrium $P^{*}$.

To calculate time lag $\tau$, we will determine an upper bound $\xi_{+}$ on $\xi_{0}$ that is independent from $\tau$. To do this, we assume that at $\xi = \xi_{0}$, $\forall$ values of $\xi$ satisfy Eq (4.21) such that $0\leq\xi \leq \xi_{+}$.

Expression (4.20) gives

Now, we maximize the right side of (4.22) as

subject to the conditions,

Therefore, we obtain that

Hence, it can be expressed as

From (4.23), it can be verified that $\xi_{0} \leq \xi_{+}$.

Also, (4.21), gives

Thus, if $\tau = 0$, then (4.21) becomes $\xi_{0}^{2} < P_{12}+Q_{12}$, and from (4.22), $P_{11}\xi_{0}^{2} = P_{13}+Q_{13}-Q_{11}\xi_{0}^{2}$; that is, $\xi_{0}^{2} = \frac{P_{13}+Q_{13}}{P_{11}+Q_{11}}$. Therefore, we confirm that at $\tau = 0$, the equilibrium point $P^{*}$ where all the populations exist is locally asymptotically stable if $(P_{13}+Q_{13}) < (P_{11}+Q_{11})(P_{12}+Q_{12})$ holds. However, for small $\tau > 0$, (4.24) must hold.

Putting (4.22) into (4.24) and rearranging the expressions, we get

Using the bounds, we obtain

and

From (4.25), we obtain that $\Psi_{1}\tau^{2}+\Psi_{2}\tau \leq \Psi_{3}$, with

Thus, we get the expression

then the Nyquist criteria ^[37] holds for $0\leq \tau \leq \tau_{+}$. This $\tau_{+}$ is the estimated maximum length of the time lag to preserve the stability of the limit cycle.

5.   Numerical simulations

This section demonstrates a few numerical examples to verify our analytical results. Also, we will visualize the role of time delay $\tau$ on the system's stability (2.2). In order to do a numerical study of system (2.2), we will use the values of the parameters given in Table (1) that are obtained in ^[22] using Chest radiograph score data. We would like to refer to ^[22] for an explanation of the range, values of the parameters, and their units.

First, we will verify that whatever the value of time delay $\tau$, the virus-free equilibrium point is locally asymptotically stable if $\chi_{0} < 1$, i.e., $\mu < \frac{d_{2}d_{3}}{\beta E(0)}$. In order to verify this, we take the fixed parameter set as: $d_{1} = 10^{-1}$, $\beta = 0.65$, $d_{2} = 0.11$, $\mu = 0.03 < \frac{d_{2}d_{3}}{\beta E(0)}$, $d_{3} = 5.36$ and the initial condition as: $E(0) = 22.41$, $I(0) = 2.59$, $v(0) = 0.061$. From Figure (1), it is clearly observed that in the presence or absence of time delay, the virus-free equilibrium point is locally asymptotically stable if $\mu < \frac{d_{2}d_{3}}{\beta E(0)}$. Thus, whatever the value of time delay in the process of antiviral responses of uninfected epithelial cells against the virus, if the virus production rate $\mu$ is lower than a threshold value $\frac{d_{2}d_{3}}{\beta E(0)}$, then the system itself is capable of clearing the virus.

Now, we will check the stability of system (2.2) at the coexisting equilibrium point for $\tau = 0$ and $\tau \neq 0$. For this, we take the fixed parameter set as: $d_{1} = 10^{-1}$, $\beta = 0.65$, $d_{2} = 0.11$, $\mu = 0.24 > \frac{d_{2}d_{3}}{\beta E(0)}$, $d_{3} = 5.36$ and the initial condition as: $E(0) = 22.41$, $I(0) = 2.59$, $v(0) = 0.061$. For the above parameter set and the initial condition, system (2.2) has two equilibrium points: $P^{\prime}(22.41, 0, 0)$ and $P^{*}(6.26356, 14.6786, 0.65725)$. In absence of time delay, i.e., $\tau = 0$, the eigenvalues associated with virus-free equilibrium $P^{\prime}(22.41, 0, 0)$ are $-5.958$, $0.4878$, and $-0.1$, which indicates that $P^{\prime}$ is an unstable saddle-node with no oscillating behavior. Hence, around the equilibrium point $P^{\prime}$, the system (2.2) shows unstable behavior. On the other hand, corresponding to the equilibrium $P^{*}(6.265, 14.68, 0.6578)$, the eigenvalues are $-5.431$, $-0.284$, $-0.114$, which indicates that $P^{*}$ is a stable node with no oscillating behavior. Thus, system (2.2) is stable around the equilibrium $P^{*}$. It is clearly observed from Figure (2) that all the populations exist over finite time with stable natures around the equilibrium $P^{*}(6.26356, 14.6786, 0.65725)$.

After checking the stability of system (2.2) in the absence of delay, we are now interested in checking the stability by varying the time delay parameter $\tau$. For $\tau = 1$, the system (2.2) shows similar stability behavior like the case $\tau = 0$ around the coexisting equilibrium $P^{*}(6.26356, 14.6786, 0.65725)$ (see Figure (3)). However, when $\tau = 12$, the system exhibits oscillatory behavior at first. It stabilizes as time increases (see Figure (4)), which indicates that as time delay increases, i.e., the time required for uninfected epithelial cells to respond to free virus increases, the cell populations oscillate first and are then able to stabilize themselves.

We further increase the value of $\tau$ to $12.7646$, at which point we have found that the system experiences a Hopf bifurcation, i.e., for this critical value of $\tau = 12.7646$ system stability is switched to instability via a limit cycle solution. Thus, we have observed a stable limit cycle around the coexisting equilibrium $P^{*}$ (see Figure (5)). Moreover, we numerically computed the values of $p_{1} = 28.6852$, $p_{2} = -0.4165$, and $p_{3} = -0.010076$ for which Eq (4.13) has a real root, $\Omega_{0} = 0.16535$ and thus $\tau_{0} = 12.7646$. Hence, Eq (4.13) has a purely imaginary eigenvalue. So, an isolated periodic solution, i.e., a limit cycle has appeared at $\tau_{0} = 12.7646$. The effect of $\tau$ on system stability is shown in Figure (8). Also, from Figure (8) it can be confirmed that system (2.2) bifurcates from one solution to two solutions (called maximum and minimum solutions) at $\tau_{0} = 12.7646$. Further, the system collapses for an approximate value of $\tau > 45.5$. A limit cycle solution of system (2.2), i.e., periodic oscillation of the population, has been observed for $\tau = 20 > 12.7646$ (see Figure (6)). For $\tau = 30 > 12.7646$, system (2.2), shows unstable behavior with an aperiodic nature (see Figure (7)).

Biologically, if the response time of uninfected epithelial cells toward the free virus increases to $12.7646$, then the cell population can control the free virus through competition. However, with further increases in the antiviral response time of uninfected epithelial cells to the virus, the cell populations are unable to stabilize the virus infection, and the body loses its stability, which also collapses as response time increases.

6.   Conclusions

In this study, we have developed a modified delay differential model to describe the within-host viral kinetics of SARS-CoV-2. The response of uninfected epithelial cells to the virus or the interaction of the virus with uninfected epithelial cells was incorporated as a nonlinear process into the model. Furthermore, a discrete-time delay has been introduced to describe the time required for uninfected epithelial cells to activate suitable antiviral responses by generating immune cytokines and chemokines. We have found that all the solutions to the formulated model are nonnegative and bounded. The basic reproduction number, $\chi_{0}$ was computed and found as $\frac{\mu \beta E(0)}{d_{2}d_{3}}$. In the absence of time delay, the virus-free equilibrium point $P^{\prime}$ was asymptotically stable for $\chi_{0} < 1$, and the co-existing equilibrium $P^{*}$ was stable for $\chi_{0} > 1$, provided that few certain conditions hold. However, in the presence of a time delay, the virus-free equilibrium point $P^{\prime}$ is locally asymptotically stable for $\chi_{0} < 1$, and the coexisting equilibrium $P^{*}$ loses its stability via a Hopf bifurcation. In addition, we have calculated the time lag over which the system maintains its stability around the co-existing equilibrium, $P^{*}$. The main finding of this study was that the considered system was stable for faster antiviral responses of epithelial cells to the virus, i.e., quick antiviral responses of epithelial cells stabilized patients' bodies by neutralizing the virus. However, if the antiviral response time of uninfected epithelial cells lengthens, the virus infects all the uninfected epithelial cells and may spread the disease over the body, which may require other treatment strategies to neutralize the virus's efficacy.

The main limitation of this study was that the model was straightforward and considered only two cell populations along with the virus concentration. However, while the SARS-CoV-2 virus interacted with the human body, other immune cells, such as natural killer cells, CD4$^{+}$ T cells and CD8$^{+}$ T cells, were also responsible for their antiviral responses. Thus, the model can be extended to explore the immune system's role by considering the effect of these immune cells ^{[15,20,38,39]}. To measure the vaccine or any other treatment effect that can be used in the model for eliminating the virus from the body ^[17,24]; this will be carried out in our future study.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Conflict of interest

The authors declare there is no conflict of interest.