Our proposed mathematical model of Sp colonization (Figures 1A,B) is a system-level representation of the prominent interactions between Sp, the airway epithelium, and immune cells (a–j) in Figures 1A,B, that were identified based on the empirical evidence from numerous experimental in vivo and in vitro studies as detailed below.

Under homeostatic conditions, a population of commensal bacteria, Sp, resides in the lumen on the apical side of the airway epithelium, where they are contained by a competent epithelial barrier integrity (Beisswenger et al., 2007) (Figures 1Aa,Ba) and immune responses mediated by neutrophils and macrophages (Dick et al., 2008; Standish and Weiser, 2009) (Figures 1Ab,Bb). Through disrupted barrier, apically located Sp can translocate to reach the blood vessel (Beisswenger et al., 2007) (Figures 1Ac,Bc), where they are either killed by resident immune cells that circulate in the blood (Li et al., 2002; Li, 2004), or grow uncontrollably and result in invasive infection if the immune cells cannot contain the translocated Sp (Silverstein and Rabadan, 2012). The amount of the translocated Sp in the blood vessel is therefore a determinant of whether disrupted Sp-host interactions cause serious infection such as sepsis.

Translocation of Sp occurs through the airway epithelial barrier, whose integrity is regulated by the apically located bacteria load. The bacteria bind to Pattern-Recognition immune receptors, specifically Toll-like receptors (TLR2s) that are preferentially expressed on the apical side of the airway epithelial cells (Melkamu et al., 2009), and activate the TLR signaling cascade (Figures 1Ad,Bd). The activation of the TLR cascade in epithelial cells decrease the barrier integrity of the airway epithelium (Figures 1Ae,Be) by TLR-mediated activation of proteases that damage the epithelial cells (Oggioni et al., 2004; Schmeck et al., 2004; Attali et al., 2008; Tieu et al., 2009) and by reduction of the barrier recovery rate due to the increased expression of the transcriptional repressor SNAIL1, which inhibits the expression of claudin, a component of the tight junctions (Clarke et al., 2011).

Active TLR signaling also induces recruitment of neutrophils from the neutrophil pool in the blood vessel, via the release of IL-17 (Zhang et al., 2009) that activates neutrophil-attracting interleukins IL-8 (Lindén, 2001) (Figures 1Af,Bf). The recruited neutrophils trigger transmigration of macrophages to the site of infection (Zhang et al., 2009), further potentiating the immune responses to the apically located pathogens (Figures 1Ag,Bg), whereas macrophages on the lumen restrict neutrophil transmigration (Zhang et al., 2009) by releasing neutrophil-repellent anti-inflammatory cytokines (Knapp et al., 2003) (Figures 1Ah,Bh). Transmigrating neutrophils release barrier degrading proteases (Chin et al., 2008) to reduce the barrier integrity (Nash et al., 1987; Nusrat et al., 1997; Zemans et al., 2009) (Figures 1Ai,Bi). The reduced barrier integrity in turn allows more transmigration of both neutrophils and macrophages from the blood vessel to the site of infection (Nash et al., 1987) (Figures 1Aj,Bj).

The model elucidates the main control structure of the system that maintains homeostatic interactions between commensal bacteria, Sp, and the host, via a dynamic interplay between environmental stressors, epithelial barrier integrity and immune responses (Figure 1B). At the apical side of the airway epithelium, Sp load is regulated via activation of TLRs, which induce immune responses that decrease the bacterial load but also reduce the epithelial barrier integrity. While the reduced barrier integrity enables transmigration of immune cells from the blood vessel for effective killing of Sp at the apical side of the epithelium barrier, it also allows transmigration of Sp from the apical side of the epithelium to the blood vessel, potentially causing systemic infection (sepsis). The dynamic interplay between the immune responses and the epithelial barrier integrity are further modulated by their mutual inhibition.

