Dose–Response Modeling: Extrapolating From Experimental Data to Real‐World Populations

Dose–response modeling of biological agents has traditionally focused on describing laboratory‐derived experimental data. Limited consideration has been given to understanding those factors that are controlled in a laboratory, but are likely to occur in real‐world scenarios. In this study, a probabilistic framework is developed that extends Brookmeyer's competing‐risks dose–response model to allow for variation in factors such as dose‐dispersion, dose‐deposition, and other within‐host parameters. With data sets drawn from dose–response experiments of inhalational anthrax, plague, and tularemia, we illustrate how for certain cases, there is the potential for overestimation of infection numbers arising from models that consider only the experimental data in isolation.


INTRODUCTION
Mathematical models have long been used to help interpret experimental data arising from challenges of infectious disease agents to animals. Any model is an abstraction of reality and necessarily simplifies the underlying processes, often to arrive at a solution that is mathematically tractable. However, models of any particular phenomenon should be carefully constructed and the associated assumptions, where possible, fully evidence-based and clearly stated, especially if chosen for convenience in absence of such an evidence base.
One concept often used in challenge studies is that members of a population have a variety of in-herent tolerance levels that, if exceeded by exposure level, result in response. Irwin was the first to coin the expression "individual effective dose" that has been widely used since (Irwin, 1937). In contrast to this concept which assumes "maximum synergism" between bacteria (Meynell & Stocker, 1957), models may alternatively assume the independent action of organisms. In such models, the variability with respect to host is replaced with stochasticity with respect to bacteria, whereby each individual organism has a probability of alone invoking response (through survival and sufficient proliferation). This is often termed a "hit" and one or more of these could be required to initiate response (Teunis & Havelaar, 2000). A third concept that may be considered is that of a bacterial threshold, which when reached results in host response. Here, the host starts with an initial load and through bacterial reproduction and death, the number of organisms eventually reaches zero (resolution) or a threshold (response).
Recent work on Bacillus anthracis has attempted to review dose-response models (both individual effective dose and single-hit) against a set of criteria (Toth et al., 2013). These criteria are that the model 67 0272-4332/21/0100-0067$22.00/1 © 2020 The Authors. Risk Analysis published by Wiley Periodicals LLC on behalf of Society for Risk Analysis This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. is fitted to data (rather than expert elicitation), the resulting fit is consistent with the pattern of infection arising from the Sverdlovsk incident (Meselson et al., 1994;Wilkening, 2006), the model is biologically motivated, and finally that the incubation period is inferred. Models more complicated than the exponential model (one form of single-hit type) are eventually rejected in Toth et al. (2013) because the extra complexity (as given by additional parameters) is not supported by the model fits (i.e., the fit is not statistically better with additional complexity and biological rationale in the model). These are sound conditions when interpreting a data set in isolation, however, no consideration was given to variation in say, individual susceptibility or clustered dosages when extrapolating away from experimental data to a human population. The purpose of this study is then to illustrate, with plausibly motivated but parsimonious models, the potential impact of extrapolating inferences based on experimental data to wider populations. In deriving models, we use language relevant to anthrax infection but will also discuss model fits to other biological organisms.

METHODS
To model the dose-response for an individual in a population to a challenge of B. anthracis, we assume that there are two key events for each bacterium: germination of the spore, after which point infection is inevitable in the host, or the removal of the spore from the system after which infection (arising from that spore) is impossible; this is the approach taken by Brookmeyer, Johnson, and Barry (2005). There is effectively a competition between the two processes and so if the removal of the spore were to occur at time σ and germination at time τ then if σ ≤ τ infection is stopped but if σ > τ infection occurs. This is conceptually often called a competing-risks framework: occurrence of one event precludes occurrence of the other.
There are a number of natural questions arising for any probabilistic interpretation of this framework: (1) What distributions/processes govern σ and τ within an individual? (2) What impact does inhalation of multiple spores have on the response? (3) What differences are there in σ and τ (or related values) between individuals?
Here we attempt to provide a framework to address these questions for infection with B. anthracis. We will develop a general model below, then make use of three experimental data sets to illustrate the potential impact of extrapolating away from carefully calibrated experimental settings to wider scenarios. In the absence of data for these wider scenarios, we illustrate on an order of magnitude basis how changes in the additional parameters can influence inferred infection rates. Limitations in deriving and parameterizing the models are discussed.

