SARS-CoV-2 Dynamics in the Mucus Layer of the Human Upper Respiratory Tract Based on Host–Cell Dynamics

Abstract

A thorough understanding of the inhalation dynamics of infectious aerosols indoors and infection dynamics within the host by inhaled viruses such as severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) plays an important role in the assessment and control of infection risks indoors. Here, by combining computational fluid–particle dynamics (CFPD) and host–cell dynamics (HCD), SARS-CoV-2 infection dynamics in the mucus layer of the human upper airway were studied. To reproduce the diffusive and convective transport of the virus in the nasal cavity–nasopharynx by mucociliary motion, a three-dimensional (3D)-shell model with a mucus layer was developed. The initial virus concentrations for HCD calculation were estimated based on the deposition distribution of droplets with representative sizes analyzed by CFPD. To develop a new HCD model, the target-cell-limited model was integrated with the convection–diffusion equation. Additionally, the sensitivity of the infection rate β to the infection dynamics was systematically investigated. The results showed that the time series of SARS-CoV-2 concentration in the mucus layer strongly depended on diffusion, convection, and β. Although the SARS-CoV-2 dynamics obtained here have not been verified by corresponding clinical data, they can preliminarily reveal its transmission mode in the upper airway, which will contribute to the prevention and treatment of coronavirus disease 2019.

1. Introduction

The spread of the coronavirus disease 2019 (COVID-19), which has beset the people worldwide for nearly two years as of February 2022, remains an ongoing issue. To date, more than 0.44 billion transmission cases have been confirmed, including over 5.9 million deaths [1]. The severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), a novel coronavirus with a diameter of approximately 0.1 μm, plays a critical role in causing pneumonia [2,3].

To better understand the infection dynamics of SARS-CoV-2, many researchers have conducted extensive clinical case studies and evaluated its transmission and the pathogenesis of COVID-19. Among them, several studies [4,5,6,7,8] have assessed the viral load from multiple sampling sites (nasopharyngeal, sputum, blood, urine, and stool) of patients in different regions, ages, and stages for several days and analyzed the epidemiological characteristics and clinical course of COVID-19 to provide a clear understanding of its natural history. In addition, some experts [9,10,11] focused on the dynamic analysis of SARS-CoV-2 in the saliva of patients with a consistently high viral load. Therefore, saliva has the potential to be a non-invasive specimen for SARS-CoV-2 diagnosis and viral load monitoring, which can reduce the risk of nosocomial transmission. However, further investigation is needed to elucidate the correlation between salivary viral load and disease severity. Unlike previous studies, Zou et al. [12] observed the changes in the viral load in the upper respiratory tract of 17 patients using nasal and pharyngeal swabs and detected high viral loads shortly after symptom onset, with higher viral loads being detected in the nose than in the throat; however, comprehensive data were needed to determine transmission dynamics. He et al. [13] reported the temporal patterns of viral shedding in 94 confirmed COVID-19 patients and simulated the characteristics of SARS-CoV-2 infectivity in 77 infector–infectee transmission pairs. The findings indicated that the viral load was the highest in the throat swabs at the time of symptom onset and inferred that the infectious peak time was at or before symptom onset; hence, there could be a lot of potential for pre-symptomatic transmission. Similar results were reported by To et al. [14], but their findings also underscored the importance of strict infection control and early use of powerful antiviral drugs in high-risk groups. Kissler et al. [15] measured viral RNA trajectories in 68 patients. Among these, the peak viral concentration and duration of the viral proliferation and clearance phases were inferred in 46 acutely infected patients. The results showed that symptom severity did not affect the rapid peak of SARS-CoV-2 concentration. Kim et al. [16] indicated that the viral load dynamics of SARS-CoV-2 may differ from those of previously reported coronavirus infections, such as SARS-CoV, after studying the viral load dynamics in the first two confirmed patients in South Korea. Severe acute respiratory syndrome coronavirus 2 can infect multiple cell types, and then replicate in them, including epithelial cells in the respiratory tract [17]. In particular, there is evidence that unlike SARS-CoV, SARS-CoV-2 can replicate in the upper respiratory tract, as reported by Wölfel et al. [18]. They performed detailed virological analyses of nine cases, confirming that the virus was actively replicating in the larynx. These findings have important implications for the prevention of infection. In addition, the viral loads of several symptomatic and asymptomatic patients with nasopharyngeal swabs showed that nasal epithelial cells were the initial point of infection and transmission [19].

