Parameter Optimization of a Viral Dynamics Model in the Mucus Layer of the Human Nasal Cavity-Nasopharynx Based on Computational Fluid-Particle and Host-Cell Dynamics

. Respiratory diseases, such as COVID-19 (coronavirus disease 2019) caused by SARS-CoV-2 (severe acute respiratory syndrome coronavirus 2) , have posed a threat to human health. For infection control and a better understanding of the pathogenesis, this study mainly focused on elucidating the virus dynamics in the mucus layer of the human nasal cavity-nasopharynx, using coupled computational fluid-particle dynamics (CFPD) and host-cell dynamics (HCD) analyses. To reproduce virus transportation in the mucus layer by mucociliary motion, a three-dimensional-shell model was created using the data obtained from computed tomography (CT) of the human upper airway. By considering the mucus milieu, the target-cell-limited model was coupled with the convection-diffusion term to develop the HCD model. Parameter optimization has been shown to have a great impact on the accuracy of model prediction; therefore, this study proposes a method that divides the geometric model into multiple regions and uses Monolix for non-linear mixed effects modeling for pharmacometrics. The results showed that data from human inoculation challenge trials could be used to estimate the corresponding parameters. The models developed and used with optimized parameters can provide relatively accurate predictions of virus dynamics, which could contribute to the prevention and treatment of respiratory diseases.


Introduction
Respiratory diseases pose a threat to human health, more so during this age of coronavirus disease 2019 . Understanding the pathogenesis of respiratory viruses, accurately predict and grasp virus dynamics, and identify hotspots of infection for disease prevention and treatment are all crucial. Clinical and experimental data are essential for predicting and simulating virus dynamics. At the early stage of COVID-19, nearly no one understood the nosogenesis of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), and inoculation experiments could not be performed quickly enough; therefore, only limited clinical data could be obtained to ensure safety. Wölfel et al. [1] reported a detailed virological analysis of nine COVID-19 cases, providing evidence of active viral replication in the upper respiratory tissues and finding a peak viral load measured in throat swabs on the 4 th day since the onset of symptoms. Previous studies on virus dynamic simulations, including the study by Hernandez-Vargas et al. [2], were mostly implemented based on these data. A mathematical model of SARS-CoV-2 dynamics in infected patients at different onset times was proposed in their study based on a target cell-limited model. The optimal model to fit the data included immune responses, which indicated a slow immune response that peaked 5-* Corresponding author: li.hanyu@kyudai.jp 10 d after the onset of symptoms. Perelson et al. [3] highlighted the role of a simple target cell-limited model and its evolution into more complex models to capture more knowledge of virology and immunology. Nevertheless, most current studies have ignored viral dynamics in the mucus layer of the respiratory tract. The main function of the mucus layer is to remove unwanted particles, such as inhaled bacteria and viruses, through the rhythmic pulsation of the cilia on the epithelial tissues, which is called mucociliary clearance [4]. Therefore, Li et al. [5] referred to a model establishment of mucociliary clearance in the nasal cavity and related studies of intranasal drug administration by Shang et al. [6] [7] combined with a target cell-limited model and static particulate exposure in the nasal cavity, to develop a host-cell dynamics (HCD) model of virus in the mucus layer. However, the results showed that the influence of parameters should not be underestimated and requires further discussion. Later, as experts learnt more about SARS-CoV-2, Killingley et al. [8] developed a novel SARS-CoV-2 human challenge model. While ensuring its ethics and safety, an intranasal inoculation of a D614G-containing pre-alpha wild-type virus in 36 young volunteers was performed to establish the virus dynamics in the main infection process of SARS-CoV-2, which is of significance for public health recommendations and strategies affecting the transmission of SARS-CoV-2. Therefore, this study used experimental data, as suggested by previous studies, to perform the fitting optimization of the HCD model parameters and summarize the systematic method.

Methods
Based on the previous studies, a scenario was set up as briefly described: two CSPs (computationally simulated persons) were standing face-to-face in a 3×3×3-m 3 displacement ventilation room 1-m apart. To simulate the worst-case scenario of cough droplet transmission, the person next to the inlet was the infected person coughing, and the person on the other end was a healthy person with a nasal breathing cycle of 4 s. To assess the distribution of infectious droplet deposition in the respiratory tract of susceptible individuals, data obtained from computed tomography was used to create a three-dimensional (3D) upper airway model coupled with a healthy CSP. A total of 45,000 particles were generated, and their particle sizes ranged from 1 to 10 μm (approximately 70% of the droplets produced in a cough). Lagrangian tracking analysis of the particles was performed using computational fluid-particle dynamics (CFPD) unsteady simulation. Evaporation or condensation was not considered in this study. The results showed that more than 90% of droplets were deposited in the nasal cavity. As shown in Fig. 1, mucus in the respiratory tract is secreted by sinuses and goblet cells in the epithelium, and roughly contains two layers: the upper layer is a high-viscosity gel layer that is composed of 97% water and 3% other substances (mucin, salts, lipids, immunoglobulin, etc.), with a thickness of 0.5-5 μm, while the lower layer is a ciliary-liquid layer with a thickness of 7-10 μm, with a viscosity similar to that of water [9] [10]. Given the concentration of virus-laden droplets in the nasal cavity and the direction of mucus flow, the nasal cavity-nasopharynx part was selected for the study of viral dynamics to better compare with nasopharyngeal swab data from experimental or clinical cases. A 3D shell model with a mucus layer with a thickness of 15 μm was constructed for the HCD analysis.

