Does the exponential Wells–Riley model provide a good fit for human coronavirus and rhinovirus? A comparison of four dose–response models based on human challenge data

The risk assessments during the COVID‐19 pandemic were primarily based on dose–response models derived from the pooled datasets for infection of animals susceptible to SARS‐CoV. Despite similarities, differences in susceptibility between animals and humans exist for respiratory viruses. The two most commonly used dose–response models for calculating the infection risk of respiratory viruses are the exponential and the Stirling approximated β‐Poisson (BP) models. The modified version of the one‐parameter exponential model or the Wells–Riley model was almost solely used for infection risk assessments during the pandemic. Still, the two‐parameter (α and β) Stirling approximated BP model is often recommended compared to the exponential dose–response model due to its flexibility. However, the Stirling approximation restricts this model to the general rules of β ≫ 1 and α ≪ β, and these conditions are very often violated. To refrain from these requirements, we tested a novel BP model by using the Laplace approximation of the Kummer hypergeometric function instead of the conservative Stirling approximation. The datasets of human respiratory airborne viruses available in the literature for human coronavirus (HCoV‐229E) and human rhinovirus (HRV‐16 and HRV‐39) are used to compare the four dose–response models. Based on goodness‐of‐fit criteria, the exponential model was the best fitting model for the HCoV‐229E (k = 0.054) and for HRV‐39 datasets (k = 1.0), whereas the Laplace approximated BP model followed by the exact and Stirling approximated BP models are preferred for both the HRV‐16 (α = 0.152 and β = 0.021 for Laplace BP) and the HRV‐16 and HRV‐39 pooled datasets (α = 0.2247 and β = 0.0215 for Laplace BP).

respiratory virus.These models are based on two principles: the estimation of the intake dose of the infectious agent and the estimation of the probability of infection under a given intake dose.The two most commonly used dose-response models for calculating the infection risk of respiratory viruses are the exponential and β-Poisson (BP) models.Both models assume a random distribution of the number of copies in the exposed medium described by the Poisson probability distribution.The main difference between the exponential and BP models is in the probability defining the host suscepti-bility p (%).The exponential model assumes that each host has an equal probability of getting infected with a single dose defined as where p (-) is defined as the number of ingested TCID 50 doses that will cause an infection, that is the infectious dose (ID).If the exposed dose contains some number of d (TCID 50 ) doses, the total infection probability defined by the exponential model is calculated according to the following exponential equation: The modified exponential version known as Wells-Riley is the classical and widely used model to quantitatively assess airborne infection risk (Riley et al., 1978;Wells, 1955).The Wells-Riley model implicitly calculates the airborne infection risk using the concept of quantum-one quantum is defined as the number of inhaled IDs required to infect 63.2% of susceptible persons in an enclosed space (d = p).The exponential Wells-Riley model has extensively been used to evaluate the airborne infection risk of respiratory diseases (Azimi et al., 2020;Gao et al., 2016;Patterson et al., 2017;Wagner et al., 2009).Although the exponential model assumes that all people are equally susceptible to infection (i.e., p = constant), the BP function assumes that the host susceptibility p holds a value between 0 and 1 and follows the β distribution.Using the Kummer confluent hypergeometric function 1F1(,  + , − d), the final expression of the exact BP model is Due to the numerical complexity in its model specification and the difficulty in parameter estimation in the exact BP model, an approximate version of the BP model was simplified by Furmoto and Mickey (1967) using the following Stirling approximation: Although the approximate BP model provides more flexibility compared to the exponential model due to the increased number of parameters, the Stirling approximation restricts this model to the general rules of  ≫ 1 and  ≪ .However, this approximation accuracy remained largely ignored in practice as the rules were considered too general to follow.Not until exactly 50 years after the approximate model was published by Xie et al. (2017) provide specific rules for which the approximate BP model worked with high accuracy.Namely, when validated against 85 datasets, the study by Xie et al. (2017) showed that an accurate approximation was achieved under the following constraints:  > (22 ⋅ ) 0.5 for 0.02 <  < 2. While revisiting studies applying and recommending the Stirling approximated BP model for some human respiratory viruses (Jones & Su, 2015), it was observed that the majority of the approximate BP did not comply with the requirements defined by Xie et al. (2017).To refrain from these requirements while also preserving the simplicity of an approximate model without sacrificing the quality of adequacy, this study attends to exploit the Laplace approximation of Kummer hypergeometric function in the exact BP model rather than adhering to the conservative Stirling approximation.The explicit and relatively simple analytical solution to the Laplace approximation of the Kummer hypergeometric function was introduced by Buttler and Wood (2002).This study also demonstrated that the Laplace approximation generates relative errors of no more than 0.1% compared to the exact Kummer hypergeometric function, even for a wide range of conditions, including quite small values of  and  near 1 (Butler & Wood, 2005).Therefore, one of the objectives of this study is to utilize the Laplace approximation of the Kummer hypergeometric function in the exact BP model and compare this novel approach to the Stirling approximated BP model and exponential model for dose-response datasets of human respiratory viruses when the requirements set by Xie et al. (2017) are not met.The datasets of human respiratory airborne viruses available in the literature for human coronavirus (HCoV-229E) and human rhinovirus (HRV-16 and HRV-39) are used to compare the four dose-response models.For practical purposes, the outputs of the four dose-response models are also compared when coupled with a mass-balance model to simulate the long-distance airborne infection risk inside a classroom.