The Competing-Risks Framework
One may derive the single-hit models from the competing-risks framework, and this derivation is repeated here for transparency and to allow clarity in assumptions made and extensions possible. Let us assume that the germination time for each spore in each individual is modeled by some probability density function f G (t, G ) with associated parameter set G and mean μ G and similarly the removal time is given by f R (t, R ) and mean μ R , with corresponding cumulative distribution functions, F G and F R , respectively. We assume that germination and removal are independent processes.
The probability β that σ ≥ τ for any given spore in each exposed individual (and thus infection subsequently occurring) is then (1) where h G and h R represent the hazard functions of germination and removal, respectively. If the hazard function over time for germination is proportional to the hazard of removal then h R = γ h G and and so a single parameter describes the underlying process, the ratio of the removal time to germination time. In the situation where the hazards of germination and removal are constant so the removal and germination times are Exponentially distributed with means μ R and μ G , (and it follows that γ = μ R /μ G , the assumption used from this point forward in article), this is equivalent to the form for β observed by Brookmeyer et al. (2005).
For an individual challenged with a dose D of spores, some fraction K (where K ∼ Binomial (D, φ) with number of trials D, and probability of success φ) will be deposited in the lung space. From this point, the competing-risks framework (as discussed above) may be applied. Alternatively one could construct a form of h G and h R in (1) that includes deposition directly, but for simplicity we continue with the former formulation explicitly coupling this faster-acting event. The parameter φ may vary depending on activity level of the host and particle size inhaled, but will be assumed not to vary between individual spores, while D will depend on (atmospheric) dispersion of the agent and so will depend on the individual's position relative to the source. Each inhaled spore i will then be subject to a Bernoulli test with probability φβ of causing infection in the host. The probability Q of removal of all spores without infection occurring is then which can be evaluated, provided with some estimate of the parameters D, γ , and φ, for an individual receiving the challenge. At a population level these values will vary, here we propose possible options for modeling this variation.

Aerosol Dispersion
Ordinarily one assumes that the number of spores in a unit volume of inhaled air is modeled as a Poisson process (Haas, Rose, & Gerba, 2014), that is we assume that spores are randomly distributed in the atmosphere local to the inhaler and with expected dose proportional to the amount of air inhaled, D ∼ Pois(d). In this case, we find the doseresponse function takes the form of an Exponential function where α = φ/(1 + γ ). If the expected inhaled dose d is overdispersed then the uncertainty in dose may be reasonably modeled by a -Poisson process (more regularly called a Negative-Binomial distribution), and the resulting form of Q is given by where the parameter ζ ∈ (0, ∞) (arising from the negative-binomial distribution) controls the level of overdispersion. This equation can be rewritten in the form of an exponential whereα = ζ ln(1 + φ/(ζ (1 + γ ))). Given that Equations (4) and (6) have the same functional form, fitted to a single experimental data set we clearly cannot tell the amount of overdispersion from the resulting value of the superparameter, without further contextual information.

Variation in Germination and Removal
We note from Equation (2) that γ is a ratio and so if we assume that the mean germination and mean removal times, μ G and μ R , respectively, are both distributed in a particular host such that then the ratio γ = μ R /μ G can be shown to have the probability density function The parameter β X (where X represents either G or R) is the index of dispersion and we define κ = β G /β R for future notational brevity. If the index of dispersion is the same for germination and removal, so κ = 1, Equation (4) becomes where p F q represents the generalized hypergeometric function (Slater, 1966) given by and (x) n denotes a rising factorial or Pochhammer symbol. Note, the generalized hypergeometric function will be used throughout this article. In Equations (4) and (6), the impact of overdispersion compared to a Poisson-distributed dose could not be determined as the models had the same functional form, however, if overdispersion in dose is included within (7), the dose-response relationship becomes The assumption that the dispersion indices for germination and removal are equal is returned to later in this article. It is worth noting that Equation (7) is called the "Beta-Poisson" model (Teunis & Havelaar, 2000) and when ρ → ∞ but λ ∼ O(1) then This approximate dose-response relationship is commonly applied in the literature and confusingly also called the Beta-Poisson model (Haas et al., 2014;Huang, Bartrand, Haas, & Weir, 2009;Tamrakar, Haluska, Haas, & Bartrand, 2011;Teunis & Havelaar, 2000).