To accurately predict the evolution of viral load, researchers in the field of mathematical or theoretical biology have proposed reasonable methods, that is, the host–cell dynamics (HCD) model. Hernandez-Vargas et al. [17] first studied the human–host mathematical model based on a target-cell-limited model for SARS-CoV-2. The results showed that the viral replication number of SARS-CoV-2 in the host was consistent with the broad value of human influenza infection. In addition, a complex model that considered immune cell responses suggested slow immune response peaks at 5 to 10 d after symptom onset. Fatehi et al. [20] presented an intracellular dynamics model of COVID-19 based on random subjects. They tried to combine these predictions with an intercellular model of intra-host infection dynamics to fit clinical data [18], and then to obtain a profile of disease progression in recovered patients without treatment. A multi-scale semi-mechanistic model was developed by Dogra et al. [21] to predict the systemic viral distribution dynamics and viral excretion of SARS-CoV-2 by combining human cell-scale viral dynamics with physiological processes related to viral transport and innate and adaptive immune responses; the model was well calibrated with published in vivo and clinical data. However, the limitation of this model was that it could not be verified clinically with the prediction of extrapulmonary compartments. Sadria et al. [22] combined the adaptive immunity term to the dynamic model of the innate immune response to influenza A virus infection and developed a new model represented by SARS-CoV-2 specific antibody to study SARS-CoV-2 dynamics. Although the clinical relevance of the analysis and conclusions might be limited by the simplification presented in the model and the continuing lack of characterization of the nature and immune response of SARS-CoV-2, this did not prevent it from describing the complex dynamics of SARS-CoV-2 in systemic circulation, as well as the innate and adaptive immune responses of the host. In addition, some scholars have evaluated the effects of drugs on viral load dynamics, based on studies of viral dynamics and in conjunction with other disciplines. Gonçalves et al. [23] matched a mathematical model of viral dynamics with in vivo data to estimate the parameters driving viral replication. Their study combined pharmacokinetic/pharmacodynamic and viral dynamics modeling to predict the effects of lopinavir/ritonavir, hydroxychloroquine, interferon (IFN)-β-1 A, and remdesivir on viral load dynamics. Nevertheless, Czuppon et al. [24] studied SARS-CoV-2 dynamics using a stochastic analog of a standard target-cell-limited model, which considered the virus eclipse period and whether it was a noninfectious virion, and found that precise values of the basic reproduction number R₀ and the burst size N within the host were critical for predicting the outcome. In addition, they used a more sophisticated model to evaluate the effects of different drugs on the early stages of viral infection and found that the maximum delay in viral infection was achieved by reducing viral production in infected cells. Both studies agreed that early intervention with drugs was necessary to inhibit viral load growth.