Data selection
The outcome of interest was human challenge studies reporting dose-response data for infectious respiratory viruses.Publications with datasets that facilitated further analysis included the exposure doses, the total inoculated humans per dose, and the number of infected humans per dose, or the publications had sufficient data to determine these parameters.Data met the following conditions for dose-response analysis: (i) At least three dose levels were used and (ii) the number of doses with a response rate other than 0% or 100% was equal to or greater than the number of dose-response model parameters.These conditions are consistent with recommendations by Haas et al. (1999).Studies were identified through searches of the PubMed and Web of Science databases and by cross-referencing cited references.Amar Aganovic reviewed all studies to assess whether studies met the data conditions for dose-response analysis.The search procedure resulted in finding two studies reporting dose-response models derived

Laplace approximation of the Kummer hypergeometric function
The Laplace approximation of the Kummer hypergeometric function when applied to the exact BP will have the following form (Butler & Wood, 2002): where the Laplace approximation of the Kummer hypergeometric function [6]: and

Dose-response model fit and evaluation
The model parameters k, α, and β were determined by fitting compiled experimental datasets using the following maximum likelihood estimation (Holcomb et al., 1999;McCullagh & Nelder, 1989): where for each trial i, we have a dose   and a group of   volunteers of which   are infected.P i is the value obtained from the exponential and the BP models.The goodness-of-fit estimates of the model parameters are those values that minimize Y.The performance of the neural network is then compared to the better fit model by the following Akaike information criterion (AIC) (Akaike, 1974): where k is the number of parameters in the model, and l is the log-likelihood: To determine whether one model m i is statistically better than other models, we will use the relative likelihood by exponentiation of the ΔAIC = AIC i − AIC min : This allows us to calculate the relative probability of each model by normalizing each likelihood value: where M is the total number of candidate models (Burnham & Anderson, 2004).If p(m i ) > 0.95, then candidate model m i is considered significantly better than the other.The ability to pool data for the exponential and BP model is assessed via the difference between the deviance of the pooled dataset  and the sum of each individual optimized deviance (1, 2, …) and Δ is compared to a Χ 2 distribution with df = (sum of the number of parameters used in fitting individual datasets) -(the number of parameters used in fitting the pooled dataset) (Watanabe et al., 2010): Confidence intervals for parameters for each model were determined with bootstrapping methods.

