Spatial and temporal dynamics of SARS-CoV-2: Modeling, analysis and simulation.

A reaction-diffusion viral infection model is formulated to characterize the infection process of the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) in a heterogeneous environment. In the model, the viral production, infection and death rates of compartments are given by the general functions. We consider the well-posedness of the solution, derive the basic reproduction number R 0 , discuss the global stability of uninfected steady state and explore the uniform persistence for the model. We further propose a spatial diffusion SARS-CoV-2 infection model with humoral immunity and spatial independent coefficients, and analyze the global attractivity of uninfected, humoral inactivated and humoral activated equilibria which are determined by two dynamical thresholds. Numerical simulations are performed to illustrate our theoretical results which reveal that diffusion, spatial heterogeneity and incidence types have evident impact on the SARS-CoV-2 infection process which should not be neglected for experiments and clinical treatments.

Introduction

In late December 2019, the first official case of coronavirus disease 2019 (COVID-19) was reported in Wuhan City, China, with devastating consequences around the world [1]. COVID-19, the pandemic caused by SARS-CoV-2, has found an abundantly susceptible individuals with no previous immunities [2], [3]. Quite many models have been proposed to investigate the ongoing COVID-19 pandemic from the epidemiological point of view [4], [5], [6], [7], [8], [9]. In addition, more and more evidence in epidemiology shows that environmental factors and individual movement have a significant impact on the disease transmission [4], which inspired some studies on the spatial diffusion of COVID-19. Ahmed et al. [5] introduced a reaction-diffusion COVID-19 model to show that the diffusion of individuals has a dramatic influence on the transmission dynamics and steady-state stability of the disease spread which helps suggest some control strategies. Lippold et al. [6] established a spatially resolved SIQRD model to study the first and second wave outbreak dynamics of COVID-19 in Germany, which affected the local government decisions relevant to public health through policies and regulations. Barwolff [7] presented a two-dimensional approach to describe the spatial spread of COVID-19 in Germany. Zhu and Zhu [8] formulated an epidemic model with a time delay and spatiotemporal heterogeneity to investigate the effect of self-limiting treatment on the prevention and control of COVID-19. Zhang et al. [9] introduced a COVID-19 model with diffusion and Beddington-Deangelis type incidence to study the transmission dynamics of the disease and their simulation result indicates that the diffusion behavior has an evident impact on the spread of COVID-19 in spatial heterogeneous and homogeneous environments.

A number of researchers have studied the SARS-CoV-2 infection mechanism through modeling and analysis. To et al. [10] presented an observational cohort study to investigate temporal profiles of viral load in posterior oropharyngeal saliva samples and serum antibody responses during infection by SARS-CoV-2. Du and Yuan [11] proposed a model to investigate the immune responses in SARS-CoV-2 infection. Pucci et al. [12] considered a model for the molecular impact of SARS-CoV-2 infection on the Renin-Angiotensin system. Contreras et al. [13] presented an ODE viral target model on personalized virus load curves for acute viral infections. Tavares and Gomes [14] developed a model of the natural killing (NK) cells and viruses that are involved in the infection process of the SARS-CoV-2 in human body. It is revealed that NK cells have the potential to correct the delicate immune balance needed to effectively overcome the SARS-CoV-2 infection. Néant et al. [15] established a SARS-CoV-2 viral kinetics model associated with mortality in hospitalized patients from the French COVID cohort. However, the above virus models focus on dynamical analysis and numerical simulations without the impact of the viral random diffusion on SARS-CoV-2 infection. Indeed, considerable evidence has revealed that the viral random diffusion in a heterogeneous environment is a crucial factor for studying the viral infection mechanism [16], [17].

