A Multiscale Model of COVID-19 Dynamics.

COVID-19, caused by the infection of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), has been a global pandemic and created unprecedented public health challenges throughout the world. Despite significant progresses in understanding the disease pathogenesis and progression, the epidemiological triad of pathogen, host, and environment remains unclear. In this paper, we develop a multiscale model to study the coupled within-host and between-host dynamics of COVID-19. The model includes multiple transmission routes (both human-to-human and environment-to-human) and connects multiple scales (both the population and individual levels). A detailed analysis on the local and global dynamics of the fast system, slow system and full system shows that rich dynamics, including both forward and backward bifurcations, emerge with the coupling of viral infection and epidemiological models. Model fitting to both virological and epidemiological data facilitates the evaluation of the influence of a few infection characteristics and antiviral treatment on the spread of the disease. Our work underlines the potential role that the environment can play in the transmission of COVID-19. Antiviral treatment of infected individuals can delay but cannot prevent the emergence of disease outbreaks. These results highlight the implementation of comprehensive intervention measures such as social distancing and wearing masks that aim to stop airborne transmission, combined with surface disinfection and hand hygiene that can prevent environmental transmission. The model also provides a multiscale modeling framework to study other infectious diseases when the environment can serve as a reservoir of pathogens.

Introduction

In late December 2019, highly contagious pneumonia of unknown etiology was first reported in Wuhan, China (Zhu et al. 2020). A novel strain of coronavirus was isolated from patients and later named severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2). Since then, the Coronavirus Disease 2019 (COVID-19) has spread to over 210 countries and territories, creating unprecedented public health challenges throughout the world. As of January 12, 2021, more than 91 million cases and 1.9 million deaths had been reported. Despite many theoretical, experimental and clinical studies, our current knowledge on the fundamental mechanisms of COVID-19 transmission and infection remains limited.

Mathematical modeling provides a powerful theoretical means to study COVID-19. A large number of mathematical and statistical models have been proposed (e.g., Kucharski et al. 2020; He et al. 2020; Li et al. 2020; Liu et al. 2020a, c; Tang et al. 2020; Weitz et al. 2020; Wu et al. 2020a, b; Zhao and Feng 2020; Hellewell et al. 2020; Chen et al. 2020; Olabode et al. 2021; Musa et al. 2022). Almost all these models are concerned with the transmission and spread of the disease (i.e., the between-host dynamics) at the population level. On the other hand, very little modeling effort has been devoted to the within-host dynamics of SARS-CoV-2. When the coronavirus enters the human body, there are complicated interactions between the pathogen and host cells taking place, which directly shapes the disease risk and infection severity for the individual hosts and which, in a collective manner, may subsequently impact the epidemic patterns (Liu et al. 2020b; He et al. 2020). Due to the inadequate investigation of the within-host dynamics of COVID-19 and their connection to the population-level epidemics, there are several fundamental questions that remain unanswered or only partially answered at present; for example, how does the viral load change inside the human body, what are the short-term and long-term interactions between SARS-CoV-2 and the host cells, and how does the within-host pathogen development affect the population-level disease transmission and spread?

As a pilot study to address these questions, we develop a multiscale modeling framework in this paper to investigate the between-host and within-host dynamics of COVID-19 and their impact on each other. At the individual host level, our model characterizes the virus–cell interactions and their time evolution within the human body. At the population level, our model describes the disease transmission and spread through multiple transmission routes. Most of the between-host COVID-19 models published thus far are based on the susceptible-exposed-infected-recovered compartmental framework or its variants, with a focus on the direct, human-to-human transmission pathway (Chan et al. 2020). On the other hand, a recent experimental study found convincing evidences that SARS-CoV-2 was detectable in aerosols for up to 3 h, on copper for up to 4 h, on cardboard for up to 24 h, and on plastic and stainless steel for up to 3 days (van Doremalen et al. 2020). There are also consistent and strong evidences that the infection can spread via airborne transmission (Greenhalgh et al. 2021). All these indicate a significant risk of the indirect, environment-to-human transmission pathway, particularly airborne and fomite transmission, for SARS-CoV-2. Additionally, the novel coronavirus has been found in the stool of some infected individuals (Zhang et al. 2020), which may contaminate the aquatic environment through the sewage water and add another possible route of environmental transmission for COVID-19 (Yeo et al. 2020). Therefore, quantifying such indirect transmission routes in a modeling study could help us to better understand the transmission mechanisms of COVID-19.

Our between-host model explicitly includes the environment-to-human transmission by incorporating the concentration of SARS-CoV-2 in the environment, which interacts with the human hosts: susceptible individuals may be infected by contracting the environmental coronavirus, and infected individuals may shed the pathogen (through coughing, sneezing, etc.) back to the environment. In addition, we introduce the concentration of SARS-CoV-2 inside the human body as a within-host variable that interacts with the host cells and that represents the individual viral load. We then bridge the within-host and between-host dynamics via a two-way coupling. Since the coronavirus may enter the human body through the environment, we assume that the within-host viral load depends on the environmental pathogen concentration. Meanwhile, since the within-host pathogen level is directly associated with the individual symptoms and infection severity, we assume that the human-to-human transmission rates depend on the viral load within the human body.