Airborne infection risk
Typically indoor airborne transmission risk models are either created using a simplified approach by analytically solving the conservation of viral concentration equations for the contaminants under quasi-ideal and quasi-uniform assumptions while not considering the airflow dynamics inside the box (Aganovic et al., 2021;Aganovic, Bi, et al., 2022;Aganovic, Cao et al., 2022;Aganovic et al., 2023).The main assumption of this mode is that the virus-carrying aerosols uniformly and instantaneously dispersed across space.If n denotes the number of RNA viral copies in the ambient air, the unique solution of viral concentration in an indoor environment with complete mixing ventilation at time t, n(t) is calculated as where n 0 is the initial viral concentration ( RNA m 3 ) at time t = 0.The removal terms due to ventilation  vent (h −1 ), the aerosol deposition rate  dep (h −1 ), the viral inactivation rate  inact (h −1 ), and respiratory absorption  abs (h −1 ) are summed into one term called the infectious virus removal rate for complete mixing ∑ (h −1 ): ) is the source emission rate of aerosolized RNA copies.For the sake of brevity, the equations for calculating the removal terms and the source emission rate are not repeated here.All the equations can be found in recent publications (Aganovic et al., 2023).In the absence of data on the relationship between TCID 50 and the number of RNA copies for HCoV-229E and HRV 39 and 16, we used data for the SARS-CoV-2 wild strain for which 1 TCID 50 ≈ 10 4 RNA copies (Sender et al., 2021) for human coronavirus and rhinovirus FEB 3 for human rhinovirus for which 1 TCID 50 = 53797 RNA copies (Parker et al., 2015).Thus, the airborne dose (TCID 50 ) after a time interval t(s) can be calculated as d(t) = 10 −4 ⋅ n(t) for human coronavirus and d(t) = IR is the inhalation rate of the exposed subject, which was assumed to be the inhalation rate for resting and standing averaged at 0.52 m 3 /h (Adams, 1993).
To assess the relative impact of different models on airborne transmission risk, a simple case study was investigated with the same dimension characteristics and the number of persons present for each case considered.The simple scenario consisted of a classroom with an area of 64 m 2 and 3 m height with 1 infected and 20 susceptible persons present.As only long-distance airborne transmission risk is considered, all the persons were distanced 1.5 m.The time exposure considered was 180 min.The indoor temperature was in the range of 20-25 • C (Aganovic et al., 2023).

Dose-response model for human coronavirus HCoV-229E
The exponential model was the best fitting model for the HCoV-229E datasets from the study by Bradburne et al. (1967) and significantly better (p = 0.99) fit than the second-

Dose-response model for the pooled datasets of human rhinovirus HRV-39 and HRV-16
The exponential model was the best fitting model for the variant HRV-39 (Hendley et al., 1972), but it was not significantly a better fit than the rest of the BP models (p < 0.95), as shown in Table 3.The Stirling approximated BP model did not comply with the specific rules provided by Xie et al. (2017) for accurate approximation.On the contrary, for both the antibody and low antibody HRV-16 datasets, the goodness-of-fit criteria did not approve the exponential model (Y = 41.27> χ 2 = 11.07),whereas the Laplace approximated BP model was the best fitting model.For both the antibody-free and low body HRV-16 datasets, there was no significant difference between the Laplace approximated BP model and Stirling approximated BP model; however, in the Stirling approximated BP  was larger than .As the difference between the pooled model's deviance and the sum of the individual models' deviances was larger than the critical χ 2 for the exponential model, both the HRV-39 and HRV-16 datasets could not be pooled for the exponential model.For the rest of the three BP models, this was not the case as there was no statistically significant evidence to reject the null hypothesis.
The dose-response models for the pooled datasets of human rhinovirus are shown in Figure 2. The Laplace BP model was the best fitting model for the pooled datasets with a deviance of Y = 2.0.There was no significant difference between the Laplace and Stirling approximated BP models (p < 0.95).The Stirling approximated model did not again meet the criteria for an accurate approximation ( > ).The fit for the exact BP model was rejected for the pooled dataset.Therefore, the Laplace approximate is the final recommended model for the dose-response assessment of human rhinovirus.

Airborne infection risk
We compared the different dose-response models for evaluating airborne infection risk for the case scenario explained in Section 2.3.Figure 3 shows the four different dose-response models for two different viral loads of human coronavirus: (i) c v = 10 9 RNA∕mL, corresponding to doses in the range d = 0.02 − 0.17 TCID 50 and (ii) c v = 10 11 RNA∕mL corresponding to doses in the range d = 2 − 17 TCID 50 , respectively.As may be seen from Figure 2, for the lower viral load c v = 10 9 RNA∕mL, the differences between the three BP models are relatively large, whereas for higher viral loads c v = 10 9 RNA∕mL, these differences almost become negligible.On the other hand, regardless of the dose range, the exponential model that did not comply with the goodnessof-fit criteria for the pooled data underestimated the airborne cross-infection risk substantially compared to the BP models.
The reason for this may be explained by the differences in the steepness of the gradients of the dose-response curves for human coronavirus at lower doses d = 0.02 − 0.17 at c v = 10 9 RNA∕mL versus the gradient steepness at higher doses d = 2 − 17 TCID 50 at c v = 10 11 RNA∕mL.The differences in the gradient steepness of the four models at low doses can be observed in Figure 1.
Figure 4 shows the four different models for human rhinovirus for two different viral loads: (i) c v = 10 9 RNA∕mL, corresponding to doses in the range d = 0.0038 − 0.0303 TCID 50 and (ii) c v = 10 11 RNA∕mL corresponding to doses in the range d = 0.38 − 3.03 TCID 50 , respectively.Similarly, as for the human coronavirus, the lower viral load c v = 10 9 RNA∕mL, the differences between the three BP models are relatively large, whereas for higher viral loads c v = 10 9 RNA∕mL, these differences almost become negligible (Figure 4).Although rejected for the pooling criteria, we compared the exponential model to the BP models, and as in the case of the human coronavirus, it underestimated the airborne cross-infection risk substantially when compared to the BP models.
Again, the reason for this may be explained by the differences in the steepness of the gradients of the dose-response curves for human rhinovirus at lower doses versus the gradient steepness at higher doses, as shown in Figure 2.