Deposition
The parameter φ models the between-host variability in the dose that actually is retained following exposure. As this is a probability, it is sensible to assume this follows a beta-distribution (with associated parameters α φ and β φ ). Indeed, Heppell, Egan, and Hall (2017) find outputs from the multiple path particle dosimetry (MPPD) model (Price, Asgharian, Miller, Cassee, & de Winter-Sorkina, 2002) can be approximated well with a beta-distribution. Incorporation of this additional source of variability can considerably increase the complexity of the doseresponse model, depending on what other factors are considered and upon what restrictions are placed on these assumptions. With an overdispersed dose, and variation in germination and removal rates (albeit with equal index of dispersion), the dose-response relationship is defined by The parameters have been rewritten as ϕ = α φ /(α φ + β φ ) which is the mean of the given beta-distribution and ω = α φ + β φ which controls variation (note, the true variance can then be written as ϕ(1 − ϕ)/(1 + ω)). A slightly simpler relation is arrived at if the dose is assumed to be Poisson-distributed (i.e., not overdispersed):

Additional Models
All previous models have assumed that the index of dispersion is equal in germination and removal rates. If we remove this constraint, the models can quickly become intractable. However, where the dose is overdispersed and variation in deposition is not permitted (as in Equation (9)), allowing for an index of dispersion not equal to unity gives where F 3 is the Third Appell hypergeometric function (Appell, 1925) and we have defined ρ = α R + α G and λ = α R /(α G κ ) to be the ratio of the underlying distribution mean values.
Allowing for a beta-distribution to describe the variation in deposition modifies the dose-response function from Equation (13) to where ϕ is the mean of the beta-distribution that models φ and ω = α φ + β φ controls the associated variance of the distribution. This can be rewritten in the notation of a Kampé de Fériet function (Weisstein, Weisstein) given by This formula captures variability around individuals at roughly the same geographic location, though some of these parameters could be expected to change over greater distances. Note, both Equations (13) and (15) can be challenging to evaluate numerically.
A further two cases are included in this section for completeness. These are the dose-response models generated by assuming that no variation is present in removal and germination rates, though allowing for variation in deposition rates. These models are given by for Poisson-distributed, and -Poisson-distributed doses, respectively. There is little reason to expect these a priori, though it is interesting to note that they take the same functional form as Equations (7) and (9), albeit with a different interpretation of the parameters.

Summary and Implementation
A complete list of dose-response models discussed is given in Table I and associated parameters are given in Table II. In this article, the standard R functions are used for simpler models (Equations (4) and (6)), the gsl package in R is used for evaluating 1 F 1 and 2 F 1 functions (Equations (7), (9), (16), and (17)). When φ > ζ the final argument is larger than one and convergence may be an issue, in this situation formula 15.3.5 from Ratio of mean rates of removal and germination Superparameter-ratio of mean rates Ratio of dispersion indices for germination and removal κ = β G /β R Abramowitz and Stegun (1972) are used to improve convergence and compared with formula 15.3.4 and most stable result for dose regime considered is chosen. There are no further identified issue with numerical convergence of these functions for parameterizations inferred once one of these transforms is made. We call a Python library for 2F2 function (12) which suffers convergence for large parameter choices, in this study this is an issue for the anthrax data set below. We do not numerically evaluate the more complex functions derived for additional models given by (11), (13), or (15) in this study.

