Spatio-Temporal Modeling of Immune Response to SARS-CoV-2 Infection

Abstract

COVID-19 is a disease occurring as a result of infection by a novel coronavirus called SARS-CoV-2. Since the WHO announced COVID-19 as a global pandemic, mathematical works have taken place to simulate infection scenarios at different scales even though the majority of these models only consider the temporal dynamics of SARS-COV-2. In this paper, we present a new spatio-temporal within-host mathematical model of COVID-19, accounting for the coupled dynamics of healthy cells, infected cells, SARS-CoV-2 molecules, chemokine concentration, effector T cells, regulatory T cells, B-lymphocytes cells and antibodies. We develop a computational framework involving discretisation schemes for diffusion and chemotaxis terms using central differences and midpoint approximations within two dimensional space combined with a predict–evaluate–correct mode for time marching. Then, we numerically investigate the model performance using a list of values simulating the baseline scenario for viral infection at a cellular scale. Moreover, we explore the model sensitivity via applying certain conditions to observe the model validity in a comparison with clinical outcomes collected from recent studies. In this computational investigation, we have a numerical range of 10⁴ to 10⁸ for the viral load peak, which is equivalent to what has been obtained from throat swab samples for many patients.

1. Introduction

In December 2019, the coronavirus disease 2019 (COVID-19) was first identified in Wuhan and rapidly invaded the whole of China and the world [1]. This disease was officially named COVID-19 by the World Health Organization (WHO) on 11 January 2020; then it was announced as a global health emergency on 11 March 2020 [2,3]. Since then, many studies have been done to investigate this global pandemic. As a part of these efforts, we here aim to develop a mathematical model to describe and simulate the interactions between SARS-CoV-2 molecules and the immune system within a two dimensional cellular space. In the following, we present a brief overview of SARS-CoV-2 in terms of its spreading properties, viral replication and mechanism of infection.

Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) is a novel virus that has a genetic similarity to SARS-CoV, which is responsible for SARS disease, while the new virus (SARS-CoV-2) causes coronavirus disease 19 (COVID-19) [4]. Both SARS-CoV and SARS-CoV-2 have a common link in terms of mechanism of infection where both of them use the spike proteins (S-protein) to bind Angiotensin-converting enzyme 2 (ACE2), which in turn helps them to reach and attack the cells [4]. Despite this similarity, it has been found that SARS-CoV-2 binds to ACE2 slower than the old virus (SARS-CoV) where this slow binding causes a longer incubation period. Hence, this time delay of onset symptoms of COVID-19 leads to a faster spread of the pandemic around the world [3].

Although there is a lack of enough experimental studies about immune responses for people who have symptomatic/asymptomatic COVID-19 infections, there are some studies which investigate the interactions between the virus and many components of the immune system including the innate and the adaptive immune responses [5,6,7]. Once coronavirus-2 gains access inside the host cell using its spike protein (S) and the target cell receptor (ACE2) (in the same manner as SARS-CoV), it begins uncovering its genetic materials and producing some essential proteins such as envelope (E), nucleocapsid (N) and also spike (S) for motivating the interactions with target cells as well as assembling and releasing a new virus generation [8,9]. The ability of viral replication drives to having a local inflammatory response which motivates the immune system to recognize the virus and start the primary line of defense, that is, the activation of neutrophils, natural killer cells, naive T and B cells and the secretion of chemokine [10,11]. As long as the infection persists for a certain period of time, T cells pass through many stages including their proliferation, activation, migration, and differentiation, ultimately they come in the form of the effector T cells involving both helper (CD4+) and cytotoxic T cells (CD8+), and the regulatory T cells (Treg) [10]. Besides the cellular immunity, B cells play a key role in humoral immunity, especially after their transformation into plasma B cells (via the assistance of the armed/follicular helper T cells [10,12,13]). These types of 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 [4,8,14]. Antibodies specific to SARS-CoV-2 are mostly a very essential tool for neutralizing the virus via decreasing its binding ability to ACE2, but we should take into account its potential threat of enhancing the immunopathogenesis [15]. However, there is biological evidence that shows that the antibody levels are going to be decreased over time through infection with SARS-CoV-2, while particular components of the immune system, such as the memory T/B cells, will maintain at a high level for almost 6 to 8 months [6].

One major obstacle that creates difficulties for dealing with SARS-CoV-2 infection is the quick activation and migration of T cells to the site of viral replication, which in turn leads to releasing a huge amount of inflammatory cytokines such as IL-1, IL-6, and IFN-$\gamma$. This storm of cytokines causes the cytokine release syndrome resulting in serious damage to vital organs such as lungs as well as causing viral sepsis [11,16]. On other hand, there is no clue so far suggesting that the evolution of B cells over time will participate in this pathogenesis. However, we cannot generalize this scenario because it has been found that some of infected patients with SARS-CoV-2 who were treated in the ICU or who died have low levels of T cells as well as B cells [7].

Over the past decades, mathematical modeling of immune responses to viral infections has developed in terms of either realizing the biological concepts or improving the mathematical techniques to obtain better outputs [17,18,19,20]. Since the outbreak of COVID-19, many mathematical research projects have been performed, but the majority of these mathematical models focus on the temporal dynamics only (ODE models) which either use the so-called SIR-models (compartmental models) to predict the population progress between compartments such as [21,22,23,24,25], or describe the within host dynamics to investigate immune responses to SARS-CoV-2 particles within the cellular scale as proposed in [22,26,27]. However, these within host models were being built on previous studies for similar infectious diseases such as influenza [18,28] and hepatitis B [29]. To our knowledge, we could not find any spatio-temporal model for describing the spatial distribution of the immune responses to SARS-CoV-2 within a two dimensional cellular space. Nevertheless, there are a few preprint mathematical studies using PDEs to describe the spread of the population compartments with respect to time and their physical locations at different regions which means these mathematical approaches are still a type of the compartmental class models. Therefore, in Section 2, we present a first attempt for deriving a new mathematical model, which is based on a coupled system of partial differential equations to obtain a reasonable description for the interaction between coronavirus-2 and its response to the components of immune system.

The rest of this paper is summarized as follows. Section 3 includes model calibration and developing a computational framework to compute the coupled reaction diffusion system within two spatial dimensions. Section 4 involves a comparison of our numerical results with a few cases collected from recent clinical studies. Lastly, in Section 5 we present a discussion on a few key points including hints at future works.

2. Modeling Framework

Building on some biological observations and mathematical models (temporal approaches) [7,8,10,11,12,13,16,22,28,30,31,32,33,34,35,36], we develop a new spatio-temporal mathematical model to explore the interactions between the most essential components of the immune system and SARS-CoV-2 particles within a cellular perspective. This dynamical interaction occurs on a spatial domain $\aleph \in {\mathbb{R}}^2$ over a time interval $[0,T]$ and it involves the healthy host cells $h(t,x)$, viral-infected cells $i(t,x)$, SARS-CoV-2 molecules $v(t,x)$, chemokines $c(t,x)$, the effector T cells $f(t,x)$, the regulatory T cells $r(t,x)$, B-lymphocytes cells $\beta (t,x)$ and antibodies $\alpha (t,x)$. In a similar manner to the modeling derivation method under a continuum approach and the conservation law as detailed in [37], we assume that $\mathcal{L}$ is a fixed (random) volume surrounded entirely by a smooth surface $\mathcal{Z}$. Then, the appearance or disappearance of any component within $\mathcal{L}$ is expressed by the flux of that component through the boundary $\mathcal{Z}$ and the proliferation and/or degradation of that component. In other words, we suppose that $\mathcal{C}$ represents the concentration of a component at a fixed time t and a position $\mathrm{x}=(x,y)$, $\mathcal{J}$ denotes the net flow of a component and $\mathcal{P}$ stands for the proliferation and/or the degradation of a component. Hence, we have the following conservation equation:

$$\begin{array}{c}\hfill \frac{\mathrm{d}}{\mathrm{dt}}{\int }_{\mathcal{L}}^{}\mathcal{C}(\mathrm{x},t)\mathrm{dx}={\int }_{\mathcal{Z}}^{}-\mathcal{J}(\mathrm{x},t)\mathrm{d}\mathcal{Z}+{\int }_{\mathcal{L}}^{}\mathcal{P}\mathrm{d}\mathcal{L},\end{array}$$

where $\mathcal{Z}=\partial \mathcal{L}$. Now, using the divergence theorem, we can differentiate to get the following conservation equation for each component of a system, namely,

$$\begin{array}{c}\hfill \frac{\partial \mathcal{C}}{\partial t}=-\nabla ·\mathcal{J}+\mathcal{P}.\end{array}$$

The above method has been used to model many biological applications such as [34,38,39,40]. Thus, we can denote the resulting 8D-reaction-diffusion-chemotaxis operator $\mathcal{K}$ by the following expression:

$$\mathcal{K}{(h,i,v,c,f,r,\beta ,\alpha )}^T=0,$$

(1)

where the differential operator $\left(\mathcal{K}\right)$ of order 2 will be explained in detail in the following context.

Host $h(t,x)$ & Infected cells $i(t,x)$ Following biological observations and spatio-temporal mathematical approach that have been taken from various experimental studies for similar infectious diseases such as [31,32], we suppose that both type of cells are exercising a random motility on the spatial domain ℵ at rates $D_{{}_h}>0$ and $D_{{}_i}>0$, respectively. The density of the host cells is assumed to be produced logistically at a rate ${\mu }_{{}_h}$, while the decay of its density comes from many reasons at a rate ${\delta }_{{}_{11}}$. At the same time, the infected cells are produced when SARS-CoV-2 virus binds the host cells by the S protein at a rate ${\mu }_{{}_i}$ through the targeted receptor, namely, ACE2 [11]. Moreover, the density of infected cells is decreased over time due to two main factors, namely, effector T cells $f(t,x)$ at a rate ${\delta }_{{}_{21}}$[10,33] and natural deaths at a rate ${\delta }_{{}_{22}}$. These assumptions can be formalised mathematically through the following expressions:

$$\left\{\begin{array}{c}\frac{\partial h}{\partial t}=D_{{}_h}{\nabla }^{{}^2}h+{\mu }_{{}_h}h(1-\frac{h}{h_{{}_0}})-{\delta }_{{}_{11}}h,\hfill \\ \frac{\partial i}{\partial t}=D_{{}_i}{\nabla }^{{}^2}i+{\mu }_{{}_i}hv-{\delta }_{{}_{21}}if-{\delta }_{{}_{22}}i.\hfill \end{array}\right.$$

(2)

SARS-CoV-2$v(t,x)$ As mentioned in [8], SARS-CoV-2 has the ability to motile inside the respiratory system, therefore, we assume that coronavirus-2 diffuses randomly through ℵ at a random motility rate $D_{{}_v}$. Furthermore, the interactions between immune cells and viruses increase the secretion level of chemokines, which in turn impacts virus movement [8]. Adopting this assumption, we suppose that SARS-CoV-2 migration is directed towards a higher concentration of a chemical gradient at a rate ${\chi }_{{}_v}$. Moreover, coronaviruses are capable of viral replication actively within throat for at least 5 days (depending on the duration of the onset of symptoms) as detailed in the clinical studies [41,42], where this viral replication occurs due to the existence of the gene (ORF1) [36]. Hence, we take into account its replication within infected cells at a rate ${\mu }_{{}_v}$. However, besides the natural decay of SARS-CoV-2 density at a rate ${\delta }_{{}_{32}}$, the antibodies (produced by B cells) play a key role in restricting or preventing virus entry inside the host cells [10], hence SARS-CoV-2 particles are assumed to be decreased due to antibodies increasing level at a rate ${\delta }_{{}_{31}}$. These words can be mathematically expressed via the following equation:

$$\left\{\begin{array}{c}\frac{\partial v}{\partial t}=D_{{}_v}{\nabla }^{{}^2}v-{\chi }_{{}_v}\nabla ·(v\nabla c)+{\mu }_{{}_v}i-{\delta }_{{}_{31}}\alpha v-{\delta }_{{}_{32}}v.\hfill \end{array}\right.$$

(3)

Chemokine Concentration $c(t,x)$ Chemokines are small proteins secreted by many cells where they have the ability to direct the surrounding cells chemotaxically towards the source of the chemokines [43]. Besides the homeostatic function of chemokines [43], chemokines have been found in many COVID-19 patients [8,10] for attracting immune cells towards the inflammatory sites [8]. Thus, it is assumed to be produced proportionally to the viral-infected cells and T cells at a rate ${\mu }_{{}_c}$, with taking into consideration its natural decay at a rate ${\delta }_{{}_{41}}$. Moreover, we suppose that chemokines can move randomly through the tissue at a random motility coefficient rate $D_{{}_c}$. Therefore, we have the following equation:

$$\left\{\begin{array}{c}\frac{\partial c}{\partial t}=D_{{}_c}{\nabla }^{{}^2}c+{\mu }_{{}_c}if-{\delta }_{{}_{41}}c.\hfill \end{array}\right.$$

(4)

The effector $f(t,x)$ and the regulatory T cells $r(t,x)$ In general, T cells are one of the most important components of immune system un which they play a key role in the initiation of the adaptive immune responses during the immunity cycle [10]. Regardless of the complexity of T cells functions and multi-stage processes including naive, helper, cytotoxic and memory T cells, it ultimately turns into the effector T cells $\left(f\right(t,x\left)\right)$[10]. Additionally, we also account for the dynamics of the regulatory T cells $\left(r\right(t,x\left)\right)$ due to its significant role in immune homeostasis and in controlling the proliferation of effector T cells [30]. Biological evidence shows that the efficiency of T cells depends on its motility pattern, spatial distribution and motility major obstacles [33]. Therefore, we take into account the random movement on the spatial domain for both T cells (effector and regulatory) with motility rates $D_{{}_f}$ and $D_{{}_r}$, respectively [7,34]. In addition to that, the directional migration of T cells is assumed to be in the form of chemotaxis towards gradients of the secreted chemokines (with the chemotaxic rates ${\chi }_{{}_f},{\chi }_{{}_r}$) [7,8,33,34]). Building on relevant studies [28,35], this specific study [22] proposes a model describing the proliferation form of T cells in the context of COVID-19. Thus, we assume that the proliferation of the effector T cells arises naturally at a rate ${\mu }_{{}_f}$ and its growth depends on the virus distribution with taking into consideration the maximum carrying capacity $k_{{}_f}$. Furthermore, the proliferation rate of the regulatory T cells keeps track of the production of the effector T cells at a rate ${\mu }_{{}_r}$ with a carrying capacity rate $k_{{}_r}$. However, both types of T cells are supposed to be decreased naturally at rates ${\delta }_{{}_{52}}$ and ${\delta }_{{}_{61}}$, respectively. Besides this, the regulatory T cells regulate the immune system and also suppress the proliferation of effector T cells before damage occurs [30]. Consequently, this interaction between the two types of T cells leads to a decay in the overall density of the effector T cells at a rate ${\delta }_{{}_{51}}$. Hence, the evolution of the effector and the regulatory T cells are represented by the following equations:

$$\left\{\begin{array}{c}\frac{\partial f}{\partial t}=D_{{}_f}{\nabla }^{{}^2}f-{\chi }_{{}_f}\nabla ·(f\nabla c)+{\mu }_{{}_f}\frac{v}{v+k_{{}_f}}-{\delta }_{{}_{51}}rf-{\delta }_{{}_{52}}f,\hfill \\ \frac{\partial r}{\partial t}=D_{{}_r}{\nabla }^{{}^2}r-{\chi }_{{}_r}\nabla ·(r\nabla c)+{\mu }_{{}_r}\frac{f}{f+k_{{}_r}}-{\delta }_{{}_{61}}r.\hfill \end{array}\right.$$

(5)

B-lymphocytes cells $\beta (t,x)$ & Antibodies $\alpha (t,x)$ During the early stages of infection with the SARS-CoV-2 virus, T cells (armed/follicular helper T cells) and their released cytokines during activation play an important role in the activation process of B cells [12,13,16] including the early responses of B cells against the N protein [8]. Therefore, we suppose that the production level of T cells increases the activation level of B cells at a rate ${\mu }_{{}_{\tilde{\beta }}}$. Moreover, the proliferation of B cells is caused via virus growth tracking at a rate ${\mu }_{{}_{\beta }}$ with a maximum carrying capacity $k_{{}_{\beta }}$. Meanwhile, its density decay is naturally occurring at a rate ${\delta }_{{}_{71}}$. Furthermore, we account for the produced antibodies by B cells at a rate ${\mu }_{{}_{\alpha }}$[7,10]. However, the density of antibodies is assumed to be degraded due to inhibition induced via binding with SARS-CoV-2 as well as by natural reasons at a rate ${\delta }_{{}_{81}}$[7,10]. Both B cells and antibodies are assumed to exercise a random motility on the spatial domain at diffusion coefficients $D_{{}_{\beta }}$ and $D_{{}_{\alpha }}$, respectively. We also take into account the chemotaxis directional movement of B cells towards inflammatory concentration [8,10] at a rate ${\chi }_{{}_{\beta }}$. The above assumptions are formalised mathematically by the following equations:

$$\left\{\begin{array}{c}\frac{\partial \beta }{\partial t}=D_{{}_{\beta }}{\nabla }^{{}^2}\beta -{\chi }_{{}_{\beta }}\nabla ·(\beta \nabla c)+{\mu }_{{}_{\beta }}\frac{v}{v+k_{{}_{\beta }}}+{\mu }_{{}_{\tilde{\beta }}}f-{\delta }_{{}_{71}}\beta ,\hfill \\ \frac{\partial \alpha }{\partial t}=D_{{}_{\alpha }}{\nabla }^{{}^2}\alpha +{\mu }_{{}_{\alpha }}\beta -{\delta }_{{}_{81}}\alpha .\hfill \end{array}\right.$$

(6)

Considering the above modeling hypotheses, the full model is written as follows:

$$\left\{\begin{array}{c}\frac{\partial h}{\partial t}=D_{{}_h}{\nabla }^{{}^2}h+{\mu }_{{}_h}h(1-\frac{h}{h_{{}_0}})-{\delta }_{{}_{11}}h,\hfill \\ \frac{\partial i}{\partial t}=D_{{}_i}{\nabla }^{{}^2}i+{\mu }_{{}_i}hv-{\delta }_{{}_{21}}if-{\delta }_{{}_{22}}i,\hfill \\ \frac{\partial v}{\partial t}=D_{{}_v}{\nabla }^{{}^2}v-{\chi }_{{}_v}\nabla ·(v\nabla c)+{\mu }_{{}_v}i-{\delta }_{{}_{31}}\alpha v-{\delta }_{{}_{32}}v,\hfill \\ \frac{\partial c}{\partial t}=D_{{}_c}{\nabla }^{{}^2}c+{\mu }_{{}_c}if-{\delta }_{{}_{41}}c,\hfill \\ \frac{\partial f}{\partial t}=D_{{}_f}{\nabla }^{{}^2}f-{\chi }_{{}_f}\nabla ·(f\nabla c)+{\mu }_{{}_f}\frac{v}{v+k_{{}_f}}-{\delta }_{{}_{51}}rf-{\delta }_{{}_{52}}f,\hfill \\ \frac{\partial r}{\partial t}=D_{{}_r}{\nabla }^{{}^2}r-{\chi }_{{}_r}\nabla ·(r\nabla c)+{\mu }_{{}_r}\frac{f}{f+k_{{}_r}}-{\delta }_{{}_{61}}r,\hfill \\ \frac{\partial \beta }{\partial t}=D_{{}_{\beta }}{\nabla }^{{}^2}\beta -{\chi }_{{}_{\beta }}\nabla ·(\beta \nabla c)+{\mu }_{{}_{\beta }}\frac{v}{v+k_{{}_{\beta }}}+{\mu }_{{}_{\tilde{\beta }}}f-{\delta }_{{}_{71}}\beta ,\hfill \\ \frac{\partial \alpha }{\partial t}=D_{{}_{\alpha }}{\nabla }^{{}^2}\alpha +{\mu }_{{}_{\alpha }}\beta -{\delta }_{{}_{81}}\alpha .\hfill \end{array}\right.$$

(7)

Each of the above coupled dynamics are followed by initial and boundary conditions which will be given in detail in Section 3.2.

3. Model Calibration and Computational Approach

In this section, we aim to prepare our dynamical system to be solved numerically and to provide a realistic simulation in the hope that it will coincide with a range of clinical findings. To this done, we first rewrite model (7) in a dimensionless form, then we select the initial seeds for system components involving shapes and quantities. Then, we calibrate the model and develop a numerical scheme to represent a scenario for immune responses to SARS-CoV-2 infection based on experimental data.

3.1. Non-Dimensionalization

To explore the performance of our proposed spatio-temporal model (7), we begin with non-dimensionalising the system in order to remove the physical dimensions of the variables to solve the system numerically. To that end, we rescale the space variable $x:=(x,y)$ to be related to the length of the considered region, namely, $L=0.1\mathrm{cm}$ in parallel with rescaling time with a reasonable scaling parameter $\tau =\frac{L^2}{D}$ where $D={10}^{-6}\frac{c{m}^2}{s}$ is the spatial chemical diffusion rate as referenced in [44]. Based on that, we have $\tau ={10}^4$ s which is equivalent to about 2.78 h. Now, we rewrite the dependent variables as follows: $\overline{h}(\overline{t},\overline{x},\overline{y})=\frac{h(t,x,y)}{h_0}$, $\overline{i}(\overline{t},\overline{x},\overline{y})=\frac{i(t,x,y)}{i_0}$, $\overline{v}(\overline{t},\overline{x},\overline{y})=\frac{v(t,x,y)}{v_0}$, $\overline{c}(\overline{t},\overline{x},\overline{y})=\frac{c(t,x,y)}{c_{{}_0}}$, $\overline{f}(\overline{t},\overline{x},\overline{y})=\frac{f(t,x,y)}{f_{{}_0}}$, $\overline{r}(\overline{t},\overline{x},\overline{y})=\frac{r(t,x,y)}{r_{{}_0}}$, $\overline{\beta }(\overline{t},\overline{x},\overline{y})=\frac{\beta (t,x,y)}{{\beta }_{{}_0}}$ and $\overline{\alpha }(\overline{t},\overline{x},\overline{y})=\frac{\alpha (t,x,y)}{{\alpha }_{{}_0}}$ where $\overline{t}=\frac{t}{\tau }$, $\overline{x}=\frac{x}{L}$, $\overline{y}=\frac{y}{L}$, while $h_{{}_0},i_{{}_0},v_{{}_0},c_{{}_0},f_{{}_0},r_{{}_0},{\beta }_{{}_0}$ and ${\alpha }_{{}_0}$ are proper fixed parameters. Furthermore, for the remaining variables, we select the following key parameters: ${\overline{D}}_{{}_h}=\frac{D_{{}_h}\tau }{L^2}=\frac{D_{{}_h}}{D}$, ${\overline{\mu }}_{{}_h}={\mu }_{{}_h}\tau$, ${\overline{\delta }}_{{}_{11}}={\delta }_{{}_{11}}\tau$, ${\overline{D}}_{{}_i}=\frac{D_{{}_i}\tau }{L^2}=\frac{D_{{}_i}}{D}$, $\overline{{\mu }_{{}_i}}=\frac{{\mu }_{{}_i}\tau h_0{v}_0}{i_{{}_0}}$, ${\overline{\delta }}_{{}_{21}}={\delta }_{{}_{21}}\tau f_{{}_0}$, ${\overline{\delta }}_{{}_{22}}={\delta }_{{}_{22}}\tau$, ${\overline{D}}_{{}_v}=\frac{D_{{}_v}\tau }{L^2}=\frac{D_{{}_v}}{D}$, ${\overline{\chi }}_{{}_v}=\frac{{\chi }_{{}_v}{c}_{{}_0}\tau }{L^2}$, ${\overline{\mu }}_{{}_v}=\frac{{\mu }_{{}_v}\tau i_{{}_0}}{v_{{}_0}}$, ${\overline{\delta }}_{{}_{31}}={\delta }_{{}_{31}}\tau {\alpha }_{{}_0}$, ${\overline{\delta }}_{{}_{32}}={\delta }_{{}_{32}}\tau$, ${\overline{D}}_{{}_c}=\frac{D_{{}_c}\tau }{L^2}=\frac{D_{{}_c}}{D}$, ${\overline{\mu }}_{{}_c}=\frac{{\mu }_{{}_c}\tau i_{{}_0}{f}_{{}_0}}{c_{{}_0}}$, ${\overline{\delta }}_{{}_{41}}={\delta }_{{}_{41}}\tau$, ${\overline{D}}_{{}_f}=\frac{D_{{}_f}\tau }{L^2}=\frac{D_{{}_f}}{D}$, ${\overline{\chi }}_{{}_f}=\frac{{\chi }_{{}_f}{c}_{{}_0}\tau }{L^2}$, ${\overline{\mu }}_{{}_f}=\frac{{\mu }_{{}_f}\tau }{f_{{}_0}}$, ${\overline{k}}_{{}_f}=\frac{k_{{}_f}}{v_{{}_0}}$, ${\overline{\delta }}_{{}_{51}}={\delta }_{{}_{51}}\tau r_{{}_0}$, ${\overline{\delta }}_{{}_{52}}={\delta }_{{}_{52}}\tau$, ${\overline{D}}_{{}_r}=\frac{D_{{}_r}\tau }{L^2}=\frac{D_{{}_r}}{D}$, ${\overline{\chi }}_{{}_r}=\frac{{\chi }_{{}_r}{c}_{{}_0}\tau }{L^2}$, ${\overline{\mu }}_{{}_r}=\frac{{\mu }_{{}_r}\tau }{r_{{}_0}}$, ${\overline{k}}_{{}_r}=\frac{k_{{}_r}}{f_{{}_0}}$, ${\overline{\delta }}_{{}_{61}}={\delta }_{{}_{61}}\tau$, ${\overline{D}}_{{}_{\beta }}=\frac{D_{{}_{\beta }}\tau }{L^2}=\frac{D_{{}_{\beta }}}{D}$, ${\overline{\chi }}_{{}_{\beta }}=\frac{{\chi }_{{}_{\beta }}{c}_{{}_0}\tau }{L^2}$, ${\overline{\mu }}_{{}_{\beta }}=\frac{{\mu }_{{}_{\beta }}\tau }{{\beta }_{{}_0}}$, ${\overline{k}}_{{}_{\beta }}=\frac{k_{{}_{\beta }}}{v_{{}_0}}$, ${\overline{\mu }}_{{}_{\tilde{\beta }}}=\frac{{\mu }_{{}_{\tilde{\beta }}}\tau f_{{}_0}}{{\beta }_{{}_0}}$, ${\overline{\delta }}_{{}_{71}}={\delta }_{{}_{71}}\tau$, ${\overline{D}}_{{}_{\alpha }}=\frac{D_{{}_{\alpha }}\tau }{L^2}=\frac{D_{{}_{\alpha }}}{D}$, ${\overline{\mu }}_{{}_{\alpha }}=\frac{{\mu }_{{}_{\alpha }}\tau {\beta }_{{}_0}}{{\alpha }_{{}_0}}$, ${\overline{\delta }}_{{}_{81}}={\delta }_{{}_{81}}\tau$. Subsequently, we substitute the above rescaled variables into the coupled reaction diffusion system (7) to obtain the following dimensionless version with removing the over-lines for easier readability:

$$\left\{\begin{array}{c}\frac{\partial h}{\partial t}=D_{{}_h}{\nabla }^{{}^2}h+{\mu }_{{}_h}h(1-h)-{\delta }_{{}_{11}}h,\hfill \\ \frac{\partial i}{\partial t}=D_{{}_i}{\nabla }^{{}^2}i+{\mu }_{{}_i}hv-{\delta }_{{}_{21}}if-{\delta }_{{}_{22}}i,\hfill \\ \frac{\partial v}{\partial t}=D_{{}_v}{\nabla }^{{}^2}v-{\chi }_{{}_v}\nabla ·(v\nabla c)+{\mu }_{{}_v}i-{\delta }_{{}_{31}}\alpha v-{\delta }_{{}_{32}}v,\hfill \\ \frac{\partial c}{\partial t}=D_{{}_c}{\nabla }^{{}^2}c+{\mu }_{{}_c}if-{\delta }_{{}_{41}}c,\hfill \\ \frac{\partial f}{\partial t}=D_{{}_f}{\nabla }^{{}^2}f-{\chi }_{{}_f}\nabla ·(f\nabla c)+{\mu }_{{}_f}\frac{v}{v+k_{{}_f}}-{\delta }_{{}_{51}}rf-{\delta }_{{}_{52}}f,\hfill \\ \frac{\partial r}{\partial t}=D_{{}_r}{\nabla }^{{}^2}r-{\chi }_{{}_r}\nabla ·(r\nabla c)+{\mu }_{{}_r}\frac{f}{f+k_{{}_r}}-{\delta }_{{}_{61}}r,\hfill \\ \frac{\partial \beta }{\partial t}=D_{{}_{\beta }}{\nabla }^{{}^2}\beta -{\chi }_{{}_{\beta }}\nabla ·(\beta \nabla c)+{\mu }_{{}_{\beta }}\frac{v}{v+k_{{}_{\beta }}}+{\mu }_{{}_{\tilde{\beta }}}f-{\delta }_{{}_{71}}\beta ,\hfill \\ \frac{\partial \alpha }{\partial t}=D_{{}_{\alpha }}{\nabla }^{{}^2}\alpha +{\mu }_{{}_{\alpha }}\beta -{\delta }_{{}_{81}}\alpha .\hfill \end{array}\right.$$

(8)

Thus, the above form of system is used during the computation procedures with taking into consideration the key parameters for conversion purposes.

3.2. Initial and Boundary Conditions

The initial spatial distribution of the host cells is supposed to be localised naturally everywhere (heterogeneous) in the computational domain $\aleph =[0,2]\times [0,2]$ (see Figure 1a), namely:

$$\begin{array}{c}\hfill h(0,x)=\frac{3}{4}+\frac{sin\left(5\pi {∥x∥}_{{}_2}\right)+sin\left(5\pi {∥(1,0)-x∥}_{{}_2}\right)}{8},x\in \aleph .\end{array}$$

Regarding the initial spatial distribution of SARS-CoV-2 particles, we assume initially the presence of the virus takes place within ℵ at multiple sites near by high concentrations of the initial spatial distribution of target cells as shown in Figure 1b, that is:

$$\begin{array}{c}\hfill v(0,x)=\sum _{i,j\in \{0,\pm 0.2,\pm 0.4\}}^{}{k}_{{}_c}\left(exp(-k_{{}_w}\frac{{∥x-(1+i,1+j)∥}_2^2}{\sqrt{\Delta x\Delta y}})\right)\phi \left(x\right),x\in \aleph ,\end{array}$$

where $\phi \left(x\right)$ is a smooth function used to approximate values to zero outside the ball $\mathbf{B}\left(\right(1+i,1+j),0.25)$, while $k_{{}_c}$ and $k_{{}_w}$ are positive parameters used to control the amount of viral concentration and viral spatial distribution width in the maximal reference region. Moreover, since we suppose at $t=0$ no pre-existing infected cells, chemokines, the effector T cells, the regulatory T cells, B cells and antibodies, hence we have the following initial conditions:

$$\left\{\begin{array}{c}i(0,x)=0,x\in \aleph ,\hfill \\ c(0,x)=0,x\in \aleph ,\hfill \\ f(0,x)=0,x\in \aleph ,\hfill \\ r(0,x)=0,x\in \aleph ,\hfill \\ \beta (0,x)=0,x\in \aleph ,\hfill \\ \alpha (0,x)=0,x\in \aleph .\hfill \end{array}\right.$$

(9)

Finally, our system is assumed to be solved within a fixed computational domain $\aleph =[0,2]\times [0,2]$ with no molecular transport across the domain boundary $\partial \aleph$, thus we impose Neumann zero flux boundary conditions for all components, namely:

$$\begin{array}{c}\hfill \langle \nabla k,\nu \rangle =0\forall k\in \{h,i,v,c,f,r,\beta ,\alpha \}\mathrm{on}(0,\infty )\times \partial \aleph ,\end{array}$$

where $\nu$ indicates the outward unit normal for vector field to the computational domain boundary $\partial \aleph$.

3.3. Computational Implementation

For computational purposes, we set the maximum spatial domain to be considered within a square region, namely, $\aleph :=[0,2]\times [0,2]$ which equals to a physical domain of length $[0,0.1]\times [0,0.1]\mathrm{cm}$. We discretise $\aleph \subset {\mathbb{R}}^2$ by an equal (uniform) spatial mesh in both directions of length $\Delta x=\Delta y:=1\times {10}^{-2}$ which is equivalent to $1\times {10}^{-3}\mathrm{cm}$ in a physical dimensional domain. For time, we discretise the time interval $[t_0,t_0+\Delta t]$ by a uniformly time step of size $\widehat{\delta t}:=\frac{\Delta t}{\mathcal{N}}$ with time nodes $\mathcal{N}+1>1$. Thus, at any spatially discretized point $(p\Delta x,q\Delta y)$, $\forall p,q:=0,...$ and at any discretized time node $t_{{}_n}=t_0+n\widehat{\delta t}$, $\forall n=0,1,...,\mathcal{N}$, we denote the concentrations for all components of the system (8) by ${\mathcal{M}}_{{}_k}(t_{{}_n},(x_{{}_p},y_{{}_q}))$$\forall k\in \{h,i,v,c,f,r,\beta ,\alpha \}$. Further, we express the spatial operators for each component of the system (8) by ${\mathcal{S}}_{{}_k}(t_{{}_n},(x_{{}_p},y_{{}_q})):={\nabla }^{{}^2}{\mathcal{M}}_{{}_k}(t_{{}_n},(x_{{}_p},y_{{}_q}))$ and ${\widehat{\mathcal{S}}}_{\widehat{k}}(t_{{}_n},(x_{{}_p},y_{{}_q})):=\nabla {\mathcal{M}}_{\widehat{k}}(t_{{}_n},(x_{{}_p},y_{{}_q}))$$\forall \widehat{k}\in \{v,f,r,\beta \}$. Hence, at any discretized spatial temporal node $(t_{{}_n},(x_{{}_p},y_{{}_q}))$, we approximate the diffusion terms in (8) by a central difference scheme using a fictitious grid points of the form ${\left({\mathcal{M}}_{{}_k}\right)}_{p\pm \frac{1}{2},q}^n,{\left({\mathcal{M}}_{{}_k}\right)}_{p,q\pm \frac{1}{2}}^n$, namely:

$$\begin{array}{c}\hfill {\left({\mathcal{S}}_{{}_k}\right)}_{p,q}^n\approx \frac{{\left({\left({\mathcal{M}}_{{}_k}\right)}_{{}_x}\right)}_{p+\frac{1}{2},q}^n-{\left({\left({\mathcal{M}}_{{}_k}\right)}_{{}_x}\right)}_{p-\frac{1}{2},q}^n}{\Delta x}+\frac{{\left({\left({\mathcal{M}}_{{}_k}\right)}_{{}_y}\right)}_{p,q+\frac{1}{2}}^n-{\left({\left({\mathcal{M}}_{{}_k}\right)}_{{}_y}\right)}_{p,q-\frac{1}{2}}^n}{\Delta y}.\end{array}$$

Additionally, we approximate the chemotaxis terms in (8) by the following expression:

$$\begin{array}{c}\hfill {\left({\widehat{\mathcal{S}}}_{\widehat{k}}{\mathcal{S}}_c\right)}_{p,q}^n\approx \frac{{\left({\mathcal{M}}_{{}_{\widehat{k}}}\right)}_{p+\frac{1}{2},q}^n{\left({\left({\mathcal{M}}_{{}_c}\right)}_{{}_x}\right)}_{p+\frac{1}{2},q}^n-{\left({\mathcal{M}}_{{}_{\widehat{k}}}\right)}_{p-\frac{1}{2},q}^n{\left({\left({\mathcal{M}}_{{}_c}\right)}_{{}_x}\right)}_{p-\frac{1}{2},q}^n}{\Delta x}\\ \hfill +\frac{{\left({\mathcal{M}}_{{}_{\widehat{k}}}\right)}_{p,q+\frac{1}{2}}^n{\left({\left({\mathcal{M}}_{{}_c}\right)}_{{}_y}\right)}_{p,q+\frac{1}{2},}^n-{\left({\mathcal{M}}_{{}_{\widehat{k}}}\right)}_{p,q-\frac{1}{2}}^n{\left({\left({\mathcal{M}}_{{}_c}\right)}_{{}_y}\right)}_{p,q-\frac{1}{2}}^n}{\Delta y},\end{array}$$

where

$$\begin{array}{c}\hfill {\left({\left({\mathcal{M}}_{{}_k}\right)}_{{}_x}\right)}_{p+\frac{1}{2},q}^n:=\frac{{\left({\mathcal{M}}_{{}_k}\right)}_{p+1,q}^n-{\left({\mathcal{M}}_{{}_k}\right)}_{p,q}^n}{\Delta x},\\ \hfill {\left({\left({\mathcal{M}}_{{}_k}\right)}_{{}_x}\right)}_{p-\frac{1}{2},q}^n:=\frac{{\left({\mathcal{M}}_{{}_k}\right)}_{p,q}^n-{\left({\mathcal{M}}_{{}_k}\right)}_{p-1,q}^n}{\Delta x},\\ \hfill {\left({\left({\mathcal{M}}_{{}_k}\right)}_{{}_y}\right)}_{p,q+\frac{1}{2}}^n:=\frac{{\left({\mathcal{M}}_{{}_k}\right)}_{p,q+1}^n-{\left({\mathcal{M}}_{{}_k}\right)}_{p,q}^n}{\Delta x},\\ \hfill {\left({\left({\mathcal{M}}_{{}_k}\right)}_{{}_y}\right)}_{p,q-\frac{1}{2}}^n:=\frac{{\left({\mathcal{M}}_{{}_k}\right)}_{p,q}^n-{\left({\mathcal{M}}_{{}_k}\right)}_{p,q-1}^n}{\Delta x},\end{array}$$

and

$$\begin{array}{cc}\hfill {\left({\mathcal{M}}_{{}_{\widehat{k}}}\right)}_{p\pm \frac{1}{2},q}^n& :=\frac{{\left({\mathcal{M}}_{{}_{\widehat{k}}}\right)}_{p,q}^n+{\left({\mathcal{M}}_{{}_{\widehat{k}}}\right)}_{p\pm 1,q}^n}{2},\\ \hfill {\left({\mathcal{M}}_{{}_{\widehat{k}}}\right)}_{p,q\pm \frac{1}{2}}^n& :=\frac{{\left({\mathcal{M}}_{{}_{\widehat{k}}}\right)}_{p,q}^n+{\left({\mathcal{M}}_{{}_{\widehat{k}}}\right)}_{p,q\pm 1}^n}{2}.\end{array}$$

For time discretisation, we use the predictor corrector method to approximate the solutions i.e., on the interval $[t_n,t_{n+1}]$, we predict the density for each component $\left({\mathcal{M}}_{{}_k}\right)$ by the explicit Euler method, namely:

$$\begin{array}{c}\hfill \begin{array}{c}{\left({\tilde{\mathcal{M}}}_{{}_k}\right)}_{p,q}^{n+1}={\left({\mathcal{M}}_{{}_k}\right)}_{p,q}^n+\delta t\left(\mathcal{F}(h_{p,q}^n,i_{p,q}^n,v_{p,q}^n,c_{p,q}^n,f_{p,q}^n,r_{p,q}^n,{\beta }_{p,q}^n,{\alpha }_{p,q}^n)\right),\hfill \end{array}\end{array}$$

where the function $\mathcal{F}(·)$ denotes all terms of the right hand side of the system (8) including the diffusion $\left({\mathcal{S}}_k\right)$ and chemotaxis $\left({\widehat{\mathcal{S}}}_{\widehat{k}}\right)$ operators. Then, we correct the predicted solutions via the following Trapezoidal approximation rule as well as using predict-evaluate-correct mode to obtain a better approximation of the solution, namely,

$$\begin{array}{c}\hfill \begin{array}{c}{\left({\widehat{\mathcal{M}}}_{{}_k}\right)}_{p,q}^{n+1}={\left({\mathcal{M}}_{{}_k}\right)}_{p,q}^n+\frac{1}{2}\delta t(\mathcal{F}(h_{p,q}^n,i_{p,q}^n,v_{p,q}^n,c_{p,q}^n,f_{p,q}^n,r_{p,q}^n,{\beta }_{p,q}^n,{\alpha }_{p,q}^n)+\hfill \\ \mathcal{F}({\tilde{h}}_{p,q}^n,{\tilde{i}}_{p,q}^n,{\tilde{v}}_{p,q}^n,{\tilde{c}}_{p,q}^n,{\tilde{f}}_{p,q}^n,{\tilde{r}}_{p,q}^n,{\tilde{\beta }}_{p,q}^n,{\tilde{\alpha }}_{p,q}^n)).\hfill \\ {\left({\mathcal{M}}_{{}_k}\right)}_{p,q}^{n+1}={\left({\mathcal{M}}_{{}_k}\right)}_{p,q}^n+\frac{1}{2}\delta t(\mathcal{F}(h_{p,q}^n,i_{p,q}^n,v_{p,q}^n,c_{p,q}^n,f_{p,q}^n,r_{p,q}^n,{\beta }_{p,q}^n,{\alpha }_{p,q}^n)+\hfill \\ \mathcal{F}({\widehat{h}}_{p,q}^n,{\widehat{i}}_{p,q}^n,{\widehat{v}}_{p,q}^n,{\widehat{c}}_{p,q}^n,{\widehat{f}}_{p,q}^n,{\widehat{r}}_{p,q}^n,{\widehat{\beta }}_{p,q}^n,{\widehat{\alpha }}_{p,q}^n)).\hfill \end{array}\end{array}$$

Thus, following the same procedures, we obtain solutions on the computational domain ℵ for each component of (8), that is, for all ${\left({\mathcal{M}}_{{}_k}\right)}_{p,q}^n$$\forall k\in \{h,i,v,c,f,r,\beta ,\alpha \}$.

3.4. Data Estimation

In general, parameter estimation in computational immunology is still a major problem due to either the lack of experimental data at a cellular level or the differences between results from study to study. The estimations of parameters are impacted by many factors such as cells shape, size and type, methods and laboratory equipments as well as the interpretation of the laboratory data [45,46]. However, simulations may provide an efficient manner for immunological exploration and prediction. Therefore, the model validation could be investigated without knowing the data where this common methodology is used during the exploratory steps [47]. Although the shortage of experimental studies about COVID-19, we select a list of values for model parameters to simulate the baseline case for which they are considered as a starting point for the numerical investigation. The baseline values (as shown in

Table 1. Baseline parameters values used in numerical simulations for the spatio-temporal model of COVID-19. ${(}^{*})$ indicates the modified values.

Parameter	Non-Dimensional Value	Dimensional Value	References
$D_{{}_h}$	$9\times {10}^{-4}$	$9\times {10}^{-10}$ cm${}^2$ s${}^{-1}$	[44]
$D_{{}_i}$	$1\times {10}^{-4}$	$1\times {10}^{-10}$ cm${}^2$ s${}^{-1}$	Est. $(D_{{}_i}=\frac{D_{{}_h}}{9})$
$D_{{}_v}$	$5\times {10}^{-4}$	$5\times {10}^{-10}$ cm${}^2$ s${}^{-1}$	[44]
$D_{{}_c}$	$1.06\times {10}^{-3}$	$1.06\times {10}^{-9}$ cm${}^2$ s${}^{-1}$	[34]
$D_{{}_f}$	$1.83\times {10}^{-3}$	$1.83\times {10}^{-9}$ cm${}^2$ s${}^{-1}$	[34]
$D_{{}_r}$	$1\times {10}^{-3}$	$1\times {10}^{-9}$ cm${}^2$ s${}^{-1}$	[34]
$D_{{}_{\beta }}$	$1.83\times {10}^{-3}$	$1.83\times {10}^{-9}$ cm${}^2$ s${}^{-1}$	Est. $(D_{{}_{\beta }}=D_{{}_f})$
$D_{{}_{\alpha }}$	$1.83\times {10}^{-3}$	$1.83\times {10}^{-9}$ cm${}^2$ s${}^{-1}$	Est. $(D_{{}_{\alpha }}=D_{{}_f})$
${\chi }_{{}_v}$	$1.67\times {10}^{-4}$	$1.66\times {10}^{-9}$ cm${}^2$ s${}^{-1}$	Est. $({\chi }_{{}_v}={\chi }_{{}_f})$
${\chi }_{{}_f}$	$1.67\times {10}^{-4}$	$1.66\times {10}^{-9}$ cm${}^2$ s${}^{-1}$	[34]
${\chi }_{{}_r}$	$1.67\times {10}^{-4}$	$1.66\times {10}^{-9}$ cm${}^2$ s${}^{-1}$	[34]
${\chi }_{{}_{\beta }}$	$1.67\times {10}^{-4}$	$1.66\times {10}^{-9}$ cm${}^2$ s${}^{-1}$	Est. $({\chi }_{{}_{\beta }}={\chi }_{{}_f})$
${\mu }_{{}_h}$	$0.462{}^{*}$	$0.462\times {10}^{-1}$ cells mL${}^{-1}$ s${}^{-1}$	[26,27]
${\mu }_{{}_i}$	$0.826$	$1.736\times {10}^{-20}$ mL (RNA copies)${}^{-1}$ s${}^{-1}$	[22,27]
${\mu }_{{}_v}$	$0.255$	$6\times {10}^{-3}$${}^{*}$ s${}^{-1}$	[27]
${\mu }_{{}_c}$	$0.174$	$0.115\times {10}^{-4}$ s${}^{-1}$	[26,27]
${\mu }_{{}_f}$	$0.694$	$0.347\times {10}^{-1}$${}^{*}$ mL cells${}^{-1}$ s${}^{-1}$	[19,27]
${\mu }_{{}_r}$	$0.077$	$0.11\times {10}^{-2}$ mL cells${}^{-1}$ s${}^{-1}$	Est. $({\mu }_{{}_r}=\frac{{\mu }_{{}_f}}{30})$
${\mu }_{{}_{\beta }}$	$0.0116$	$0.11\times {10}^{-3}$${}^{*}$ mL cells${}^{-1}$ s${}^{-1}$	[19,27]
${\mu }_{{}_{\tilde{\beta }}}$	$0.579$	$0.11\times {10}^{-4}$ mL cells${}^{-1}$ s${}^{-1}$	Est. $({\mu }_{{}_{\tilde{\beta }}}=\frac{{\mu }_{{}_{\beta }}}{10})$
${\mu }_{{}_{\alpha }}$	$0.694$	$0.81\times {10}^{-7}$ s${}^{-1}$	[19,27]
${\delta }_{{}_{11}}$	$0.0162$	$1.62\times {10}^{-6}$ s${}^{-1}$	[19,27]
${\delta }_{{}_{21}}$	$0.0579$	$1.15\times {10}^{-8}$ mL cells${}^{-1}$ s${}^{-1}$	[27]
${\delta }_{{}_{22}}$	$0.0011$	$1.15\times {10}^{-7}$ s${}^{-1}$	[19,27]
${\delta }_{{}_{31}}$	$0.116$	$1.15\times {10}^{-4}$ ml molecules${}^{-1}$ s${}^{-1}$	[19,27]
${\delta }_{{}_{32}}$	$0.04$	$4.05\times {10}^{-6}$${}^{*}$ s${}^{-1}$	[27]
${\delta }_{{}_{41}}$	$0.0016$	$1.66\times {10}^{-7}$ s${}^{-1}$	[34]
${\delta }_{{}_{51}}$	$0.52$	$3.47\times {10}^{-7}$ mL cells${}^{-1}$ s${}^{-1}$	[34]
${\delta }_{{}_{52}}$	$0.0034$	$3.47\times {10}^{-7}$ s${}^{-1}$	[34]
${\delta }_{{}_{61}}$	$0.0034$	$3.47\times {10}^{-7}$ s${}^{-1}$	[34]
${\delta }_{{}_{71}}$	$0.023$	$2.31\times {10}^{-6}$ s${}^{-1}$	[19,26]
${\delta }_{{}_{81}}$	$0.0081$	$8.1\times {10}^{-7}$ s${}^{-1}$	[19,27]
$k_{{}_f}$	$8.14\times {10}^{-1}$	≈$2.9\times {10}^6$	[22]
$k_{{}_r}$	$5.81\times {10}^3$	≈$2.9\times {10}^6$	[22]
$k_{{}_{\beta }}$	$8.14\times {10}^{-1}$	≈$2.9\times {10}^6$	[22]
$\Delta x$, $\Delta y$	$1\times {10}^{-2}$	$1\times {10}^{-3}$ cm	[44]
$\Delta t$	$1\times {10}^{-3}$	10 s	−

) have been built on clinical and mathematical results, namely [19,22,26,27,34,44]. We estimate diffusion rates building on the size and the shape of cells and particles as detailed in [44] with taking into consideration the unit of diffusion as the area per unit time where we mainly use centimeters as units for measuring length and seconds for time. Further, we pick out the other motility rates and chemotaxis coefficients for immune system components based on biological ranges specified in [34].

For production and decay rates, we select the baseline values depending on few recent mathematical studies [22,26,27] used temporal (within host) modelling approaches for COVID-19. However, most of these data have been collected from previous studies for similar infectious diseases such as dengue infection, HIV infection and AIDS. Therefore, we still face a lack of experimental data for infection with SARS-CoV-2, especially at the cellular scale.

4. A Comparison of Computational and Clinical Results

In this section, we aim to investigate numerically the dynamics of host cells $h(t,x)$, infected cells $i(t,x)$, SARS-CoV-2 particles $v(t,x)$, chemokines $c(t,x)$, the effector T cells $f(t,x)$, the regulatory T cells $r(t,x)$, B-lymphocytes cells $\beta (t,x)$ and antibodies $\alpha (t,x)$. These dynamics are computed at each discretized spatio-temporal node $(t_{{}_n},(x_{{}_p},y_{{}_q}))$ within a computational domain of size $\aleph =[0,2]\times [0,2]$ which equals to a physical domain $[0,0.1]\times [0,0.1]\mathrm{cm}$. For time steps, we run the numerical simulation for 120 stages (120,000$\Delta t$), which equals approximately 333.34 h. However, SARS-CoV-2 particles may remain in the body for more than 60 days or sometimes reaching 150 days as discussed in [5,19,48], which corresponds to 259,200$\Delta t-$650,000$\Delta t$ (in non-dimensionalised time step). Moreover, we choose the following initial values for the variables those defined in Section 3.1 based on some biological and mathematical observations [19,20,27], $h_{{}_0}=4\times {10}^5,i_{{}_0}=3\times {10}^{-4},v_{{}_0}=357\times {10}^4,c_{{}_0}=0.1,f_{{}_0}=500,r_{{}_0}=150,{\beta }_{{}_0}=100,{\alpha }_{{}_0}=0.1$. In the following context, we discuss the validation of the suggested model (8) with regards to experimental observations in terms of the viral load of SARS-CoV-2 during time. As mentioned in Section 1, there are a bunch of studies that discuss a within host modeling of SARS-CoV-2 interactions with immune system components, but these studies consider only the temporal dynamics of this interaction. Hence, there is a lack of both computational and experimental data for the spatial distribution of SARS-CoV-2 within the cellular scale. However, we still wish that our suggested spatio-temporal model (8) will stimulate experimental works to validate our numerical findings.

Now, in order to check the validation of our modeling assumptions, we first select a list of values for the model parameters as presented in Table 1 which mainly based on the following published studies [19,22,26,27,34,44]. Using this baseline values, we show the spatial distribution for all dynamics as illustrated in Figure 2 where we present (for three different time periods, namely, ≈2.5 days, 7 days, 14 days): healthy cells (a), infected cells (b), SARS-CoV-2 particles (c), chemokines (d), the effector T cells (e), the regulatory T cells (f), B-lymphocytes cells (g) and antibodies (h). It is obvious to see that the spatial distribution of immune components follows the same distribution patterns of virus. Moreover, for the parameter regimes that we have selected as a starting point for the computational investigation, we find that the total viral load of SARS-CoV-2 increases since day one of infection. This increase of virus occurs due to the slow growth of T and B cells over time (see Figure 3a,b). However, comparing our computational result with patients data for the whole clinical course that have been considered in [49], we conclude the following, in the clinical data, the average virus load in swab test was about $3.44\times {10}^5$ copies per mL on the day 5 which equals to $9.64\times {10}^{-2}$ in non-dimensional values, while numerically, the viral load on the day 5 is almost equal to $9.58\times {10}^{-2}$, that is, $3.42\times {10}^5$ (see Figure 4a–i). Nevertheless, the baseline scenario of the numerical simulation represents a failure of immune responses to infection with SARS-CoV-2 because it is clearly to see that the density of infected cells increases from day one of infection until 14 days. This immunity failure probably occurred due to the late response of T cells (especially for asymptomatic cases) as presented in Figure 5 where we compare the effector T cells evolution at early time stages for two cases of immune responses to SARS-CoV-2, namely, baseline scenario and a better case of immunity responses (it will be explained in detail in the next context). Further, the slow response of T cells is negatively reflected on B cells growth and the activation process of antibodies.

In the following, we investigate the impact of variation of parameters on the level of SARS-CoV-2 viral load over a certain period of time. To that end, we first denote the virus domain ${\aleph }^{*}\left(t_{{}_0}\right)$ to be a sub-domain located entirely within our computational domain, that is, ${\aleph }^{*}\left(t_{{}_0}\right)\subset \aleph \left(t_{{}_0}\right)$. Then, we define a new function $\psi (·,·)$ which enables increasing or decreasing the level of all system components exclusively at each grid point containing a non-zero solution of virus i.e., at every ${\left(v\right)}_{p,q}^n{|}_{v\ne 0},\forall p,q\in {\aleph }^{*}\left(t_{{}_0}\right)$. Accordingly, $\psi (·,·)$ is expressed as follows:

$${\psi }_{{}_{\pm }}(p,q):=\left\{\begin{array}{c}\zeta \mathrm{if}(x_{{}_p},y_{{}_q})\in \overline{{\aleph }^{*}\left(t_{{}_0}\right)},\hfill \\ 1\mathrm{if}(x_{{}_p},y_{{}_q})\notin \overline{{\aleph }^{*}\left(t_{{}_0}\right)},\hfill \end{array}\right.$$

(10)

where $\overline{{\aleph }^{*}\left(t_{{}_0}\right)}$ is the topological closure of ${\aleph }^{*}\left(t_{{}_0}\right)$ and the percentage variable $\zeta$ is given by the following choices:

$$\zeta :=\{100\%,150\%,200\%\}.$$

(11)

Thus, the system (8) can be rewritten in the following form:

$$\left\{\begin{array}{c}\frac{\partial h}{\partial t}=D_{{}_h}{\nabla }^{{}^2}h+{\mu }_{{}_h}h(1-h)-{\delta }_{{}_{11}}h,\hfill \\ \frac{\partial i}{\partial t}=D_{{}_i}{\nabla }^{{}^2}i+{\psi }_{{}_{-}}{\mu }_{{}_i}hv-{\delta }_{{}_{21}}if-{\delta }_{{}_{22}}i,\hfill \\ \frac{\partial v}{\partial t}=D_{{}_v}{\nabla }^{{}^2}v-{\chi }_{{}_v}\nabla ·(v\nabla c)+{\psi }_{{}_{-}}{\mu }_{{}_v}i-{\delta }_{{}_{31}}\alpha v-{\delta }_{{}_{32}}v,\hfill \\ \frac{\partial c}{\partial t}=D_{{}_c}{\nabla }^{{}^2}c+{\mu }_{{}_c}if-{\delta }_{{}_{41}}c,\hfill \\ \frac{\partial f}{\partial t}=D_{{}_f}{\nabla }^{{}^2}f-{\chi }_{{}_f}\nabla ·(f\nabla c)+{\psi }_{{}_{+}}{\mu }_{{}_f}\frac{v}{v+k_{{}_f}}-{\delta }_{{}_{51}}rf-{\delta }_{{}_{52}}f,\hfill \\ \frac{\partial r}{\partial t}=D_{{}_r}{\nabla }^{{}^2}r-{\chi }_{{}_r}\nabla ·(r\nabla c)+{\psi }_{{}_{-}}{\mu }_{{}_r}\frac{f}{f+k_{{}_r}}-{\delta }_{{}_{61}}r,\hfill \\ \frac{\partial \beta }{\partial t}=D_{{}_{\beta }}{\nabla }^{{}^2}\beta -{\chi }_{{}_{\beta }}\nabla ·(\beta \nabla c)+{\psi }_{{}_{+}}{\mu }_{{}_{\beta }}\frac{v}{v+k_{{}_{\beta }}}+{\psi }_{{}_{+}}{\mu }_{{}_{\tilde{\beta }}}f-{\delta }_{{}_{71}}\beta ,\hfill \\ \frac{\partial \alpha }{\partial t}=D_{{}_{\alpha }}{\nabla }^{{}^2}\alpha +{\psi }_{{}_{+}}{\mu }_{{}_{\alpha }}\beta -{\delta }_{{}_{81}}\alpha .\hfill \end{array}\right.$$

(12)

Just to note that the system (12) is only computed by using one of the functions $({\psi }_{{}_{+}},{\psi }_{{}_{-}})$, for example, if we plug the function $\left({\psi }_{{}_{+}}\right)$ into system (12), the other function $\left({\psi }_{{}_{-}}\right)$ will be implemented as a matrix of ones during computations. Based on that, we can divide the parameters into two main groups, namely, (A) parameters that, if increased, would decrease the overall viral load over time (motivation occurs at a rate ${\psi }_{{}_{+}}$) and (B) parameters that, if increased, would increase the level of viral load (that occurs at a rate ${\psi }_{{}_{-}}$).

Now, we start the sensitivity analysis using the percentage variable $\zeta$ to explore the impact of parameters defined in group (A) upon the overall result. In general, increasing this group of parameters leads to a proportional decay on the overall viral load of SARS-CoV-2 over time. In other words, increasing the percentage variable $\left(\zeta \right)$ brings a clear influence towards controlling the infection with SARS-CoV-2. This result is obvious in Figure 4a, where we show the total viral load of SARS-CoV-2 over 120 stages which corresponds to approximately 14 days. The black curve represents the basic scenario using the baseline values presented in Table 1, while the blue, green and red curves illustrate the same case, but with employing the percentage variable $\zeta :=100\%$, $\zeta :=150\%$ and $\zeta :=200\%$ respectively. The experimental study [41] shows the amount of viral load of SARS-CoV-2 for more than 80 patients over 12 to 15 days. The peak of viral load occurs between 5 to 6 or 7 days post symptom onset. However, the median time of symptom onset was varying from 2 to 7 days according to [22]. Here, we mainly focus on two patients from Beijing where their daily clinical samples were collected from throat swabs, urine, stool and sputum as presented in [41]. The first patient has a peak of viral load in throat swap in almost 7 days, while the second patient has that peak by the day 5. Hence, in order to compare our computational findings with the above samples results, we find that in the case when $\zeta :=200\%$, simulations show that the peak of viral load takes place at stage 43 which corresponds to the day 5 where this result coincides with clinical samples collected from the throat swabs for the second patient (see Figure 4a(iv)). Moreover, the viral load copies (per mL) at the peak point is ranging from ${10}^4$ to ${10}^8$ as collected from clinical data [41,49], while in our numerical results, when $\zeta :=200\%$, the viral load (on the day 5) is equal to $9.69\times {10}^{-2}$ copies per ml (≈$3.46\times {10}^5$—dimensional value). However, in the case when $\zeta :=100\%,150\%$, the peak of viral load appears at stage 80 and 54 (which corresponds to the day 9.25 and the day 6.25) with concentrations of $9.99\times {10}^{-2}$ and $9.79\times {10}^{-2}$ copies per ml (which equals to $3.565\times {10}^5$ and $3.495\times {10}^5$), respectively. Consequently, for both cases of $\zeta :=150\%,200\%$, the viral load tends to peak within the same time range compared to data collected from clinical trials. Further, in Figure 6a, we show the spatial distribution of virus particles at the final stage (14 days) for the three cases of $\zeta$, namely, $\zeta :=$ (i) $100\%$, (ii) $150\%$ and (iii) $200\%$. Thus, this result confirms that increasing the parameter values specified in group (A) leads to a better scenario with regards to infection with SARS-CoV-2.

On other hand, to observe the impact of the parameters defined in group (B) on the overall virus spatial expansion and the total viral load during interactions between virus and immune components, we first let ${\psi }_{{}_{-}}(p,q)$ be denoted by the formula defined in (10), while substituting the function ${\psi }_{{}_{+}}(p,q)$ by one into model (12) $\forall p,q\in {\aleph }^{*}\left(t_{{}_0}\right)$. Then, using parameters values presented in Table 1 and changing $\zeta$ by percentages determined in (11), we show a comparison of viral load of SARS-CoV-2 for the above scenario starting from 2.78 h (stage 1) to 14 days (stage 120). Numerically, we obtain a higher viral load and more expansion of spatial distribution on the computational domain whenever $\zeta$ is increased (see Figure 4b and Figure 6b). This undesirable growth of virus can be explained by the fact that both rates of infection and viral replication have been considered within group (B), therefore increasing these rates by ${\psi }_{{}_{-}}$ will mostly cause more viral load over time. Following similar scenarios from clinical findings for died patients as presented in [41], the viral load on the day 8 was at least equal to ${10}^6$ for nasal swab and $1.34\times {10}^{11}$ for sputum sample. Hence, in a comparison to our simulations, the viral load on the day 8 (at stage 69) is greater than ${10}^6$ (copies per mL) in dimensional values for all cases of $\zeta$, namely, when $\zeta :=100\%$: viral load is $0.376$ (non-dim) which equals to ≈$1.34\times {10}^6$ (dim); when $\zeta :=150\%$: viral load is $0.687$ (non-dim) ≈$2.45\times {10}^6$ (dim); when $\zeta :=200\%$: viral load is $1.08$ (non-dim) ≈$3.85\times {10}^6$ (dim). Thus, rising the parameters defined in group (B) drives to a serious situation of infection with SARS-CoV-2.

Another factor that affects SARS-CoV-2 infection is the variation of random motility rates for cells, virus and chemokines on the spatial domain. To check this, we decide to duplicate, triplicate and quadruplicate the diffusion rates for all system components using best outcomes in regards to viral load control, namely, substituting the function ${\psi }_{{}_{+}}$ as defined in (10) into system (12) with $\zeta :=\{150\%,200\%\}$, while ${\psi }_{{}_{-}}:=1$ (these cases are illustrated in Figure 4a(iii,iv)). Accordingly, in Figure 7 we show the evolution of viral load of SARS-CoV-2 over 14 days using three values of random motility rates for all system components i.e., $2D_{{}_k},3D_{{}_k},4D_{{}_k}\forall k\in \{h,i,v,c,f,r,\beta ,\alpha \}$ where Figure 7a represents the viral load evolution in the case when $\zeta :=\{150\%\}$, while Figure 7b shows the same case but when $\zeta :=\{200\%\}$. In general, numerical results show that there is a viral load growth whenever diffusion coefficients are increased. This result could be explained as follows, fast spatial diffusion of system components (especially infected cells and virus particles) within the tissue (our computational domain ℵ) could expand the area of the dynamical interaction as well as transfer the virus to respiratory system which that may drive to have a late and a huge immune responses causing the cytokine release syndrome and leading to serious damage of vital organs and death in the worst cases [8,11,16].

5. Discussion

In contrast to the temporal mathematical models (either between-host or within-host models) as presented in [21,22,23,24,25,26], we here suggested a new spatio-temporal (within-host) mathematical model to describe the interactions at the cellular scale between healthy cells, viral infected cells, SARS-CoV-2 particles and the main parts of immune system, namely, effector and regulatory T cells, B-lymphocytes cells and antibodies in the presence of chemokines concentration. Modeling assumptions and settings were carefully selected building on recent biomedical and mathematical studies as detailed in Section 2. The numerical simulations provided reasonable outcomes for the within-host infection scenario which harmonized with data collected from clinical samples. However, we have faced many challenges regarding data availability as well as investigating the model outputs in a comparison with observations collected from clinical studies.

To investigate the model validation, we have improved a part of the computational framework introduced in [50] and then extended in [39]. Besides the implementation process of the finite difference method including approximations of the diffusion and chemotaxis terms, we have inserted additional strategies involving predict-evaluate-correct mode (see Section 3.3) with taking into account a multi-site of infection on the computational domain and assuming a heterogeneous initial condition for host cells (see Section 3.2) as well as applying certain conditions such as increasing the density of system components exclusively on the virus domain ${\aleph }^{*}\left(t_{{}_0}\right)$. We investigated numerically multiple scenarios compared to the baseline parameter regimes presented in Table 1 where the majority of these values are selected from [19,22,26,27,34,44]. However, simulation results using these values showed a failure case of immune responses to SARS-CoV-2 as presented in Figure 4a(i). Hence, we decided to define the functions ${\psi }_{{}_{\pm }}(·,·)$ (which exercise their role exclusively at ${\aleph }^{*}\left(t_{{}_0}\right)\subset \aleph \left(t_{{}_0}\right)$) in order to study the sensitivity of solutions using three choices of the percentage variable $\left(\zeta \right)$. Based on that, we conclude the following, increasing the production rates of some immune components, namely, (${\mu }_{{}_f},{\mu }_{{}_{\beta }},{\mu }_{{}_{\tilde{\beta }}},{\mu }_{{}_{\alpha }}$) by a size of ${\psi }_{{}_{+}}(·,·)$ leads to a clear decay in the total viral load during time, conversely, rising the parameters (${\mu }_{{}_i},{\mu }_{{}_v},{\mu }_{{}_r}$) by a quantity of ${\psi }_{{}_{-}}(·,·)$ encouraged the viral replication and reduced T cells production which in turn caused a growth of the total viral load. Additionally, in the case when ${\psi }_{{}_{+}}(·,·)$ is considered into model (12), we have found that staying at the peak of viral load differs as the percentage variable changes, for example, when $\zeta :=100\%$, the viral load reaches the peak after passing 9.25 days and remains on peak for another 36 h and 6.6 min, then decreases on the day 10.5 (for the other cases of $\left(\zeta \right)$, see

Table 2. Duration and amount of viral load on the peak for each case of $\zeta$ in the case when ${\psi }_{{}_{+}}(·,·)$ is considered into model (12) and ${\psi }_{{}_{-}}:=1$.

Case	Duration of Viral Load Peak	Viral Load (Non-Dim)	Viral Load (Dim)
$\zeta :=100\%$	36.11 h (stage 80–92)	$9.99\times {10}^{-2}$	$3.565\times {10}^5$
$\zeta :=150\%$	25 h (stage 54–62)	$9.79\times {10}^{-2}$	$3.495\times {10}^5$
$\zeta :=200\%$	25 h (stage 40–48)	$9.69\times {10}^{-2}$	$3.46\times {10}^5$

). Furthermore, increasing the motility rates for system components would play a central role on the growth of the viral load over time. In general, for all cases of ($\zeta$), our numerical results almost fit clinical outcomes for some patients data collected from [41,49] in terms of the amount of viral load on peaks (in dimensional values—copies per ml) as well as the time of the peak of viral load.

In this work, we have mainly presented a mathematical framework for modeling the immune responses to SARS-CoV-2 using a system of partial differential equations. We have developed a computational approach in order to simulate some scenarios of viral infection at a cellular level. Numerically, we have obtained a better control of SARS-CoV-2 infection in the case when the production level of T cells increases, and B cells and antibodies are exclusively in the spatial domain of the virus. From an immunological perspective, the ability to generate a strong immune response varies between people due to multiple factors. Therefore, we hope that the assumptions formulated in this work will stimulate suitable experimental studies that would investigate the production level of immune components at the infection site. However, in the future, we aim to provide further insight into the other factors such as the resulting pathological damages caused by either virus interactions with surrounding components or the cytokine storm caused by immune system. Moreover, questions regarding the qualitative behavior of the model, the local existence and uniqueness remain an open problem which could be investigated in a later work. Further, topics of COVID-19 are one of the most important research priorities nowadays; therefore, mathematical works should keep track of these efforts which require further improvement and future work. Thus, the proposed modeling framework and the computational technique can be extended to follow new findings from the forthcoming experimental studies.