A Model for In-vivo Delivered Dose Estimation for Inhaled Bacillus

ARTICLE pubs.acs.org/est

A Model for In-vivo Delivered Dose Estimation for Inhaled Bacillus anthracis Spores in Humans with Interspecies Extrapolation Mark H. Weir*,† and Charles N. Haas‡ †

Department of Fisheries and Wildlife, Michigan State University, 303 Manly Miles Building, East Lansing, Michigan 48824, United States ‡ Department of Civil Architectural and Environmental Engineering, Drexel University, Philadelphia, Pennsylvania, United States

bS Supporting Information ABSTRACT: The quantitative yardstick for quantitative microbial risk assessment (QMRA) is the dose response assessment phase. In this phase of the QMRA paradigm a mathematical model is used to describe the relationship between host response (infection, disease, etc.) and pathogen dose. There are, however, key uncertainties which if addressed can expand our understanding of the dose response relationship and improve its accuracy. The dose response models most frequently used in this phase of QMRA are based on the average exposed dose (i.e., inhaled, ingested, etc.). However once inhaled, spores are considered infectious after being transported to a specific region of the lungs (alveoli), therefore, average exposed dose does not account for this required spore transport through the respiratory system. It is the aim of this manuscript to develop a model for the in vivo delivered dose to the alveolated region of the lungs that accounts for losses of spores through the respiratory system. A stochastic system is used to account for the physics in the respiratory system that account for the various sinks during respiration. This stochastic system is then integrated into the exponential and beta Poisson dose response models. The stochastic model is also then expanded to the respiratory systems of guinea pigs and rhesus macaques as these are common animal models. This work develops a framework for a new class of dose response models accounting for host physiology, making progress to understanding dose response heterogeneity among hosts.

1. INTRODUCTION Quantitative microbial risk assessment (QMRA) has proven itself an invaluable tool in decision and policy making. Within QMRA, the dose response assessment is the quantitative yardstick with which a risk level that can be characterized is estimated. The dose response assessment is best described as the fitting or development of a mathematical model describing the expected response from a host when exposed to a dose of a pathogen. Two dose response models have been used extensively in order to describe the dose
response relationship of various pathogens.1 Equation 1 shows the exponential model, with a single parameter (k) and eq 2 shows the beta Poisson model with two parameters (R and N50). The parameters for the models are optimized using animal or human controlled dosing trials. Therefore the optimization of the parameters is based on the average exposed dose to the host, which does not account for potential losses in the respiratory system (or gastrointestinal system for water or foodborne pathogens). Consequently there is the potential to address heterogeneity observed in the modeled dose response relationship by accounting for these losses in the system. This amounts to developing a means of estimating the in vivo delivered dose the host experiences by accounting for the physics r 2011 American Chemical Society

of the respiratory system and effect on the exposed dose. Exponential : pðResponseÞ ¼ 1
expð
k 3 doseÞ

ð1Þ

Beta Poisson : pðResponseÞ "

#
R dose 3 21=R
1 ¼ 1
1þ N50

ð2Þ

Dose response models for Bacillus anthracis (B. anthracis), the causative agent of anthrax, have been developed previously.2 Specifically, Bartrand et al. 2 developed a quantitative depiction of the effect of spore diameter on the dose response of B. anthracis for multiple host species. With the understanding of the effect of spore size this led to the need to develop a model for effective dose based on transport physics of spores through the respiratory system. The scope of this work is to develop an easily adaptable Received: March 17, 2011 Accepted: May 31, 2011 Revised: May 24, 2011 Published: June 13, 2011 5828

dx.doi.org/10.1021/es200901e | Environ. Sci. Technol. 2011, 45, 5828–5833

Environmental Science & Technology model of the fate and transport through the respiratory system of mammals, specifically humans. The structure of the model framework is desired to be easily modified so as to expand its use to additional host species or exposure pathways in the future. Therefore a stochastic model that simulates the physics of spore fate and transport through the human respiratory system was chosen. In conjunction with the modeling of the physical systems, we have developed a correction factor which converts the average exposed dose to an in vivo delivered dose, that can be directly included in the dose response models. As this dose correction factor depicts the physics of respiration and deposition of inhaled spores, this establishes a framework for physiologically based dose response models (PDRM) in QMRA. This also marks a step in understanding the heterogeneity of the dose response between and within host groups.

