Bitensorial formulation of the singularity method for Stokes flows

1.   Introduction

COVID-19 continues to devastate countries and regions throughout the world. As of August 2022, over 600 million cases and nearly 6.5 million deaths have been reported worldwide ^[1]. The ongoing pandemic and the persistence of the virus in human populations underscore the gap between our current knowledge and intervention strategies for COVID-19 and the complex mechanisms for the transmission and spread of the disease.

It is commonly accepted that the principal mode of COVID-19 transmission is the deposition of respiratory droplets through direct contact between infectious and susceptible hosts, referred to as the human-to-human transmission. On the other hand, there is increasing evidence that SARS-CoV-2 can be transmitted indirectly from the contaminated environment, referred to as the environment-to-human transmission, which includes both the airborne transmission and fomite transmission ^[1,2]. When infectious individuals speak, cough, sneeze or exhale, they release respiratory particles that contain the virus. Those particles of larger size, such as droplets, would quickly fall on nearby surfaces and objects. Those smallest particles (typically smaller than 5 $\mu$m), in contrast, could suspend and float in the air as aerosols for an extended period of time, which can be breathed in when other people pass by. In addition, individuals could catch the virus by contacting the contaminated surfaces or objects (i.e., fomites) and then touching their eyes, noses or mouths.

Environmental survival and viability of many types of coronaviruses have been investigated ^[3,4,5]. Particularly, an experimental study published in 2020 found that SARS-CoV-2 was detectable in aerosols after 3 hours, on copper after 4 hours, on cardboard after 24 hours, and on plastic and stainless steel after 3 days ^[6]. In another experimental investigation, SARS-CoV-2 retained its infectivity for a period of 15–16 hours when suspended in aerosols ^[7]. A number of biological and physical studies further confirmed the viability and persistence of this virus in the air and on inanimate surfaces ^[8,9,10,11]. These findings indicate the possibility of airborne and fomite transmission for SARS-CoV-2.

The current paper will focus on the modeling and analysis of the airborne transmission. A Lancet article provided strong arguments for the airborne transmission of SARS-CoV-2 ^[12]. A Science paper reviewed recent advances in understanding airborne transmission of several respiratory pathogens, including SARS-CoV-2 ^[13]. More recently, a Nature article reflected the change of views at WHO regarding the feasibility and risk of airborne transmission related to COVID-19 ^[14]. An experimental study using air and surface samples estimated that the probability of COVID-19 infection is about 1 per 100 exposures to SARS-CoV-2-laden aerosols through inhalation and about 1 per 100,000 exposures from contacting contaminated surfaces ^[15]. In other words, the risk of airborne transmission is about 1000 times higher than that of fomite transmission for COVID-19. Similarly, it was found through an animal study that SARS-CoV-2 transmission via fomites is markedly less efficient than airborne transmission ^[16].

At present, however, there are divided opinions and views on how often the virus spreads through the airborne transmission route, what proportion of COVID-19 infections are acquired through airborne transmission, and how much it contributes to the pandemic. As an effort toward addressing these issues, the current paper aims to use mathematical modeling to quantify the role of airborne transmission in the spread of COVID-19.

A large number of mathematical and computational models have been published to investigate the transmission, spread and progression of COVID-19 (e.g., ^{[17,18,19,20,21,22,23,24,25,26]}) as well as the effects of COVID-19 vaccination (e.g., ^[27,28,29]). Some of these models have considered the viral dynamics in the environment and the environment-to-human transmission route ^[24,25], though not specifically for the airborne transmission. A limitation of these prior modeling studies is the absence of data for the virus in the environment, as such data were not available at the time of these publications. Consequently, only human host data (cases, mortalities, vaccination coverage, etc.) were used in fitting those models, and the accuracy of their prediction on the environmental viral dynamics was unclear. Data of SARS-CoV-2 concentration in the air have been available recently, as several experimental studies have been conducted to investigate the airborne transmission path for SARS-CoV-2 and measure its concentration using various sampling, filtering, and analytical techniques ^[30,31,32]. Particularly, a longitudinal surveillance study published in 2022 ^[15] reported the concentration of SARS-CoV-2 in the air at different times at the University of Michigan (UM) and its correlation with the COVID-19 cases on the campus, providing longitudinal data for both the human host and the viral concentration in the same environment. Our present paper is mainly based on these datasets for the UM campus.