To the best of our knowledge, the investigation of viral diffusion with spatial heterogeneity in SARS-CoV-2 infection is still at the preliminary stage. Some research shows that SARS-CoV-2 can spread from the respiratory tract to the heart, brain, and many other organs in the body within days of a patient’s infection, and can persist for months [18], which suggests that diffusion of the virus happens soon after infection and changes over time. For example, in the early stage, the virus concentration in respiratory tissues was significantly higher, but as the disease progressed, the virus concentration in pulmonary and extrapulmonary tissue categories gradually approached a similar level, which may be related to a weak immune response outside of the respiratory tract [18]. In addition, from a clinical point of view, alveolar epithelial cells migration refers to the movement of the cells after receiving a signal or sensing a concentration gradient of certain substances. During the process of migration, a cell repeats the cycle of extending the protruding foot forward and then pulling the cell body [19], [20]. In the present work, we make an effort to propose a spatially heterogeneous reaction-diffusion SARS-CoV-2 infection model with six compartments (target cells, non-productive infected cells, productive infected cells, non-infectious virus, infectious free virus and natural killing cells) to study the infection dynamics of SARS-CoV-2 within the host. Specifically, the production, infection and death rates in the model are described by more general functions.

The remainder of this paper is structured as follows. We present our model in Section 2 and prove the well-posedness in Section 3. In Section 4, we derive the basic reproduction number by using the next generation operator, and study the global stability of uninfected steady state and the uniform persistence. We discuss the existence and global attractivity of three space-independent equilibria in Section 5. Finally, numerical simulations and a conclusion are presented in Sections 6 and 7, respectively.

Model formulation