DISCUSSION
To control a pandemic outbreak, it is necessary to understand how the virus transmits to a susceptible individual and eventually infects the community.To accurately quantify transmission risk, a dose-response function should be based on experimental data that represent this process.In the case of coronaviruses and rhinoviruses, transmission occurs through exposure to respiratory fluids carrying the infectious virus (Leung, 2021).The virus-carrying respiratory droplets and aerosols can be produced through all expiratory activities, including breathing, talking, coughing, and sneezing, from both symptomatic and asymptomatic individuals.These infectious aerosols and droplets may come into direct contact with susceptible individuals by being inhaled from the surrounding air or by indirect contact when the susceptible individual touches a surface contaminated by infectious respiratory fluid (Wang et al., 2021).
Once recognized as the main route of COVID-19 spread, identifying the relative importance of different engineering controls targeting the spread of COVID-19 in indoor environments requires accurate prediction of the transmission risk.In this context, there is a need for predictive risk assessment tools for better understanding when planning effective strategies to minimize risks associated with airborne transmission.The concept behind the mathematical tools used so far for modeling airborne transmission risk is based on coupling dose-response models with a box model containing a source and sink of contaminants (Aganovic et al., 2023).The infection risk inside the box can be modeled with a simplified approach by analytically solving the conservation of mass equations for the contaminants under quasi-ideal and quasiuniform assumptions/conditions, as shown in Section 2.3.On the other hand, the two most extensively used dose-response models for calculating the infection risk of respiratory viruses are exponential and the Stirling approximated BP models (Watanabe et al., 2010;To & Chao, 2010).
It is important to emphasize the difference between the exponential and BP models at higher doses as the exponential-based Wells-Riley model was almost exclusively used during the pandemic for both prospective and retrospective risk assessments.The exponential Wells-Riley model does not take into account the variation in host susceptibility; rather it assumes that all individuals exposed to a respiratory virus respond equally to the same amount of the intake dose, which is not the case in reality (van Slujis et al., 2017).Unlike the exponential model, the BP model accounts for interindividual variability of the immune status and the host's sensitivity to the pathogen.
The only existing dose-response model for SARS-CoV far originated from datasets for mice and not for humans (Watanabe et al., 2010).In the absence of human challenge data for SARS-CoV-2 during the pandemic, the exponential dose-response model derived by Watanabe et al. (2010) served as a base for the many risk assessments based on the exponential Wells-Riley models (Buonanno, Morawska et al., 2020;Bazant & Bush, 2021;Buonanno, Stabile et al., 2020;Cortellessa et al., 2021;Schijven et al., 2021) However, In studies with animal models, there are additional host-specific concerns.Despite similarities differences in susceptibility between mice and humans exist for respiratory viruses (Heykers et al., 2019;Zhang et al., 2021), including SARS-CoV-2.Therefore, a fundamental doseresponse relation is missing for a more realistic evaluation of the human-to-human transmission risk of coronavirus.On the other hand, rhinovirus dose-response models based on datasets from human challenge studies have been developed (Patterson et al., 2017).However, the datasets were tested for an approximate BP dose-response model that did not comply with the Stirling approximation that restricts this model to the general rules of  ≫ 1 and  ≪ .Furthermore, the existing human rhinovirus dose-response models (Jones & Su, 2015) did not even follow the requirements for an accurate approximation under the following requirements:  > (22 ⋅ ) 0.5 for 0.02 <  < 2 as defined by Xie et al. (2017).To refrain from these requirements, we developed and tested a novel BP model by using the Laplace approximation of the Kummer hypergeometric function instead of the conservative Stirling approximation.The models were additionally compared to the BP model based on an exact calculation of the hypergeometric function.Based on currently available data, and based on AIC, the BP exponential model was the best fitting model for the HCoV-229E and for HRV-39 datasets, whereas the Laplace approximated BP model followed by the exact and Stirling approximated BP models are preferred for both the HRV-16 and the HRV-16 and HRV-39 pooled datasets.Still, as the Stirling approximated dose-response model did not even follow the requirements for an accurate approximation, the novelty of applying the relatively simple Laplace approximation has the added advantage of not being subject to such restrictions.If it is supported or favored based on AIC or some other objective criterion, the Laplace approximated BP should be utilized to assess infection risk.However, all three BP models varied substantially in between for very low doses.Although this was not the case at higher doses, the relatively low doses based due to assumed lower viral loads present a more realistic scenario compared to the high dose super-spreading events.On the other hand, the conventionally used exponential dose-response model generated lower infection risk compared to the BP models by a substantial amount for all datasets.
A methodological issue with the study is the low number of human datasets for corona-and rhinovirus available in the literature for modeling.However, human challenge studies are rarely very large, and despite the small group numbers, the dose-response relations are well defined.This study does by no means offer a conclusive dose-response model for either of the two viruses considered but rather amplifies the importance of the differences between dose-response models, especially for low doses.Another limitation of this study is that the human volunteers from the datasets used in this study were administered through intranasal instillation.Still, given sufficient fluid volume, intranasally instilled materials have been shown to move to the lungs (Southam et al., 2002;van Slujis et al., 2017).In studies with human volunteers, this means the exposure conditions and doses administered should reflect natural transmission conditions.
The case study presented in this analysis also highlights some of the knowledge gaps that still need to be addressed in terms of human corona-and rhinovirus exposure risks.Future human challenge data studies should focus more on low-dose data due to the following two reasons: (i) It is in low-dose ranges that the dose-response models differ; (ii) low doses present a more realistic scenario compared to high doses for risk assessment of airborne cross-infection risks.