Hence, our modeling framework incorporates multiple transmission routes (both human-to-human and environment-to-human pathways) and connects multiple scales (both the population and individual levels). The within-host interactions normally occur on the time scale of hours to days, which is referred to as the fast dynamics. In contrast, the between-host transmission and spread is on the scale of weeks, months to years, referred to as the slow dynamics. To facilitate the mathematical analysis, we will first separate the two time scales so that a thorough investigation can be conducted to each of the fast and slow systems. Then we will combine the between-host and within-host systems and study the coupled dynamics. Our detailed analysis on the local and global dynamics of the fast system, the slow system and the coupled system shows that rich dynamics, including both forward and backward bifurcations, emerge with the coupling of the within-host and between-host models.

The remainder of the paper is organized as follows. Section 2 is devoted to the formulation of the COVID-19 model linking within-host and between-host dynamics, and the derivation of the basic reproduction number for the model. We analyze the model by using the bifurcation theory and fast-slow analysis in Sect. 3, and conduct numerical simulation and model validation in Sect. 4. Finally, we conclude the paper in Sect. 5 with some discussions.

Model

Let S, E, I, R denote the number of susceptible, pre-symptomatic infected, symptomatic infected and recovered host individuals. Let Z and V be the concentration of coronavirus in the environment and within the host, respectively. T and denote the concentration of target cells and infected target cells within the host, respectively. To link the within-host and between-host interactions, we extend the population-level model proposed in Yang and Wang (2020) with the inclusion of the within-host dynamics (Wang et al. 2020), which takes the formThe symptomatic infected class has fully developed disease symptoms and can infect others. The pre-symptomatic infected class is in the incubation period; COVID-19 patients do not show symptoms but are still capable of infecting others. There are some other models that include the exposed (assumed to be non-infectious), asymptomatic and symptomatic groups explicitly (Zhao and Feng 2020; Ngonghala et al. 2020; Xue et al. 2020; Tang et al. 2020). Because we will couple it with within-host models, we keep the between-host model with a minimum number of variables.

The parameter represents the generation rate of susceptible individuals, is the natural death rate, is the period of incubation between infection and the onset of symptoms, is the disease-induced death rate, is the rate of recovery, and are the respective viral release rates to the environmental reservoir from the pre-symptomatic infected and symptomatic infected individuals, and is the viral removal rate from the environment. The functions and represent the direct, human-to-human transmission rates between pre-symptomatic infected and susceptible individuals, and between symptomatic infected and susceptible individuals, respectively. The function is the indirect, environment-to-human transmission rate. We assume that , and are all non-increasing functions of E, I, and Z, respectively, because higher levels of E, I and Z would motivate stronger control measures to prevent transmission. Higher viral load can result in a higher transmission rate. Thus, we assume that and increase as V increases. For example, in the simple case, , and can be constant, or they may take the nonlinear form (44) as shown in numerical simulations. For mathematical analysis of the multiscale model, we avoid the complexity of keeping track of the dynamics within each individual and only formulate the within-host model for a conceptual “average” individual. The within-host model will be fitted to the average viral load of a group of patients. The parameter b is the generation rate of target cells, is the viral infection rate, d is the death rate of target cells, q is the death rate of infected cells, p is the viral production rate, and c is the viral clearance rate. The environment can transmit virus to the individual at a rate . Below we summarize the assumptions on the functions. The coupled dynamics have two distinct time scales. The within-host interactions normally occur on the time scale of hours to days, referred to as the fast dynamics. In contrast, the between-host transmission and spread is on the scale of weeks, months to years, referred to as the slow dynamics. These two scales are coupled by a small constant .

, and are positive functions.

, for , and .

and are both concave downward (i.e., the Hessian matrices of and are negative semidefinite).

, and for , and for .

Model (1) has a disease-free equilibrium solution (DES), given by The new infection matrix and transmission matrix areandBy the next-generation matrix method (van den Driessche and Watmough 2002), the basic reproduction number is defined as the spectral radius of , i.e.,whereLet . Here and are the between-host and within-host threshold parameters, respectively. More specifically, is the sum of three terms, where (resp. , ) measures the contribution to the between-host basic reproduction number from direct pre-symptomatic-infected-to-susceptible human-to-human transmission (resp. direct symptomatic-infected-to-susceptible human-to-human transmission, indirect environment-to-human transmission).

Analysis

Slow System

For the slow time scale t, our system can be studied using the full model (1) by setting , which is referred to as the slow system. In analyzing the slow system, the fast system is treated at its quasi-steady state.

Since , the within-host dynamics can be studied byLetting in (5) leads toand V is determined bywhereIf , it leads to the virus-infection-free ES or the coexistence state of virus and infected target cellsIf , it follows from (H4) that , and hence and . In this case, equation (7) has two zeros with , and is the only biological feasible equilibrium solution. The infected target cells and viral loads both persist and the corresponding ES isDifferentiating (7) with respect to Z along yieldsThis implies that along as and Inequality (11) shows that V is a strictly increasing function of Z. Similarly, differentiating (10) with respect to Z along and using give usIt follows from (6), (11) and (H4) that , and . This implies that .

