Mathematical modelling of COVID-19 disease dynamics: Interaction between immune system and SARS-CoV-2 within host

1.   Introduction

Human civilization has been facing different infectious pathogen-caused epidemics for a long time. Severe Acute Respiratory Syndrome (SARS-CoV), a coronavirus class, was first identified in 2003 in Asia, with 8422 cases and an $11\%$ fatality rate ^[1]. Another class of coronavirus, Middle East Respiratory Syndrome (MERS-CoV) was reported with approximately 1572 human cases in 2012 in Arabian Peninsula ^[2]. SARS-CoV-2, the new class of novel coronavirus has now become a serious global threat due to its first advent in Wuhan, Hubei province, China in the year 2019 ^[3]. In spite of taking different social interventions and pharmaceutical measures by different countries, SARS-CoV-2 spreads over the globe within a very short time and tossed the biggest challenge ever to global health ^[4]. Every morning the world gets introduced with record-breaking cases in both infection and death caused by SARS-CoV-2. To pull down the race of the SARS-CoV-2 pandemic, researchers are attempting their best to discover strategies to explore the virus-protein structure of the SARS-CoV-2 ^[5], phylodynamics of the SARS-CoV-2 ^[6], and build up the antibody against SARS-CoV-2 infection by developing vaccines ^[7].

Mathematical models are often useful in characterizing the infection dynamics and forecast disease severity. Many pathogens were studied using a simple model including target cell, infected cell, and a virus. The target cell limited model was used by many researchers to study the in-host dynamics of HIV ^[8,9,10], Hepatitis ^[11], Ebola ^[12,13], and Influenza viruses ^[14,15]. Many mathematical models ^{[16,17,18,19]} focused on the transmission dynamics of the SARS-CoV-2 virus and involved different pharmacological, non-pharmacological interventions to reduce the spread of the COVID-19 disease. However, only a few models explained the replication cycle of SARS-CoV-2 within the host and interactions between this virus and the host immune system ^[20,21,22]. A brief study ^[22] on SARS-CoV-2 virus dynamics within-host discussed the target cell model with an eclipse phase, and secondary infection in presence of lymphocytes. The study found a rapid rise in viral load to reach the peak at the initial stage pointing to a plateau phase caused by lymphocytes. Also, an adaptive immune response results in decreasing the viral load at the last stage. Hernandez-Vargas et al. ^[21] investigated the interaction between the SARS-CoV-2 virus and the immune system within-host. Though the authors discussed a couple of models explaining the SARS-CoV-2 dynamic within the host, the immune cell (T-cell) response model performed better than others to fit the data indicating a slow immuno-response reaching the peak in 5–10 days after symptom onset.

Some researchers ^{[23,24,25,26,27]} studied the COVID-19 disease severity and mortality based on the data collected from the categorized patients. The viral load of SARS-CoV-2 in a group of patients with divergent SARS-CoV-2 serotypes was measured in the study of Fajnzylber et al. ^[23]. The authors stated that the spectrum of SARS-CoV-2 viremia is linked to inflammation markers and disease severity, such as elevated C-reactive protein (CRP), IL-6, and a low lymphocyte count. Through a virological study, Wolfel et al. ^[26] found sequence-distinct virus replication in the throat, lung tissue, and upper respiratory. The pulmonary cell affinity with the COVID-19 disease intensity was identified by computer tomography (CT-scan) and radio-graph of the patients' chest ^[25]. Later, Li et al. ^[24] used chest radiograph scores to determine the lung cell destruction due to the SARS-CoV-2 virus infection. They proposed a model describing the kinetics of the interactions between the virus particles and the lung epithelial cells within-host. Using the chest radiograph data, the authors determined the model parameters and found the basic reproduction number around 3.79 for the SARS-CoV-2 growth within the host.