The COVID-19 pandemic caused by SARS-CoV-2 is responsible for a horrible health devastation with profoundly harmful consequences for the economic, social and political activities of communities around the world [21]. Medical studies have shown that SARS-CoV-2 predominantly infects airway, alveolar epithelial cells (AECs) and macrophages [22], [23], [24]. Different from the earlier works on mathematical modeling of SARS-CoV-2 infection, we consider the basic and core elements: target cells (AECs), infected cells evolved into an unproductive state , infected cells evolved into a productive state (, the internal virus in alveolar infected cells (, and free virus in the plasma coming from the virus in alveolar infected cells. The schematic diagram of the infection process of SARS-CoV-2 is described in Fig. 1 .

Fig. 1

Schematic diagram of the infection process of SARS-CoV-2 within the host.

To investigate the spatial diffusion of SARS-CoV-2 in a heterogeneous environment, we denote by the distribution of alveolar epithelial cells, non-productive infected epithelial cells, productive infected epithelial cells, internal virus particles in cells, and free virus particles in plasma at position and at time , respectively. Based on the schematic diagram shown in Fig. 1 we propose the following model:with the no-flux boundary condition and initial conditionwhere is a spatial bounded domain with the smooth boundary and is the unit normal vector on . Here, the term represents the spatial diffusion of and is the space-dependent diffusion rate of , . We denote by the net rate of alveolar epithelial cells and denote by the density of the newly non-productive infected cells, where and the function represent the space-dependent infection rate and general incidence function, respectively. The non-productive infected epithelial cells transform to productive infected cells at a rate . The term represents the bursting rate of . The internal viruses enter plasma at a rate . The meanings of other parameters are listed in Table 1 .

Table 1

Description of the parameters in model (2.1).

Parameter	Description
D_i(x)	The diffusion rate of the ith compartment at location x, i=1,2,3,4,5
β(x)	The infection rate of V_f(x,t) at location x
δ_I(x)	The natural death rate of A_I(x,t) at location x
δ_P(x)	The natural death rate of A_P(x,t) at location x
η(x)	The removal rate of A_I(x,t) at location x
q(x)	The burst size of A_P(x,t) at location x
μ_1(x)	The natural death rate of V_int(x,t) at location x
p(x)	The removal rate of V_int(x,t) at location x
μ_f(x)	The natural death rate of V_f(x,t) at location x

Throughout the paper, we make the following hypotheses:

and for and . There is a constant such that for and for . This implies that there also admit and such that .

Suppose that and the partial derivatives and are positive for and iff . Suppose that for and . There exists a constant such that .

Suppose that for , respectively. Suppose that hold iff , respectively, and for , for , for , for . There exist constants and such that for , for , for , and for .

Well‐posedness of system (2.1)

In this section, we consider the global existence and uniqueness of the solution of system (2.1). We denote

Consider the parabolic systemFrom Assumption 2.1 , we see that , where and .

, Lemma 3.2] System(3.1)admits a unique positive steady statewhich is globally attractive in.

Let be equipped with the supremum norm . So is an ordered Banach space. Set equipped with norm , where , , . It is known that for and , the semigroup denoted bywhere , are the Green functions associated with , , , , and , respectively. In view of the no-flux boundary condition, it follows from [25, Corollary 7.2.3] that are compact and strongly positive for . We have that there exist positive constants such that for , where are the principle eigenvalue of subject to the no-flux boundary condition.

Define , wherefor and . Let be the solution of system (2.1) with the initial values . Then system (2.1) can be rewritten as

Following [26, Corollary 4], we have

For system(2.1)with the initial value, there admits a unique non-negative mild solution. In addition,is a classical solution of system(2.1).

For system(2.1)with the initial value, letbe the solution of system(2.1)for. Then,is ultimately bounded.

Let be a non-negative solution system (2.1) for in Lemma 3.2. Suppose that . Then we have . From Assumption 2.1 , we haveFrom Lemma 3.2 and the comparison principle of parabolic equations, there exists a constant such that for . Then, it follows from Assumption 2.1 that

Consider a comparison systemIt is clear that the eigenvalue problem associated with system (3.2) admits one principle with respect to a strongly positive eigenfunction . Thus, system (3.2) has a solution , where and . Then there exists an such that , which leads to a contradiction with . Hence, . Then we obtain the global existence of . It suffices to show the ultimately boundedness of . From Lemma 3.2, it follows that is ultimately bounded for , i.e. .

SetThen, we haveAccording to the comparison principle, there are and such that for . Denote by the eigenvalue of subject to the non-flux boundary condition. Then . It follows from [27, Theorem 2.4.7] that there is a positive constant such that for . Let . From Assumption 2.1 , for we havewhich implies .

Similarly, we can obtain that and are ultimately bounded. □

Theorem 3.3 indicates that the solution semiflow is point dissipative on . From [27, Theorem 2.6] and [28, Theorem 3.4.8], we can immediately obtain

The solution semiflow of system(2.1)admits a compact global attractor.

Basic reproduction number and global stability

It follows from Lemma 3.1 that system (2.1) has an uninfected steady state , where satisfies system (3.1) and globally attractive in . Linearizing system (2.1) at the uninfected steady state , for infectious compartments and we have

Let . Then system (4.1) can be rewritten asClearly, system (4.2) is a cooperative system. Then system (4.2) has a unique principle eigenvalue with a strongly positive eigenfunction . Let be the solution semigroup associated with system (4.1) and defineThus, the distribution of total new infectious cells isBased on the definition of the next generation operator, we can obtain the basic reproduction number of system (2.1) as follows

, Theorem 3.1] The principal eigenvalueof system(4.1)has the same sign as. Moreover, the uninfected steady stateis locally asymptotically stable when. Otherwise, it is unstable.

The uninfected steady stateof system(2.1)is globally asymptotically stable when.

By the comparison principle of parabolic equations, we know that uniformly for . Without loss of generality, we suppose that , . By Lemma 3.1, there admits a such that . Then, from Assumption 2.1 we haveAssume that , where is the eigenfunction associated with the principle eigenvalue . Using the comparison principle yieldswhich implies . Furthermore, follows from the theory of asymptotically autonomous semiflows [29]. Hence, according to Lemma 4.1, we obtain the global asymptotic stability of . □