2. MODELING METHODS 2.1. Fundamentals of Infection and Physiology. The infection process can be best described as a two-stage process.3 The two-stage approach depicts both the transport and pathogenesis in a combined model, thus allowing for the entire infection process to be captured by the model. For B. anthracis the inhaled spores must reach the alveoli in order to germinate and replicate, thus being capable of developing an infection.4,5 The model in this paper depicts the first of the two stages, the fate and transport of the inhaled spores. Modeling the full complexity of the respiratory system, accounting for all of the diameter changes and intermediary bifurcations is considered computationally infeasible.6 Prior work has attempted this for the upper respiratory tract, a portion of the bronchioles and to a limited extent the individual alveoli. However when extended to terminal bronchioles, all alveoli or attempts to combine these models into a single model for the entire respiratory system are considered computationally infeasible.7
9 Therefore system simplification is required to model the entire respiratory system in a computationally tractable manner. A simplification of the respiratory system must be based on how the respiratory system has evolved as a compartmentalized system. A three-region model (Table 1) of the human respiratory system was developed for advanced radiological dosimetry.9 This three-region model was developed based on physiological compartmentalization, where region 1 (R1) is from the nares (nostrils) to the back of the throat, region 2 (R2) is from the back of the throat to just beyond the main bronchi, and region 3 (R3) is from the beginning of the bronchioles to the alveoli. The model is built as a Markov chain model where each of the regions are the primary states of the model and additional states are constructed based on flow through and deposition in the respiratory system, and survival of the spores. This model simulates transport of spores to a delivery location of the alveolated region of the lungs. The second stage model coupled to the first, the subject of further work, models the transport of spores through the respiratory system, transport of spores into the alveoli, and the pathogenesis of inhaled B. anthracis. 2.2. Markov Chain Model. A Markov chain model is a stochastic means of modeling physical systems. A schematic of the Markov chain can be seen in Figure 1, where the Greek letter λ signifies the loss rate from the associated states. Loss rates are mathematical functions that describe the physics of spores being removed from the states. The Markov chain is constructed by

ARTICLE

Table 1. Three-Region Model of the Respiratory System9 included anatomical

region

region

region

structures

number

symbol

name

1

R1

nasopharynx

2

R2

trachiobronchial

3

R3

pulmonary

nares, nasal cavity, oral cavity nasopharynx, larygopharynx, larynx trachea, bronchi, first branches of bronchioles respiratory bronchioles, acinus, alveolar sacs

filling the Markov matrix with transition probabilities that estimate the probability of a spore transitioning from one physical compartment (state) to another, i.e., depositing from air to a surface. This Markov chain model consists of 8 states in order, starting with state 1: bulk fluid in R1, respiratory system surface of R1, bulk fluid in R2, respiratory system surface of R2, bulk fluid of R3, respiratory system surface of R3, inactivation of B. anthracis due to its natural survival rate in air, and exhalation from the respiratory system. As can be seen in figure 1, States 1, 3, and 5 represent the bulk fluid in each of the regions, and States 2, 4, and 6 represent the surface of the respiratory system. When particles (or spores) deposit on the respiratory system surface they are prevented from resuspension by the mucociliary escalator, and eventually expectorated within 12 h post deposition.7 Thus States 2, 4, and 6 can be considered absorbing states, therefore, only the natural inactivation rate of spores in air is the loss mechanism from these states. The transition probabilities that are the functioning portion of the Markov chain model are derived from first order loss rates associated with the individual states. Loss rates are first order since inactivation of spores is first order. Let us denote the probability of being in state i and remaining in state i during a single time step (Δt) as pii. Given the sum total of losses from state i (λi), the exponential survival probability informs pii as shown in eq 3.10,11 pii ¼ expð
λi 3 ΔtÞ

ð3Þ

In turn an unconditional probability can be used to estimate the probability of being in state i and transitioning to state j (pij) during Δt. This is shown in eq 4, where the compliment of pii is multiplied by the ratio of the specific loss rate associate with transitioning from state i to state j (λj) and λi.11 pij ¼

λij ½1
pij  λi 3

ð4Þ