Against the SARS-CoV-2 virus infection, no therapeutic and prophylactic medications are yet approved. Patients with acute symptoms are administered oxygen therapy, ventilation, antibiotics, and corticosteroids to support the respiratory functions and to control inflammatory reactions. These instructions allow time for activation of the adaptive immune response ^[28,29]. Different possible medical schemes include the use of antiviral drugs inhibiting the replications of the virus cycle. Repurposing existing anti-viral drugs can be an effective drug strategy against the SARS-CoV-2 that could reduce the time and limit the cost compared to the de novo drug discovery. In a study of drug re-purpose designed for SARS-CoV-2, Zhou et al. ^[30] proposed sixteen repurposable drugs based on network proximity analyses (NPA) of drug targets within the human host. The authors also proposed three prospective drug combinations such as mercaptopurine and melatonin, sirolimus and dactinomycin, toremifene and emodin, against the SARS-CoV-2 virus infection.

In any infections within the host, natural killer (NK) cells serve as first-line protection. NK cell and T-cell are necessary for a strong confrontation and clearance of the virus from the host. Among various anti-viral activities, these cells are accountable for directly killing the target cells and the release of immunomodulatory cytokines ^[31]. After conducting a thorough literature review, it is found that only a few research ^[21,22,24] have been carried out reporting the SARS-CoV-2 viral growth in the host under the impact of T-cells. However, no studies have been found which analyses the interaction between the immune system and SARS-CoV-2 within the host considering T -cells and NK-cells together.

A deterministic mathematical model based on physiologically plausible assumptions performs better in forecasting disease severity over time in the absence of relevant experimental data. A deterministic model reflecting the short-term expected antiviral effect during SARS-CoV-2 development within the host is required in addition to forecasting illness severity. As a result, developing a mathematical model under the effect of both T -cells and NK-cells becomes essential for qualitative analysis of the viral load within the host. To comprehend the SARS-CoV-2 growth dynamics within the host under immune surveillance, a mathematical model is therefore developed in the present study by analysing the interaction between the immune system and SARS-CoV-2 within the host integrating the function of NK cells and T-cells. In essence, the model looks into ways to control the pace of SARS-CoV-2 viral replication within a host, because a quicker viral replication rate inside the host is directly related to community transmission. This research will help physicians, scientists, and policymakers to understand how patients' immunity interact against COVID-19 disease, as well as to devise disease treatment strategies.

2.   Model formulation

The mathematical framework of the proposed model is based on the interaction between virus particles and different lymphocytes within-host including the most common lymphocytes, NK cells, and T-cells. T-cells are in control of both adaptive and innate immune responses while NK cells, which are part of the innate immune system, are the first immune effectors against virus infection before any specific immunity arrives. Furthermore, in the absence of specific antibodies, the NK cell recognizes and kills the distressed cell ^[32].

The target cell of the SARS-CoV-2 virus in the respiratory system is not completely known. Though there is a dispute in current literature about the primary compartments of SARS-CoV-2 infection, it is usually believed that the primary infection succeeds in respiratory epithelial cells ^[33]. Based on the results came out from the studies ^[34,35], epithelial cell in the respiratory system has been considered in the present model as a primary target of the SARS-CoV-2 virus. The infected epithelial cells are assumed to be capable of producing new virus particles at a constant rate. Taking the interactions of the virus-immune system and virus-epithelial cells, the proposed model is developed as follows:

where $E, E^i, V, N, L$ represent the number of susceptible epithelial cells, infected epithelial cells, viral load, natural killer cells, and T-lymphocytes respectively; $a_3$ is the infection rate of the virus; $c_1$ is the production rate; and the terms $a_2$, $b_1$, $c_2$, $d_2$, and $e_2$, are the natural death rates corresponding to $E, E^i, V, N,$ and $L$ respectively.

