**DOI:**10.1128/AEM.01190-10

## ABSTRACT

*Francisella tularensis* can be disseminated via aerosols, and once inhaled, only a few microorganisms may result in tularemia pneumonia. Effective responses to this threat depend on a thorough understanding of the disease development and pathogenesis. In this study, a class of time-dose-response models was expanded to describe quantitatively the relationship between the temporal probability distribution of the host response and the *in vivo* bacterial kinetics. An extensive literature search was conducted to locate both the dose-dependent survival data and the *in vivo* bacterial count data of monkeys exposed to aerosolized *F. tularensis*. One study reporting responses of monkeys to four different sizes of aerosol particles (2.1, 7.5, 12.5, and 24.0 μm) of the SCHU S4 strain and three studies involving five *in vivo* growth curves of various strains (SCHU S4, 425, and live vaccine strains) initially delivered to hosts in aerosol form (1 to 5 μm) were found. The candidate models exhibited statistically acceptable fits to the time- and dose-dependent host response and provided estimates for the bacterial growth distribution. The variation pattern of such estimates with aerosol size was found to be consistent with the reported pathophysiological and clinical observations. The predicted growth curve for 2.1-μm aerosolized bacteria was highly consistent with the available bacterial count data. This is the first instance in which the relationship between the *in vivo* growth of *F*. *tularensis* and the host response can be quantified by mechanistic mathematical models.

*Francisella tularensis* is the causative agent of tularemia. It is an intracellular pathogenic species of Gram-negative bacteria, replicating mainly in macrophages, and has also been found in amoebae (29). Interest in this pathogen was raised due to its high infectivity, ease of dissemination, and consequently potential use as a biological weapon (5, 18, 33). It can be easily disseminated via aerosols that once inhaled may result in tularemia pneumonia, a severe form of disease with high mortality if untreated (27). Two primary subspecies, type A and type B, have been classified. For type A *F. tularensis*, known as one of the most infectious pathogens, only a few organisms may cause infection (24, 25). The U.S. Centers for Disease Control and Prevention have classified *F. tularensis* as a category A bioterrorism agent for public health preparedness. Given the rapid-progression characteristics of inhalation tularemia, effective responses to this threat depend on a thorough understanding of its likelihood and temporal course. Since the survival and growth of the pathogens are believed to be the cause of the disease, the quantification of the relationship between host response and bacterial kinetics becomes critical. Previously, to quantify microbial growth in culture media and food, the modified Gompertz model and the Baranyi model were widely studied (1, 2, 34). However, these prior studies only described the bacterial growth curve *in vitro* and did not quantitatively associate the *in vivo* microorganism growth with the probability of a host response. Quantitative microbial risk assessment that estimates the risk of adverse consequences from exposure to certain doses of pathogens has been developed by researchers (4, 6, 12, 28). Huang et al. (15, 16) suggested a class of time-dose-response (TDR) models by incorporating time dependencies that model the *in vivo* bacterial kinetics into the classical dose-response models. The resulting models demonstrated adequate flexibility in quantifying the time to response caused by disparate pathogens. However, the authors did not provide a direct verification of the biological validity of the models with *in vivo* pathogen growth. This study aimed at such verification with data from monkey, which has been considered the animal model that most closely mimics inhalational diseases in humans and which was used in past dose-response studies for category A agents (3, 17). Open literature was searched for survival dose-response data and bacterial viable count data for monkeys infected by *F. tularensis* via the inhalation route. The suggested models by Huang et al. were further developed to model quantitatively the hypothetical relationship between the host response and instantaneous microorganism numbers *in vivo*. The resulting models were fit to survival dose-response data, and the estimates of bacterial dynamics for different aerosol sizes were inferred. The estimates were then compared with the data for bacterial growth and the clinical and pathological findings for verification. Variation in the median lethal dose (LD_{50}) with aerosol size and exposure route was also discussed.

## MATERIALS AND METHODS

Dose-dependent survival data.Day and Berendt (7) exposed 4- to 5-kg monkeys (*Macaca mulatta*) to the SCHU S4 strain of type A *F. tularensis.* The bacteria were administered in the forms of aerosol particles of different sizes. The aerosols were produced with a University of Chicago Toxicity Laboratories (UCTL)-type atomizer (22) fitted with an Aerotec tube (13) (for 2.1-μm particles) or a spinning top (20) (for other particle sizes). The number of bacteria inhaled by each animal was controlled by varying the exposure time. The survival data are shown in Table 1 .