Accordingly, the slow system can be written as

Equilibrium Solutions

We will analyze the equilibrium solutions (ESs) of the slow system (13). Note that the associated ES satisfiesSolving (14), we obtainwithWhen , it leads to the disease-free ES .

When , substituting (15) into the second equation of (14) and solving S as a function of I at the equilibrium, we havewhere withStraightforward calculation yieldsandBy (H3), the Hessian matrix of is negative semidefinite, and hence the first three terms in the bracket of the right-hand side of (18) are non-positive. By (H2) and (12), and . Thus, Likewise, one can verify that and . Hence, , and this implies thatThis shows that the intersections of the curves and in determine the nontrivial equilibrium solutions of the slow system (13), as the nontrivial equilibrium solution is of the formwhere is a positive intersection of these two curves.

Since and , we haveAdditionally, is a strictly decreasing line and is concave upward. To study the existence of the nontrivial equilibrium of (13), we consider three cases. This result implies that system (13) has at most three biologically feasible equilibria, i.e., if , this system has two equilibria: disease-free equilibrium (DFE) and endemic equilibrium (EE); If , it can have one, two or three equilibria (and one of them is the DFE) depending on the parameter values. This indicates that the system (13) may undergo a backward bifurcation.

. In this case, , these two curves have a unique intersection in , denoted as .

. In this case, . In view of (16), . When , these two curves do not intersect; when , there exists a unique intersection .

. In this case, and there are three possibilities for the intersection of these two curves in . (i) There is no intersection if for ; (ii) They have one intersection in if there exists such that and ; (c) Otherwise, there are two intersections in .

Existence of Forward and Backward Bifurcations

In this section, we study the bifurcation in terms of I with as a bifurcation parameter. Note that is a multiple of (i.e., with ). We treat as an independent parameter and the rest parameters fixed. It suffices to analyze the bifurcation diagram of I with as a bifurcation parameter. If , then the corresponding . At the equilibrium of system (13), when , it leads to the disease-free equilibrium (DFE) and it follows from van den Driessche and Watmough (2002)[Theorem 2] that the DFE is locally asymptotically stable (resp. unstable) when , i.e., (resp. , i.e., ).

In the remainder of this section, we assume that and intersect in the first quadrant (i.e., system (13) has equilibrium solutions with ). Thus, at the equilibrium,To compute , we calculate . Differentiating (21) with respect to I yieldsBy the implicit function theorem, if . Assume that as . Integrating (22) yieldsIn view of (19), we know as a function of I is positive and strictly concave upward.

Case I: If (i.e., ), then and hence and is a positive and strictly concave downward function for . In this case, system (13) exhibits a transcritical bifurcation at (i.e., ).

Case II: If (i.e., ), there exists a unique such that , when and when . Denote and the corresponding by . Clearly, and . Moreover, has two branches, denote by and when (i.e., ), where (a) , (b) and , (c) and for , and (d) . Specifically, is a strictly decreasing and concave downward function, whereas is a strictly increasing and concave downward function as . On the other hand, has only one branch and is well-defined, denote as , when (i.e., ). Besides, , when and for . So we rename by . In this case, system (13) exhibits a backward bifurcation.

The result for the occurrence of transcritical and backward bifurcations is summarized in the following lemma and illustrated in Fig. 1.

Fig. 1

Illustration of bifurcation diagram of I as a function of

Lemma 1

System (13) undergoes a transcritical bifurcation at (i.e., ), which further induces a backward bifurcation if , and a forward bifurcation if .

Stability Analysis for Forward Bifurcation

Suppose . As , , and , this yields a biologically feasible domain of system (13), which is given byIt is easy to see that is positively invariant for system (13).

Assume that . By Lemma 1, system (13) exhibits a forward bifurcation. Together with the result in Sect. 3.1.1 leads to the following: In what follows, a global stability analysis is performed to study the forward bifurcation. The associated result is established in the following two theorems.

If , system (13) admits a unique equilibrium solution, the DFE , that is locally asymptomatically stable when .

If , system (13) has two equilibria: the DFE and the endemic equilibrium (EE) . Moreover, the DFE is unstable.

Theorem 1

Suppose that

Assume that either or is satisfied. If , then the DFE of system (13) is globally asymptomatically stable in .

If , the disease is uniformly persistent in the interior of .

Proof

When and , the DFE is the unique equilibrium of (13) and it is locally asymptotically stable. If either or holds, by (H1), (H3) and (12), we havewhere when , which is obtained from (11). Let . This implies thatwhereOne can verify that and , where . Consider a Lyapunov functionDifferentiating L along solutions of system (13) leads toIf , . and if and only if .1 This implies . Using the first and the fourth equations of (13) yields and . Thus, the largest invariant set where contains only one point, i.e., the DFE . In view of LaSalle’s Invariant Principle (Salle 1976), the DFE is globally asymptotically stable in when .

If , we know that is unstable. It follows from the persistence theory (Thieme 1993) and a similar argument as in the proof of Gao and Ruan (2011)[Proposition 3.3] that instability of implies uniform persistence of system (13) in the interior of .