Our model further assumes two switches, an R-switch for TLR activation and an S_v-switch for the growth of the transmigrated bacteria in the blood vessel, based on the experimental evidence described below. The R-switch for Sp-mediated activation of TLRs reflects the observations that low concentrations of Sp do not cause activation of TLR signaling, while high concentrations lead to a sharp increase in TLR activity with hysteresis (He et al., 2009; Shalek et al., 2013; Sung et al., 2014). We model the R-switch by a perfect switch, which is a phenomenological representation of the bistable switch (Sung et al., 2014), and is described by the off- and on-states (R = R_off and R_on) with the activation (S⁺) and inactivation (S⁻) thresholds for the critical concentrations of apically located Sp that abruptly and sharply turn on-or-off TLR activity (Figure 1C and Equation 2). The S_v-switch reflects the observations that transmigrated bacteria in the blood vessel (S_v) either overgrow (Benton et al., 1997) or are contained by resident immune cells depending on the bacterial concentration (Supplementary Figure 3B). We model the S_v-switch with a switching threshold of Sv*, above which the infiltrated bacteria in the blood vessel grow exponentially (Figure 1D).

The resulting model is described by a hybrid system of five ODEs (Equation 1). The nominal values of the 24 model parameters (Table 1) were derived by fitting the model outcome to datasets from 11 independent studies, namely three in vivo studies (Benton et al., 1997; Zhang et al., 2009 and our own experiment) and eight in vitro studies (Nash et al., 1987; Coyne et al., 2002; Lagrou et al., 2003; Attali et al., 2008; Chin et al., 2008; Komori et al., 2011; Hathaway et al., 2012; Kwok et al., 2012), as detailed in the Supplementary Material. Our model therefore provides a coherent mathematical framework to explain both in vivo and in vitro data.

One of the dataset used for the parameter estimation was obtained from in vivo studies in Zhang et al. (2009), where the mice were challenged with 10⁷ CFU of Sp and recovered their healthy state, which is characterized by nonzero apical commensal bacterial load that does not trigger host responses. Our model was fitted to reproduce the experimental measurement in Zhang et al. (2009) for the apical bacterial load (S_a) and the concentrations of neutrophils (N) and macrophages (M) (Figure 2A). Both the experimental data and our model simulation demonstrate that the transient Sp challenge (increase of S_a) triggers a transient increase in N and a subsequent increase in M. These immune responses can bring S_a down to a homeostatic level, when the saturation limit for S_a is not high enough, which enabled the recovery of the mice from the bacterial challenge within 7 days without demonstrating invasive infection.

The simulation of our model with the data-calibrated nominal parameters further predicts the dynamics of three variables that were not measured in this experiment, TLR activity (R), barrier integrity (B) and infiltrated Sp in the blood vessel (S_v), thereby explaining the underlying mechanism of the healthy recovery from a bacterial challenge. Upon a Sp challenge, the apically located bacterial load becomes high enough (Sa(0)>S+) to activate TLRs (R = R_on), which trigger recruitment of immune cells to the site of infection. These immune responses bring the initially high S_a down to below S⁻, where the R-switch turns off (R = R_off) and stays off as S_a remains below S⁺, as suggested by the focal point analysis (see Methods). The epithelial barrier integrity (B) continuously decreases while the R-switch is on (R = R_on), allowing bacteria to invade the blood vessel, as demonstrated by a rise in S_v. However, a healthy clearance of S_v is achieved without causing sepsis, since the peak of S_v remains below the threshold, Sv*, of the S_v-switch. Note that both the experiments and our model simulation demonstrate that M stays high while B is kept high after N goes to zero, suggesting the importance of M as an immune mediator that does not compromise the barrier integrity.

The healthy recovery behavior described above is characterized in our model by convergence to the off-state of TLR activity (R = R_off) accompanied by the containment of S_v (Sv<Sv*), and is robustly observed under perturbations to the parameter values. Among 10,000 simulations conducted by randomly sampling parameter values from an uniform distribution over two orders of magnitude around the nominal values, 83% of the simulations demonstrated a healthy recovery from a transient Sp challenge (Figures 3A,B). In 98% of the healthy recovery cases computationally observed, S_v reached its peak while R = R_on (Figure 3E), suggesting that the appropriate host responses via TLR activation are responsible for containing S_v. The appropriate level of the barrier damage by active TLRs enables effective recruitment of immune cells that can reduce S_a, but prevents excessive transmigration of S_a to S_v, keeping the S_v lower than the threshold, Sv*. The robust appearance of the healthy recovery in our model simulations confirms that that our model can coherently explain the mechanism behind the healthy recovery of the host from pneumococcal colonization, which can be effectively cleared by the natural host responses without any treatments.