*In vivo* viable count data.To quantify the biological responses of animals, measurements of the bacterial burden in various organs are typically conducted after sacrificing animals and aseptically removing their organs. Since *F. tularensis* predominantly infects lung cells when administered via the inhalation route, in this study only the *in vivo* bacterial count data for entire lungs of monkeys post-respiratory exposure were analyzed. While the SCHU S4 strain is the preferable pathogen model for consistency with the available animal survival data, the data for other strains of *F. tularensis* were also located for comparison purpose. Through an extensive literature search, three studies with adequate data were found and were summarized as follows.

White et al. (31) exposed rhesus monkeys (*Macaca mulatta*) with an average weight of 3.5 kg to 1- or 8-μm particles of the aerosolized SCHU S4 strain. These aerosol particles were generated with a modified Henderson apparatus by either a Collison spray head (14) (for 1-μm particles) or a vibrating reed (32) (for 8-μm particles). Only the data of viable counts for the 1-μm particles were presented in their work.

Eigelsbach et al. (9) infected monkeys (*Macaca irus*) weighing 1.8 to 4.9 kg to aerosols of the SCHU S4 strain and an attenuated live vaccine strain (LVS), respectively. The particles of aerosols in a range between 1 and 5 μm in diameter were disseminated with a nebulizer.

Schricker et al. (26) exposed rhesus monkeys (*Macaca mulatta*) weighing 2 to 3.5 kg to the aerosolized 425 strain of type B *F. tularensis*, with particles ranging from 1 to 5 μm, produced by a nebulizer.