To facilitate a detailed analysis and assessment for the airborne transmission of SARS-CoV-2, we will employ an epidemic model of the simplest form, named SIRC, which includes both the human-to-human and environment-to-human transmission pathways. We will conduct mathematical analysis and numerical simulation, and will measure the contributions of the two different transmission routes to the spread of COVID-19. Furthermore, we will formulate the model in a multi-scale setting that includes the fast and slow time scales, and will carefully compare our model with the standard SIR model in both analytical and numerical means. We will utilize the datasets reported from the UM campus for our model validation and comparison.

The remainder of this paper is organized as follows. In Section 2, the mathematical formulation of our SIRC model with dual transmission routes is presented and its main dynamical properties are summarized (with details provided in the Appendices). An asymptotic analysis is then conducted to compare the SIRC model with the standard SIR model. In Section 3, model parameters are estimated using real data, and numerical simulation is conducted to verify analytical predictions and to evaluate the contribution of each transmission route. Subsequently, the relationship between the SIRC and SIR models is further analyzed and quantified, and various scenarios involving different time periods for the viral viability are numerically explored. Finally, conclusions are drawn and some additional discussion is made in Section 4.

2.   Model formulation and analysis

We utilize mathematical modeling based on differential equations to investigate the transmission dynamics of COVID-19. We employ a parsimonious model that incorporates both the human-to-human and environment-to-human (i.e., airborne) transmission routes. Our model involves only four compartments: the susceptible individuals (denoted by $S$), the infected and infectious individuals (denoted by $I$), the recovered individuals (denoted by $R$), and the coronavirus concentration in the air (denoted by $C$).

We refer to our model as the SIRC model, which is described as follows:

The parameter $\Lambda$ is the population influx rate, $\beta_I$ is the direct (human-to-human) transmission rate, $\beta_C$ is the indirect (environment-to-human) transmission rate, $\mu$ is the natural death rate for the human hosts, $\gamma$ is the rate of recovery from the infection, $w$ is the disease induced death rate, $\phi$ is the virus shedding rate by infectious individuals, and $\delta$ is the rate at which SARS-CoV-2 loses its viability in the air. We assume that all these model parameters are positive constants.

Below we summarize the standard dynamical properties of system (2.1), with details provided in the Appendices.

Basic reproduction number. The basic reproduction number $\mathcal{R}_0$ of our SIRC model is given by

where the first part in the brackets represents the contribution from the direct, human-to-human transmission route, and the second part represents the contribution from the indirect, airborne transmission route. See Appendix A for the derivation of $\mathcal{R}_0$.

Equilibrium analysis. When $\mathcal{R}_0 < 1$, the disease-free equilibrium (DFE) is the only equilibrium of system (2.1), and it is globally asymptotically stable. When $\mathcal{R}_0 > 1$, the DFE becomes unstable, and there is a unique endemic equilibrium which is globally asymptotically stable. See Appendix B for the relevant theorems and proofs.

These results show that the equilibrium dynamics of our SIRC model follow a regular forward transcritical bifurcation that is characterized by a sharp threshold at $\mathcal{R}_0 = 1$. From an epidemiological point of view, this indicates that disease control measures would have to reduce the basic reproduction number below unity in order to contain the infection. Possible reduction of $\mathcal{R}_0$ could be achieved by targeting either the human-to-human transmission route or the environment-to-human transmission route, or both.

Next, we consider the classical SIR model with population influx and vital dynamics for a comparison with our SIRC model:

System (2.3) can be obtained from system (2.1) by setting $\beta_C = 0$ and removing the virus concentration compartment. The basic reproduction number of this SIR model is given by

Clearly, when $\beta_C = 0$, $\mathcal{R}_0 = \mathcal{R}_0^{SIR}$. On the other hand, we see that when $\delta \to \infty$, $\mathcal{R}_0 \to \mathcal{R}_0^{SIR}$. The parameter $\delta$ is crucial in determining the dynamics of the coronavirus in the air and the relevant environment-to-human disease transmission. Particularly, we are interested in the scenario when $\delta$ is very large, representing a very short period of viability for the virus in the air. Mathematically, the limit $\delta \to \infty$ gives a singular perturbation to system (2.1).

To analyze this scenario, we let

be the length of the viability period for SARS-CoV-2 in the air, and examine the solution for a small $\varepsilon > 0$ using techniques of singular perturbations ^[33,34,35].