For the solutionof system(2.1)under the initial conditionand, if, then there is a constantsuch thatIn addition, system(2.1)admits at least one positive steady state.

LetThe process of proof is divided into three steps

Step 1: . If , then or or or . Without loss of generality, we suppose that . Then holds from the third equation of system (2.1) and Assumption 2.1 , and holds from the first and second equations of system (2.1) and Assumption 2.1 . Furthermore, follows from the last equation of system (2.1) and Assumption . Then, uniformly for by the theory of asymptotically autonomous semiflows.

Step 2: we show that if , then there is a constant such thatIt follows from Lemma 3.1 that there exists a such that . Suppose that it is not true, then there admits a such that . It implies that there is such that for . Then we haveLet be the eigenfunction associated with the principal eigenvalue . Assume that such that . Thenwhich leads to . This is a contradiction.

Step 3: If , then there admits a such that . Define a continuous function by . It is easy to see that and or . Then . Thus, we can verify that of in converges to and . Moreover, it is easy to see that does not form any cycle in and is isolated in . Based on Theorem 3.4 and [30, Theorem 3], there exists such thatDue to the ultimate boundedness of the solution of system (2.1) and Assumption 2.1, we can obtain that . Thus, . Set . The uniform persistence of system (2.1) follows immediately. In view of the uniform persistence of system (2.1), we can obtain that system (2.1) has at least one positive steady state according to [16, Theorem 3]. □

Reaction‐diffusion COVID‐19 infection model with humoral immunity

To study the effect of humoral immunity on the infection dynamics of COVID-19 within the host, following [14] we formulate the following system which includes active natural killer cellswith the following no-flux boundary condition and initial conditionwhere represents the distribution of active natural killer cells at position and at time , represents the killing rate of by active NK cells (), and and represents the activation rate and natural death rate of NK cells, respectively. All other model parameters are the same as given in Table 1.

for , there exists a positive constant such that for and for .

Assume that Assumptions2.1and5.1hold. Then for the initial value, system(5.1)has a nonnegativity and ultimately bounded solution for.

Assume that Assumptions2.1and5.1hold. Then there admit two threshold valuesandwithsuch that

(1) If, then system(5.1)has an infection-free equilibrium.

(2) If, then system(5.1)has a humoral inactivated equilibriumbesides.

(3) If, then system(5.1)has a humoral activated equilibriumbesidesand.

Let be a positive equilibrium. Then

From (5.2f), we have , i.e. or . If , it follows from Assumption 5.1 that . From (5.2a)-(5.2e), we have and

Assumptions 2.1 and 5.1 imply that exist, strictly increasing with . DefineThen, .

From Assumption 2.1, we have for and . Since , it follows thatwhere . Eq. (5.3) has a solution corresponding to the uninfected-free equilibrium . Denote . Based on Assumptions 2.1 and 5.1, we get

Note that . ThenSinceit follows thatFrom Assumption 2.1, we have . Hence, when and there is such that . Thus, and . Namely, the humoral inactivated equilibrium exists when .

If , then . By replacing by in Eq. (5.2a), we find . Set . Assumptions 2.1 and 5.1 imply that is strictly decreasing with and . Hence, there exists a unique such that . Then, and . From , it follows that . ThenHence, a humeral activated equilibrium exists if . In view of (5.4) and (5.5), we have  □

Let Assumptions2.1and5.1hold and. Then bothandhave the same sign as.

From Assumptions 2.1 and 5.1, we find

Suppose that . From the condition of the existence of and , we haveFrom (5.6a)–(5.6c), we see that have the same sign as , which is a contradiction. Hence, we can verify that . Note thatFrom (5.6b)–(5.6d), we know that has the same sign as . □

Suppose that Assumptions2.1and5.1hold. (1) If, then the uninfected equilibriumis globally attractive. (2) If, then the humeral inactivated equilibriumis globally attractive. (3) If, then the humeral activated equilibriumis globally attractive.