The regeneration of the epithelial cells is considered as $a_1 = a_2 \times E(0)$ ^[21], where $a_2$ is the natural death rate of the epithelial cells and $E(0)$ represents the number of initial epithelial cells. The term $a_3EV$ in the first equation denotes virus particles infecting a healthy epithelial cell at an infection rate of $a_3$. The growth of infected epithelial cells is denoted by $a_3EV$ in the second equation. The first term $c_1E^i$ in the third equation implies the reproduction of new SARS-CoV-2 virus from an infected cell at a rate $c_1$. Also, $c_3VN$ and $c_4V L$ represent the local interaction dynamics of the virus $(V)$ with natural killer cells $(N)$ and $T$ - lymphocytes $(L)$. Due to these interactions, virus particles are reduced at the rates $c_3$ and $c_4$ respectively. Similarly, $d_1$ implies the constant external source of the natural killer cells as expressed in the fourth equation of the model. The term $e_3$ is the natural recruitment rate of T-lymphocytes and $e_1LV / (n + V)$ is the proliferation in response to the presence of virus particles. This functional proliferation form is employed here as the model assumes the recruitment of $T$-lymphocytes due to the signals to release different cytokines (prompting molecule that regulates and mediates the immune system). The half-saturation constant $n$ is the$T$-lymphocytes proliferation by a virus which is estimated as $1.26 \times 10^5$ for a sigmoidal function to the best data fit ^[21].

The model is considered under the initial conditions $E(t) \geq 0, E i(t) \geq 0, V (t) \geq 0, N(t) \geq 0$ and $L(t) \geq 0$ with non-negative value of the model parameters. It is obvious that the right-hand parts of the system of Eq (2.1) are continuous and satisfy the Lipschitz condition ^[36]. The uniqueness and existence criterion of the ordinary differential equation confirms that model (2.1) has a unique solution with the non-negative initial conditions ^[36]. The first equation of the model (2.1), is considered to explore the non-negativity property of the solution curves as follows:

The solution trajectory of inequality (2.2) is obtained as

The inequality (2.3) gives $E(t) \geq 0$ for $t \rightarrow \infty$ which indicates that for every $t \in [0, \infty)$ the solution trajectory $E(t)$ will be entirely in the positive orthant for the given non-negative initial condition and parameter values. In a similar approach, it is easy to show $E^i(t)\geq 0, V(t)\geq 0, N(t) \geq 0, L(t) \geq 0$ for $t \rightarrow \infty$.

The first equation of the model (2.1) also reveals that

i.e., for long term behaviour $t\rightarrow \infty, E(t)$ is bounded by ${\frac{a_1}{a_2} = \Pi}_1(say)\in \mathbb{R}.$

Moreover, the solution trajectories of $E^i(t), V (t)$, and $L(t)$ of the system (2.1) give

where $p(0), q_1(0), q_2(0),$ and $r_1(0)$ represent the initial values of the corresponding functions in terms of state variables $E, V, N$ and $L$. Since $E$ and $V$ remain constant in the second equation of the model (2.1), the supremum of $E^i(t)$ is given by

Defining a function $F(t) = E(t) + Ei(t) + V (t) + N(t) + L(t)$ and its derivative along the trajectories yields,

where $H = \beta + \mu + c_1\Pi_2$, $h = \min \; \{a_2, b_1, c_2, d_2, e_2\}$, $\beta = (a_1 + d_1 + e_3)$, and $\mu = e_1LV /(n + V)$.

Inequality (2.6) reveals that $F(t)$ is bounded by $H/h$, whereas $E(t)$ and $E^i(t)$ have bounds $\Pi_1$ and $\Pi_2$ respectively. Therefore, $N(t), V(t)$ and $L(t)$ are eventually bounded with an other bound $\Pi \in \mathbb{R}$. Denoting $M = \max\; \{\Pi, \Pi_1, \Pi_2\}$, for sufficiently large time scale $t$, it reveals that $E(t) \leq M$, $E^i(t)\leq M$, $V (t)\leq M$, $N(t)\leq M$ and $L(t) \leq M$.

Therefore, under the non-negative initial conditions with the positive parameter values, the solution space of the system (2.1) gives

3.   Model analysis

Evaluation of the basic reproduction number $R_0$ at equilibrium points is the simplest way to investigate the disease dynamics described by a mathematical model. Setting the right-hand side of the model equations to zero, the first equilibrium point $\eta (E_0, E_0^i, V_0, N_0, L_0)$ in absence of the SARS-CoV-2 virus of the model (2.1) is obtained as