We rewrite system (2.1) as

Note that the last equation in system (2.5) is on a faster time scale than that of the first three equations, due to the multiplication of a small $\varepsilon$ to the time derivative. To reflect this difference, we introduce a fast time variable $\tau$ by

Then system (2.5) can be transformed to

We refer to system (2.5) as the slow system, since it is formulated by the slow time variable $t$. On the other hand, we refer to system (2.6) as the fast system, since it is based on the fast time variable $\tau$. These two systems have the same phase portraits for $\varepsilon > 0$, but they have different limit behaviors at $\varepsilon = 0$: the limit of the fast system describes the dynamics for a small neighborhood of time $0$, whereas the limit of the slow system describes the dynamics on a larger time interval away from $0$. We provide details below for the asymptotic analysis.

Setting $\varepsilon = 0$ in system (2.5), we obtain the slow limit system:

System (2.7) is equivalent to the SIR model (2.3). Hence, when $\varepsilon \to 0$, the dynamical behavior of the SIRC model will approach that of the SIR model as long as the time $t$ is not very close to $0$. This scenario has practical relevance as the focuses in many epidemiological applications are how an epidemic would spread after its onset and what would be its long-term progression, and those are concerned with relatively large time.

This result is consistent with our previous observation that $\mathcal{R}_0 \to \mathcal{R}^{SIR}_0$ when $\delta \to \infty$. The biological interpretation is that when the viability period of the virus in the environment is sufficiently short, the impact of the environment-to-human transmission could be disregarded for large time $t$. Consequently, the long-term dynamics of the disease transmission could be described by the standard SIR model without incorporating the $C$ compartment.

When $\varepsilon \to 0$, the dynamical behavior of the SIRC model near time $0$ is described by the limit of the fast system (2.6). Setting $\varepsilon = 0$ in (2.6) leads to the fast limit system

System (2.8) yields

Hence, in a small neighborhood of time $0$, the variables $S$, $I$ and $R$ can be treated as constants. In contrast, the solution for $C(\tau)$ exhibits a rapid change, as $C(\tau)$ quickly decreases to $0$ when $\tau$ is increasing. This is analogous to a boundary layer solution in fluid dynamics ^[36,37].

Practically, however, the value of $\delta$ is a finite number, and its quantitative effects on the disease transmission dynamics needs to be carefully evaluated. In the next section, we will compare numerical simulation results with these asymptotic predictions, and quantify the impact of the length of the viral viability period on the dynamics of the SIRC model using real data.

3.   Numerical simulation

3.1.   Data fitting

In contrast to the large body of COVID-19 infection data for human populations that have been reported by WHO, CDC and numerous other sources, data for the environmental concentration of SARS-CoV-2 are scarce and only available recently. In particular, an experimental study published in 2022 ^[15] provides longitudinal data for the virus concentration in the air at the University of Michigan (UM), with correlation to the weekly data of campus-wide COVID-19 cases reported by the UM COVID-19 Dashboard ^[38]. Hence, the work supplies longitudinal data for both the human cases and the virus concentration for a single location, the UM campus, in a relatively closed setting with a stable population. We will thus utilize the UM datasets as our main data source for fitting and simulation in the present study.

We use the COVID-19 cases data at UM from September 2020 to May 2021 for 35 weeks total. The UM COVID-19 Dashboard reported new infections on a weekly basis ^[38]. With a reasonable assumption of the average recovery period, we can then construct weekly data for the $S$, $I$ and $R$ compartments in our model by

where $t$ denotes the week number, $t = 1, \, 2, \, \cdots, \, 35$, and the constant $N$ denotes the population size at the University of Michigan. We estimate the total campus population to be around $N = 50,000$, using the enrollment statistics available on the UM website. Meanwhile, four data points for the viral concentration at UM, in terms of genome copies per liter of air (gc/L), are available from ^[15] for the week number of 5, 7, 10 and 22. We will use both the COVID-19 cases and the viral concentration data.

In order to implement our model and obtain practically meaningful results, we need to determine the values of all the model parameters. Fitting the model to real data is an important approach to estimate parameter values. However, when the number of parameters to be fitted is large, there are several potential issues for data fitting: the accuracy of fitting might be low, the results may not be biologically feasible, and parameter identifiability could be a serious concern ^[39,40]. To overcome this challenge, we apply the following multi-step procedure to parameterize our model, with the goal of minimizing the number of parameters to be fitted by data.