The transmission mechanism of this virus is changeable and intricate, and aerosol transmission is commonly considered as an important and hot topic in the current discussion. Droplets carrying the virus spread when a patient sneezes and coughs. Of these, larger particles settle naturally, while small particles can remain suspended in the air for a long time. Therefore, five detailed case studies and an airborne physical model were used by Prentiss et al. [25] to estimate the characteristic number of SARS-CoV-2 required to induce infection in each case and to provide a basis for estimating the risk of daily activities. If these virus-laden droplets are inhaled by others, some will be deposited in the upper respiratory tract. Similar to the influenza A virus (IAV) studied by Haghnegahdar et al. [26], airborne transmission of infectious materials from infected people is often considered the decisive mechanism of transmission. They developed and used an experimentally validated computational fluid and particle dynamics (CFPD-HCD) model to simulate the transport and immune system response of inhaled IAV-loaded droplets in the human respiratory system. In addition to considering the droplet interaction with water vapor, a one-way coupled Euler–Lagrange-format CFPD model was used to simulate the aerodynamics of virus-carrying droplets in the respiratory system with transient breathing patterns. The simulation results of the local deposition patterns could be transferred to the HCD model as an input to determine the regional host dynamic response. The results showed that most of the droplets in the case of nasal inhalation were concentrated in the nasal cavity, with a few concentrated in the pharynx, but in the case of oral inhalation, the droplets were mainly deposited in the oral cavity and pharynx. Accordingly, the activated dendritic cells had a similar distribution.

Nevertheless, it is important to note that when virus-laden droplets enter the upper respiratory tract, they are deposited in the mucus layer on the airway surface. A series of activities, such as convection, diffusion, and apoptosis in the mucus layer before the virus reaches the target cells, are often ignored in viral load calculations. It is impractical and unethical to experiment with viruses in the human respiratory tract. At this point, the computational fluid dynamics (CFD) simulations should be considered. Given that the mucus layer is thin compared to the size of the airway structure, modeling and calculating it is not easy. Some scholars have tried a few methods to simulate mucus motion, especially in the field of intranasal drug delivery, which is worth learning from. For the nasal cavity, Rygg et al. [27] proposed a method to transform the nasal cavity into a surface-based two-dimensional (2D) model, introducing mucus uniformly into the surface computational domain, delivering it to the nasopharynx, and performing CFD simulations. On this basis, Shang et al. [28] developed two new models: the unwrapped-surface model and three-dimensional (3D)-shell surface model using the surface unwrapping technique that preserves the nasal topology. Comparing the two models, the shell model was recommended because the results of the calculated mucus flow were more consistent with the experimental data.

As mentioned above, in this study, a combination of CFPD and HCD with a 3D shell human airway model was proposed to better characterize the infection dynamics of SARS-CoV-2 along the mucus layer of the upper respiratory tract.

2. Materials and Methods

2.1. SARS-CoV-2 Dynamics Prediction Based on the Target-Cell-Limited Model

This study was based on a 3D respiratory tract model created from the CT data of a healthy human by using Mimics (Materialise NV), as shown in Figure 1. The geometric model was divided into six sites, including the nasal vestibule, central nasal passage, nasopharynx, trachea, mouth cavity, and lower respiratory tract. Moreover, its surface mesh and volume mesh were created by using ANSYS/Spaceclaim for CFPD model. Assuming the worst-case scenario, in which a healthy person who has a breathing rate of 7.5 L/min (gentle breathing state) is in a close conversation with an infected person, 10,000 droplets coughed up by the infected person are assumed to be directly inhaled by the other person through nasal inhalation. Therefore, the mouth cavity and lower respiratory tract were neglected in this study.

Based on previously reported studies [29,30], the distribution of droplets released from a human cough is summarized and plotted in Figure 2. It can be clearly seen that more than 70% of the droplets were less than 10μm in diameter; therefore, five particle sizes (1, 2.5, 5, 7.5, and 10 μm) were selected as representatives consistent with the previous study [31] for the Lagrangian particle tracking analysis. Evaporation or condensation was not considered in this study. The regional and total deposition results were calculated according to the proportion in Figure 2 and detailed information is shown in

Table 1. Basic initial conditions for the HCD model.

Regions		Vestibule	Central Nasal Passages	Nasopharynx-Larynx	Trachea	Total
Surface area (m2)		3.71 × 10−3	1.55 × 10−2	9.88 × 10−3	4.72 × 10−3	3.38 × 10−2
Vmucus (mL)		0.06	0.23	0.15	0.07	0.51
T(0) [cells]		9.27 × 107	3.88 × 108	2.47 × 108	1.18 × 108	8.46 × 108
Particles	1 μm	0	8	5	1	14
	2.5 μm	0	16	10	2	28
	5 μm	0	14	14	18	46
	7.5 μm	0	42	49	72	163
	10 μm	0	88	71	84	243