Application of Theory to Dose-Response Data Sets
In the previous section, equations were derived from the competing-risks framework that allow for variation in dose-dispersion, dose-deposition, and in the germination and removal processes. In typical dose-response experimental data, some of these factors may naturally vary and hence in attempting to describe the data, models should include them. However, other factors are controlled for and any extrapolation to real-world scenarios necessitates an understanding of these parameters. In this section, we investigate these additional factors.
The language used throughout this article has focused on the sporulating bacterium B. anthracis. We make use of data from two additional bacteria: Yersinia pestis, the causative agent of plague; and Francisella tularensis, the causative agent of tularemia. Descriptions of the mode of infection and the suitability of both agents for applying the modeling framework discussed in this study are given below. Data sets are identified for all three agents.
The pathogenesis of Y. pestis infection (by inhalation) is complicated and not fully understood. Several factors have been elucidated and center on the ability of Y. pestis to suppress and evade the host immune system in an initial anti-inflammatory phase (Li & Yang, 2008). To initiate infection, Y. pestis is initially phagocytozed by host leucocytes. If taken up by neutrophils bacteria are eliminated, whereas in macrophages they may resist cell-mediated destruction, express a number of virulence factors, then escape host cells with a capacity to proliferate with reduced risk of further phagocytosis (Cavanaugh & Randall, 1959;Lukaszewski et al., 2005). Thus success or failure in surviving macrophage destruction may determine whether infection persists or the bacterium is eliminated. This model of infection is analogous to that in the competing-risks model discussed throughout this article, and is thus an appropriate mathematical model to describe inhalational plague infection. Here we make use of two data sets that, to our knowledge, have not been used in published dose-response models before. Both data sets involve inhalational exposures of CO92 strain Y. pestis to Cynomolgus macaques. These data sets were obtained from Van Andel et al. (2008) and Warren et al. (2011); from this point forward we refer to the combined data sets simply as the plague data set. Fit to an exponential model (Equation (4)), the ID 50 (dose required to cause infection in 50% of the exposed, susceptible population) for this combined data set is 88 organisms. Data were collected by exposing 40 participants to individual doses.
In a previous study, Wood, Egan, and Hall (2014) demonstrate that a mechanistic birth-death-survival model for inhalational tularemia infection is approximated well by an exponential model and can be interpreted as a hit model. The hit could be the failed elimination of a bacterium following phagocytosis: an unsuccessful attempt to kill allows the bacterium to multiply to number several hundred within the macrophage, before rupture and release into the extracellular environment for the process to be repeated. Similarly to Y. pestis, this describes a competing-risks type framework (clearance vs. survival and proliferation) and thus we can use the models developed in this article. The data sets featured in Wood et al. are human inhalational experiments using the virulent SCHU-S4 strain bacteria (Saslaw, Eigelsbach, Prior, Wilson, and Carhart, 1961;Sawyer, Jemski, et al., 1966). Fit to an exponential model, the ID 50 for this combined data set is 10 organisms. These data only include a relatively small range of doses and so are unlikely to support more complicated models when making inference directly but in this context the effect of structure in the models is being shown. The Saslaw study exposed 20 volunteers to a range of doses, with doses of 10, 20, 23, and 46 organisms being received by two participants each (Saslaw et al., 1961).
A third comparison data set is obtained from the experiment of Altboum et al. (2002). Here, guinea pigs are exposed inhalationally to Vollum strain B. anthracis. Fit to an exponential model, the ID 50 for this data set is 1.5 × 10 5 organisms. In total, 51 Guinea pigs were used in study with 4, 8, 12, 12, 8, and 7 receiving each dose in increasing order.