In the first step, we make several assumptions to simplify some of the terms in system (2.1). Since the university is a relatively closed environment not subject to large fluctuations in population within an academic year, we assume that $\Lambda = 0$. Since the time frame of our concern is relatively short, we ignore the natural death and set $\mu = 0$. Additionally, since college students are generally not facing a high risk of disease-induced deaths, we assume that $w = 0$.

In the second step, we find the values of the parameters $\gamma$ and $\delta$ from the literature. A commonly accepted value for the average recovery period of COVID-19 is two weeks ^[1,23,24,41], which leads to $\gamma = 1/2$ per week. Meanwhile, it is reported in a recent experimental study ^[6] that SARS-CoV-2 was viable in aerosols for 3 hours (or, 1/56 week), which leads to $\delta = 56$ per week.

In the third step, we conduct preliminary fitting of our SIRC model by using only the human cases data from UM, without considering the viral concentration data at all. Consequently, the model could produce reasonable results for the COVID-19 cases at UM, but generate poor predictions for the viral concentration. Nevertheless, we find that in different test scenarios the values of the indirect transmission rate $\beta_C$ and the viral shedding rate $\phi$ vary significantly, whereas the value of the direct transmission rate $\beta_I$ remains almost a constant $1.0 \times 10^{-5}$ per person per week. We also note that this value of the direct transmission rate is consistent with those in a few previous COVID modeling studies ^[20,24] after a scaling by the population size. Hence, we fix $\beta_I = 1.0 \times 10^{-5}$ per person per week, and use more careful data fitting in the next step (i.e., the fourth step) to fine tune the two remaining parameters $\beta_C$ and $\phi$.

In the fourth step, we conduct detailed data fitting to estimate $\beta_C$ and $\phi$, using both the human infection data and the viral concentration data at UM. The initial conditions are provided by $I(0) = 400$, $R(0) = 12$, $S(0) = N - I(0)-R(0)$, and $C(0) = 0.06$ gc/L, based on the UM datasets ^[15,38]. The Matlab function lsqnonlin, which employs a least-squares procedure, is used to fit the SIRC model to the data. Specifically, we fit our numerical results to the cumulative number of cases per week, calculated by $I(t)+R(t)$, and the viral concentration, calculated by $C(t)$. We set the time step size as $\Delta t = 0.01$ week in our numerical simulation. The parameter values and their 95% confidence intervals are presented in Table 1. The fitting curve for the cumulative cases is shown in Figure 1. The normalized root mean square error (NRMSE) for the fitting is $0.088$, where the NRMSE for the cumulative cases is calculated as

where $n$ is the total number of weeks, $x_{i}$ is the number of cumulative cases by week $i$, and $y_{i}$ is the model prediction of cumulative cases by week $i$.

Meanwhile, we compare the reported data with the model predictions based on the fitted parameter values for the four state variables: the numbers of the susceptible (S), infected (I) and recovered (R) individuals and the concentration (C) of SARS-CoV-2, in Figure 2. We observe generally good agreement between the simulation results and the data for the susceptible and recovered numbers. For the infected number, the data exhibit an oscillatory pattern, and our simulation curve appears to correctly catch the average behavior of the data. Generally, models with time-dependent transmission rates could better match the oscillatory infection data (see, e.g., ^[42,43]). The oscillation also takes place in the viral concentration data, though we note that the insufficiency of the concentration data, where only four data points are available, has further limited the accuracy of our simulation curve.

Next, we apply the standard SIR model (without vital dynamics and disease-induced mortality) to the same data, in order to compare with our SIRC model. To that end, we set $C = 0$ in the SIRC model; or, equivalently, set $\Lambda = w = \mu = 0$ in system (2.3), to obtain the SIR model. This model has only two parameters, $\beta_I$ and $\gamma$, where we again set $\gamma = 1/2$ per week. Then we fit the model to the number of cumulative cases to estimate $\beta_I$. The viral concentration data are not relevant to this model. The estimated value for $\beta_I$ and its 95% confidence interval are presented in Table 2. The NRMSE for this fitting is $0.088$, the same (up to three decimal places) as that for the SIRC model fitting. The fitting and simulation curves for the SIR model are presented in Figure 3. We observe a very similar pattern as that from the SIRC model predictions, except that the SIR model does not output the viral concentration.