Percentage		0.00%	1.68%	1.49%	1.77%	4.94%
V(0) (copies/mL)		0	5.64 × 10−3	4.90 × 10−3	6.11 × 10−3	1.66 × 10−2

.

The deposition of SARS-CoV-2-laden droplets in the upper respiratory tract is shown in Figure 3. This visual result was obtained by scaling and simply summing up the deposition data for the five particle sizes referring to their proportion of droplet size distribution. These deposition distribution results of SARS-CoV-2-laden droplets in the upper respiratory tract shown in Figure 3 are based on our previous reported CFPD analysis [31]. These droplets were mainly concentrated around the nasal cavity and throat, which may be responsible for nasal congestion and a sore throat in the early stages of infection. As can be seen from the percentage of deposited particles in Table 1, more than 90% of the fine particles entering the lower portion may cause serious symptoms.

In this study, SARS-CoV-2 dynamics were preliminarily calculated and described using the target-cell-limited model in HCD [17], which has been verified to fit the trend of clinical data well. The formulas used are as follows:

$$\frac{dT}{dt}=-\beta T{V}^{ }$$

(1)

$$\frac{dI}{dt}=\beta TV-\delta I$$

(2)

$$\frac{dV}{dt}=\frac{p^{\prime }}{V_{mucus}}I-cV$$

(3)

where T, I, and V represent the number of susceptible target cells, infected cells, and viral load, respectively. Referring to the initial conditions and parameter values of the validated model, infection occurred 3 d before symptom onset was assumed. For the initial conditions, the initial number of target cells, T(0), was estimated according to the regional surface area of the airway and the surface area of each epithelial cell that was selected as 4 × 10⁻¹¹ m²/cell [32]. Moreover, the initial virus load, V(0), was calculated from the droplet sizes and the percentage of deposition, as well as the amount of virus carried by each droplet, assumed to be 10⁵ copies/mL, which depended greatly on the infection stages of the individuals. For the parameters, according to the mean values of the fitting parameters in Esteban’s study [17], which more closely fit the clinical data, that is, the onset of infection is 3 d before the onset of symptoms, the infection rate of the target cells, β, was assumed to be 4.71 × 10⁻⁸ (copies/mL) day⁻¹, and the apoptosis of infected cells occurred at the rate of δ = 1.07 day⁻¹. Meanwhile, the virus can be cleared at a rate of c = 2.4 day⁻¹. Nevertheless, Equation (3) was modified slightly from the original model, and the rate of virus production in each part of the respiratory tract was calculated separately to replace the uniform rate p of virus production ((copies/mL)/day/cell), making the results closer to the actual situation. Hence, the virus production rate was set as p’ = 0.74 copies/day/cell, and the volume of mucus, V_mucus, was calculated as the product of surface area and thickness, which was set as 15 μm [28]. The results of this model were compared with the clinical data of nine patients, and Wölfel [18] verified its feasibility.

2.2. Coupling of Target-Cell Limited Model and Convection–Diffusion Model

The upper airway is the first entrance for inhaled infectious droplets to enter the human body. The deposition fraction of droplets in the upper airway, including the nasal cavity–nasopharynx, is significant. The SARS-CoV-2 dynamics in the human airway are closely related to the mucus layer that attaches to epithelial cells. Mucus is composed of two parts: the upper layer is a high-viscosity gel layer composed of 97% water and 3% other substances (mucin, salts, lipids, and immunoglobulins, etc.), with a thickness of 0.5~5 μm, and the lower layer is a 7~10-μm thick periciliary fluid layer infiltrating the cilia with viscosity similar to that of water [11,33,34]. The primary function of the mucus layer is to remove unwanted particles such as inhaled bacteria and viruses by rhythmic beating of cilia on the epithelium, which is called mucociliary clearance [34]. To highlight viral changes, as shown in Figure 4, the mucus layer was simplified to a single layer.