When the host inhales spore-contaminated air, the spores are transported into the respiratory system with a volumetric flow rate Q. The air is then exhaled with a volumetric flow rate B which with an assumption of constant volume complete mixing, is the same as Q12 using a value of 125 cm3/min. While inhalation and exhalation flow rates are very variable based on gender, health condition, age, physical exertion, and a number of other factors, this value of 125 cm3/min was used by Heyder et al.13 for their lung cast experiments, values of which are used in the model. The loss rates associated with being inhaled deeper into the respiratory system can be generalized in eq 5, where VRx is the volume of the current region (i.e., if beginning in R1 then VRx is the volume of R1, and in turn λxy is the loss of a spore in the higher region x transitioning to region y deeper 5829

dx.doi.org/10.1021/es200901e |Environ. Sci. Technol. 2011, 45, 5828–5833

Environmental Science & Technology

ARTICLE

Figure 1. Computational schematic depicting the structure of the Markov chain model used for the first stage model.

in the respiratory system). λxy ¼

Q þB V 9Rx

ð5Þ

Then spores in a lower region of the respiratory system may be lost to a higher region of the respiratory system, via exhalation. This loss from the lower regions of the respiratory system to higher ones can be generalized in eq 6. λyx ¼

B V 9Ry

ð6Þ

In addition to spores being lost due to bulk fluid transport from one region to another, deposition can occur. For deposition sedimentation, impaction, and diffusion are accounted for as the main loss mechanisms. Sedimentation is considered, since although counterintuitive to most, the majority of the bronchioles and terminal bronchioles are perpendicular to the gravity vector when in a sitting or standing position. The sedimentation rate is based on the terminal settling velocity (vts) shown in eq 7, which is an approximation based on particle diameter (dp) and is accurate for particles up to 50 μm in diameter in air. The quotient of vts and mean diameter of the region in which the spore is currently suspended estimates the sedimentation rate for the spore in that region. " # 0:166 2 ð7Þ vts ¼ 0:0018 3 dp 3 1 þ dp By taking the sum of the sedimentation rate and rates of deposition due to impaction and diffusion for the region in which the spore is currently suspended (DIRx) obtained from Heyder et al.,13 the loss of spores from suspension in bulk fluid in region x can be estimated (eq 8). vts þ DIRx ð8Þ λdeposition ¼ dRx The physiological parameters for human, guinea pig, and rhesus macaques, DI values at a dp of 1 μm, and the functional forms and estimated values for all loss rates in the Markov chain model can be seen in Supporting Information. Inactivation of spores is assumed constant and a value from Sinclair et al.14 for ambient air is used. 2.3. Extrapolation to Additional Host Species. Due to the desired framework of the model, changes such as anatomical dimensions are straightforward, so long as the anatomical structure

Figure 2. Results from Markov chain model after 60 min of human respiration were simulated.

is similar (i.e., not moving from mammalian to amphibian lungs). The differences between humans and animals can be compared using host specific volumes and mean diameters of the three regions for humans,12 guinea pigs,13 and rhesus macaques,14 values of which are summarized in the Supporting Information. Because the same dynamics affect transport and deposition of spores; the loss rate functional forms remain the same, only the parameter values for the loss rates are altered. 2.4. Code and Model Verification. The stochastic system developed for this model simulates the mass balance of spores transported through the respiratory systems. A 2-fold process of verification was performed to ensure that the mass balance and code embodying it operated properly. Code verification was enacted by having a researcher familiar with coding in MATLAB checking for errors that may have been missed by the original author. The code is then executed line-by-line, or loop-by-loop, monitoring for the expected or required output (i.e., vectors are filled with correct information). The code passed this verification step. Verification of the mass balance is a more involved process, which is divided into two main steps. First the states are decoupled from the entire Markov chain and a controlled number, 100 spores, is introduced to the individual states. Then each state’s output is monitored to account for all of the 100 spores, and monitored to ensure that only the loss rates affecting that state are causing a loss of spores from the state. For example, 5830

dx.doi.org/10.1021/es200901e |Environ. Sci. Technol. 2011, 45, 5828–5833

Environmental Science & Technology

ARTICLE

Table 3. Parameter Alterations for the Sensitivity Analysis of Physiological Parameters parameter alteration control