Let denote the unique endemic equilibrium when . To establish the global stability of system (13) in this case, we need to make some assumptions: for . By the similar argument as that in Yang and Wang (2020)[Theorem 2.2], we obtain the following result. A poof is given in Appendix A for completeness.

;

;

Theorem 2

Suppose that and (C1)-(C3) holds. If , the endemic equilibrium is globally asymptotically stable in the interior of .

Stability Analysis for Backward Bifurcation

In this section, we assume that . In this case, slow system (13) may exhibit a backward bifurcation according to Lemma 1. By the similar argument as that in Sect. 3.1.3, we have the following result. As shown in Fig. 1(a), the bifurcation diagram of I as a function of consists of three branches: the top, middle and bottom branches, where the bottom (resp. middle, top) branch is composed of the DFE (resp. the EE , the EE ). Thus, it remains to investigate the stability of the middle and top branches in the case where .

If and exists, system (13) admits up to three equilibria. The DFE always exists and it is locally asymptomatically stable. Two EEs, and , coexist when , and they merge into one when . There is no positive EE when .

If , system (13) has two equilibria: the DFE and the EE . Moreover, the DFE is unstable and the EE is globally as asymptomatically stable in the interior of when and holds, where the global stability of the EE is established by the same method as that in Theorem 2.

Denote

Theorem 3

The middle branch of the equilibrium solutions of system (13) in the backward bifurcation diagram is unstable.

As the fourth equation (i.e., the equation for R) is decoupled from the rest of system (13), it suffices to use the following system to prove the instability of the EE Linearizing system (25) at leads to the corresponding Jacobian matrixIt follows from the direct calculation thatwhere m, and are defined in (17). Note that, at , for andDifferentiating (28) with respect to I yieldsBy (28), (29) and (16), can be written asClearly the sign of is determined by the sign ofOn the other hand, at the EE, , i.e.,which implies thatDifferentiating both sides of (32) in terms of and doing a simple algebraic manipulation, we findSince and at , . In view of (30) and (31), at . This shows that is unstable.

Theorem 4

Suppose (i) is sufficiently small, and (ii) at . Then the top branch of the equilibrium solutions of system (13) in the backward bifurcation is locally asymptotically stable.

First, we assume that . In this case, it is easy to verify that the characteristic equation for eigenvalues of the linearized system of (25) at (i.e., an equilibrium solution on the top branch) iswhere is the identity matrix, is defined in (26) andOne eigenvalue is . To prove the remaining three eigenvalues have negative real parts, by the Routh–Hurwitz criterion, it suffices to show thatIt is obvious that and Additionally, note that , andThus, at Hence . This implies .

Along the top branch of the equilibrium solutions of system (13), . Using the same argument as in the proof of Theorem 3, one can verify that at . By , we have .

On the other hand, it follows from direct calculation thatIn view of , equation (34) and the positiveness of all the model parameters involved, we have , which implies that as . This shows that is locally asymptotically stable when .

By the continuous dependence of the spectrum on the model parameters, the local stability of the top branch of the equilibrium solutions of system (13) remains to be true when is sufficiently small.

Fast System

In the fast time scale , the full system (1) can be written asLetting yieldswhich is referred to the fast system.

If , the fast system (36) exhibits a forward bifurcation (Smith and Leenheer 2003). It admits at most two biologically feasible equilibria. The virus-infection-free equilibrium always exists, and the coexistence of virus and infected target cells is determined by a threshold parameter . More specifically, If , recalling (9) and the analysis in Sect. 3.1, the fast system (36) has a unique endemic equilibriumAccordingly, the global stability of the endemic equilibrium is summarized as follows.

if , there is no coexistence state, and the virus-infection-free equilibrium is globally asymptotically stable in .

if , the coexistence state is globally asymptotically stable in .

Theorem 5

For each , the endemic equilibrium of the fast system (36) is globally asymptotically stable in .

Proof

Consider a Lyapunov functionwhere the constant is to be determined. It is easy to verify that for and if and only if .

For simplicity, we write by . Differentiating along the solution of the fast system (36) givesUse the equilibrium equations of (36) , and Equation (38) can be simplified asSetting and using yieldIt follows immediately from the arithmetic-geometric mean inequality that and if and only if . Thus, is globally asymptotically stable.

The Full System

The time scale separation allows us to conduct a fast-slow analysis to analyze the dynamics of the slow and fast system, which is useful to gain insights into the dynamics of the full system.

It is easy to see from the results in Sects. 3.1 and 3.2 that the full system (1) admits at most four biologically reasonable solutions. Let denote the DFE, which always exists. Letdenote the boundary equilibrium solutions (bES).

In the case of (i.e., the slow system exhibits a forward bifurcation), the slow system has a unique EE . Hence, the full system admits the bES and the endemic equilibrium (EE)with

In the case of (i.e., the slow system exhibits a backward bifurcation), the slow system (13) allows up to two EEs, denoted aswith and t representing the middle and top branches of equilibrium solutions of the slow system, respectively. Thus, the full system admits the bES and at most two EEsfor