Further, at endemic equilibrium, none of the state variables is identically zero in presence of the virus within-host. In that case, the second equilibrium point $\eta^* (E_0^*, E_0^{i*}, V_0^*, N_0^*, L_0^*)$ for model (2.1) is given by:

The infection ability of a pathogen within-host is quantified by determining $R_0$ for the corresponding model ^[37]. This number is employed to evaluate the average secondary infections resulting from an infected cell. Mathematically, $R_0 < 1$ implies that an infected cell (pathogen) might infect less than one cell, resulting in an invasion of the disease from the host. In contrast, $R_0 > 1$ suggests a suppression of the host cells resulting in the progress of the disease within-host.

The concept of the next-generation matrix $(NGM)$ ^[38] is used to determine the basic reproduction number $R_0$. The system (2.1) has three infected states $E, E^i$ and $V$, and two uninfected states, $N$ and $L$. Linearising the Eqs (1)–(3) of the system (2.1) as a subsystem of infection, the following transmission matrix $(T)$ and transition matrix $(Q)$ are obtained:

where $\alpha = c_2+\dfrac{c_3d_1}{d_2}+\dfrac{c_4e_3}{e_2}$.

Since the first two columns of the transmission matrix $(T)$ are zero, the next generation matrix ^[38] will be $K = E^TTQ^{-1}E$ where the axillary matrix $E$ is given by:

Using the next-generation matrix $K$, the $R_0$ of the model (2.1) can be written as

3.1.   Steady state analysis

The stability of model (2.1) at the virus-free equilibrium point $\eta$, is determined by the following lemma:

Lemma 1. The disease-free equilibrium of model (2.1) is locally asymptotically stable for $R_0 < 1$ and unstable for $R_0 > 1$.

Proof. The Jacobian matrix of the model at disease-free equilibrium point $\eta$ can be written as

The characteristic equation of this Jacobean matrix is given by

where $k_1 = (b_1 + \alpha) > 0$ and $k_2 = b_1\alpha - (a_1a_3c_1)/(a_2)$.

It is clear that the first three eigen-values obtained from Eq (3.4) are negative as $a_2 > 0$, $d_2 > 0$ and $e_2 > 0$. As per the Descarte's sign rule ^[39], the rest of the eigen-values are negative for $k_1 > 0$ and $k_2 > 0$. Rewriting the condition $k_2 > 0$ in $R_0$, it gives $k_2 = 1 - R_0 > 0$, that means, $k_2 > 0$ if only if $R_0 < 1$. Thus all of the eigenvalues are negative. On the contrary, $k_2$ becomes negative if $R_0 > 1$ which provides at least one of the eigenvalues is positive. Hence, Lemma 1 follows.

The analytical computation of the characteristic equation at endemic equilibrium point $\eta^\ast$ produced a large expression and therefore, the determination of the sign of eigenvalues at $\eta^\ast$ becomes very complicated. Rather, the qualitative behaviour of the solution trajectories at $\eta^\ast$ was investigated through a series of numerical explorations.

4.   Numerical simulation and discussion

A series of numerical simulations were performed using Python (V. 3.7) SciPy.integrate package and odeint library to solve the ordinary non-linear differential equations of the present model. The model equations comprised 14 parameters and assumed 5 initial conditions. The parameter values were estimated based on some homeostatic process, and available statistical data of SARS-CoV-2 infected patients ^[21,24]. Considering the physical significance of the parameters, numerical investigations corresponding to the effects of external NK cell influx, anti-viral therapy, and their combined effect over the SARS-CoV-2 growth projection have been analyzed under immune surveillance in section 4.1–4.3 respectively.