2.3 (10 ) 1.0 (10
8)

1.7 (10
3) 5.05 (10
8)

2 (Q)

2.01 (10-8)

5.05 (10-8)

-4

2.98 (10 )

host animal

[confidence interval]

[confidence interval] a

0.556; [0.0644, 1.323]

guinea pig

147,411; [17,400, 357,700]

rhesus macaque

92,000; [29,400, 154,560] a

2.392; [0.764, 4.019]

humans

2,500
55,000 b

2.085
45.87

a

Developed in ref 2. b Confidence Interval not possible as range only given in ref 18 that the LD50dd is based on.

if spores are in R1 the only losses are deposition to R1 surfaces, bulk fluid transport to R2, exhalation, and inactivation of the spores. This ensures that spores are being neither created nor destroyed in the states and are accounted by the loss mechanisms. Then the states are recoupled into the full Markov chain and the same process of introducing a controlled number of spores, 20,000 in this case, and the losses for the system are monitored for again. This verification step also passed.

3. RESULTS The model is executed so as to simulate 60 min of breathing once the host is exposed to air contaminated with a known number of spores. As can be seen in Figure 2 after an exposed dose of 100,000 spores, there is a steady increase in the number of spores in the third region, which is considered the delivery location for this model, the alveolated region of the lungs (where the air sacs are located in the terminal bronchioles). Entry of the spores into the alveoli is modeled with the second stage model which is out of the scope of this manuscript as this second stage is a separate modeling effort and focuses on B. anthracis pathogenesis. From the first stage model, delivered doses can be computed using dose response data for 1 μm spores 17 shown in Supporting Information. As can be seen in Figure 3 a linear relationship between exposed and delivered dose is exhibited. This allows for a constant correction factor for humans of 0.000834 for humans (referred to as ηddhm). In other words, only 0.0823% of inhaled spores in the exposed dose will reach the alveolated region of the lungs after 60 min of respiration. The same linear relationship between exposed and delivered dose

1.9 (10
5)


6

2.5 (10
5)


5

2.0 (10
4)

5.3 (10 )

2 (Q) and 2(V) and 2 (d)

LD50dd (spores)

1.3 (10
4)


6

1.4 (10 )

2 (Q) and 2(V)

1.7 (10-3)


5

1.3 (10 )

0.5 (d)

LD50 (spores)

2.6 (10
5)


4

2 (V) 0.5 (V)

2 (d)

Table 2. Alteration in LD50 with and without Considering Delivered Dose Correction for Guinea Pig and Rhesus Macaques

rhesus macaque ηdd


6

3.7 (10 )

0.5 (Q)

Figure 3. Linear dose correction factor for humans, guinea pigs, and rhesus macaques showed the same relationship as can be seen in Supporting Information.

guinea pig ηdd

1.5 (10 )

holds for guinea pigs and rhesus macaques (plots in Supporting Information). Constant correction factors of ηddgp of 3.77(10
6) for guinea pigs and ηddrm of 2.6(10
5) for rhesus macaques were observed as well. This difference in delivered dose correction allows for an extrapolation of the median lethal dose (LD50) to a median lethal delivered dose (LD50dd), which can be seen in Table 2 for humans, guinea pigs, and rhesus macaques (LD50 from ref 18 for humans, and 2 for guinea pigs and rhesus macaques). This lethal dose correction helps to explain why guinea pigs that have a greater clinical susceptibility to B. anthracis compared to other hosts that have a greater LD50. An interesting result, as can be seen in Table 2, is that the LD50 for guinea pigs and rhesus macaques based on delivered dose are now estimated to be approximately 1 and 4 spores, respectively, and reduced to between 2 and 46 spores for humans. This is a large decrease from the unaltered LD50 which did not account for an effective delivered dose to the susceptible region of the respiratory system. One of the reasons for such a dramatic decrease in the LD50 for these host species is the fact that only removal mechanisms are considered for this first stage model. This highlights the great efficiency with which respiratory systems protect the health of potential hosts. The coupled model accounts for spore germination and bacterial growth in the respiratory system, therefore it is highly likely that the LD50 specific to a pathogen burden will not be such a great decrease from the exposed dose LD50. The differences in correction factor between host species can likely be attributed to the much higher respiration rate as well as the regional volume difference. Human adults have a respiration rate of 13 breaths min
1 which is less than half that of the rhesus macaques at 29
30 breaths min
1, with guinea pigs even greater at an average of 60 breaths min
1. A sensitivity analysis can be seen in Table 3. The methods summarized in Bartrand et al.2 have been charted and shown in Supporting Information showing that changes in volumetric flow rate have the greatest effect on resulting correction factor, closely followed by overall lung volume. It is interesting to note that while the dose correction factor for both guinea pigs and rhesus macaques was changed most significantly by both volumetric flow rate and total lung volume, the dose correction factor for rhesus macaques showed a greater sensitivity overall to lung volume, as opposed to guinea pigs where volumetric flow rate was the overall driver of sensitivity. This constant correction factor can be directly included into the current dose response models, resulting in a dose
response relationship using delivered dose as a metric. Adapting the exponential and beta Poisson dose response models into the 5831