3.2.   Model comparison

The data fitting results in Section 3.1 show that the SIRC and SIR models have very similar predictions for the numbers of susceptible, infected and recovered people, indicating that the presence of the viral aerosols in the air may be insignificant and the airborne transmission route may only play a minor role in the transmission and spread of COVID-19. We now conduct more detailed analysis and comparison of the models to verify this observation.

Our assumption of $\Lambda = w = \mu = 0$ leads to the reduced SIRC model:

The basic reproduction number $\mathcal{R}_0$ given in Eq (2.2) is not well defined for system (3.1) where $\Lambda = \mu = 0$. Nevertheless, if we consider the infected sub-system that consists of the $I$ equation (i.e., the second equation) and the $C$ equation (i.e., the last equation), and evaluate its Jacobian matrix at the initial time with $S(t) = S(0)$, we obtain

Simple algebra yields that when $S(0) < \frac{\delta \gamma}{\delta \beta_I + \phi \beta_C}$, both eigenvalues of the Jacobian matrix have negative real parts, which implies that the disease will die out. In contrast, when $S(0) > \frac{\delta \gamma}{\delta \beta_I + \phi \beta_C}$, one of the eigenvalues becomes positive, implying a disease outbreak.

Hence, the basic reproduction number for system (3.1) can be defined as

where $\mathcal{R}_I$ and $\mathcal{R}_C$ measure the contributions from the human-to-human and environment-to-human (i.e., airborne) transmission routes, respectively. Substituting the parameter values from Section 3.1, we obtain

and

The fact that $\mathcal{R}_0 > 1$ could explain the COVID-19 outbreak at the UM campus (see Appendix B for the mathematical theory). It shows that both $\mathcal{R}_I$ and $\mathcal{R}_C$ contribute to the disease outbreak, as each is lower than $1$ in magnitude but their combination exceeds $1$. Nevertheless, we clearly observe that $\mathcal{R}_I$ takes a dominant role in shaping the overall disease risk, whereas $\mathcal{R}_C$ only contributes about 2.3% to the disease risk.

Next, we quantify the incidence (i.e., force of infection) and the contribution to the infection size for each of these two transmission routes. To do that, we plot the two incidence terms $\beta_I SI$ and $\beta_C SC$ in Figure 4, where we see that the human-to-human incidence is significantly higher than the environment-to-human incidence. Meanwhile, we reduce the original model (3.1) to two simpler models: one from setting $\beta_C = 0$ and the other from setting $\beta_I = 0$, each with a single transmission mode. We then run each simpler model and generate the infected number $I(t)$ and the cumulative number $(I+R)(t)$. The results are presented in Figures 5 and 6, where we observe that the airborne transmission route is effective only for a short period of time and its contribution quickly approaches 0 after the first few weeks. This observation agrees with the asymptotic prediction in Section 2. Overall, it is clear from Figures 5 and 6 that the vast majority of the infection numbers are contributed by the human-to-human transmission route. These results are consistent with the estimates of the reproduction numbers $\mathcal{R}_I$ and $\mathcal{R}_C$ for the two transmission routes.

Thus, based on the current parameter setting, our numerical findings seem to suggest that the environmental virus concentration (i.e., the C compartment in our SIRC model) may be reasonably ignored and that the SIRC model may be replaced by an SIR model. We proceed to verify this by an asymptotic analysis of system (3.1), assuming that the value of $\delta$ is large but finite. Similar to our singular perturbation analysis in Section 2, the last equation (i.e., the equation for the environmental virus) of system (3.1) is regarded as on a fast time scale, while the other three equations (for the human hosts) are on a slow time scale. We then separate the two time scales, where a slow variable can be treated as a constant in the fast scale, and a fast variable can be approximated by its quasi-steady state in the slow scale.

Since $\delta$ is a large but finite number, the last equation implies that the virus concentration $C$ would quickly converge to its quasi-steady state:

where the slow variable $I$ is treated as a constant in this fast process. Substitute Eq (3.3) into the first two equations of system (3.1), we obtain a slow dynamics system

System (3.4) is equivalent to a standard SIR model (without vital dynamics and disease-induced mortality), where the effective transmission rate is given by

Based on the parameter values of the SIRC model (see Table 1), we obtain