To prove Part (1), we define the Lyapunov functionalwhere . Then, we have

TakeIn view of , we have

According to the divergence theorem associated with the no-flux boundary condition, we getandFrom Assumptions 2.1 and 5.1, we can verify that . Hence, if and only if and the equality sign holds if and only if . Thus, we obtain that is globally attractive by the LaSalle invariance principle.

To prove Part (2), we define the Lyapunov functionalwhereand are given by (5.7). Then, we have

Note thatA straightforward calculation yieldsNote thatwhere for . Thus, we can verify that if and only if , and the equality holds if and only ifClearly, we can obtain the global attractivity of from the LaSalle invariance principle.

We are left to Part (3). To this end, we definewhereWe know that satisfiesThen, we haveHence, we further deriveFrom Assumption 5.1, we know that if and only if , and the equality holds if and only if and . Then, by applying the LaSalle invariance principle, we can immediately obtain that the humoral activated equilibrium is globally attractive. □

Numerical simulations

In this section, we are devoted to conducting some numerical simulations to illustrate our theoretical results. For simplicity, we set Following [16], we choose .

Dynamical evolution of system (2.1)

Take the spatial region . The values of parameters are listed in Table 3 [14] and we take the initial value of system (2.1) as follows:

Table 3

The parameters and numerical values of system (2.1).

Parameters	Value	Unit	Parameters	Value	Unit
Λ(x)	1.2×10^4	L^−1day^−1	η(x)	8	day^−1
β(x)	β(1+0.5sin(2πx))	L^−1day^−1	q(x)	400	virions
d(x)	0.14	day^−1	p(x)	0.3	day^−1
δ_I(x)	4	day^−1	μ_f(x)	10	day^−1
δ_P(x)	2	day^−1	D_2(x)	0.06	mm^2 day^−1
μ_1	10	day^−1	D_3(x)	0.05	mm^2 day^−1
D_1(x)	0.07	mm^2 day^−1	D_4(x)(D_5(x))	0.03(0.024)	mm^2 day^−1

We choose in Table 3 and calculate the value of in (4.3) by MATLAB. From Theorem 4.2, we know that the uninfected steady state is globally asymptotically stable when . In fact, it can be observed from Fig. 2 that the solution of system (2.1) converges to as time goes to infinity. We choose , then . In Fig. 3, we can observe that the virus persists within the host, which agrees well with the theoretical result described in Theorem 4.3.

Fig. 2

The dynamical evolution of system (2.1) when . It can be observed that the solutions of system (2.1) converge to as time goes to infinity, which agrees well with Theorem 4.2.

Fig. 3

The dynamical evolution of system (2.1) when . We can observe that the virus persists within the host, which agrees well with the theoretical result described in Theorem 4.3.

Impact of model parameters on for system (2.1)

To study the impact of the spatial factor on , we first set and the value of other parameters are the same as given in Table 3. In Fig. 4 (a), it can be observed that is an increasing function with respect to the parameter . In Fig. 4(b), we set with . It is shown that the influence of the parameter on the basic reproduction number is quite significant. The above simulation results indicate that the impact of spatial heterogeneity on SARS-CoV-2 infections can not be ignored during clinic treatment and medical experiment. On the other hand, we conduct some simulations to investigate the impact of diffusion rates of four infected compartments on the value of . From Fig. 4(c–f), we can see that , and have tremendous influence on the value of , among which the diffusion rate has the greatest impact on . This indicates that the random diffusion of the free virus in the plasma is a core factor in determining whether the infection outbreaks or not. However, we notice that the random diffusion does not affect the value of . This may be because this type of virus, in the productive infected cells, does not enter the plasma, so the diffusion of this type of virus has almost no effect on the infection process of SARS-CoV-2 within the host.

Fig. 4