The remaining 17% of the simulations with parameters perturbed from their nominal values demonstrated systems dynamics that correspond to serious infection or inflammation. They are classified into three disease phenotypes, depending on the states of the S_v- and R-switches (Figure 3A). The state of the S_v-switch determines whether sepsis occurs due to invasive infection of Sv(>Sv*), and that of the R-switch determines whether immunological scarring occurs due to persistent host responses caused by R = R_on. Immunological scarring refers to the cumulative and long-term effects of immune response to pathogens, including tissue remodeling and altered immune responses to new pathogenic challenges, that persist after the pathogenic organism has been cleared (Fonseca et al., 2015). In our simulations, 13% demonstrated sepsis without immunological scarring (Sv>Sv* and R = R_off) and the other two disease phenotypes, immunological scarring (R = R_on) with and without sepsis, were observed 2% each (Figure 3A).

Immunological scarring is characterized by a persistent on-state of the R-switch due to S_a staying above S⁻ (Figure 3C). The R-switch triggers persistent host responses leading to sustained immune responses which are not strong enough to decrease S_a below S⁻ but cause persistent barrier damage. Note that the peak of N is much lower in the immunological scarring phenotype than that in the healthy recovery and sepsis phenotypes with R = R_off (Figures 3B,D), resulting in weak immune responses that are not sufficient to decrease S_a. As a result, the host becomes vulnerable to a second bacterial attack due to the damaged barrier and the sub-threshold concentration of S_v (<Sv*), which stay as silent remainders of the first pathogenic challenge.

Sepsis is characterized by outgrowth of S_v once it surpasses the threshold Sv* (Figure 3C). In 99.7% of the sepsis phenotypes simulated by our model, the onset of sepsis occurs (when Sv=Sv* is achieved) while R is on (Figure 3G), suggesting that whether sepsis occurs or not is determined by the dynamics of S_v while R is on. It is similar with the healthy recovery case, where S_v reaches its peak below the threshold Sv*, while R is on. Moreover, the duration of R = R_on is much longer for the sepsis phenotype compared to the healthy recovery phenotype (Figure 3H), suggesting that persistent host response may allow excessive transmigration of Sp into the blood vessel above Sv*.

When both the S_v- and R-switches are on, sepsis is accompanied by immunological scarring (Figure 3E), where the barrier is severely damaged and S_v continues increasing above Sv*, while S_a remains above S⁻.

The four phenotypes, including a healthy phenotype and three disease phenotypes, correspond to different patient cohorts observed in the clinic. Healthy recovery from colonization is the most common outcome of host-pneumococcal interactions (Austrian, 1986) and corresponds to patients who can clear their symptoms from transient infection without any antibiotics treatment. Sepsis corresponds to patients who would develop systemic infection as a consequence of dysregulated transepithelial crossing of bacteria if no treatment is applied (Clarke et al., 2011; Siegel and Weiser, 2015). Immunological scarring corresponds to tissue-damaging inflammation that prevails even after clearance of the pathogens (Periselneris et al., 2015). The long-term deleterious consequences of such sterile inflammation and the associated tissue restructuring/damage are considered to underlie many diseases, including pulmonary fibrosis associated to previous Sp infections (Knippenberg et al., 2015), chronic obstructive pulmonary disease (Garcha et al., 2012) and cancer (Elinav et al., 2013; Pradere et al., 2014). A sustained activation TLR is recognized to be an important molecular player responsible for this tissue damage (Pradere et al., 2014), as in our model. Sepsis with immunological scarring corresponds to patient cohorts who would develop a severe infection with long-term deleterious effects in absence of treatment (Leibovici, 2013).

To identify the model parameters that affect the states of the R- and S_v-switches thereby determine the four phenotypes, we conducted the global parameter sensitivity analysis of our model with respect to R and S_v, respectively, using both Sobol and eFAST sensitivity indices (Marino et al., 2008; Cannavó, 2012).