dx.doi.org/10.1021/es200901e |Environ. Sci. Technol. 2011, 45, 5828–5833

Environmental Science & Technology

ARTICLE

forms shown in eqs 9 and 10, respectively, allow for delivered dose specific parameters. Thus for the exponential model the parameter k has been altered to a delivered dose specific parameter (k*), and in turn R to R* and N50 to N50* for the beta Poisson model. 

Exponential : pðResponseÞ ¼ 1
expð
k 3 dose 3 ηdd Þ

ð9Þ

Beta Poisson : pðResponseÞ "



dose 3 ηdd 3 21=a
1 ¼ 1
1þ  N50

#
a ð10Þ

The fit of the models did not change due to the correction being a constant factor affecting the dose and given the structure of the data, only the grouped parameters ηdd 3 k* (for the exponential) or ηdd/N*50 (for the beta Poisson) can be estimated. Therefore ηdd cannot be independently optimized using animal model data. This is due to the delivered dose correction not affecting the observed probability of response (eq 11) to which the models must be optimized.1 Where in eq 11 πoi is the observed probability of response, pi is the observed positive responses and ni is the total number of hosts in the dosing group. πoi ¼

pi ni

ð11Þ

This model for delivered dose is then coupled with the second stage, where the transport into the alveoli and pathogenesis of the spores is modeled. The combination of the two models provides an understanding of the in vivo body burden as can be estimated from the exposed dose. However this first stage can be separated from the coupled model so as to assess the effectiveness of this stage, since the data requirements are less intensive than the second stage. By coupling the two stages this develops a means of describing the main infection process for inhaled B. anthracis spores. The modeling framework developed allows for not only the separation of the two stages but also the ease of adaptation of the model. Since it is established as a stochastic model, adaptation to other hosts species as demonstrated is straightforward, but also the general framework can easily be adapted to additional exposure routes beyond inhalation. The correction factors developed also allow for a great understanding of how physiology of the hosts affects the dose that is delivered to the susceptible location of the lungs. This leads to an understanding of host heterogeneity for the dose response relationship based on host physiology thereby also developing an initial framework for PDRM. Also this work leads us to a greater understanding of the mitigation and building rehabilitation goals given a greater comprehension of the efficiency of protection that the respiratory system presents and its affect on the dose response of inhaled B. anthracis spores.

’ ASSOCIATED CONTENT

bS

Supporting Information. This material is available free of charge via the Internet at http://pubs.acs.org.

’ AUTHOR INFORMATION Corresponding Author

*E-mail: [email protected].