The impact of model parameters on for system (2.1). Other model parameters are the same as given in the case of Fig. 3. (a) It can be observed that is an increasing function with respect to the parameter . (b) We set with . It is shown that the influence of the parameter on the basic reproduction number is quite significant. (c-f) We can see that , and have tremendous influence on the value of , among which the diffusion rate has the greatest impact on .

Impact of different incidence types on the viral peak time and value

Given that the incidence functions in our model are more general, in this section we discuss the impact of several types of incidence functions that satisfy Assumption 2.1 (A2) on the peak viral load and the time to reach the peak. For this purpose, we choose three incidence types, bilinear, saturated, and Beddington-DeAngelis, which are widely applied in the literature [9], [16], [17], [31]. We can clearly see from Fig. 5 (c-d) and Table 4 that, the viral peak time of ranges from 9.99 days to 11.41 days and the viral peak time of ranges from 11.26 days to 12.78 days. In fact, the average incubation period of COVID-19 is 3–7 days [33], and the peak arrival time of the virus is 5–6 days after the onset of COVID-19 symptoms [34]. This means that the viral load peaks at 8–13 days after being infected with SARS-CoV-2. Hence, our simulation results sounds reasonable. Moreover, compared with the saturation incidence and Beddington-DeAngelis incidence, the estimated peak value for the bilinear incidence is the largest and the estimated peak time for the bilinear incidence is the smallest. This indicates that the incidence type plays a crucial role in the process of SARS-CoV-2 infections.

Fig. 5

The impact of incidence type on viral peak value and time. Other model parameters are the same as given in the case of Fig. 3. Compared with the saturation incidence and Beddington-DeAngelis incidence, the estimated peak value for the bilinear incidence is the largest and the estimated peak time for the bilinear incidence is the smallest. This indicates that the incidence type plays a crucial role in the process of SARS-CoV-2 infections.

Table 4

The impact of different incidence types on the viral peak time and value. Other parameters are the same as shown in Table 3.

Incidence type	Forms	Parameter values	Peak time (days) and value (virions/mL)
			V_int(π/2,t)	V_f(π/2,t)
Bilinear incidence	βAV_f	β=2.36×10^−6	(9.99, 3.667×10^8)	(11.26, 1.091×10^8)
Saturated incidence I	βAV_f1+m_AA	β=2.36×10^−6,m_A=2×10^−6[9]	(10.73, 3.45×10^8)	(11.73, 1.026×10^8)
Saturated incidence II	βAV_f1+m_fV_f	β=2.36×10^−6,m_f=4×10^−7[9]	(10.94, 3.143×10^8)	(12.05, 9.346×10^7)
Beddington-DeAngelis	βAV_f1+n_1A+n_2V_f	n_1=1.5×10^−6,n_2=3×10^−7[9]	(11.41, 2.958×10^8)	(12.78, 8.80×10^7)

Dynamical evolution of system (5.1)

Since SARS-CoV-2 predominantly infects alveolar epithelial cells (AECs) in the lung and target cell specificity, we concentrate on the 2-dimensional region , and denote as a 2-dimensional circular disc. The initial values of system (5.1) are given byThe initial values and are the same as given in Table 3.

We choose , , , , and other parameters are the same as given in Table 3. After calculations, we obtain . From Theorem 5.5 (1), we know that the uninfected equilibrium is globally attractive. In fact, we can observe that the infected compartments of system (5.1) converge to 0 as the time goes to from Figs. 6 and 7 . We choose , , , , and other parameters are the same as given in Table 3. Then we have . We can see that system (5.1) admits a humoral inactivated equilibrium which is globally attractive in Figs. 8 and 9 . This agrees well with Theorems 5.3 (2) and 5.5 (2). We choose , , , , and other parameters are the same as given in Table 3. Then we have . From Figs. 10 and 11 , it can be observed that the solution of system (5.1) is globally attracted to the humoral activated equilibrium , which is in agreement with Theorems 5.3 (3) and 5.5 (3).

Fig. 6