Basic HCD model
When virus-laden droplets are deposited on the mucus layer of the respiratory tract along with inhaled air, a series of dynamic viral activities occur. For instance, the virus transports and diffuses with the mucus flow and combines with the receptors of target cells to infect them. The proliferation, release, and apoptosis of viruses occur in the infected cells. In the present study, the target-cell limited model was considered as the basic formula, and the HCD model required for this study was updated by combining the convection-diffusion term with it, as follows: where T, I, and V represent the number of susceptible target cells, number of infected cells, and viral load, respectively. The initial target cell number, T(0), was estimated based on the surface area of the airway and surface area of each epithelial cell. The initial viral load, V(0), was calculated by investigating the droplets of different sizes deposited, and the amount of virus carried by each droplet. In terms of the parameters, β represents the rate at which the virus infects target cells in copies/mL/day. Infected cells were cleared at a δ (/day) rate owing to cytopathic viral effects and immune responses. The virus was cleared at a rate of c (per day). The infected cells in various parts of the respiratory tract produced and released viruses at the rate of p' (copies/day/cell). The volume of mucus, Vmucus, was the product of the surface area and thickness of the mucus. The parameter values have been shown to have a great impact on the model, and to better understand the virus dynamics in the mucus layer and visualize it, we conducted experiments with regards to parameter optimization.

Parameter optimization
Here, we proposed a method for parameter fitting optimization. First, the geometric model was divided into k compartments with equal surface areas. Three partitioned modes (k = 1, 10, and 20) were preliminarily considered, as shown in Fig. 2. According to the simulated droplet deposition data and proportion of viruses in each droplet, the total initial virus titer was calculated as V(0) = 0.43 copies/mL, and the initial virus concentration Vk(0) of each compartment is set based on the particle deposition location. Under the initial conditions, the initial number of target cells, T(0), was estimated according to the surface area of the airway and the surface area of each epithelial cell, and 4 × 10 -11 m 2 /cell was selected [11]. According to a rough estimate [12], approximately 20% of epithelial cells in the respiratory tract express receptors and coreceptors, such as ACE2 and TMPRSS2 [12]. This estimate was used to set the initial target cell number Tk(0) in each compartment. The initial number of infected cells, Ik(0), was assumed to be 0. Subsequently, combined with the virus human challenge vaccination data, the value was substituted into Monolix for parameter fitting.
In addition, the influence of the number of layers according to the properties of the mucus layer on the fitting results was considered. The density ρ and viscosity μ of mucus were set to 1000 kg/m 3 and 12 Pa·s, respectively [6]. The nasopharynx is located at the back of the nose and is thought to be the outlet for mucus produced in the nasal cavity.

One-layer multi-compartment model
According to the existing assumptions [6], the mucus layer can be simplified into one layer, as shown in Fig.  3.
where Tk, Ik, VLk, Vmucus, and k represent the number of susceptible target cells, number of infected cells, viral load, and mucus volume in each compartment, respectively. Q represents the flux of mucus flow between adjacent compartments, which was estimated as the product of the outlet area and velocity. As mentioned by Shang et al. [6], the outlet velocity, which is the posterior nasal mucus flow rate of 10 mm/min, was selected here.

Two-layer multi-compartment model
In addition, the manner in which mucus retains two layers was taken into account, as shown in Fig. 4.

Fig. 4. Viral dynamics of the two-layer multi-compartment model
In the gel layer, convection is caused by mucus flow, whereas the liquid in the sol layer is relatively stationary, owing to the existence and swinging of cilia. Therefore, the infected cells are assumed to release the virus into the sol layer only by diffusion. After diffusion to the gel layer, convection and diffusion occur along with mucus flow; thus, the developed model can be expressed as follows: where VLsol,k, VLgel,k, Vsol,k, and Vgel,k are the viral load and mucus volume of the sol and gel layers in each compartment, respectively. D represents the virus diffusion coefficient in the mucus (m 2 /s), while Δx is the distance between the sol and gel layer centers and Ak is the area of the interface between the sol and gel layers in each compartment. The other parameters are the same as those defined in section 2.2.1.

Results and discussions
We performed a parametric fitting of viral load data from 12 untreated participants of previous studies [8] by Monolix and then 12 parameter combinations were substituted into the HCD model, and CFD was used for simulation. The best parameter combinations fitting the experimental data for diffusion or convection-diffusion cases were obtained, which were shown in Fig. 5.  For all cases, it could be known that β was concentrated on the order of magnitude of ~×10 -8 , whereas this value was slightly lower for the case of one compartment. The convection-diffusion case yielded slightly higher values than the diffusion case. δ was almost distributed in the range of 1-2 /day. The distributions of c and p' were similar: for the diffusion cases, the values of c were concentrated within 5/day, while the values of p' were distributed below 15 copies/day/cell; for the convectiondiffusion cases, the values of c spanned four orders of magnitude, and the value of p' was relatively large, wherein the two-layer results are relatively stable. This finding could be related to the convection created by mucus flow. According to the results in Fig. 5, the best fitting case was the convection-diffusion case with two layers-20 compartments, which was visualized in Fig.6.

Conclusions
This study provides a practical method for parameter optimization in a mathematical model based on a 3Dshell model that coupled CFPD and HCD models to predict virus dynamics in the mucus layer of the human airway. The number of mucus layers, partitioning method of the geometric model, and addition of convection-diffusion have a certain degree of influence on the parameter fitting. The parameter combination obtained in this study had a good effect on the specific optimization of virus dynamic prediction. Based on these results, it is possible that virus remains in the front of the nasal cavity even if the viral load in the nasopharynx is reduced to a lower level after a long time. This information will help prevent and treat respiratory diseases.