The local stability of the DFE and bES is an immediate consequence of the results from Sects. 3.1 and 3.2, as the full system would be decoupled in those cases. The global stability of the DFE and EE is established in Theorems 6 and 7. The associated existence and stability of equilibrium solutions are summarized in Tables 1 and 2, for which DNE means “does not exist,” l.a.s. and g.a.s. are locally asymptotically stable and globally asymptotically stable, respectively, and superscript asterisk means that the result holds under certain condition and superscript n indicates that the result is verified numerically.

Table 1

Steady-state stability classification of the full system when

Region	DFE E_0	bES E_b	EE E_e
R_0=max{R_0w,R_0b}<1	l.a.s.	DNE	DNE
R_0b<1<R_0w	Unstable	l.a.s.	DNE
R_0w<1<R_0b	Unstable	DNE	g.a.s.^∗
min{R_0w,R_0b}>1	Unstable	Unstable	g.a.s.^∗

Table 2

Steady-state stability classification of the full system when

Region	DFE E_0	bES E_b	EE E_em	EE E_et
R_0b<R0bc,R_0w<1	l.a.s.	DNE	DNE	DNE
R0bc<R_0b<1,R_0w<1	l.a.s.	DNE	Unstable	l.a.s.^n
R_0w<1<R_0b	Unstable	DNE	DNE	g.a.s.^∗
R_0b<R0bc,R_0w>1	Unstable	l.a.s.	DNE	DNE
R0bc<R_0b<1,R_0w>1	Unstable	l.a.s.	Unstable	l.a.s.^∗
min{R_0w,R_0b}>1	Unstable	Unstable	DNE	g.a.s.^∗

Let which exists by (H4), and . It is easy to verify that is a positive invariant region for the full system (1). The global stability of the DFE of our full system is established in the following result.

Theorem 6

Suppose that . If , then the DFE of the full system (1) is globally asymptotically stable in .

Note that, for all and , where the first and second inequalities are obtained from (H2) and (H3), respectively, and the last equality follows from our assumption . Similarly, one can verify that for and . Additionally, using , we have for .

Denote . It follows from direct calculation that, for all ,where and are defined in (2) and (3), respectively. Since both and are nonnegative, it follows from the Perron–Frobenius Theorem that the non-negative matrix has a non-negative left eigenvector associated with the eigenvalue .

Let us consider the Lyapunov function Differentiate along solutions of system (1), we haveBy , . Moreover, if and only if . Since at least one entry of w is positive, if , then The second equation in (40) implies that . Using the third equation in (40), we have . By , and hence the fourth and fifth equations imply the only solution to is . In view of the first, fourth and sixth equations of the full system (1), , and . Similarly, if for some , one can verify that implies , and . This shows that the largest invariant set where is the singleton . By LaSalle Invariance Principle, the DFE is globally asymptotically stable in when .

To establish the global stability of the EE , denoted by , we further assume that, for ,

;

;

, which is the same as (C3);

.

Theorem 7

Let denote the interior of and suppose that assumptions (A1)-(A4) hold. If , then the EE of the full system is globally asymptotically stable in .

For simplicity, we denote , , , and . Consider the following Lyapunov functionwith , , , and , where () are positive constants to be be determined later and . It is easy to see that and if and only if for .

Differentiating (42) along the solution of the full system (1) yieldsLet for and . By the similar arguments as that in the proof of Theorems 2 and 5, we find thatBy (A1)-(A4), (H2) and the arithmetic-geometric mean inequality,Setit follows from (11) that with the chosen , and , the right-hand side of (43) is zero. This shows that in . Additionally, it is easy to verify that if and only if for , which implies that the largest invariant set where is the singleton . It follows from the LaSalle Invariance Principle (Salle 1976) that is globally asymptotically stable in .

Model Validation and Simulation