which is almost identical to the fitted transmission rate of the SIR model; i.e., $1.0230 \times 10^{-5}$ (see Table 2). This implies that, in a longer term, the SIRC model may be well approximated by an SIR model, with the effective human-to-human transmission rate modified by Eq (3.5) to incorporate the impact of the airborne transmission rate.

This connection between the SIRC and SIR models is also reflected by the relationship between their respective basic reproduction numbers. The basic reproduction number for the standard SIR model with a transmission rate $\beta$ is well known and given by

Using Eq (3.5) and replacing $\beta$ by $\beta_I + \frac{\phi}{\delta} \beta_C$, we obtain Eq (3.2) for the basic reproduction number of the SIRC model (3.1).

The observations thus far can be summarized as follows: (1) The SIRC and SIR models yield similar predictions for the epidemic characteristics associated with the human population. (2) In the SIRC model, the human-to-human transmission route is dominant in shaping the disease risk, in determining the strength of the incidence, and in contributing to the infection size; the airborne transmission route appears to only play a minor role. (3) Practically, the SIRC model can be replaced by an SIR model with a modified transmission rate given by Eq (3.5).

We emphasize, however, that these conclusions are based on a specific parameter setting, particularly for the length of viability of SARS-CoV-2 in the air. Our parameter value $\delta = 56$ per week currently in use is based on an experimental measurement ^[6] that the virus is viable in aerosols for 3 hours. At present, there is no well-accepted value for the viability period of SARS-CoV-2 in aerosols, which is likely to have large variations. Hence, in the next section, we will conduct numerical simulation to investigate the impact of different values of $\delta$ on the model outcome, treating $\delta = 56$ as the base value.

3.3.   Simulated scenarios

We first consider a hypothetical scenario where $\delta$ is increased by 10-fold from the base value; i.e., $\delta = 560$. This means that the time for the persistence of SARS-CoV-2 in the air is reduced by a factor of 10. Keeping other parameter values unchanged, we conduct simulation to this scenario and present the results in Figure 7. As can be naturally expected, the infection level $I(t)$ and the viral concentration $C(t)$ are both decreased with the increased value of $\delta$. The reduction of the infection level, though, is inconsequential, given that the parameter $\delta$ impacts $I(t)$ in an indirect manner and that the base value $\delta = 56$ is already sufficiently large. On the other hand, the reduction of the viral concentration is more significant since $C(t)$ directly depends on the parameter $\delta$.

Next, we consider scenarios where $\delta$ is decreased from the base value. In a different experimental study ^[7], it is reported that SARS-CoV-2 can remain infectious in aerosols for 15 hours (or, 1/11.2 week). This would lead to $\delta = 11.2$ per week, which represents a reduction from the base value by a factor of 5. The simulation results for this scenario are presented in Figure 8. We observe a substantial growth for the viral concentration $C(t)$, where the peak value of the curve is increased by 7 times in comparison to that in the base scenario. Meanwhile, the infection level $I(t)$ also exhibits a significant rise, with the peak value increased from around 400 in the base scenario to around 680 in the current scenario. Moreover, simple calculations yield the reproduction numbers associated with the current setting:

Thus, airborne transmission would account for more than 10% of the total transmission risk. These results show that with the reduced value of $\delta = 11.2$, the airborne transmission route becomes significantly more important in shaping the transmission risk, the infection size and other epidemic characteristics, suggesting that the viral concentration compartment $C(t)$ may become indispensable toward accurate predictions of the epidemic. They also indicate that the SIR model may not be a good approximation to the SIRC model any more. This can be seen from the perturbation analysis in Section 3.2, which relies on the assumption that $\delta$ is a large number. When $\delta$ becomes much smaller, the accuracy of the asymptotic properties derived in Eqs (3.3)–(3.5) would be considerably reduced.

Biologically, such a finding is feasible, as the longer the virus remains infectious in the air, the higher the risk is for susceptible individuals to contract the infection through the airborne transmission route. With a sufficiently long viability period (corresponding to a sufficiently small $\delta$), the airborne transmission route could even play a major role in the transmission and spread of the disease. Figure 9 displays the simulation results for another scenario where $\delta$ is decreased by a factor of 10 from the base value; i.e., $\delta = 5.6$. We clearly see more substantial increases in the infection level $I(t)$ and the viral concentration $C(t)$, compared to those in Figure 8.

4.   Discussion