The analysis identified the three most sensitive parameters for the propensity to turn on both the R-switch (to develop immunological scarring) and the S_v-switch (to develop sepsis) (Figure 4): the rate of bacterial transmigration through the barrier (θ_S), the bacterial carrying capacity (μ_S), and the killing rate of bacteria by macrophages (ϕ_MS) further confirming the importance of macrophages. Simulations with systematic variations of these three parameters further suggest that they affect the occurrence of sepsis and of immunological scarring, as well as how quickly these occur after the Sp challenge (Figure 5).

These three sensitive parameters have a direct correspondence with risk factors for disease triggered by Sp that have been reported in the experimental literature. For example, increase in θ_S can be caused by co-infection, which damages the barrier directly or by having triggered previous immune responses (McCullers, 2014). Increase in μ_S is caused by previous infections, for example by influenza virus, that damage the tissue, increase nutrient contents (Siegel et al., 2014), or shift the microbiome composition affecting the dynamics of the different bacterial populations (McCullers, 2014). ϕ_MS can be affected for example by severe asthma (Liang et al., 2014).

Other parameters that are also affected by co-infection were not identified to be very sensitive for the propensity to develop sepsis (increase in S_v) or unresolved host responses (increase in R) (Figure 4). These parameters include S⁺ which can increase as a consequence of TLR2 desensitization caused by a previous influenza virus infection (Didierlaurent et al., 2008), M that may increase as a consequence of previous infectious events (La Gruta et al., 2007; Yin et al., 2013), and the size of the neutrophil pool (N_v) which may decrease by chemotherapy or severe infections (Dick et al., 2008). The unsensitivity to the initial conditions can be partially explained by the existence of a unique stable steady state corresponding to the healthy recovery.

In the septic behavior observed in 15% of the simulations (Figures 3C,D), the sepsis occurred (S_v increases above Sv*) within 2 days post Sp challenge in 79% of the cases (Figure 6). When sepsis is accompanied with immunological scarring (R stays on and S_v increases above Sv*, Figure 3D), the time to sepsis is longer than when it is not (Figure 6). The computationally predicted rapid onset of the sepsis is consistent with experimental observations in Andonegui et al. (2009) that the mice either survived or died within 36 h upon Sp challenge applied directly into the lumen of the lungs. The results suggest that rapid treatments within 36 h are crucial to prevent the onset of sepsis.

Increasing the immune activity, for example by activation of adaptive immune responses, could be an effective way to decrease the risk of sepsis onset, as it elevates the switching threshold, Sv*, which depends on the strength of the resident immune cells. While the adaptive immunity could be activated naturally, the time for activation of the adaptive immune responses (which involves the de novo differentiation of naive T cells into mature T cells) by infiltrated pathogens was experimentally evaluated to be more than 2 days in mice (Zheng and Flavell, 1997). Such slow activation of the adaptive immunity therefore cannot prevent the onset of sepsis within 36 h.

These results suggest that prophylactic activation of the adaptive immune responses, for example by vaccinations, could be an effective strategy to prevent the incidence of sepsis, as demonstrated by the protective effects of Sp vaccination in mice (Cao et al., 2013). It is also consistent with the clinical suggestions to use vaccines as a preventive measurement against transient bacterial challenges in the high-risk patients (World Health Organization, 2012).

Using the proposed model, we investigate optimal treatment regimens and determine the minimal strength and duration of antibiotics treatment that are required to prevent or revert the pathological consequences of a transient Sp challenge. The minimal use of antibiotics is important for tackling the problem of antibiotics resistance (Nuorti et al., 1998), since the emergence of antibiotic-resistant Sp strains has been associated to the excessive use of antibiotics (Schrag et al., 2001; Prina et al., 2015b). We consider two different types of bactricidal antibiotics treatments in our modeling framework: apical application of antibiotics in the luminal side of the mucosa that decreases S_a and can thereby turn off the R-switch and stop the immunological scarring, and systemic application of antibiotics in the blood vessel that directly decreases S_v to prevent the onset of invasive infection (described in the Methods Section 4.3).