In this section, we fit the multiscale model (1) to both the within-host and between-host data of COVID-19. On the basis of the fitting and parameter estimation, we evaluate the influence of a few infection characteristics and treatment on the disease dynamics. The viral load data we used to fit are the mean viral load measured by reverse transcription quantitative polymerase chain reaction (RT-qPCR) from samples of posterior oropharyngeal saliva of 23 patients admitted to the Princess Margaret hospital and Queen Mary hospital in Hong Kong between Jan 22 and Feb 12, 2020 (To et al. 2020). The data of COVID-19 cases are from the Johns Hopkins University Coronavirus Resource Center, which tracks COVID-19 cases through a map-based dashboard and is updated multiple times per day. We used the data of confirmed cases in Florida between March 13th and August 17th, 2020 (https://coronavirus.jhu.edu/region/us/florida).

Data Fitting

In model (1), we chose the functions , , and as below. The choice of these functions was inspired by Yang and Wang (2020). We assumed that they have a base transmission rate, increase as the viral load increases, and decrease as the level of infection increases.The initial conditions of the variables are set to S(0)=, , , , , cells/ml (Wang et al. 2020), , copies/ml. We estimated parameters , , , , , p, c, by fitting the model to the viral load data and COVID-19 case data simultaneously. The other parameters are fixed according to either the literature or our best estimates. The viral shedding rate by pre-symptomatic infected persons was estimated to be per day (Yang and Wang 2020). The viral load is higher during the early stage of infection (Wang et al. 2020). Thus, we assumed the viral shedding rate by symptomatic infected persons to be smaller and chose per day. The values of and are small (Yang and Wang 2020) and chosen to be on the order of to . Because the turnover of the target cells of SARS-CoV-2 infection is slow (Wang et al. 2020), we assumed that both the generation and death rate of target cells (b and d) are 0.01 per day. The constant in the transmission function is fixed to and the influence of its variation on the dynamics will be investigated later.

We fit the model to both the viral load and confirmed case data simultaneously using the R programming language. The root mean square (RMS) between the model prediction and the data is calculated and parameter estimates are based on the best fitting that achieves the minimum RMS. Figure 2 shows that the model provides a good fit to both the within-host and between-host data. The fitted parameters and other fixed parameters are listed in Table 3. It should be noted that the estimate of the scale-adjusting parameter might be different for variants such as delta and omicron, compared with the wild type. However, the difference in the time scales between viral strains should be much smaller than that between the within-host and between-host scales. We also note that some non-pharmaceutical control measures have been implemented during the period of data collection. These measures were not explicitly included in the model but their effect may have been reflected in the transmission rate between hosts.

Fig. 2

Best fit of the multiscale model to the mean viral load data (left panel) of Hong Kong patients (To et al. 2020) and COVID-19 cases in Florida (right panel) between Mar 13th and August 17th, 2020 (https://coronavirus.jhu.edu/region/us/florida) (color figure online)

Table 3

Model parameters and values

Parameter	Description	Value	Refs
C_0	Constant in the transmission β_E and β_I	10^-2	Fitting
C_1	Constant in the transmission β_E	10^7	Fitting
C_2	Constant in the transmission β_I	74	Fitting
η_0	Viral transmission rate from	5.28×10^-6 day^-1	Fitting
	environment to host
κ	Viral infection rate	2×10^-5 day^-1	Fitting
p	Viral production rate	12 day^-1	Fitting
c	Viral clearance rate	0.3 day^-1	Fitting
ϵ	Scale parameter	5.28×10^-4	Fitting
β_E_0	Constant in the transmission β_E	4×10^-7	See text
β_I_0	Constant in the transmission β_I	10^-6	See text
β_Z_0	Constant in the transmission β_Z	10^-6	See text
Λ	Population growth rate (in Florida)	557.13 day^-1	http://edr.state.fl.us/content/population-demographics/reports/econographicnews-2018-v2.pdf
μ	Natural death rate (in Florida)	2.78565×10^-5 day^-1	http://edr.state.fl.us/content/population-demographics/reports/econographicnews-2018-v2.pdf
α	Transition rate from pre-symptomatic infected to symptomatic infected	1/7 day^-1	Spencer et al. (2020)
ω	Disease-induced death rate	0.016 day^-1	https://www.worldometers.info/coronavirus/usa/florida
γ	Recovery rate	0.096 day^-1	https://www.worldometers.info/coronavirus/usa/florida
ξ_E	Viral shedding by pre-symptomatic infected persons	2.3 day^-1	Yang and Wang (2020)
ξ_I	Viral shedding by symptomatic infected persons	1.5 day^-1	See text
C_3	Constant in the transmission β_Z	10^3	See text
δ	Rate of viral clearance in environment	1 day^-1	Geller et al. (2012)
d	Death rate of target cells	0.01 day^-1	See text
b	Generation rate of target cells	0.01 day^-1	See text
q	Death rate of infected cells	2 day^-1	Wang et al. (2020)

The Influence of Environment, Incubation and Treatment

On the basis of the data fitting and parameter estimates, we performed numerical simulations to investigate the influence of a few factors, including the environment, incubation, and potential antiviral treatment, on the infection dynamics within individuals and the spread of the disease at the population level. Recall that represents the contribution to the viral load within symptomatic infected individuals from the environment. We simulated the model with different values of . The predicted dynamics are very sensitive to the change of (Fig. 3). For example, when increases slightly from 0.01 to 0.012, the peaks of symptomatic infected I, pre-symptomatic infected E and environmental virus Z increase significantly and the time to reach the peaks is largely shortened (Fig. 3). In all the cases, the disease persists at the population level and the viral load within symptomatic infected individuals declines to a very low level (i.e., ) but does not go to zero due to the contribution from environmental transmission. Without environmental transmission (i.e., ), the viral load will decline to zero quickly ( Fig. 7).

Fig. 3

Simulation of the coupled model with different values of . The other parameters are from Table 3. In the right bottom figure of the viral load V, the blue dash line represents the detection limit (color figure online)

Fig. 7

Simulation of within-host viral dynamics with different values of . When , there is no environmental transmission. All the other parameters are the same as those in Fig. 3 (color figure online)

It should be noted that the final states shown in Fig. 3 are not the steady states. The simulation over a longer time interval reaching the steady state is shown in Appendix Fig. 8. We further simulated the model with different values of (the constant in the environmental transmission function ). As increases to (i.e., the value we fixed in the data fitting Fig. 2 and simulation in Fig. 3), the time for the system to stabilize increases (Appendix Fig. 9). This explains why under the parameter estimates from data fitting it takes a long time for the variables to reach the steady states after some oscillations, as shown in Appendix Fig. 8.

Fig. 8

Simulation of the coupled model with different values of over a longer time interval. The simulation over a shorter time interval is shown in Fig. 3 (color figure online)

Fig. 9

Model simulation with different values of (the constant in the environmental transmission function ). We fixed at 0.01 and the other parameters are the same as those in Fig. 3 (color figure online)

The duration of incubation is an important clinical characteristic in symptomatic infected disease surveillance, prevention and control (Lessler et al. 2009). Some studies showed that the incubation duration ranges from 1.8 to 12.8 days for COVID-19 (Leung 2020; Ki 2020; Jiang et al. 2020). In Fig. 4, we simulated model (1) using different values of to investigate how the incubation duration (i.e., ) influences the dynamics of COVID-19. As the duration of incubation increases, more infections are observed and the time to reach the epidemic peak is also shortened. This is not surprising as COVID-19 patients can transmit SARS-CoV-2 during the incubation period. This result highlights the need of identifying infected but still asymptomatic people in the early stage of infection.

Fig. 4

Model simulation with different duration of incubation (i.e., ). The other parameters are from Table 3 (color figure online)

Several drugs that target different aspects of COVID-19 pathogenesis have been proposed. Some have been approved by the FDA while others are currently being tested in clinical trials of different stages (Sanders et al. 2020). In our model, we assumed that when COVID-19 patients are treated with antiviral therapies, the viral production rate p is reduced to , where is the treatment efficacy. The simulation in Fig. 5 shows that antiviral treatment can significantly delay the spread of the disease among population. However, it does not change the final steady state. This is because the viral load declines to a very low level (see Appendix Fig. 8) even without treatment (i.e., the within-host basic reproduction number with the parameter values obtained from data fitting). “Perfect” therapy that can completely block viral replication cannot reduce the duration of viral persistence within symptomatic infected individuals or diminish the time to recovery (see in Appendix Fig. 10). This is different from the prediction by a model that only considered the within-host dynamics (Wang et al. 2020). The discrepancy is due to the viral transfer from the environment to the host in the multiscale model. If the symptomatic infected individual is isolated after infection, then the individual does not have a chance to contract environmental virus. In this scenario, antiviral treatment will accelerate viral eradication and diminish the recovery time in the symptomatic infected individual (Wang et al. 2020). In summary, our simulation shows that antiviral treatment only delays the occurrence of the disease outbreak. It would not control the disease if the environmental transmission to hosts cannot be prevented.

Fig. 5

Simulation with different drug efficacy that blocks viral production within symptomatic infected individuals. The other parameters are from Table 3 (color figure online)

Fig. 10

Simulation of viral load dynamics with different drug efficacies . The other parameters are the same as those in Fig. 5 (color figure online)

Bistability

Table 2 shows that the model has two steady states, and , when , and . Although it is challenging to perform a formal mathematical analysis on the stability of the steady state in this case, numerical simulation with different initial conditions suggests that is locally asymptotically stable (Appendix Fig. 11). Furthermore, using initial conditions of the susceptible and symptomatic infected as an example, we identified the region in which the model converges to the disease-free steady state or the endemic steady state (Fig. 6). This shows that the multiscale model exhibits a phenomenon of bistability. With lower levels of initial susceptible and symptomatic infected people, the disease is predicted to die out; when either the number of susceptible or symptomatic infected people exceeds a threshold value, the disease will persist.

Fig. 11

Model convergence to the endemic steady state under differential initial conditions. a The steady state of susceptible S is 565.31, consistent with the calculation ; b the steady state of infected I is 2449.68; c the steady state of E is 2140.61, consistent with the calculation ; d the steady state of R is 422.21, consistent with the calculation ; e the steady state of Z is 8383.85, consistent with the calculation ; f the steady state of viral load is 2794.67, consistent with the calculation ; g the steady state of T is 0.31 cells/ml, consistent with the calculation ; h the steady state of infected cells is 0.0035 cells/ml, consistent with the calculation (color figure online)

Fig. 6

Bistability region of the initial susceptible and symptomatic infected in which the model converges to either the steady state or . The parameters are from Table 3 (color figure online)

Discussion

A large number of mathematical models have been developed to study the dynamics of COVID-19 (for example, see Kucharski et al. 2020; Li et al. 2020; Liu et al. 2020c; Tang et al. 2020; Weitz et al. 2020; Wu et al. 2020a; Zhao and Feng 2020). However, most of the models are focused on the disease transmission at the population level. Several models have been used to study the within-host virus dynamics (Wang et al. 2020; Hernandez-Vargas and Velasco-Hernandez 2020; Gonçalves et al. 2020; Sanche et al. 2020; Goyal et al. 2020). We have not seen any models developed to investigate the coupled dynamics of COVID-19. In this paper, we developed a multiscale model to study the dynamics of COVID-19 at both the within-host and between-host levels. In addition to the assumption that the transmission rate between hosts depends on the viral load within symptomatic infected individuals in some coupled models for other infectious diseases (Tuncer et al. 2016; Childs et al. 2019; Martcheva et al. 2015; Barfield et al. 2018; Nikin-Beers et al. 2018; Dorratoltaj et al. 2017; Numfor et al. 2014; Feng et al. 2015, 2012, 2013; Cen et al. 2014), we introduce a function into the equation of viral load to represent the inhalation/ingestion rate of the coronavirus from the environment into the human body. This represents the contribution of the environmental reservoirs to the growth of the viral load within individual hosts. Therefore, our model provides a two-way connection between the individual viral kinetics and the population-level disease transmission. In comparison with the above-mentioned coupled models, the transmission rate between hosts in this model was assumed to also depend on the level of infection because a higher level of infection would lead to stronger control measures and consequently reduce the transmission rate. The model is also different from that in Feng et al. (2015, 2013), in which the authors studied a coupled model specifically for an environmentally-driven infectious disease. There was no infection between infected and susceptible individuals considered in that model.

Our analysis approach for the subsystem (i.e., slow or fast subsystem) is similar to that in previous papers including (Yang and Wang 2020). Mathematical analysis of the coupled model is more challenging. We have analyzed the model using the bifurcation theory and a fast-slow analysis, where the latter approach is built upon the time scale separation for the dynamical processes associated with the individual host and the environment and human population. Our detailed analysis on the local and global dynamics of the fast system, the slow system and the coupled system shows that rich dynamics, including both forward and backward bifurcations, emerge with the coupling of the within-host and between-host models. Numerical simulation illustrated the bifurcations (Fig. 1). We have also numerically identified a region in which the disease either persists or dies out (Fig. 6). The emergence of bistability further suggests that the initial populations of susceptible or symptomatic infected can dictate the fate of disease transmission. These bifurcation and bistability results indicate the significance of the proposed coupled modeling approach and highlights the challenges in the prevention and control of the ongoing pandemic.

Both mathematical analysis and the simulation in Fig. 3 suggest that the environmental transmission can contribute to the COVID-19 outbreak and persistence. Potential antiviral treatment can delay the disease outbreak but cannot prevent its emergence. Even when treatment significantly suppresses the viral production within the infected individual, the virus can persist in the host (Fig. 5). This may sound bizarre for a specific individual but it actually makes sense for the conceptual “average” individual considered in the model. If the disease persists at the population level, then the virus will also persist within the average individual due to environmental transmission. Although environmental transmission can contribute to viral persistence within the average host when the disease persists at the population level, it does not mean that environmental transmission plays a dominant role in the disease spread. It depends on the relative magnitude of the transmission terms , , and . Most spread of COVID-19 might go through airborne human-to-human transmission (Zhang et al. 2020). Therefore, comprehensive intervention measures such as social distancing and wearing masks, combined with surface disinfection and hand hygiene that can prevent environmental transmission, should be implemented to mitigate the spread of COVID-19 (Pradhan et al. 2020). Our study provides a modeling framework that can be used to evaluate the potential influence of environmental transmission on the disease dynamics.

This is the first attempt to link within-host and between-host dynamics of COVID-19. The validation of multiscale models is usually very challenging as it requires data from all the scales for the same individual or the same region. In this paper, we used the mean viral load of the patients from Hong Kong to represent the within-host dynamics of symptomatic infected individuals and the accumulated cases from Florida to represent the between-host dynamics. Although patients in Florida may have different within-host dynamics, we can use the above as an example to study the mutual impact of the dynamics at one level on the other. Another limitation of our multiscale model is that it assumes individual hosts have the same internal states at a time (i.e., the within-host model is for a conceptual “average” individual). To overcome this limitation, one can couple the within-host and between-host dynamics via the infection age of individuals using the nested modeling approach (Gilchrist and Sasaki 2002; Martcheva et al. 2015), which allows incorporating the staged progression nature of the disease at the population level. There are also some other potential ways to account for the individual within-host heterogeneity in the model, such as adding a subscript i on the within-host variables T, , and V for each individual i or discretizing the population with the level of the viral load as in a recent study (Lin et al. 2021). A more complex model incorporating all the individual heterogeneity will be extremely challenging, if not impossible, for mathematical analysis of the model. Lastly, we used a simple model at each scale to make it analytically more trackable. Immune response was shown to play a significant role in controlling viral replication within infected individuals (Wang et al. 2020; Goyal et al. 2020). The within-host model can be extended to include immune responses but we speculate that it may not have a significant impact on the disease spread at the population level (note that the viral load is suppressed to a very low level in our simulations). A between-host model including more compartments, spatial heterogeneity, stochasticity, and vaccination status would be more realistic to study the disease spread. In a recent work (Lin et al. 2021), Lin et al. used a region-specific model to study the COVID dynamics in a number of populous metropolitan statistical areas in the USA.

To summarize, we have developed a multiscale model to study the interaction between within-host viral replication, the environment, and the disease spread at the population level. The analysis has generated some insights that would not be obtained if the model only considers within-host or between-host dynamics. This work also provides a modeling framework for studying other infectious diseases when the environment can serve as a reservoir of pathogens.