There were no time-to-response data reported for these studies. The *in vivo* viable count data are normalized by the inhaled dose and shown in Table 2. For the data given only in graphical format, a digitizing software program, Engauge 4.1 (http://digitizer.sourceforge.net), was used to extract the information. The relative errors for extracted data points are generally less than 1% based on the absolute errors given by the digitizing software.

Other data.A set of traditional dose-response data for *F. tularensis* infection via oral exposure in the work of Quan et al. (21) are used in Discussion for investigating the effects of the inoculation route on the response. In the original study, albino mice were infected orally with drinking water contaminated with 10^{4} to 10^{8} organisms of a highly virulent strain, Aa. Mortalities were reported based on groups of 22 mice, except that 11 animals were administered only 10^{6} cells. Data points are plotted in Fig. 5.

All the data used in this study are based on test animals that were neither immunized nor administered prophylaxis, so the resulting models should be valid only for previously unexposed populations.

TDR model.A class of time-dose-response (TDR) models incorporating the time postinoculation into the classical dose-response models for microbial infection has been suggested by prior studies (15, 16). The parameter *k* in the exponential dose-response model and the parameter *N*_{50} in the beta-Poisson model (see below) were set equal to functions of time that represent *in vivo* bacterial kinetics.

The exponential TDR model can be given as
$$mathtex$$\[F(d,t){=}1{-}e^{{-}G(t;{\theta},{\ldots})k_{0}d}\]$$mathtex$$(1) The beta-Poisson TDR model can be given as
$$mathtex$$\[F(d,t){=}1{-}\ \left[1{+}\frac{d}{\frac{N_{50}}{G(t;{\theta},{\ldots})}}{\times}(2^{\frac{1}{{\alpha}}}{-}1)\right]^{{-}{\alpha}}\]$$mathtex$$(2) *F*(*d*, *t*) is the probability of a positive response (death). *d* is the average exposed dose. *t* is the time postinoculation. *G*(*t*; θ,…) is a cumulative distribution function that reflects the *in vivo* kinetics of a single microorganism; the parameter *θ* is the time parameter. *k _{0}* is the dose parameter for the exponential model, describing the probability that a microorganism can remain, survive, and multiply in the host body in order to produce an infectious focus (response);

*G*(

*t*;

*θ*,…)

*k*

_{0}is therefore the cumulative probability that an organism survives and multiplies in order to initiate a response at time

*t*postexposure. For the beta-Poisson model,

*α*is the exponential fitting parameter;

*N*

_{50}, equivalent to LD

_{50}, is the dose required to produce a response in 50% of the subjects.

*N*

_{50}/

*G*(

*t*;

*θ*,…) is the dose required to produce a 50% response at or before time

*t*. As

*t*approaches infinity, the TDR models become the classical dose-response models (equation 1 or 2 with

*G*(

*t*;

*θ*,…) equal to 1).

Modeling bacterial growth.In the current study, we expanded these models to examine theoretically the relationship between *in vivo* bacterial growth, the host response, and the exposed dose. Denoting the first partial derivatives of equations 1 and 2 with respect to *t*, which represent the probability density function of the time to response as *f*(*d*, *t*), a hazard rate function, *λ*(*d*, *t*), which presents the conditional response probability of an unresponsive subject(s) at time *t*, can be given as
$$mathtex$$\[{\lambda}(d,t){=}\frac{f(d,t)}{[1{-}F(d,t)]}\]$$mathtex$$(3) where *f*(*d*,*t*) = ∂*F*(*d*,*t*)/∂*t*.

For the exponential TDR model, we solve *λ*(*d*,*t*) based on equations 1 and 3 and obtained
$$mathtex$$\[{\lambda}(d,t){=}(k_{0}d)g(t;{\theta},{\ldots})\]$$mathtex$$(4) For the beta-Poisson TDR model,
$$mathtex$$\[{\lambda}(d,t){=}\ \left[1{+}\frac{G(t;{\theta},{\ldots})d}{N_{50}}{\times}(2^{\frac{1}{{\alpha}}}{-}1)\right]^{{-}1}\left[{\alpha}(2^{\frac{1}{{\alpha}}}{-}1)\frac{d}{N_{50}}\right]g(t;{\theta},{\ldots})\]$$mathtex$$(5) *g*(*t*; *θ*,…) is the probability density function of *G*(*t*; *θ*,…). Assuming the conditional response probability is proportional to the instantaneous number of microorganisms *in vivo*,
$$mathtex$$\[{\lambda}(d,t){=}kN\]$$mathtex$$(6) where *N* is the instantaneous number of microorganisms at a given time, and *k* is the coefficient of host-pathogen interaction.

Based on equations 4 to 6, we solve for the normalized *in vivo* number *N*/*d*, which reflects bacterial growth, as follows.

For the exponential TDR model,
$$mathtex$$\[\frac{N}{d}{=}\frac{k_{0}}{k}g(t;{\theta},{\ldots})\]$$mathtex$$(7) For the beta-Poisson TDR model,
$$mathtex$$\[\frac{N}{d}{=}\ \left[1{+}\frac{G(t;{\theta},{\ldots})d}{N_{50}}{\times}(2^{\frac{1}{{\alpha}}}{-}1)\right]^{{-}1}\left[{\alpha}(2^{\frac{1}{{\alpha}}}{-}1)\frac{1}{N_{50}k}\right]g(t;{\theta},{\ldots})\]$$mathtex$$(8) For the exponential model, it is shown by equation 7 that both the magnitude and temporal distribution of the normalized number are independent of the initial dose, and the distribution is determined solely by *g*(*t*; θ,…), the kinetics of a single microorganism. This is consistent with the assumption of the classical exponential dose-response model that each ingested organism acts independently of the others and has an equal probability of producing an infectious focus (12).

Compared with the exponential model, the trend of bacterial growth described by the beta-Poisson model (equation 8) is more complicated due to the term $$mathtex$$\({\{}1{+}[G(t;{\theta},{\ldots})d/N_{50}]{\times}(2^{\frac{1}{{\alpha}}}{-}1){\}}^{{-}1}\)$$mathtex$$. When the dose is sufficiently low to make $$mathtex$$\([G(t;{\theta},{\ldots})d/N_{50}]{\times}(2^{\frac{1}{{\alpha}}}{-}1){\ll}1\)$$mathtex$$, equation 8 can be simplified to yield $$mathtex$$\[\frac{N}{d}{=}\ \left[{\alpha}(2^{\frac{1}{{\alpha}}}{-}1)\frac{1}{N_{50}k}\right]g(t;{\theta},{\ldots})\]$$mathtex$$(9) When the dose becomes very large, $$mathtex$$\([G(t;{\theta},{\ldots})d/N_{50}]{\times}(2^{\frac{1}{{\alpha}}}{-}1){\gg}1\)$$mathtex$$, equation 8 can be simplified as $$mathtex$$\[N{=}\frac{{\alpha}}{k}\ \frac{g(t;{\theta},{\ldots})}{G(t;{\theta},{\ldots})}\]$$mathtex$$(10) Equations 8 to 10 for the beta-Poisson TDR model demonstrate that at a low dose (equation 9), the bacterial growth curve follows the kinetics of a single inoculated microorganism. With increasing exposed doses, the growth tends to diminish, indicating a decreasing average growth capacity of microorganisms. When the dose becomes sufficiently large, both the magnitude and the distribution of bacterial numbers become independent of the initial dose as shown by equation 10, a monotone decreasing function. These characteristics indicate a competing interference between microorganisms, which is consistent with the assumption of variable pathogen-host survival probability for the classical beta-Poisson dose-response model (12).

Statistical methods.In this study, the TDR models (equations 1 and 2), with *G*(*t*; *θ*,…) equal to candidate functions including the inverse-Weibull (16), inverse-exponential (16), the lognormal (19, 23), Weibull (10) and gamma (15, 28) distributions, were fit to the data from Table 1. Model parameters were optimized, and the deviances were determined via the method of maximum-likelihood estimation (MLE) (12), employing a binomial distribution as described in previous studies (15, 16). The MLE was performed using the R programming language (www.r-project.org, accessed on 8 August 2007). Models are considered to exhibit a statistically acceptable fit if the model minimized deviances were less than the 95% confidence value of the *χ*^{2} distribution (*χ*^{2}_{0.95,} * _{df}*;

*df*=

*m*−

*n*, where

*df*is the degrees of freedom,

*m*is the number of doses, and

*n*is the number of parameters). For models with the same number of parameters, the best-fit model is the one that yielded the lowest deviance. A model with more parameters is consider to be a better model than one with fewer parameters if the reduction of deviance by applying the additional parameters was greater than

*χ*

^{2}

_{0.95,}

_{Δ}

*, where Δ*

_{df}*df*is the difference of degrees of freedom between the two models. Confidence intervals for best-fit models can be determined using bootstrap analysis in R with samples drawn from the data sets.

Comparison of predicted bacterial kinetics with *in vivo* data.Of the parameters on the right side of equations 7 and 8, all can be determined by fitting the TDR models to the host survival data except the host-pathogen interaction coefficient *k*. Since the parameter *k* affects only the magnitude but not the distribution of the given equations, it is feasible that the distribution of bacterial growth can be inferred based on the host response. For verification, we compared the observed normalized viable counts to the estimate for 2.1-μm aerosols, given the similar sizes of inhaled aerosols.

## RESULTS

General fitting results.Equations 1 and 2 with candidate *G*(*t*; *θ*,…) were fit to the data listed in Table 1. In all cases, the beta-Poisson TDR model did not provide a statistically significant improvement in fit over the exponential. The deviances and parameter estimates of the best-fit models are shown in Table 3. Statistically acceptable fits were exhibited for all the individual data sets. The exponential TDR model incorporating the Weibull distribution provided the best fit to the responses caused by 2.1- and 7.5-μm particles, while it exhibited the best fit to those initiated by 12.5- and 24.0-μm particles by incorporating the gamma and inverse-Weibull distributions, respectively. In Fig. 1, the probability density function of the best-fit model for 2.1-μm particles is compared with the observed densities of deaths at different dose levels. It can be seen that the TDR curves are well consistent with the data. The cumulative function of the incorporated Weibull distribution is given by
$$mathtex$$\[G(t;{\theta}_{1},{\theta}_{2}){=}1{-}e^{{-}\ \left(\frac{t}{{\theta}_{2}}\right){\theta}_{1}}\]$$mathtex$$(11) where *θ*_{1} is the shape parameter and *θ*_{2} is the scale parameter.

Size effect.For the exponential TDR models, the temporal distributions of conditional response probability and bacterial growth are conceptually determined solely by *g*(*t*; *θ*,…), as shown by equations 4 and 7. The *g*(*t*; *θ*,…) value estimated for the aerosol exposures to the four different-size particles (see Table 3 for the parameters) were plotted in Fig. 2. It is noticeable that the estimate for the 2.1-μm aerosols is identical with that for the 7.5-μm aerosols, while the *g*(*t*; *θ*,…) values for the larger aerosol particles have a trend of right shift and widening that indicates a less-acute disease development. This parallels the pathophysiological and clinical observations in the original study: for animals exposed to the 2.1- and 7.5-μm particles, the infection occurred primarily in the lower respiratory tract (respiratory bronchiole) and the symptoms of animals were found to be similar, whereas for those exposed to the aerosols of 12.5 μm and 24.0 μm, the major infection site was the upper respiratory tract and milder clinical symptoms were seen at onset than were observed for the small particles (7). The clinical and pathological similarities in animals inhaling aerosols of *F. tularensis* smaller than 8 μm were also observed in other prior studies. Goodlow and Leonard (11) and White et al. (31) both reported similar pathological observations and infection site (bronchiole) in rhesus monkeys exposed to 1-μm and 8-μm aerosolized strain SCHU S4. White et al. (30) found that the respiratory bronchiole was also the initial site of infection caused by the aerosol particles of the live vaccine strain with a diameter of 5 μm or less.

The LD_{50}s for the four particle sizes are shown in Table 3. It can be seen that the LD_{50} value increases with the size of the particles. The variation in the LD_{50} may result from particle deposition and retention capabilities that differ with size as well as various infection sites.

Comparison of estimated bacterial kinetics with *in vivo* data.Taking the log_{10} of both sides of equation 7 yields
$$mathtex$$\[\mathrm{log}\ \left[\frac{N}{d}\right]{=}\mathrm{log}[g(t;{\theta},{\ldots})]{+}\mathrm{log}\left(\frac{k_{0}}{k}\right)\]$$mathtex$$(12) In Fig. 3, we plot the observed log (*N*/*d*) value for the three bacteria strains versus the log [*g*(*t*; *θ*_{1}, *θ*_{2})] value, where *g*(*t*; *θ*_{1}, *θ*_{2}) is the density of the Weibull distribution for the 2.1-μm aerosols of strain SCHU S4. The scatter points for the SCHU S4 strain after the 12th hour postinoculation exhibit a strong linear association with the slope of the trend line by linear regression equal to nearly 1 (*r*^{2} = 0.92), which is markedly consistent with the prediction by equation 12 and consequently supports the view that the time dependency of the model is driven by bacterial kinetics. Since the *y* intercept of the trend line is conceptually equal to log (*k*_{0}/*k*), the *k* value for the exponential growth phase can be extrapolated. Noticeably, this consistency is not seen within 12 h postexposure, during which the observed bacterial numbers are considerably larger than expected. This discrepancy may result from an inconstant host-pathogen interaction coefficient *k* during the lag phase. It is well known that shortly after entering the host body, the pathogens first need to survive by accommodating to the environment before growing and producing any deleterious effects, so the coefficient *k* may be significantly lower than at the latter stages. Hence, with the assumption of a constant *k* value, the bacterial number during the lag phase may be underestimated. It can also be noted in Fig. 3 that the more virulent strain apparently has a higher growth rate, which demonstrates that the virulence of various strains may directly relate to their abilities to grow *in vivo*.

## DISCUSSION

This is the first report showing that the observed *in vivo* growth of *F*. *tularensis* and the host response are associated consistently in a form predicted by a class of biologically plausible models. This preliminary success in developing mechanistic models for tularemia yields information on the mechanism by which the infection develops and potential strategies for controlling it.

Aerosol particle size has also been identified as one of the most critical factors influencing the infectivity of *Bacillus anthracis*, one of the other category A bioterrorism agents. Druett et al. (8) exposed guinea pigs and rhesus monkeys to the Vollum strain of *B. anthracis* at aerosol sizes ranging from 1 to 12 μm and recorded their mortalities. Bartrand et al. (3) analyzed the dose-response data from Druett's study and estimated the LD_{50}s for the different-size aerosols. Variations of the normalized LD_{50} with aerosol size for anthrax (Bartrand's study) and for tularemia (this study) are illustrated in Fig. 4. Noticeably, for both diseases, the most significant variation in the LD_{50} with size occurs in the range of 4.5 to 8 μm, indicating that the deposition and retention abilities of aerosols in the lung may drop dramatically within this size range. In addition, it can be seen that the size dependency in respiratory tularemia is stronger than that in anthrax, which possibly illuminates the essentiality of the infection site and tissue susceptibility to the pathogenicity of *F. tularensis*. This can be further tested by a comparison with infections via other exposure routes. By fitting the classical exponential dose-response model [equation 1, with *G*(*t*; *θ*,…) = 1] to the data for oral infection with the Aa strain (Fig. 5), a LD_{50} of 5.2 × 10^{6} was estimated. Noting an LD_{50} of 3 to 5 organisms that has been previously reported for intraperitoneal injection using the same strain (15), as well as the LD_{50} for respiratory infection determined in the current study, a substantial variation in virulence with infection site is manifestly present. Further research is needed to reveal the virulence factors and biological mechanisms accounting for such diversity.

## ACKNOWLEDGMENTS

This research was funded through the Center for Advancing Microbial Risk Assessment, supported by the U.S. Environmental Protection Agency and Department of Homeland Security under the U.S. EPA Science to Achieve Results (STAR) grant program (grant number R83236201).

This work does not necessarily reflect the official positions of the above-listed agencies.

## FOOTNOTES

- Received 18 May 2010.
- Accepted 15 November 2010.
- Accepted manuscript posted online 29 November 2010.

- Copyright © 2011, American Society for Microbiology