The initial values for healthy and infected epithelial cells, SARS-CoV-2 virus, and T-cell are based on COVID-19 patient data analysis ^[21,24], which may provide realistic results for numerical simulations. The model assumes the initial conditions as follows: $E^i(0) = 2.59$ (score) ^[24], $E(0) = 25 - E^i(0) = 22.41$ ^[24]), $V (0) = 0.061$ (score) ^[24], $L(0) = 106$ ^[21], $N(0) = 105$ ^[40]. The virus clearance rate $c_2$ was fixed and characterized as a process that is not linked with the immune system. The parameter $a_1$ was estimated as ^[24]: $a_1 = a_2 \times E(0) = 0.02241$, and the other simulation parameters are tabulated in Table 1. It is noted that the parameter values for epithelial cells and SARS-CoV-2 were chosen from published literature based on chest radiograph score which is an approximation, not actual infected data for pulmonary cells ^[24].

4.1.   External influx of natural killer (NK) cell

The impact of NK cell influx from an external source on the disease progression is assessed to quantify the viral dynamics within-host. The immunity system works to clear the virus from the day of symptom onset. However, the immunity strength depends on different factors including age and co-morbidness of the patients. Figure 1 shows the effect of external NK cell influx over viral load. As seen in the figure, the viral load increase very fast at the beginning of the infection due to slow immune response and reaches the peak of 0.35 (chest score) on day 10 after symptom onset where the viral load is expressed as a chest radiograph score. A small increase in the external influx of NK cell as $d_1 = d_1 + j \ast 1000\; (0\leq j \leq 5)$ results in a reduction of viral load about 0.3 to $1.5\%$ in each loop iteration. The maximum viral load (chest score) for the corresponding loop iterations are $v = 0.3501, v = 0.3492, v = 0.3475, v = 0.3450,$ and $v = 0.3417$ respectively. Because the viral loads corresponding to $d = 3200, d = 4200,$ and $d = 6200$ are so similar, Figure 1 contains plots corresponding to d = 3200, d = 9200, and d = 13200, which illustrate the intensity of viral load relative to the change in external NK cell influx.

Figure 1(d) shows the basic reproduction number $R_0$ which gradually decreases subject to the external influx of the NK cell. The large influx of NK cells strengthens the immune system, helps to restrict viral proliferation. Diminishing $R_0$ from 3.65 to 3.33 (approx.) implies the invasion of viral replication within the host during the disease. A simulation was also run with higher initial values of the SARS-CoV-2 virus to investigate the initial viral load over the disease progression. A significant rise was observed in the chest radiograph score (viral load) which indicates a symptom severity in patients.

4.2.   Anti-viral therapy

The role of possible anti-viral therapy can reduce the SARS-CoV-2 virus replication rate. A high replication rate $c_1$ of the SARS-CoV-2 implies a rapid invasion of the susceptible epithelial lung cells $(E)$. The SARS-CoV-2 replication rate on viral load dynamics is portrayed in Figure 2. The maximum viral load of 0.35 (score) occurs in Figure 2(a), while the minimum viral load of 0.003 (score) ensures in Figure 2(e). For each loop decrement $(c_1 = c_1-j \ast 0.02)$ in the value of $c_1$, viral load was decreased by approximately 14 to $31\%$. For a lower replication rate, SARS-CoV-2 growth dynamics are found less steep than for a higher replication rate (Figure 2), which indicates that an anti-viral therapy can more effectively combat the virus. In addition, the SARS-CoV-2 virus cannot grow within-host under immune surveillance if the replication rate is less than or equal to 0.04 as seen in Figure 2(e). The threshold value of $c_1$ was found 0.04 for the defined loop iterations. The actual threshold value of $c_1$ for model (2.1) was determined as 0.09 implying that SARS-CoV-2 cannot grow in the host with a replication rate less than 0.09.

Figure 2(f) depicts the graph of basic reproduction number $R_0$ relative to the SARS-CoV-2 virus replication rate. The number $R_0$ shows a positive linear trend with virus replication rate $c_1$. A higher replication rate clearly contributes to an increased symptom and disease prevalence in patients, resulting in community disease transmission and lung cell destruction. This finding suggests that prescribing appropriate antiviral drugs to a patient sooner may reduce the symptom.

4.3.   Combination of anti-viral therapy and influx of external NK cell