The model was divided into the upper and lower parts of the pharynx because mucociliary transport forms a flow from the entire airway to the pharynx and needs to be calculated separately. Here, the infection dynamics in the upper part, including the nasal cavity and nasopharynx, were studied. As shown in Figure 5, a 3D-shell model with a 15-μm thick mucus layer was developed. Meanwhile, the surface area was approximately 0.021 m², so the mucus volume was calculated as 0.32 mL. Here, a source term was added to the continuity equation for the mucus domain [35], as follows:

$$S_m=\rho \left(\nabla ·\overrightarrow{u}\right)$$

(4)

The density ρ and viscosity μ were assumed as 1000 kg/m³ and 12 Pa·s, respectively. The nasopharynx, which is the site of the posterior nose, was assumed to be the exit of mucus produced in the nasal cavity. The mucus flow velocity at the pharynx region was assumed to be 10 mm/min [28], which determined the value of S_m. The outlets of the two nostrils were set as the pressure outlet at 0 Pa. Although two layers of mucus were simplified into one layer for calculation, and the specified shear stress of the mucus surface was set as 0, the contact surface of mucus and respiratory epithelium was set as the no-slip boundary, according to the characteristics of the cilia layer. Based on the droplet deposition results, appropriate sites were selected to simulate the local viral load with HCD and CFPD for 50 d in the milieu with mucus properties. In addition, mucus was assumed to be secreted by goblet cells in the epithelial cells [36] (sinus secretion was ignored here).

The evolution of viral load (V) was calculated by coupling the target-cell-limited model with the convection and diffusion terms. The parameter values such as T(0), V(0), β, δ, c, p’, and V_mucus were set to be the same as those in Section 2.1. The improved formulas are as follows:

$$\frac{\partial \overline{V}}{\partial t}+\frac{\partial \overline{U_i}\overline{V}}{\partial x_i}=D_m\frac{{\partial }^2\overline{V}}{\partial x_i{}^2}+\frac{p^{\prime }}{V_{mucus}}{I}_i-c\overline{V}$$

(5)

$$\frac{\partial T_i}{\partial t}=-\beta T_{i}V$$

(6)

$$\frac{\partial I_i}{\partial t}=\beta T_{i}V-\delta I_i$$

(7)

Here, U_i represents the convection velocity of mucus flow. For the estimation of the virus diffusivity D_m, the formulas selected here refer to the calculation method of drug spray diffusion in the nasal mucus layer during intranasal administration [35]:

$$\frac{D_m}{D_0}=e^{-\frac{\pi }{4}{\left(\frac{r_v+r_f}{r_g+r_f}\right)}^2}$$

(8)

$$D_0=\frac{k_{B}T}{6\pi {\mu }_0{r}_v}$$

(9)

where D₀ represents the virus diffusion coefficient in water, while r_v, r_f, and r_g express the effective radius of the virus, the mucin fiber radius, and the mucin network’s effective mesh fiber spacing, with values of 5 × 10⁻⁸ m, 3.5 × 10⁻⁹ m, and 5 × 10⁻⁸ m, respectively. To calculate D₀, a constant temperature T of 310 K in the airway was assumed, and the viscosity of water μ₀ was set as 0.00071 Pa·s. In addition, k_B is the Boltzmann’s constant at 1.380649 × 10⁻²³ J/K.

In this study, we calculated two cases: Case 1 with only virus diffusion (no mucus flow) and Case 2 with both convection and diffusion (with mucus flow).

2.3. Changes in the Infection Rate of the Target Cells β