When the patients have immunological scarring without sepsis (Figure 3C), the treatment by apical application of antibiotics should aim to reduce S_a down below S⁻ to turn off the R-switch (Figure 7A). Once the R-switch is turned off, further use of antibiotics is no longer needed, as the healthy steady state with R = R_off is locally attractive (Sa<S+) for all the parameter combinations tested (over 10, 000 simulations). The minimal treatment potency of apically applied antibiotics (minimal strength × minimal duration) to bring S_a down below S⁻ depends on the severity of the phenotype measured by the deviation of the high focal point from S⁻ (R² = 0.46804).

When the patients are susceptible for sepsis (Figure 3D), the treatment by systemic application of antibiotics should aim to reduce S_v to avoid reaching Sv* and thereby causing invasive infection (Figure 7B). The minimal strength of systemically applied antibiotics allows the maximum of S_v to reach just below Sv*, and the minimal duration of the treatment with the minimal strength corresponds to the time required for R(t) to naturally turn off by S_a reaching S⁻. The minimal treatment potency of systemic antibiotics to prevent sepsis depends on the severity of the phenotype measured by the time to reach Sv* in the absence of treatments.

When the patients are susceptible to the combination of sepsis and immunological scarring (Figure 3E), they require antibiotics that are strong enough to be apically applied until S_a decreases below S⁻ to turn off the R-switch (Figure 7C). Our model simulations predicted that the apical treatment is enough to prevent invasive infection for 48% of these cases, since the reduction of S_a also reduces S_v, but the remaining 52% of the cases require additional application of comparatively small amounts of antibiotics directly in the blood vessel.

The distributions of the minimal treatment strengths and durations that we computationally predicted can be used as a guide to design safe and effective treatment options for the three patient cohorts. For example, our results suggest that antibiotics treatment for 20 days can prevent or revert most of immune scarring (Figure 7A) or invasive infection (Figure 7B), but that a much longer antibiotics treatment is needed for patients with a propensity for both sepsis and immune scarring (Figure 7C).

The proposed model for commensal bacterial infection describes the dynamics of bacterial load on the surface of the airway epithelium barrier, S_a(t)[CFU/ml], infiltrated bacterial load, S_v(t)[CFU/ml], concentrations of neutrophils and macrophages on the surface of the mucosal barrier, N(t) and M(t) [cells/ml], and the strength of barrier integrity, B(t) relative to the maximum strength, by

The variable R(S_a(t)) denotes the S_a-dependent TLR activation level described by a perfect switch,

where t⁻ is a time slightly before the time t. The dynamics of the TLR activity stabilizes within hours (Filewod et al., 2009; Witt et al., 2009; Hoffman et al., 2015).

The growth of the bacterial load on the apical side of the barrier, S_a, is modeled by a logistic equation (Smith et al., 2011), where the growth rate is limited by a carrying capacity (saturation term μ_S) that reflects the limited availability of nutrients in the epithelial lumen (Burnaugh et al., 2008). S_a is eradicated by immune cells, N and M, and transmigrates to the basal side of the epithelial barrier. The transmigrated bacteria, S_v, is assumed to grow exponentially in the blood vessel with abundant nutrients, but is killed by resident immune cells. The capacity to contain S_v is described by the saturated degradation of S_v, leading to the complete decay of S_v if it is below the threshold Sv*, which corresponds to the unstable steady state of the ODE for S_v when Sa=S-.

Recruitment of neutrophils (N) and macrophages (M) to the site of infection from their respective pool in the blood vessel (N_v) and the airway tissues (M_v) is inhibited by the epithelium barrier integrity (B), and is enhanced by the TLR activation (R) and the recruited neutrophils, respectively. Recruitment of N is further inhibited by M. N and M decay with a respective constant decay rate, as de novo production of the immune cells does not occur outside the bone marrow (Tak et al., 2013) and they do not divide in the epithelial tissue.