’ ACKNOWLEDGMENT This research was performed under the appointment to a dissertation award program sponsored by the Department of Homeland Security (DHS), administered by the Oak Ridge Institute for Science and Education (ORISE) through an interagency agreement between the U.S. Department of Energy (DOE) and DHS. ORISE is managed by Oak Ridge Associated Universities (ORAU) under DOE contract DE-Ac05-06OR23100. All opinions expressed in this paper are the authors and do not necessarily reflect the policies and views of DHS, DOE, and ORAU/ORISE. Thanks to the Center for Advancing Microbial Risk Assessment (CAMRA) jointly funded through the DHS and Environmental Protection Agency under the Science to Achieve Results (STAR) grant program grant number R83236201, Koerner family for generous support through the Koerner Family Fellowship, and the Department of Education through Graduate Assistantships in Areas of National Need (GAANN). ’ REFERENCES (1) Haas, C. N.; Rose, J. B.; Gerba, C. P. Quantitative Microbial Risk Assessment; John Wiley and Sons Co.: New York, 1999. (2) Bartrand, T. A.; Weir, M. H.; Haas, C. N. Dose-Response Models for Inhalation of Bacillus anthracis Spores: Interspecies Comparisons. Risk Anal. 2008, 28 (4), 1115–1124. (3) Meynell, G. G.; Maw, J. Evidence for a two-stage model of microbial infection. J. Hygiene 1968, 66 (2), 273–280. (4) Titball, R. W.; Manchee, R. J. Factors Affecting the Germination of Spores of Bacillus anthracis. J. Appl. Bacteriol. 1987, 62, 269–273. (5) Guidi-Rontani, C. The Alveolar Macrophage: The Trojan Horse of Bacillus anthracis. Trends Microbiol. 2002, 10 (9), 405–409. (6) Haber, S.; Yitzhak, D.; Tsuda, A. Gravitational Deposition in a Rhythmically Expanding and Contracting Alveolus. J. Appl. Physiol. 2003, 95, 657–671. (7) Koblinger, L.; Hofmann, W. Analysis of Human Lung Morphometric Data for Stochastic Deposition Calculations. Phys. Med. Biol. 1985, 30 (6), 541–556. (8) Yeh, H.-C.; Schum, G. M. Models of Human Lung Airways and Their Applications to Inhaled Particle Deposition. Bull. Math. Biol. 1980, 42, 461–480. (9) Morrow, P. E.; Bates, D. V.; Hatch, T. F.; Mercer, T. T. International Commission of Radiological Protection Task Group on Lung Dynamics, Deposition and Retention Models for Internal Dosimetry of the Human Respiratory Tract. Health Phys. 1966, 12 (10) Nicas, M.; Sun, G. An Integrated Model of Infection Risk in a Health Care Environment. Risk Anal. 2006, 26 (4), 1085–1096. (11) Ross, S. M. Introduction to Probability Models, 9th ed.; Academic Press: Burlington, MA, 2007. (12) Weibel, E. R. Morphometry of the Human Lung; Academic Press Inc.: New York, 1963. (13) Heyder, J.; Gebhart, J.; Rudolf, G.; Schiller, C. F.; Stahlhofen, W. Deposition of Particles in the Human Respiratory Tract in the Size Range 0.005
15 μm. J. Aerosol Sci. 1986, 17 (5), 81–85. (14) Sinclair, R.; Boone, S. A.; Greenberg, D.; Keim, P.; Gerba, C. P. Persistence of Category A Select Agents in the Environment. Appl. Environ. Microbiol. 2008, 74 (3), 555–563. (15) Schreider, J. P.; Hutchens, J. O. Morphology of the Guinea Pig Respiratory Tract. Anat. Rec. 2005, 196, 313–321. (16) Martonen, T. B.; Katz, I. M.; Musante, C. J. A Nonhuman Primate Aerosol Deposition Model for Toxicological and Pharmaceutical Studies. Inhalat. Toxicol. 2001, 13, 307–324. (17) Druett, H. A.; Henderson, D. W.; Packman, L.; Peacock, S. Studies on Respiratory Infection I: The Influence of Particle Size on 5832

dx.doi.org/10.1021/es200901e |Environ. Sci. Technol. 2011, 45, 5828–5833

Environmental Science & Technology

ARTICLE

Respiratory Infection with Anthrax Spores. J. Hygiene 1953, 51, 359–371. (18) Inglesby, T. V.; O’Toole, T.; Henderson, D. A.; Bartlett, J. G.; Ascher, M. S.; Eitzen, E.; Friedlander, A. M.; Hauer, J.; McDade, J.; Osterholm, M. T.; Parker, G.; Perl, T. M.; Russel, P. K. Anthrax as a Biological Weapon, 2002: Updated Recommendation for Management. J. Am. Med. Assoc. 2002, 287, 2236–2252.

5833

dx.doi.org/10.1021/es200901e |Environ. Sci. Technol. 2011, 45, 5828–5833