Owing to the limited clinical data currently available, parameters fitted by different studies based on different patient data varied greatly. Here, the initial number of target cells, T(0), in the nasal cavity–nasopharynx could be determined by the methods and values used in previous calculations in Section 2.1. According to Equations (5)–(7), the infection rate of the target cells, β, had the main effect on both the peak viral load and the time to reach the peak viral load, which should be focused on. Among the references [17,20,37,38], β has a wide range, varying from 1.9 × 10⁻⁴ to 3.86 × 10⁻¹⁴ ((copies/mL) day⁻¹). Based on this range, this study selected a larger general range to include these values and divided them into seven cases to assess their effects on SARS-CoV-2 dynamics prediction. The detailed values are listed in

Table 2. The values of β varied among the cases.

No.	1	2	3	4	5	6	7
β ((copies/mL) day−1)	1 × 10−2	1 × 10−3	1 × 10−4	1 × 10−5	1 × 10−6	1 × 10−7	1 × 10−8

. The first two methods and conditions in Section 2.1 and Section 2.2 were used to calculate these cases simultaneously, and the influences of the mucus layer on the calculation results were compared.

3. Results

3.1. SARS-CoV-2 Dynamics Prediction of the Upper Respiratory Tract Based on the Target-Cell-Limited Model

According to the droplet deposition in Figure 3 and the initial conditions in Table 1 mentioned in Section 2.1, Figure 6 shows that the results obtained by using the target-cell-limited model to directly predict the viral load were similar to the changing trend of the clinical data of nine patients (a–i). The virus peaked at more than 10⁸ copies/mL, which was also consistent with the range of the peak value of the viral load (10⁶–10⁹ copies/mL) in their clinical data in the early stages of COVID-19 symptoms. [18] Furthermore, the viral load peaked around 5 d after infection and decreased below the detection limit at approximately 18 d after infection.

In Figure 7, the virus load in the vestibule remained at zero because there was almost no droplet left there. Meanwhile, the peak value of the viral load and the time to reach the peak at other sites were almost identical to the overall viral load, with only minor differences. Because only the initial value of the virus had a small difference, the effect on the change in the viral load was small.

3.2. SARS-CoV-2 Dynamics Prediction of the Nasal Cavity–Nasopharynx with the Mucus Layer

Among the nasal cavity–nasopharynx parts, four sites were selected for subsequent HCD analysis because of the abundance of droplets and the high initial viral concentration, as shown in Figure 8. These four hot spots were the right agger nasi, left agger nasi, left inferior nasal concha, and nasopharynx, respectively. Detailed information, including the center coordinates and initial viral load, is summarized in

Table 3. Details of the representative deposit sites.

Virus No.	0	1	2	3
Sites	Right Agger Nasi	Left Agger Nasi	Left Inferior Nasal Concha	Nasopharynx
Initial Virus Counting (copies)	2.71 × 10−4	1.17 × 10−4	2.21 × 10−5	2.95 × 10−5

.

For Case 1 with virus diffusion in the mucus environment, the change in the viral load and the number of target and infected cells within 50 d of infection (i.e., when the virus-laden droplets were deposited in the mucus layer) are shown in Figure 9 and Figure 10.

Because the actual virus detection limit was 100 copies/mL under existing technology [18], the viral load could begin to be detected 2 d after infection, with the first peak of viral load around 15 d after infection, and it remained for approximately 6 d. It then declined, and it did not fall below the virus limit until about 44 d, which is quite different from the data of clinical cases depicted in Section 3.1. Therefore, certain parameters in the HCD model should also be appropriately altered when considering virus diffusion in the mucus environment. In this study, the influence of the main parameter β on the prediction is considered in Section 3.3.

Meanwhile, the infected cells and the viral load maintained similar trends, while the majority of target cells were infected at a later stage. In addition, Figure 11 was illustrated for a more intuitive view of the SARS-CoV-2 dynamics with diffusion in the nasal cavity–nasopharynx model.