Considering Overdispersion in Dose
In typical laboratory-controlled dose-response experiments (including the data sets we have chosen to use), we may assume that dosages are unstructured, i.e., they are Poisson-distributed. It is also reasonable to assume that deposition rates show little variation between hosts since activity levels will be similar. On this basis, the dose-response relationship given by Equation (7) may be an appropriate choice to describe our data. We thus begin by fitting our three experimental data sets to this model. Best-fit parameters obtained for each data set and deposition value are given in Table III. To investigate the impact of overdispersion in the dose, we substitute these best-fit parameters into the model given by Equation (9), and vary the unknown overdispersion parameter ζ , on an order of magnitude basis (recalling that ζ ∈ (0, ∞)). It is important to notice that changing ζ in this way shows the impact of such overdispersion rather than fitting the overdispersed model to the data, having the mechanistic structures allows prediction to be made of the effect of overdispersion.
Three values of ζ are considered (ζ = 10, 1, 0.1). Fig. 1 illustrates for the three featured data sets (left to right: anthrax, plague, and tularemia) and a deposition rate of 0.1, the effects of varying the dispersion parameter ζ . Figures for the other depositiondispersion combinations are near identical to their parent models and hence omitted for brevity, however Table IV lists the response rates at the ID 50 of the original (parent) model for all parameter combinations.
For the anthrax data set, there is no discernible difference between dose-response relationships that allow the overdispersion parameter ζ to vary; this is true regardless of the level of deposition assumed. For plague and tularemia (which have lower ID 50 s), allowing for overdispersion yields dose-response relationships that predict lower infection rates across all doses though the magnitude of the effect for most cases is modest (see Table IV). The impact of the overdispersion appears to increase as each of the dispersion parameter, deposition rate, and ID 50 of the data set shrinks. Only in the extreme, i.e., for the tularemia data set with a deposition rate of 0.1, and dispersion parameter ζ = 0.1 do we see a noteworthy difference; the response rate at the previous ID 50 has reduced to 38.64% (rounded to four significant figures).
In most cases then, the added complexity of including overdispersion is not warranted, and the simpler parent model given by Equation (7) is valid, however when dealing with a highly infective organism, or a very low deposition rate, it may be important to consider the level of overdispersion.
It is interesting to note that in graph C (deposition value of 0.1, tularemia data set) the response rate appears not to tend to unity as the dose increases. This is not a numerical artifact but instead an artifact that the maximum-likelihood estimate for ρ 1 in which case (7) tends to λ/(1 + λ). Exploration of the likelihood function in the parameter-space suggests    that the given parameter values do indeed maximize the likelihood, though an alternative parameter with smaller likelihood tend to unity in a reasonable dosedependent manner.

Considering Variation in Deposition
The impact of variation on deposition can also be considered independent of overdispersion. Equation (12) was generated by considering a beta-distribution on the deposition parameter φ. The parameter ϕ is fixed by the mean deposition rate, however the remaining parameter ω that controls the variance, is free to vary.
The variation in deposition parameter ω is considered for three values: 2, 20, and 200, where increasing values of ω indicate reduced variability in the deposition rate. Note that these values of ω were chosen to be at least 2 to prevent the underlying betadistribution from being bimodal with local maxima at zero and unity; this would not be expected biologically for these causative agents, but in generality is this framework is applied to other infections where some inherent protection is extant in population the limit on ω may be relaxed. Fig. 2 illustrates the impact of incorporating variation into the deposition parameter on the doseresponse relationship, and Table V details    Note, only the data sets for plague and tularemia are present since the parameterization for the anthrax data set meant that all implementations of the hypergeometric function 2 F 2 could not converge numerically.
The differences observed when considering uncertainty on the deposition parameter are visually evident in both data sets from a deposition level of 0.5; they are however still relatively modest, for instance, the new response rate for the plague data set with a deposition rate of 0.5 and ω = 2 is 45.12%. A significant reduction in response is observed for a deposition rate of 0.1 and ω = 2 for both the plague and tularemia data sets with response rates of 27.43% and 27.83%, respectively. These results suggest that as the deposition rate decreases, the uncertainty around that deposition plays a larger role in determining the probability of infection.
For the purposes of this article, we have provided examples for deposition and dispersion independently. It is possible to consider the combined effects of variation using Equation (11), however the visualization becomes increasingly complex as the dimensionality increases.