Combinations of different therapeutic treatments have become popular in many viral disease management protocols. For instance, monoclonal antibodies with other anti-viral medicine showed an important role in Ebola virus infection ^[12]. In this section, the combined impact of anti-viral therapy and the eternal influx of NK over the SARS-CoV-2 virus growth projection has been demonstrated. The combined effect of possible anti-viral drugs and influx of external NK cells is shown in Figure 3. The influx of NK cells strengthens the patient's immune system, while antiviral therapy reduces virus replication. It is found that administering this combined therapy results in a quick reduction of viral load from the initial peak. Such a substantial reduction in viral load causes symptom relief and disease invasion in the patient.

The basic reproduction number plot Figure 3(f) shows a sharp fall due to the combined effect of NK cell influx and anti-viral drug. This reflection implies that the combined effect is more effective in suppressing the viral load within the host under immune surveillance.

5.   Sensitivity analysis

The sensitivity of a specific parameter identifies how a small change in numeric value causes a large impact on the model output. This analysis may assist to establish the possible treatment protocol for the SARS-CoV-2 disease. The sensitivity indices of the basic reproduction number $R_0$ relative to the parameters is defined as ^[12]:

where $r$ is the model parameter. The estimated sensitivity indices of the model (2.1) are summarized in Table 2. From the sensitivity indices, it is obvious that $a_1, a_3, c_1, d_2$, and $e_2$ have a positive impact on $R_0$ while the rest of the parameters shows a negative impact. The negative impact implies the invasion of viral load whereas the positive impact indicates a progression of virus load within-host. The indices $\Theta_{a_1}, \Theta_{a_3}$, and $\Theta_{c_1}$ are being equal to 1 (most sensitive) reveals that a $10\%$ increase in the corresponding parameter values causes a $10\%$ growth projection in $R_0$. Similarly, $a_2$ and $b_1$ show reverse correlation to $R_0$, i.e., $10\%$ increase in $a_2$ and $b_1$ (least sensitive) results in decreasing the value of $R_0$ by $10\%$. These results fairly agree with the numerical results (Figures 1 and 2). The constant regeneration $a_1$ shows positive sensitivity to $R_0$ as in Table 2. It means that a host with higher cells of pulmonary epithelial is more likely to become infected, and provides clarification why especially babies having lower epithelial cells are not most expected to get infected.

6.   Conclusions

A mathematical model was developed to explore the interaction between the immune system and SARS-CoV-2 within the host by integrating the function of NK cells and T-cells. A series of simulations were performed to investigate the qualitative behaviour of the model. The numerical results complied with the analytical results of the model where the solution trajectories were found bounded in the positive quadrant. A sensitivity analysis was also carried out for basic reproduction number $R_0$ to analyse the impact of the parameters on the disease progression. The model is found to be advantageous in assessing the antiviral medication efficacy and NK cell influx into the host body, and with the value of $R_0$ across the time of infection period, the model can be used in forecasting disease severity. The numerical analysis showed that the use of anti-viral therapy and influx of NK cells causes a fall in the value of $R_0$. In addition, SARS-CoV-2 cannot grow within-host if the virus replication rate $(c_1)$ is under a threshold limit of 0.09. The combined use of anti-viral therapy and influx of external NK cells together played an important role to control the SARS-CoV-2 viral growth within-host. If the proposed anti-viral combinations are capable to control or reduce the SARS-CoV-2 replication rate, it may be a good attempt to use those from the very beginning of the disease. The present model can be used further to investigate the role of antibody development and cross-reaction. Furthermore, the present model can be extended to study the innate immune response to different SARS-CoV-2 mutations. Each new variety must have its own replication rate, which has an impact on the severity of the disease and its propagation in the population. This idea would help us to develop a new mathematical model in the future study.

Acknowledgements

The authors extend their appreciation to the Deanship of Scientific Research at King Khalid University for funding this work through the research groups program under grant number RGP.1/327/42.

Conflict of interest

The authors declare that they have no conflict of interest.