For Case 2 with virus diffusion and convection, Figure 12 is a vector diagram of the stabilized mucus velocity, indicating the mucus flow in the nasal cavity–nasopharynx part. The average flow rate of mucus in the whole model was approximately 3 mm/min, and the velocity around the central nasal passages was slower than that at other sites. However, this was not consistent with the summaries of Rygg et al. [27] and Shang et al. [28], who noted that the mucociliary clearance averaged 5–6 mm/min in a typical healthy adult. Moreover, in addition to the epithelium goblet cells, they found that part of the mucus produced in the maxillary sinus was located in the central nasal passages, which we had ignored. In addition, the goblet cells secreting mucus were unevenly distributed in the nasal cavity–nasopharynx part, especially in the vestibule, and there were few goblet cells. Thus, little mucus generally enters the vestibule and flows out of the nostril due to gravity, but its low flow rate (approximately 12 mm/h) could be neglected. In this study, we assumed that mucus was uniformly secreted by the whole nasal epithelium, and the flow rate of posterior nasal mucus measured in the experiment (8–10 mm/min) was used as the boundary condition for the CFPD simulation, which might be the reason for the differences.

As shown in Figure 13, the viral load dropped from the virus-laden droplet deposition was so rapid that it decreased to close to 0 within 2 h with the mucus flow because the mucus carrying the virus flowed much faster than the virus could spread. This result indicates that the mucus flow in the human respiratory tract has a great influence on the distribution of SARS-CoV-2.

As the convective velocity generated by mucus flow was so fast, the reasonable diffusion coefficient D_m calculated and the infection rate β selected were both small, leading to the virus flowing to the lower part of the upper respiratory tract with mucus before it had time to infect the target cells, as shown in Figure 14. Likewise, Figure 15 shows the entire process of SARS-CoV-2 dynamics with diffusion and convection for 30 d.

Only the diffusion and convection of SARS-CoV-2 in the mucus layer were considered here, so the choice of parameter values was crucial for this case.

3.3. Effects of the Infection Rate β on SARS-CoV-2 Dynamics Prediction

With the different values of β summarized in Table 2, seven cases were calculated and compared with the original case with β0 = 4.71 × 10⁻⁸ ((copies/mL) day⁻¹). For the direct calculation of the HCD model with Equations (1)–(3), excessive values of β1–β3 led to high viral load and errors in calculation, which were unreasonable and not shown in Figure 16. The results showed that β had a great influence on the time when the virus reached its peak value but not on the peak value and duration.

Nevertheless, it seemed that β had a more significant impact on the viral load prediction for the cases that combined HCD with CFD coupled to the mucus environment, as shown in Figure 17. For Case 1, (a) indicated that β not only affected the peak time and duration of the virus peak but also affected the peak value to a certain extent. The larger the β value, the faster the SARS-CoV-2 infection of the target cells, the faster the viral load reached the peak, the higher the peak, and the shorter the duration. When other conditions remain unchanged, among these different magnitudes of β values, β4 was more consistent with the previous clinical data in sputum. To make it more persuasive, the viral load at the nasopharyngeal hotspot was selected as the representative site for comparison with the pharyngeal swab data, as shown in (b), to verify the feasibility of β4 in this case.

Case 2 had a general consequence with Case 1, as shown in (c), whereas the conflict between the infection rate represented by β and the mucus flow rate determined the peak viral load that could be reached in this case, where the diffusion rate was too small to affect the results. In this calculation, the large value of β1 led to the instability of the calculation and an invalid result. The low rates of Β0, β6, and β7 resulted in a rapid outflow of the virus along with the mucus, contrary to clinical data. Therefore, these lines are neglected. The findings revealed that no infection in Case 2 caused the peak viral load to reach the clinical peak. They typically peaked in one day or approximately one day after infection, and then decayed. Interestingly, almost all of these prediction lines showed a slight two-peak character, the two peaks were more than 10 d apart, and the peak values were in the range of 10⁴–10⁷ copies/mL, which were similar to the predictions mentioned in Fatehi’s study [20]. However, they coupled the terms of immune system responses with the HCD model but without mucus flow. Considering peak times and bimodal phenomena alone, β5 was the optimal value for this simulation. Nevertheless, as seen from (d), the variation in viral load that it represents was quite different from the data at the pharyngeal swab sampling site. Hence, when convection was considered, other complex virus–epithelial interactions were not considered, which require further in-depth consideration.