We have presented a parsimonious mathematical model, named SIRC, to study the airborne transmission of COVID-19. The model combines the human-to-human and environment-to-human transmission routes in a simple form, which makes it convenient to evaluate the relative contribution of each transmission pathway to the overall disease risk. We have also conducted an asymptotic analysis by separating the fast and slow time scales, and compared the dynamics with those of the SIR model. Our model validation, simulation and comparison are based on the datasets ^[15,38] reported at the UM campus that provide longitudinal data for both the human cases and the virus concentration in the air.

Our findings indicate that the role of airborne transmission of SARS-CoV-2 strongly depends on the length of time the virus remains infectious in aerosols. Mathematically, the dynamics of the SIRC model would approach that of the classical SIR model when the parameter $\delta \to \infty$; i.e., when the length of the viral viability period approaches zero. In a practical sense, however, $\delta$ can only take a finite value. Based on a single estimate of 3 hours for the period of viral viability ^[6], our numerical results show that the airborne transmission route plays a minor role, contributing to approximately $2.3\%$ of the total transmission risk. Similarly, the airborne transmission route appears to be insignificant in shaping the strength of the incidence and the size of the infection. Consistent with this, our simulation shows that the contribution of the airborne transmission pathway mainly happens during an initial period of time that does not last long. Consequently, the main difference between the SIRC and SIR models is only reflected during a short temporal interval, away from which the outcomes of the two models may not be distinguishable. This has been justified through a perturbation analysis in Section 3.2. From the data fitting point of view, the estimated transmission rate in the SIR model has already incorporated both the actual human-to-human transmission rate and the (minor) impact of the airborne transmission rate from the SIRC model. Thus, under this parameter setting with $\delta = 56$ per week, the SIRC model may be practically replaced by the standard SIR model without degrading the accuracy of the prediction.

On the other hand, the length of time that SARS-CoV-2 remains stable and infectious in the air may depend on several factors such as temperature, humidity, ultraviolet radiation, and air flow conditions ^[13,44], and may vary considerably from place to place or from time to time. As such, there is no well-accepted value for $\delta$ at present. If the potency of the virus lasts longer than a few hours, then airborne transmission could become substantially more significant, as demonstrated in our simulation scenarios in Section 3.3. In those scenarios, the SIRC and SIR models would exhibit significant differences in their outcome, and the viral concentration compartment would become indispensable in order to accurately predict the spread and progression of COVID-19.

We remark that our SIRC model, in some sense, can be regarded as a special case of some more complex models incorporating multiple transmission pathways, such as those proposed in ^[24,25]. Our simple model, however, allows a detailed mathematical analysis that completely resolves the local and global stabilities and the multi-scale dynamics, and enables a careful comparison with the standard SIR model to quantify the role of the airborne transmission for COVID-19. These findings may not possible for those more sophisticated models.

We believe that our incorporation of the SARS-CoV-2 concentration data into model fitting and simulation is new and has not been conducted before in the modeling of the transmission dynamics of COVID-19. This effort helps to quantify the role of airborne transmission in a realistic setting and promote more holistic understanding of the multiple transmission pathways for COVID-19. On the other hand, we acknowledge that the inadequacy of the concentration data has limited the accuracy of our curve fitting for the viral concentration function $C(t)$, which might subsequently impact the accuracy of some other model predictions. Future modeling efforts for COVID-19 would benefit from the availability of a sufficient amount of data for both the human population and the viral concentration in the same environment. The modeling efforts would also benefit from a deeper understanding in experimental studies regarding the viability of SARS-CoV-2 in the air.

Our modeling study shows that the probability of contracting SARS-CoV-2 from aerosols could be considerably high if the value of $\delta$ is low. Even with a large value of $\delta$, there is still a certain level of risk for airborne transmission of COVID-19, especially during a relatively short period of time immediately following the shedding of the virus into the air by infectious individuals. These results support public health guidance focused on preventions and interventions to reduce airborne transmission of SARS-CoV-2, including social distancing, mask wearing, and increasing indoor ventilation and air filtration.

Acknowledgments

We acknowledge partial support from the NSF REU program under grant number 1852288 which provided an opportunity for the first three authors (MH, BM and WY) to conduct undergraduate research at UTC. We also acknowledge partial support for JW from NSF DMS under grant number 1951345. We thank the two anonymous reviewers for their helpful comments that have improved the quality of the original manuscript.