The self-recovery of the mucosal barrier to its homeostatic level (Nusrat et al., 1997; Coyne et al., 2002; Heijink et al., 2012) is modeled in a phenomenological manner, with the recovery rate being compromised by a decreased gene expression of epithelial cell differentiation markers (Clarke et al., 2011) induced by TLR activation. The barrier is directly damaged by transmigrating neutrophils and by proteases that are activated via TLR signaling (Chun and Prince, 2009). The switch-like activation of the TLR signaling is triggered by apically located bacteria, and is modeled by a phenomenological representation (Mochan et al., 2014; Domínguez-Hüttinger et al., in press). For simplicity, the inhibition by x is modeled phenomenologically by 11 + x.

All the numerical model analysis was conducted using MATLAB version R2014a (The MathWorks, Inc., Natick, MA, USA). Numerical integration was conducted by ode15s from the initial conditions corresponding to a transient Sp challenges with Sa(0)=107 CFU/ml, S_v(0) = 0 CFU/ml, N(0) = 0 cells/ml, M(0) = 10 cells/ml, and B(0) = 1, with R(0) = 1, as in the experiments in Zhang et al. (2009). The switch-dependent governing equation was chosen by the event-location function.

We model the effects of bactericidal antibiotics, such as penicillin, ceftriaxone and amoxilin, which are commonly prescribed to treat pneumococcal infection (Mandell et al., 2007; Prina et al., 2015b) by dSvdt=-VSv(t) (systemic application of antibiotics) and dSadt=-ASa(t) (apical application of antibiotics), where V and A represent a constant strength of antibiotics that kill S_v(t) and S_a(t), respectively. The strength of the antibiotics (with a unit of 1/h) is described by V = E_VD_V or A = E_AD_A, where E_V and E_A is the antibiotics killing efficacy and D_V and D_A is the amount applied, and can be chosen in the clinic by either selecting an antibiotic with a particular killing efficacy and/or adjusting the dose administrated. The killing efficacy of an antibiotic over a specific bacterial strain is commonly evaluated by the Minimal Inhibitory Concentration (MIC), the lowest concentration of antibiotics that will inhibit the visible growth of a bacterium after overnight incubation in a kinetic growth assay. From such experimental information, the antibiotics efficacy in our model, E_A (and E_V) can be calculated as EA=κSμS(1-SMIC(t))×1DAMIC, where S_MIC(t) is the concentration of Sp exposed to an antibiotics dose of DAMIC = MIC, and hence does not increase further. This expression is obtained from the steady state equation dSMIC(t)dt=κSμS(1-SMIC(t))-EADAMICSMIC(t)=0, which holds for the value of S_MIC(1day) = S_a in these experiments.

Minimal strength of antibiotics treatment to achieve remission was determined by checking whether Sv(t)<Sv* or Sa(t)<S- is achieved while gradually increasing V or A, respectively, by an increment of 0.01. The minimal duration of treatments corresponds to the time required to achieve Sa=S-. The minimal antibiotics treatment regimens to revert immune scarring without invasive infection (Figure 7A) were calculated under the assumptions that the treatment starts when S_a reached its steady state, after a transient S_a challenge that was modeled with initial conditions of Sa(0)=107 CFU/ml, S_v(0) = 0 CFU/ml, N(0) = 0 cells/ml, M(0) = 10 cells/ml, and B(0) = 1, with R(0) = 1. For the regimens to prevent invasive infection (Figure 7B), as well as to reverse immunological scarring and to prevent sepsis (Figure 7C), we assumed that the treatment starts at the time of the S_a challenge.

We varied all the model parameters over one order of magnitude (0.1 - 10 times) around the nominal values, and the initial conditions N(0), M(0) and B(0) within the intervals [0 1,000], [0 100] and [0 1], respectively. Robustness of the healthy behavior was tested by simultaneously varying all the parameters which were sampled from uniform distributions for 10, 000 iterations. The global parameter sensitivity was evaluated by 10, 000 iterations using the Global Sensitivity Analysis Toolbox for Matlab (Cannavó, 2012), with respect to the final concentrations of S_v and R at 7 days post Sp challenge with Sa(0)=107 CFU/ml.

We conducted a focal point analysis to determine the long-term behavior of the model in absence of a Sp challenge. Following the methodology in Oyarzún et al. (2012), we considered two subsystems for Equation (1), defined by fixing R to either R = R_off or to R = R_on, and evaluated the local stability of their steady states.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.