4. Discussion

As mentioned above, some improvements were made on the results of previous studies, such as the adding convective–diffusion term to the target-cell-limited model, and combining CFPD model to analyze the SARS-CoV-2 dynamics in the mucus milieu based on the 3D-shell model. However, owing to ethical problems, the simulation results obtained in this study are difficult to verify experimentally and accurately. Therefore, as much clinical data as possible need to be collected for comparison and verification.

In this study, the distribution of the droplet’s deposition was simulated by the steady CFPD calculation, so there are many possibilities in which the transient CFPD calculation of mixed-laden particles of different sizes should be performed. As for the mucus properties, uneven secretion and distribution of mucus will also cause differences in the mucus flow field, which should also be optimized. In addition, the basic target-cell-limited model was used to calculate SARS-CoV-2 dynamics. Hence, the HCD model should also consider other parameters such as δ, the immune response, and the patient vaccination status, in addition to integration with CFPD analysis and convection–diffusion terms. The clinical data of the specific classification, such as mild, ordinary, and severe patients, should be sorted to select parameters in the HCD model.

In Case 1, considering only diffusion, is essentially similar to the virus culture experiment in a laboratory Petri dish. It is possible that the results can be verified by summarizing the laboratory data, but the initial conditions and parameters still need to be considered according to the differences between individuals. In Case 2 of convection and diffusion, the nonlinear coupling of each parameter needs to be considered in addition to the influence of β.

5. Conclusions

This study attempted to combine CFPD with HCD by coupling to the convection-diffusion term to predict the SARS-CoV-2 dynamics in the mucus milieu under different conditions. The simple simulation of droplet deposition at different parts of the upper respiratory tract, the change in the viral load, the number of cells, and the effects of various infection rates β were calculated and discussed in detail.

In the case of droplet deposition distribution in the upper respiratory tract calculated by CFPD steady-state before, for SARS-CoV-2, the trend in variation of the viral load at different sites was consistent by directly using the target-cell-limited model. Under the numerical conditions in this study, the viral load peaked around day 5, between 1 × 10⁻⁸ and 1 × 10⁻⁹ (copies/mL), which was reasonable compared to the clinical data.

When considering the SARS-CoV-2 dynamics in the mucus environment, it can be seen from the results of Case 1 and Case 2 that when the convection–diffusion term was added to the HCD model, the peak time, the peak value, and its duration were changed accordingly, which was quite different from the clinical data. This difference was particularly pronounced in Case 2 with convection caused by mucus flow. At this time, the corresponding parameters needed to be adjusted.

Then, the influence of infection rate β on the prediction of Case 1 and the original HCD model was analyzed. β had a significant impact on the three types of calculations. However, under the conditions of Case 1, various β values not only affected the peak time but also affected the peak value and duration. In this case, the viral load calculated by the magnitude represented by β4, which was 1 × 10⁻⁵ ((copies/mL) day⁻¹), was more consistent with the trend of clinical data. For Case 2, there was a slight bimodal phenomenon, which is worth investigating in conjunction with immune responses.

Taken together, although these simulation results have not been well verified experimentally, they will be an important application to help the understanding of infection dynamics of SARS-CoV-2 corresponding to more accurate initial conditions and physical parameters of SARS-CoV-2. Additionally, if an inhalable therapeutic drug for SARS-CoV-2 is developed, our application may greatly contribute to the development of a drug delivery system, which is a method for effectively transporting the drug to the infected site.