Conflict of interest

The authors declare that there is no conflict of interest.

Appendix

A.   Basic reproduction number

All parameters and variables of system (2.1) are nonnegative. Adding the first three equations, we obtain

Using the fact $I \leq \Lambda/\mu$, we have

Hence, the biologically feasible domain

is invariant for system (2.1).

Using the stand next-generation matrix technique ^[45], we can obtain the matrix $F$ representing the new infections and the matrix $V$ representing the transfer of infections:

The next-generation matrix can then be computed as

The basic reproduction number $\mathcal{R}_0$, defined as the spectral radius of the matrix $FV^{-1}$, is thus given by

B.   Equilibrium analysis

Through simple algebraic manipulation, we obtain that system (2.1) has two possible equilibrium points, $\textbf{x}_{0}$ and $\textbf{x}^{*}$, where $\textbf{x}_{0}$ is referred to as the disease free equilibrium (DFE) and $\textbf{x}^{*}$ the endemic equilibrium (EE). They are given below:

and

Theorem B.1. If $\mathcal{R}_{0} < 1$, then ${\bf x}_{0}$ is the only equilibrium point in $\Omega$ and is locally asymptotically stable. If $\mathcal{R}_{0} > 1$, then both ${\bf x}_{0} \in \Omega$ and ${\bf x}^{*} \in \Omega$, and ${\bf x}_{0}$ is unstable.

Proof. For any $\mathcal{R}_{0} > 0$, it is clear that $\textbf{x}_{0} \in \Omega$. Meanwhile, we may re-write the EE as

Obviously, $\textbf{x}^{*} > 0$ if and only if $\mathcal{R}_{0} > 1$. Now suppose $\mathcal{R}_{0} > 1$. We have

Since $\mathcal{R}_{0} > 1$, we obtain that $\mu \delta (\gamma + w + \mu) < \Lambda (\delta \beta_{I}+\phi \beta_{C})$. Thus,

This shows that $\textbf{x}^{*} \in \Omega$ when $\mathcal{R}_{0} > 1$. In addition, the local asymptotic stability of $\textbf{x}_{0}$ when $\mathcal{R}_{0} < 1$ and instability when $\mathcal{R}_{0} > 1$ follow directly from Theorem 2 of ^[45].

Next, we analyze the global stability properties of the DFE and the EE.

Theorem B.2. When $\mathcal{R}_{0} < 1$, the DFE ${\bf x}_{0}$ of system $(2.1)$ is globally asymptotically stable in $\Omega$.

Proof. Clearly, in system (2.1), $S$ and $R$ are uninfected compartments, and $I$ and $C$ are infected compartments. Define the two subsystems below:

Note that $\mathcal{G}(0, 0, S, R) = 0$. Based on the global stability result by Castillo-Chavez et al. ^[46], we only need to verify that the following two conditions hold when $\mathcal{R}_{0} < 1$:

$({\rm{H1}})$ For the subsystem $\mathcal{F}(0, 0, S, R), \; \, \textbf{x}_0$ is globally asymptotically stable;

$({\rm{H2}})$ $\mathcal{G}(\textbf{x}) = A[I, \, C]^T-\hat{G}(\textbf{x})$, where the off-diagonal elements of the Jacobian matrix $A = D\mathcal{G}(0, 0, S_0, R_0)$ are non-negative, and $\hat{G}(\textbf{x})\geq0$ for $\textbf{x}\in \Omega$.

Since

condition (H1) can be easily verified. Next, we have

which yields

and

Hence, condition (H2) is also satisfied. This would establish the global asymptotic stability of the DFE ${\bf x}_0$.

Theorem B.3. When $\mathcal{R}_{0} > 1$, the EE ${\bf x}^{*}$ of system $(2.1)$ is globally asymptotically stable in $\Omega$.

Proof. Define a Lyapunov function of conventional construction ^[47]

which has gradient

Taking the dot product of $\nabla V$ along the trajectories of system (2.1), we obtain

Note that the endemic equilibrium ${\bf x}^{*}$ satisfies the system below:

Using these equations, we can write

By collecting all like terms, we obtain

due to the inequality of arithmetic and geometric means, and we have equality only when $(I, C, S, R) = (I^{*}, C^{*}, S^{*}, R^{*})$. Hence, the endemic equilibrium is globally asymptotically stable.