DISCUSSION
A central paradigm in mathematical modelling is to build the simplest model that describes the data. While this may generally be true, if the simplest doseresponse models fit solely to laboratory-controlled data are used to make inferences regarding realworld exposures, there is the possibility for an overestimation of the number of infected cases.
In this work, we have extended Brookmeyer's competing-risks model of inhalational anthrax infection to consider how variation or uncertainty may impact upon the predicted response rate. Candidate models have been generated by considering general distributions for prescribed sources of variation. Parent models have been parameterized from real inhalational data sets for a variety of organisms, and the parameters used to inform the more complex child model. From this, we have been able to speculatively explore the impacts that population heterogeneity or parameter uncertainty may have on infection probabilities.
In the examples we have shown, allowing for variation in overdispersion has had minimal impact on the dose-response relationship in all but the most extreme case. The results suggest that the impact increases as the ID 50 for the data set, the deposition rate, and the overdispersion parameter ω decrease. This could be regarded as a high-infectivity effect, and thus, consideration of whether this model provides any added benefit over the parent model depends upon the organism being modelled. In contrast to dispersion, the impact of variation on deposition appears to be important at higher deposition levels, and potentially independent of the ID 50 of the data set. As the deposition rate decreases, increasing variability in the deposition rate results in a smaller response rate. Again, assessing whether this model should be used may be dependent on the deposition level, and hence may be organismand host-dependent.
While the framework initially described is formulated for inhalational anthrax infection, we have shown examples in tularemia and plague for which the same model could be applied. The central concept binding all three examples is that in the case of independent action, organisms may behave in a manner such that a single bacterium has two possible fates: removal (from the system, by any means) or survival with subsequent proliferation resulting in infection. The focus on inhalation as route of infection in this study is not necessarily a barrier to translation to other modes of transmission, such as ingestion or cutaneous infection. The methods may have a broader applicability worthy of further study.
To our knowledge, there are no published doseresponse models with application to human inhalational plague (i.e., applied to data of inhalational or intranasal exposure of human or nonhuman primate hosts). Using the same concept as the competingrisks model for inhalational anthrax infection, a model for inhalational plague parameterized from primate dose-response data has been posited. Further work is required to validate the modeling assumptions made, however in the absence of other models, we recommend this model be used for inhalational plague risk assessment. Note that the estimates made by this model are subject to the same concerns regarding all models that attempt to explain only the experimental data; they may overpredict infection estimates when variation is not taken into account.
Throughout the model generation process, different realizations of variation have yielded models with the same functional form. For example, assuming Poisson-distributed dispersion of organisms, a beta-distribution on the deposition parameter would produce a model with the same functional form as a model allowing for gamma distributions on the variation in germination and removal rates (with equal dispersion index). Furthermore, some authors (Teunis & Havelaar, 2000) consider an exponential model for dose-response where the unknown parameter is simply the probability of an individual bacterium causing infection, and a beta-distribution on this parameter allows for population heterogeneity at the level of either host or bacterium. Each interpretation produces models with the same functional form, and hence model fit alone is insufficient in providing evidence that a specified set of assumptions gives rise to observed behavior. Clearly, some understanding of the experimental conditions (and therefore model selection) is required if model users wish to make inferences about the processes that dictate infection events.
Mathematical tractability has been a problem throughout this article. For instance, the hypergeometric functions are often computed as a power series, and subject to computational overflows when certain parameters are sufficiently large. One measure of whether 2 F 2 may be simulated for a data set is in the ID 50 . For large ID 50 s, the infection generating event rate must be small (relative to the removal rate), resulting in at least one large parameter. This was the case for the anthrax data set, and hence the model was omitted in Fig. 2. There is therefore a tradeoff in terms of model complexity versus practicality. This may not be the case if alternative descriptions of variation are used that do not result in a Gaussian hypergeometric function and numerical and mathematical representation is subject to further work.
The modeling assumptions we have made to generate a set of candidate models require further work to validate their appropriateness and so this analysis can be considered as illustrative only. However, the assumptions made have used general functions and so serve some purpose beyond simple illustration. Further work is required to consider if the distributions we have considered capture real-world variability, or if alternatives need to be investigated.
In the context of emergency response and planning, often the most risk-averse scenario is considered. Here, this equates to selecting the competingrisks model (with the same functional form as an exponential model), as it gives a higher number of infections than other candidate models discussed in this article. However, the uncertainty needs to be considered when assigning reasonability to the worstcase scenario. If the assumption that deposition rates can be described by a nonbimodal beta-distribution is valid, then (as previously discussed) simulations for which ω = 2 may provide a bound on this uncertainty. In addition, the analysis presented here suggests when uncertainty and variation should be considered, that is, for highly infective organisms, and for low deposition rates. The work presented here examining uncertainty and variation in the absence of data therefore contributes to the overall risk assessment.

CONCLUSION
The competing-risks model of bacterial infection can be extended to incorporate sources of variation that are often controlled for in laboratories yet likely to arise in real-world scenarios. Model realizations involving additional variational parameters can be partially parameterized from their simpler parent models fitted to real data. Exploration of remaining unknown parameters has shown that model predictions may overestimate infection numbers when simpler models fit to laboratory-controlled data are used in isolation.