CONCLUSION
The objective of this study was to fit dose-response models for selected viruses capable of infecting the respiratory tract to facilitate infection risk estimates.The datasets of human respiratory airborne viruses available from the literature for human coronavirus (HCoV-229E) and human rhinovirus (HRV-16, and HRV-39) are used to compare the four dose-response models.BP.Based on goodnessof-fit criteria, the exponential model was the best fitting model for the HCoV-229E and for HRV-39 datasets, whereas the Laplace approximated BP model followed by the exact and Stirling approximated BP models are preferred for both the HRV-16 and the HRV-16 and HRV-39 pooled datasets.Despite negligible differences between the BP models at higher doses, there were considerable differences at low doses (n < 5.0 TCID 50 for human coronavirus and n < 1.0 TCID 50 for human rhinovirus).As low doses present a more objective scenario compared to the high dose super-spreading events, future human challenge data studies should focus more on low inoculation doses.

C O N F L I C T O F I N T E R E S T S TAT E M E N T
We have no conflict of interest to disclose.

F
The dose-response model fitted to the dataset for human coronavirus.BP, β-Poisson.bestranked model-the Stirling approximated BP model (Table2).The Stirling approximated BP model complied with the specific rules provided byXie et al. (2017) for accurate approximation.There was no significant difference among the exact, Stirling, and Laplace approximated BP models for this dataset (p < 0.95), as shown by the dose-response lines overlapped in Figure1.The dose-response models for the HCoV-229E datasets from the study byBradburne et al. (1967) are shown in Figure1.

F
The dose-response model fitted to the pooled datasets for human rhinovirus.BP, β-Poisson.F I G U R E 3 Long-distance airborne cross-infection risk for human coronavirus calculated based on four different dose-response models at two different viral loads.BP, β-Poisson.

F
I G U R E 4 Long-distance airborne cross-infection risk for human rhinovirus calculated based on four different dose-response models at two different viral loads.BP, β-Poisson.

TA B L E 1
Dose-response datasets selected from the literature.

Type (Reference) Age Dose (TCID 50 ) Tested Positive Negative
Comparison of estimated dose-response models based on goodness-of-fit to the observed human coronavirus.
Comparison of estimated dose-response models based on goodness-of-fit to the observed human rhinovirus.
Xie et al. (2017)with rules provided byXie et al. (2017)for accurate approximation of the Stirling function.b Antibody-free volunteers.c Low antibody volunteers.Abbreviations: AIC Akaike information criterion; BP, β-Poisson.