In silico modeling of spore inhalation reveals fungal persistence following low dose exposure

The human lung is constantly exposed to spores of the environmental mould Aspergillus fumigatus, a major opportunistic pathogen. The spectrum of resultant disease is the outcome of complex host-pathogen interactions, an integrated, quantitative understanding of which lies beyond the ethical and technical reach permitted by animal studies. Here we construct a mathematical model of spore inhalation and clearance by concerted actions of macrophages and neutrophils, and use it to derive a mechanistic understanding of pathogen clearance by the healthy, immunocompetent host. In particular, we investigated the impact of inoculum size upon outcomes of single-dose fungal exposure by simulated titrations of inoculation dose, from 106 to 102 spores. Simulated low-dose (102) spore exposure, an everyday occurrence for humans, revealed a counter-intuitive prediction of fungal persistence (>3 days). The model predictions were reflected in the short-term dynamics of experimental murine exposure to fungal spores, thereby highlighting the potential of mathematical modelling for studying relevant behaviours in experimental models of fungal disease. Our model suggests that infectious outcomes can be highly dependent upon short-term dynamics of fungal exposure, which may govern occurrence of cyclic or persistent subclinical fungal colonisation of the lung following low dose spore inhalation in non-neutropenic hosts.


Direct derivation of the parameters
The first 8 parameters (M, N v , δ N , δ C , δ F , k ND , k NF , d MF and β) listed in Table 1 were directly derived from the experimental data as described below. We use 2 significant digits by rounding the third digit.
Number of macrophages M. Hohl et al. [3] measured the absolute number of macrophages in murine bronchoalveolar lavage fluid (BAL), at 24 hours post-inoculation with an intranasal inoculum of 10 7 live A. fumigatus conidia. In our model we assume that the number of macrophages remains constant and we use the observed value of M = 0.3 × 10 6 cells ( Fig. 1 in [3]).
Number of available neutrophils N v . Assuming total murine blood volume to be 58.5ml/kg, the circulating volume of blood in a mouse with a typical weight of 25 g can be estimated as 58.5 × 25 × 10 −3 = 1.4625 ml. If the total number of white blood cells average, in mice 7.45 × 10 3 cells/µl ( Fig. 3 in [4]), of which 14.63% are neutrophils [4,5], the circulating neutrophil titre can be calculated as 7.45 × 10 3 × 0.1463 × 1.4625 × 10 3 = 1.59 × 10 6 cells. Since circulating neutrophil count is estimated as accounting for 1 − 2% of the total body neutrophil count [5], we obtain N v = 150 × 10 6 cells as a number of the total neutrophils available.
N-mediated F killing rate k NF . In vitro incubation of A. fumigatus spores (10 5 ml −1 ) with an equal number of neutrophils lead to killing of 30.6% of the spores after 2-3 hours (Table1 in [6]). Assuming the rate of change of the number of spores F in this setting is described by  [7]). Assuming the rate of change of fungal burden F in this setting is described by Ex vivo killing assay found 47 % of A. fumigatus swollen conidia were killed after 6 hours of incubation with 0.5 × 10 6 murine alveolar macrophages ( Fig. 4 in [8]). In this study alveolar macrophages were harvested from mouse lungs and seeded at a conidium/AM ratio of 1:10. We therefore derived We therefore chose the nominal value of d MF = 0.32 × 10 −6 cells −1 h −1 by taking the mean of these two values.
F proliferation rate β. Vallor et al. [9] measured the temporal change of A. fumigatus fungal burden in a neutropenic guinea pig model of invasive pulmonary aspergillosis. We assume the fungal proliferation rate is approximate to that in a mouse. From their data (Table 1) F(1h) = 10 5.9 conidial equivalents (CEs) and F(24h) = 10 7.7 CEs.
The model equation Note also that the in vitro growth rates of 0.2349 − 0.3375 h −1 measured by a microbroth method for nine A. fumigatus isolates ( Table 1 in [10]).
N-induced D efflux rate k ND . The numbers of DCs in the lung (D) and the lymph node (D lymph ) in immunocompetent mice challenged with A. fumigatus were measured over 24 hours ( Fig. 1B in [11]): D lymph (0h) = 2.5×10 4 cells and D lymph (24h) = 5×10 4 cells. Assuming that the increase in the lymph node DCs is due to migration of DCs from the lungs to the lymph nodes and that N(t) ≈ 1.5 × 10 6 cells and D(t) ≈ 0.1 × 10 6 cells for the first 24 hours post-inoculation ( Fig. 1 in [2]) for simplification, N degradation rate δ N . Basu et al. [12] measured BrdU-labeled neutrophils in murine blood at 96 and 120 hours after BrdU injection, (Fig. 5C) and calculated δ N = 0.0609 ≈ 6.1 × 10 −2 h −1 from the rate of loss of labeled neutrophils, using the formula This is the value we use in our model. The same method with prior published data of human and mouse neutrophil decay deduces similar values of 0.033 h −1 (Fig. 1 in [13]) and 0.027 h −1 (Fig. 2 in [14]). Note also that neutrophils have a similar lifetime in the circulation as in the tissue (≈ 0.015 h −1 by Fig. 2b in [15]).
C degradation rate δ C . Mehrad et al. [16] measured lung TNF-α levels in cyclophosphamide-treated (immunocompromised) mice after A. fumigatus challenge. We assume that the decrease of the lung TNF-α level after its peak at 48 hours post-inoculation is mainly due to natural degradation and derived δ C = 0.0655 ≈ 6.6 × 10 −2 h −1 using their data C(48h) = 5.8 ng/ml and C(96h) = 0.25 ng/ml The assumption of δ C C(t) >> k CD D(t), k C MF(t) in the model equation (2) was confirmed to be valid in our model simulation.