Case of It is observed that the infected compartments of system (5.1) converge to 0 as the time goes to , which is in agreement with Theorem 5.5.

Fig. 7

Case of It is observed that the infected compartments of system (5.1) converge to 0 as the time goes to , which agrees well with Theorem 5.5.

Fig. 8

Case of It can be seen that system (5.1) admits a humoral inactivated equilibrium which is globally attractive.

Fig. 9

Case of It can be seen that system (5.1) admits a humoral inactivated equilibrium which is globally attractive.

Fig. 10

Case of It can be observed that the solution of system (5.1) is globally attracted to the humoral activated equilibrium , which is in agreement with Theorems 5.3 (3) and 5.5 (3).

Fig. 11

Case of It can be observed that the solution of system (5.1) is globally attracted to the humoral activated equilibrium , which is in agreement with Theorems 5.3 (3) and 5.5 (3).

Conclusion

The COVID-19 pandemic around the world has resulted in serious influences on daily life, economics, business, industry, international transportation etc., which strongly motivates scientists and mathematicians to find more effective measures so as to help treat SARS-CoV-2 infections and control the spread of COVID-19 as soon as possible. In this paper, we proposed a diffusion SARS-CoV-2 infection model with spatial heterogeneity, in which the death rates, infection rate and viral production rate are given by more general functions. We discussed and presented the basic reproduction number , the global stability of the uninfected steady state and the uniform persistence of the system. To study the impact of humoral immunity on the infection process of SARS-CoV-2 within the host, we further formulated a reaction-diffusion SARS-CoV-2 infection model with humoral immunity and some spatially independent parameters, and established the existence and global attractivity of the uninfected equilibrium, humoral inactivated equilibrium and humoral activated equilibrium which depends on two parameters and .

Compared with some earlier related and useful studies without considering the spatial factor of SARS-CoV-2 infection [35], [36], [37], [38], our results give rise to the following implications: (1) It is clear to see an important impact of spatial heterogeneity on SARS-CoV-2 infection within the host. Based on the theoretical results, as the threshold parameter of dynamical behaviors plays a crucial role in SARS-CoV-2 infections. The diffusion rates , and all impact the value of , among which has a more significant impact on (see Fig. 3). This implies that the viral diffusion and spatial heterogeneity cannot be neglected during the clinical treatment. (2) Our simulation results reveal that incidence type is an important factor in the infection process of SARS-CoV-2, which could provide a new perspective on the disease prediction and control. This is obtained by comparing the influences of three incidence functions (two saturated incidence types and one Beddington-DeAngelis incidence type) on the viral peak time and value (see Fig. 5 and Table 4). (3) When the humoral immunity is incorporated, the dynamical analysis becomes more complex even if the model contains spatial heterogeneity only in diffusion. To visualize the dynamical evolution of the solution of the system with humoral immunity, we performed some simulations on the 2-dimensional spatial domain (see Fig. 7, Fig. 8, Fig. 9, Fig. 10, Fig. 11).

There are some limitations of our models. At this stage we do not consider the effects of treatment or vaccination on the viral infection dynamics. In addition, we only studied the case where one strain of SARS-CoV-2 persists in the host body. In order to provide more informative perspective or guidance for vaccine designing and clinical treatment, it is worth incorporating vaccination or drug efficiency as well as two competing viral strains in a diffusion model. Moreover, it is also interesting to consider the above problems within fractional settings. We will leave this as a future work.

Table 2

The initial conditions (IC’s) of system (2.1).

IC’s	Value	Range	Unit	Source
ϕ_1(x)	3×10^3	[2×10^3,4×10^4]	L^−1	[32]
ϕ_2(x)	70	-	mL^−1	[13]
ϕ_3(x)	1	-	mL^−1	Assumed
ϕ_4(x)	10	-	mL^−1	[32]
ϕ_5(x)	10^6	[600,10^6]	mL^−1	[32]