Parameter optimization
The remaining 5 model parameters in our model (α, α D [D v ], k CD , k C and δ D ) were estimated by parameter optimization using Monte Carlo Simulated Annealing [17] to fit to the in vivo experimental data, while the directly evaluated 8 parameters in the previous section were fixed. The in vivo experimental data used for the parameter optimization (Table S1) includes the dynamics of F (CFU counting) as measured in our own murine fungal challenge experiments (see Methods in the main text and Fig. 6 (a)), C (in BAL) measured at different times after a challenge by 6 × 10 6.25 CFU of A. fumigatus infection ( Fig. 1 in [1]), and N and D measured post challenge by killed Aspergillus hyphae (5 − 9 × 10 5 CFU) in mice ( Fig. 1A and 4B in [2]). Note that high inocula (order of 10 6 CEs) were used in these published experiments. We optimized the parameter values from 1000 randomly chosen initial values within the range of [0, 1] for each parameter and selected the set with the minimal error, S = ∑ 13 i=1 (R sim i − R exp i ) 2 between the experimentally measured value R exp i and its simulated value R sim i , as our nominal parameters.

Parameter sensitivity analysis
We performed sensitivity analysis to check whether our main results are robust against the perturbation in the nominal parameter values.

Global sensitivity analysis
We evaluated the global parameter sensitivity of our model, in order to identify the sensitivity of the time for N C (and the maximum fungal burden) to be reached (T c F 0 ), and the time required for clearance (to below one CE) of the initial fungal burden F 0 (τ F 0 ), to simultaneous variations in the parameters. For this purpose, we performed simulations of our model for 100, 000 sets of the parameters, where the value of each of the 13 parameters was chosen from a uniform distribution of up to 10% variation around its nominal value, and derived Sobol's global sensitivity metric [18] for both a high dose (F 0 = 10 6 CE) and a low dose (F 0 = 10 2 CE) inoculation. The parameter sets for which the model does not achieve fungal clearance within 150 hours were eliminated from the analysis (≈ 4% of the 100, 000 parameter sets investigated). Sobol's global sensitivity metrics for both F 0 = 10 6 and F 0 = 10 2 indicate that the results are most sensitive to β (the fungal growth rate), followed by M, α, Nv, k C , d MF and δ C (Fig. S1). We confirmed the convergence of the sensitivity indices with simulations of another 100000 sets of parameters. For the low dose inoculation, β has a higher relative sensitivity than for the high dose inoculation, since the earlytime dynamics is dominated by the exponential growth of the fungus (the β[F] term in our equations) that becomes exaggerated as the initial inoculum is decreased. Indeed our results demonstrate that a slight change in β for a low dose inoculation leads to qualitatively different disease outcomes (Fig. 8). For the high dose inoculation, M has a higher relative sensitivity than for the low dose infection, since high dose spores have to be killed quickly to avoid its saturation.

Persistence of low dose fungal exposure
Our model predicted (Fig 5(a)) and the experiments verified (Fig 6(b), (c)) that longer time is required to achieve N c for lower inoculum. To confirm whether this model prediction robustly holds despite parameter variations, we evaluated T c 10 2 and T c 10 6 , where T c F 0 is the time to achieve N c with the initial fungal burden F 0 , for increasing levels of variations up to 0.1, 0.2, 0.5 and 1 orders of magnitude from the nominal values for all the parameters. The results confirmed that T c 10 2 > T c 10 6 robustly holds ( Fig. S2(a)). The same analysis confirmed that τ 10 2 > τ 10 6 robustly holds despite parameter variations, when the curative neutrophil threshold N C > 0 ( Fig. S2 (b)). Note that τ 10 2 < τ 10 6 are observed if N C < 0, since the fungal burden does not grow initially (when N = 0) under this condition and thus is rapidly cleared mainly by macrophages (M), corresponding to the behaviour A in Fig. 8(b)-(c).

Existence of oscillatory behaviors
The clearance time against varied values of β for different initial fungal burdens (Fig. 8) indicates a qualitative change atF 0 = 10 4 , above which the region B disappears, suggesting that the oscillatory behaviors are observed only for initial fungal burden F 0 below this threshold burdenF 0 . Our extensive set of simulations and mathematical analysis shows that the existence of the oscillatory behaviors in low dose infection is always observed for perturbation of all the parameters for at least up to an order of magnitude from the nominal values, except for the case with N C < 0 corresponding to the behaviour A in Fig. 8(b)-(c).

Experimental data
The experimental data that we discuss in the main text is presented in Tables S1 and S2.  Figure S2: Simulation results for uniformly sampled parameter sets, each of which deviates up to 0.1 (black), 0.2 (red), 0.5 (purple/light pink) or 1 (grey/light grey) order of magnitude from its nominal value (blue circle). (a) T c 10 2 against T c 10 6 for the cases N C is reached. All results lie above the blue line suggesting that T c 10 2 > T c 10 6 . (b) τ 10 2 against τ 10 6 . Light pink and light grey points represent parameter sets where N C < 0. For N C > 0, when the neutrophils are required for clearance, all results lie above the blue line meaning that τ 10 2 > τ 10 6 .