Mutations in bacterial genes induce unanticipated changes in the relationship between bacterial pathogens in experimental otitis media

Otitis media (OM) is a common polymicrobial infection of the middle ear in children under the age of 15 years. A widely used experimental strategy to analyse roles of specific phenotypes of bacterial pathogens of OM is to study changes in co-infection kinetics of bacterial populations in animal models when a wild-type bacterial strain is replaced by a specific isogenic mutant strain in the co-inoculating mixtures. As relationships between the OM bacterial pathogens within the host are regulated by many interlinked processes, connecting the changes in the co-infection kinetics to a bacterial phenotype can be challenging. We investigated middle ear co-infections in adult chinchillas (Chinchilla lanigera) by two major OM pathogens: non-typeable Haemophilus influenzae (NTHi) and Moraxella catarrhalis (Mcat), as well as isogenic mutant strains in each bacterial species. We analysed the infection kinetic data using Lotka–Volterra population dynamics, maximum entropy inference and Akaike information criteria-(AIC)-based model selection. We found that changes in relationships between the bacterial pathogens that were not anticipated in the design of the co-infection experiments involving mutant strains are common and were strong regulators of the co-infecting bacterial populations. The framework developed here allows for a systematic analysis of host–host variations of bacterial populations and small sizes of animal cohorts in co-infection experiments to quantify the role of specific mutant strains in changing the infection kinetics. Our combined approach can be used to analyse the functional footprint of mutant strains in regulating co-infection kinetics in models of experimental OM and other polymicrobial diseases.

JD, 0000-0001-9649-4698 Otitis media (OM) is a common polymicrobial infection of the middle ear in children under the age of 15 years. A widely used experimental strategy to analyse roles of specific phenotypes of bacterial pathogens of OM is to study changes in co-infection kinetics of bacterial populations in animal models when a wild-type bacterial strain is replaced by a specific isogenic mutant strain in the co-inoculating mixtures. As relationships between the OM bacterial pathogens within the host are regulated by many interlinked processes, connecting the changes in the co-infection kinetics to a bacterial phenotype can be challenging. We investigated middle ear co-infections in adult chinchillas (Chinchilla lanigera) by two major OM pathogens: non-typeable Haemophilus influenzae (NTHi) and Moraxella catarrhalis (Mcat), as well as isogenic mutant strains in each bacterial species. We analysed the infection kinetic data using Lotka -Volterra population dynamics, maximum entropy inference and Akaike information criteria-(AIC)based model selection. We found that changes in relationships between the bacterial pathogens that were not anticipated in the design of the co-infection experiments involving mutant strains are common and were strong & 2018 The Authors. Published by the Royal Society under the terms of the Creative

Introduction
Otitis media (OM) is a common polymicrobial bacterial infection of the middle ear in children which is caused by three major bacterial pathogens: non-typeable Haemophilus influenzae (NTHi), Moraxella catarrhalis (Mcat) and Streptococcus pneumoniae (Sp) [1]. The relationships among these OM pathogens are both direct and indirect in nature. For example, quorum signals (autoinducer-2 or AI-2) secreted by NTHi help Mcat to form a biofilm and survive in the hostile middle ear environment [2]. This interaction represents a direct relationship (or an active interaction [3,4]) between NTHi and Mcat. In another case, NTHi stimulates the host immune response in the middle ear that suppresses the growth of Sp [5,6]; this interaction is an example of an indirect relationship (or a passive interaction) between NTHi and Sp. The qualitative (cooperative, competitive or neutral) and quantitative (interaction strength) nature of the active and passive interactions between the OM bacterial pathogens depend on phenotypes specific to bacterial strains and the host response [4,7]. Mechanistic understanding of how these interactions affect pathogenesis of polymicrobial diseases including OM has been a major research goal for developing vaccine candidates and other therapeutic strategies [8,9].
A common strategy to evaluate mechanistic roles of specific phenotypes of bacterial OM pathogens in vivo has been to co-infect animal models with bacterial pathogens obtained from clinical isolates and then assess changes in infection kinetics by replacing a wild-type bacterial strain with a mutant strain [2, 6,10,11]. The mutant strains are designed to produce a loss or gain of specific bacterial phenotypes of interest. However, because the bacterial phenotypes probed by a mutant strain can be tightly intertwined with the phenotypes of the other bacterial pathogens, via active and passive interactions, this task could become challenging. For example, a mutant strain of Mcat lacking the ability to receive quorum signal from NTHi conceivably results in a decrease in the cooperative interaction from NTHi to Mcat. However, the same mutation could produce unanticipated changes in other relationships such as change in cooperation/competition of Mcat to NTHi. When these unanticipated changes are strong regulators of bacterial populations in co-infection experiments, we are required to revise our mechanistic understanding regarding the role of the specific mutation in influencing co-infection kinetics.
To this end, we address the above challenge by developing a framework that provides an answer to the following question: How is it possible to assess if unanticipated changes in the relationships induced by introducing mutant strains of OM pathogens in co-infection experiments are strong or weak regulators of the bacterial populations in the experiments? We define loss or gain of phenotype(s) in a specific mutant strain as a weak regulator when the interactions between bacterial species in co-infection experiments with the mutant strain are modified according to changes in the phenotype(s) as hypothesized for the mutant strain. The change in the phenotype(s) is defined as a strong regulator when additional unanticipated interactions are altered in co-infection experiments with the mutant strain. A more precise and formal definition of the weak and strong regulators is provided in the Material and methods section. The answer to the above question will provide a quantitative way to evaluate the mechanistic role of a bacterial gene in affecting the co-infection kinetics. Our framework combines (1) in vivo bacterial load measurements in an animal model, with (2) in silico approaches comprising Lotka-Volterra (LV) population dynamic models [12,13], maximum entropy (MaxEnt) inference [14 -17] and Akaike information criterion (AIC)-based model selection [18]. The animal model we used is a Chinchilla lanigera experimental OM model [11], wherein the animals' middle ears are co-inoculated with NTHi (86-028NP) and Mcat. These strains may be wild-type or isogenic mutant strains.
Our study revealed three important findings. First, the unanticipated changes in the relationships between OM bacterial pathogens that substantially affect the co-infection kinetics are commonly present in co-infection with mutant bacterial strains in experimental OM. Second, several bacterial phenotypes are tightly correlated across co-infecting bacterial stains. Third, our combined framework provides a systematic way to deal with two common difficulties faced when analysing infection rsos.royalsocietypublishing.org R. Soc. open sci. 5: 180810 kinetic measurements in animal models: host -host variations of bacterial populations and small size of animal cohorts. Our framework can be used to design mutant strains to generate desired infection kinetics in experimental models of polymicrobial diseases such as OM and infections secondary to cystic fibrosis [3] with potential implications for therapeutic strategies [9] Figure 1. Schematic representation of our framework to determine roles of mutant bacterial strains in regulating co-infection kinetics. (a) Inter-and intra-species interactions between two bacterial species NTHi and Mcat residing within a host can be both active (solid lines) or passive (dashed lines) in nature. These interactions can be simplified and described by LV interaction parameters (fa ij g). a 11 (greater than 0) and a 22 (greater than 0) represent intra-species interactions for NTHi and Mcat, respectively. a 12 and a 21 represent the overall effect of Mcat on the growth of NTHi and NTHi on the growth of Mcat, respectively. a 12 and a 21 can be positive (competitive interaction), zero (neutral interaction) or negative (cooperative interaction). (b) Replacing a wild-type strain by a mutant strain in the co-infection experiments can change the LV interactions. These changes may be anticipated (blue 'X') or not anticipated (red 'X') based on the design of the experiment. Our framework uses data from co-infection experiments involving the wild-type strains to generate models that determine if these unanticipated changes are weak or strong regulators of bacterial kinetics. These models are compared to each other using AIC.  (see table 1 for details). Therefore, in the design of the co-infection with NTHi (wt) þ Mcat (hag), one would anticipate a lower carrying capacity for Mcat or an increase of a 22 . However, the changes in passive interactions induced by the hag mutant can also lead to changes in other LV interactions that were not anticipated, such as an increase in a 21 . In addition, the host-host variations measured in the NTHi and Mcat populations were assumed to arise from the variations of the LV interactions (fa ij g) in our models. We developed a framework (figure 1b) to address the above question. The framework is executed in two main steps.
Step 1. The wild-type NTHi (or N 1 ) and the wild-type Mcat (or N 2 ) populations in individual chinchillas were measured in co-infection experiments carried out in a cohort of n number of chinchillas at a time T (7 or 14 days) post inoculation. These data were used to generate the reference data (denoted by the subscript r) n D r . This reference dataset was composed of the mean values ð N 1 , N 2 Þ, variances ðs 2 1 , s 2 2 Þ and the covariance ðr 12 Þ calculated from the bacterial load measurements in the cohort; i.e. n D r ; f N 1 , N 2 , s 2 1 , s 2 2 , r 12 g ; f n D j r g for j ¼ 1, . . . , 5. We estimated the probability distribution functionpðN 1 , N 2 Þ of populations of wild-type NTHi (N 1 ) and wild-type Mcat (N 2 ) in the chinchilla cohort using MaxEnt (electronic supplementary material, figure S1), wherein the mean values, variances and the covariance as measured in the experiments ( n D r ) were constrained in the calculation (electronic supplemental material, §S2). Using thispðN 1 , N 2 Þ, we calculated the joint probability distribution of the interaction parametersqða 11 , a 12 , a 21 , a 22 Þ ð;qðfa ij gÞ using a MaxEntbased method (details in the Material and methods section and electronic supplementary material, figures S2 -S3). The estimated joint probability distribution functionqðfa ij gÞwas used to generate models that fall either in the weak or the strong category.
Step 2. We generated a test dataset using data from co-infection experiments where at least one of the wild-type bacterial strains was replaced by a mutant strain. The test dataset n 0 D x contains the population means, variances and covariance for the cohort which contained n 0 number of chinchillas for the same time T (7 or 14 days) post inoculation. The subscript x in n 0 D x , to be determined in our analysis, quantifies the role of the mutation in the co-infection: x ¼ weak or x ¼ strong. As described above, this distinction indicates whether the unanticipated changes in LV interactions induced by the mutant strain(s) were weak or strong regulators of the bacterial populations. To safeguard against small sizes of n 0 (approx. 10 or less), we performed bootstrapping [23] on the data, wherein we sampled n 0 data points with replacement from the original. In this way, we generated t sets and determined x in each of those t samples of n 0 D x . These t samples of n 0 D x , e.g. ð1Þ n 0 D x , . . . , ðtÞ n 0 D x , were generated by evaluating n 0 D x in t independent groups containing n 0 number of animals each. Next we evaluated which model(s) generated in step 1 for a particular co-infection experiment best described the t samples of the test dataset ð ð1Þ n 0 D x , . . . , ðtÞ n 0 D x Þ; subsequently, we assigned the category (weak or strong) of the best model to x. The best model was found (whenever possible) by comparing the Akaike information criterion (AIC) values for each model in a head-to-head pairwise manner for each of the t samples. The 'Condorcet Winner' was the model which was preferred over all others in head-to-head comparisons [24]. We chose the Condorcet winner as the best model. Further details are provided in the Material and methods section and in the electronic supplementary material.

Application of the framework on synthetic data
To test the efficacy of our method, we generated synthetic data and applied the framework developed in the previous section. We had the following three goals in mind: (i) validate the framework, (ii) determine how the strengths of the mutations and/or the host immune response affect interspecies interactions and (iii) determine the dependence of the model selection on the sample size or the number (n 0 ) of animals in the cohort. We generated the synthetic data by numerical solution of coupled ODEs that described LVtype population kinetics involving two interacting bacterial species and a host immune response (see Material and methods section). The parameters describing the inter-and intra-species bacterial interactions as well as the host immune response were drawn from uniform distribution within rsos.royalsocietypublishing.org R. Soc. open sci. 5: 180810 Table 1.  [20,21]. This mutant could lead to lower biomass or decreased carrying capacity (or higher The mcaB mutant is unable to sense the AI-2 molecules secreted by NTHi (Li et al., unpublished data). Thus, this Mcat mutant will experience less cooperation from NTHi, or increased values of The aaa mutant has the Hag programme constitutively on (Li et al., unpublished data). Therefore, it makes dense biofilms without any help from NTHi and has a larger carrying capacity or decreased value of The luxS mutant is unable to produce AI-2, produces less dense biofilms and as a result causes more acute disease with higher levels of early inflammation. Thus the mutant can produce lower cooperation to Mcat (higher a 21 ) due to the lack of quorum sensing, and lower values of carrying capacity (higher a 11 ) due to poor biofilm quality [22]. However, the increased inflammation could provide nutrients for self-growth (lower The mclR mutant does not respond to NTHi quorum sensing and forms poor biofilms (Li et al., unpublished data). Thus, Mcat will receive less cooperation from NTHi (higher a 21 ) and give less cooperation to NTHi (higher The dtgt mutant is unable to respond to quorum signals from NTHi and also has a mutation in the hag promoter. The latter effect disrupts Mcat production of this adhesion protein and thus adhesion to the substrate (presumably host epithelia but also potentially aggregation with other bacteria) ( 6 specific ranges to generate host -host variations of the co-infection kinetics (see Material and methods). We chose the parameter range for the wild-type strains such that it produced steady-state population values similar to those observed in vivo (electronic supplementary material, figures S6-S8). The average statistical variables n D r ; f N 1 , N 2 , s 2 1 , s 2 2 r 12 g, calculated from the numerical solutions at T ¼ day 7, produced the reference dataset. We generated test datasets f n 0 D x g, wherein one of the wildtype strains was replaced by a mutant strain. Specifically, we considered the following two mutant strains for species 1. The ða ðþÞ 11 Þ mutant, which is an increase in the a 11 parameter, possesses increased self-competition for species 1 compared to wild-type. The ða ðþÞ 21 Þ mutant, which is an increase in the a 21 parameter, possesses increased competition of species 1 towards species 2 as compared to wildtype. We also considered two mutant strains for species 2. The ða ðþÞ 22 Þ mutant, which is an increase in the a 22 parameter, possesses increased self-competition for species 2 compared to wild-type. The ða ðþÞ 12 Þ mutant, which is an increase in the a 12 parameter, possesses increased competition of species 2 towards species 1 as compared to wild-type. The mutants were generated by changing the ranges of the associated parameters from that of the wild-type strains (see Material and methods). For example, the range of a 11 used to generate the mutant strain a ðþÞ 11 spanned a smaller range [a 0 , b] compared that of the wild-type strain, [a, b], where, a , a 0 . Each of these mutations was performed with low, moderate and large strengths based on the relative change in the magnitude of range for the parameters. We also solved the co-infection kinetics in the presence of no, weak and strong host immune response. Thus, in total we considered 9 Â 4 different mutation experiments in silico (electronic supplementary material, figures S9 -S12). The test datasets f n 0 D x g were obtained from the co-infection kinetics involving the above mutant strains.
We show results for two mutant strains a ðþÞ 12 and a ðþÞ 22 with moderate strength mutations in the absence of any host immune response (figure 2b,c). The rest of the mutants are described in the electronic supplementary material. We followed the steps described in electronic supplementary material §S1 to determine the nature of the mutation (or x) in a test dataset n 0 D x . First, we used n D r to generate models that belonged to the weak or the strong category corresponding to the co-infection strain#1(wt)þstrain#2ða ðþÞ 22 Þ or strain#1(wt)þstrain#2ða ðþÞ 12 Þ. Next, we compared the models with the samples (t . 100) of the test dataset and evaluated the Condorcet winner model. We found that for the co-infection with strain#1(wt) þ strain#2ða strain will be small. These unanticipated changes would play a weak role in regulating the bacterial populations. By contrast, a 12 is correlated strongly with several other interaction parameters (e.g. a 21 , Corr % 20.5) (electronic supplementary material, figure S6). Thus increasing a 12 even by a moderate amount, as in the a ðþÞ 12 strain, will produce large changes in the other LV parameters. Those unanticipated changes will generate a strong effect on the bacterial populations. These expectations about the effects of the correlations among the parameters were consistent with the results obtained from our framework. Therefore, these results validate our framework. The roles of the above mutant strains in affecting co-infection kinetics for different mutation strengths change depending on the mutation strength and/or the presence of the host immune response (electronic supplementary material, figures S9 -S12). Therefore, strengths of the mutations and the host immune response are important in determining the influence of the mutations in the co-infection kinetics.
We checked the dependence of the Condorcet winner on the sample size n 0 in the test dataset. We found that even for small sample sizes (n 0 ¼ 10), the framework picked the correct Condorcet winner; however, the margin of victory increased with larger n 0 (figure 2f,g). The data were collected from both ears of the chinchillas in cohorts containing more than five animals. The bacterial counts showed large to moderate host-host variations (electronic supplementary material, figure S1). These variations could arise from the host -host differences in the physiology, anatomy and immune responses in the upper respiratory tract of the outbred population of chinchillas. The mean NTHi population was substantially larger (greater than 100 fold) than that of Mcat for the co-inoculation with the wild-type NTHi and Mcat strains (figure 3). The mean populations of the wild-type NTHi strain at day 7 were lower if co-inoculated with any mutant Mcat strain (except Mcat (mcaB)) rather than if co-inoculated with wild-type Mcat (electronic supplementary material, table S1). By contrast, the mean populations of the Mcat strains for the same co-infections  (electronic supplementary  material, table S1). In a few cases, such as NTHi (wt) þ Mcat (hag) NTHi or NTHi (wt) þ Mcat (mcaB) the covariances were positive (electronic supplementary material, table S1). Overall, the data showed a complex pattern as further explained below.

Analysis of the in vivo data
Changes in the bacterial counts due to mutations pointed towards the presence of unanticipated changes in the bacterial relationships in regulating bacterial populations. For example, compared to co-infecting with NTHi (wt) þ Mcat (wt), co-infecting with NTHi (wt) þ Mcat (hag) substantially decreases the NTHi population (almost by half ), whereas the Mcat population only decreases a small amount. As the hag mutant has lower adherence and poor biofilm formation capability compared to its wild-type counterpart, we would expect the Mcat population to decrease while having minimal consequence to the NTHi population. Because the NTHi population was substantially reduced here, it suggests a potential change in the interaction from Mcat towards NTHi. We used our framework to quantify the roles of specific NTHi and Mcat mutations in regulating the bacterial populations in coinfection experiments. As described in the previous section, the data from the co-infection experiments with the wild-type strains generated our reference dataset n D r . The models in the weak or the strong category were generated using n D r and were compared against the test datasets f n 0 D x g. The test datasets were obtained from the co-infection experiments that involved at least one mutant strain. Multiple samples of a test dataset were obtained by using bootstrapping [23]. Our analysis showed that for the majority of the cases, models with additional interactions (strong models) better described the data (at both day 7 and day 14) compared to models with no additional interactions (weak models) (figure 4). We found that the weak model described the data obtained at day 7 for the co-   figure S4) described the same data better than any weak model. Therefore, unanticipated changes in LV interactions were prevalent in co-infections with mutant bacterial strains both at early and late stages of the co-infection kinetics.

Host immune responses modulate Mcat-NTHi interactions at later stages of the infection
The mean populations of the wild-type strains of NTHi and Mcat increase as the infection progresses from day 7 to day 14 post inoculation ( figure 5). However, the covariance of the NTHi and Mcat populations becomes more negative (approx. twofold change) (electronic supplementary material, table S1). The negative correlation indicates that the populations of the two species are more mutually exclusive; that is, when one species has high abundance, then the other's is low. We hypothesized that the host immune response generated by both the pathogens could lead to a decrease in cooperation (or increase in a 12 and a 21 ) between Mcat and NTHi. We tested our hypothesis by applying our scheme with the day 7 data as the reference and day 14 data as the test. The model wherein a 12 and a 21 increase was the Condorcet winner ( figure 5). The agreement demonstrates the role of the host immune response in regulating passive inter-species interaction between Mcat and NTHi.

Discussion
Co-infection of animal models with mutant bacterial strains is a powerful tool in probing mechanisms that underlie pathogenesis of polymicrobial infections such as OM. However, the interconnected and variable nature of interactions involving bacterial pathogens within the host makes it challenging to connect specific perturbations, such as a mutation, in these experiments to mechanisms. The datadriven framework developed here provides a systematic method of addressing this challenge. The framework uses bacterial counts measured in animal hosts to quantitatively determine perturbations in the bacterial interactions induced by the replacement of a wild-type bacterial strain with a mutant strain that strongly or weakly regulates bacterial populations in the co-infection. Therefore, using this framework we are able to quantitatively assess the mechanistic role of a specific bacterial phenotype probed by an isogenic mutant strain in affecting the co-infection kinetics. Isogenic mutant strains are used for identification of bacterial determinants of colonization, persistence and virulence. Thus the quantitative information obtained from our framework will be valuable for determining specific targets for diagnostics, development of therapy and potentially vaccination. In addition, our framework also addresses the practical problem of systematically analysing bacterial count data obtained from a small size (approx. 10 animals) of animal cohorts in co-infection experiments. Application of our framework to co-infection experiments in Chinchilla lanigera co-inoculated with wild-type and mutant strains of two major OM pathogens (NTHi and Mcat) found that in a majority of the co-infections the mutant strains gave rise to unanticipated changes in the bacterial interactions, which influenced the bacterial populations substantially. The emergence of unanticipated perturbations of the bacterial interactions is probably caused by the interdependencies between the interactions and the hostile environment of the host [16]. The interdependencies can be caused due to many shared processes such as feeding on common nutrients, exchanging small molecules and the host immune response that regulates the growth of OM pathogens within the host [5,8]. Therefore, when a specific bacterial phenotype is altered in the form of a mutant strain, several other phenotypes in co-infecting OM pathogens are also altered, some of which can be non-intuitive. Our analysis of the synthetic co-infection data lends support to this speculation. We found that correlations between the LV parameters generated non-intuitive changes in several LV interactions when a specific LV interaction was perturbed in a mutant strain, and in many cases these unanticipated changes in the LV interactions were strong regulators of the bacterial populations.
Our framework required estimation of LV interactions involving the co-infecting bacterial pathogens using the measured bacterial counts. The estimated interactions between wild-type strains of NTHi and Mcat demonstrate prevalence of cooperative interspecies interactions (a 12 , 0, a 21 , 0) (electronic supplementary material, figure S3). Previous experiments noted several molecular mechanisms regarding the help of NTHi towards Mcat's growth, e.g. the quorum signal AI-2 secreted by NTHi helping Mcat to form a biofilm and thereby helping it to survive within the host [22,25]. Reciprocally, the estimated interactions also suggested cooperation of Mcat towards the growth of NTHi (or a 21 , 0). Such a cooperative effect can potentially occur through passive interactions. For example, Mcat binds and sequesters AI-2 molecules secreted by NTHi; this sequestration could help keep the AI-2 abundances at an optimal level for production of quorum signals by NTHi. Similar optimal regulation of quorum sensing has been found in mutualistic relationship between two human oral bacteria, Actinomyces naeslundii T14 V and Streptococcus oralis 34, where AI-2 secreted from the latter bacterial species help the former species to form biofilms [26].  [20,21]. Another example of cooperation involves nutrient recruitment. The inflammation caused by Mcat can produce an influx of host serum which also can provide nutrients for the growth of NTHi [27]. All three of the above sources of interaction could contribute towards generating an overall cooperative interaction from Mcat to NTHi. We found that several LV interactions (e.g. a 12 and a 21 ) involving NTHi and Mcat are tightly correlated (jCorrj . 0.5) with each other. As the MaxEnt method estimates the most spread out or uniform probability distribution phenotype that is consistent with the measured data [17], the method used here provided the most conservative estimate correlations between the LV interactions. This approach is a major departure from several methods that have been developed in recent years to evaluate interactions [28 -30] pertaining to microbiome datasets where the correlations between the interaction parameters are not analysed. Ecological models often assume LV interactions between coexisting species as uncorrelated random variables [13]. This assumption makes the calculations amenable to analytical methods. The presence of strong correlations between LV interaction parameters could have important implications in assessing general principles underlying the diversity of eco-systems [31].

Experiments
Moraxella catarrhalis persistence in the middle ear chambers of chinchillas was assessed essentially as described previously [2]. Animals were purchased on need and allowed to acclimate to the vivarium for more than 10 days before infection. No used animals showed visible sign of illness prior to infection. Chinchillas (five animals per group) were anaesthetized with isofluorane and infected via transbullar injection with both approximately 10 4 CFU of M. catarrhalis and approximately 10 3 CFU of H. influenzae. All inocula were confirmed by plate counting. Animals were euthanized at 7 and 14 days post infection, and their bullae were aseptically opened to recover possible effusion fluid. Middle ear lavage was performed using sterile PBS and also saved and combined. Bullae were then excised and homogenized in 10 ml of sterile PBS. All of fluid and homogenized samples were serially diluted and plated on brain heart infusion (BHI) agar plates to obtain viable counts of M. catarrhalis. Note that H. influenzae would not grow on BHI plate due to lack of haemin and NAD. H. influenzae bacteria were enumerated on BHI agar supplemented with 10 mg ml 21 haeme and NAD and containing 5 mg ml 21 clarithromycin, which inhibits the growth of M. catarrhalis.

Bacterial strains and growth conditions
Moraxella catarrhalis strain O35E, as a WT and parent strain in this study, is a commonly used laboratory strain [32]. Moraxella catarrhalis O35E hag::Sp containing a spectinomycin resistance cassette disrupting the hag gene is a kind gift from Dr Eric Hansen. Non-typeable H. influenzae strain 86-028NP is a nasopharyngeal isolate from a child with chronic OM [33], and its luxS::Kn mutant with a kanamycin resistance cassette disrupting the AI-2 synthase was described previously [22]. Moraxella catarrhalis strains were cultured in BHI medium (Difco), and H. influenzae strains were cultivated in BHI medium supplemented with 10 mg ml 21 of haemin chloride (MP Biomedicals) and 10 mg ml 21 of NAD (Sigma), referred as supplemented BHI (sBHI). Mixture cultures of both M. catarrhalis and H. influenzae used sBHI.

Mutant bacterial strains
M. catarrhalis mclR::spec was generated by insertional mutagenesis using the following approach. Genomic DNA was purified from M. catarrhalis O35E using the Wizard genomic DNA purification kit (Promega), essentially according to the manufacturer's instructions. Portions of the mclR (allele MCR_1062) open reading frame were amplified using the PCR and primers specific for intragenic rsos.royalsocietypublishing.org R. Soc. open sci. 5: 180810 regions (luxRUPF: CATCATGACTTGGTAACTTGCTG, luxRUPR: GCTGATCGGCAATTTGCCCCC GGGGTCGAGTGGCTTCTACACC). The resulting amplicons were cloned and a spectinomycin resistance cassette was introduced into a SmaI restriction site within the intergenic primers. This mutant allele was introduced into the parental strain by natural transformation and the resistant derivatives were confirmed by PCR and DNA sequence analysis.

Hag promoter mutants
Deletion of a tgt sequence and insertion of an aaa sequence within a predicted lux box within the hag promoter was achieved by overlap PCR; the resulting allele was introduced by natural transformation using a linked spectinomycin resistance marker (Li et al., unpublished data).  ðpÞ ij } uniquely using the above equations. Therefore, we developed a MaxEnt-based inference scheme to estimate the parameters in the chinchilla population using the measured bacterial loads. The MaxEnt-based method estimates parameters based on the measured data without any additional prior assumption. This also implies that MaxEnt estimates the 'flattest' distribution that is consistent with the measured data. A recent work [34] used maximum caliber inference, which is an extension of MaxEnt for analysing timedependent data [15,17], to estimate parameters in a gene regulatory reaction network. Parameter estimation in gene regulatory reaction networks using sparse time-dependent data represents a problem of similar spirit as the problem investigated here. First, we estimated the probability distribution function of N 1 and N 2 ,pðN 1 , N 2 Þ, in the chinchilla population using MaxEnt (electronic supplementary material, figure S1). Then, we estimated the joint probability distribution function qðfa ij gÞ in the interaction parameters fa ij g using the estimatedpðN 1 , N 2 Þ by applying MaxEnt the second time (electronic supplementary material, figures S2 -S3). The details regarding the implementation of the method for the in vivo data are provided below and in the electronic supplementary material, §S2.

Estimation ofqðfa ij gÞ
We estimatedqðfa ij gÞ usingpðN 1 , N 2 Þ by applying the MaxEnt inference the second time. We discretize the four-dimensional space spanned by a 11 , a 12 , a 21 and a 22 on a grid. a 11 and a 22 assume only positive real values, and a 21 and a 12 can assume both positive and negative real values. Specifically, we discretized a 11 from 0.027 to 2 Â 10 6 CFU and a 22 from 18.9189 to 1400 Â 10 6 CFU into 74 bins each. We discretized a 12 from 22000 to 50 Â 10 6 CFU and a 21 from 250 to 1 Â 10 6 CFU into 201 bins each. We have varied the bounds and the lattice sizes of the grids and there was no change in the qualitative results (electronic supplementary material, figure S5). In this case, the Shannon's entropy [16,17] k(fa ij g) in the above equations denotes the degeneracy factor or the number of distinct points in the a space that produce the same value of N 1 and N 2 . We calculated k(fa ij g) numerically by counting the number of lattice points in the a space that map to the same lattice point in the N-space.

Determination of the role (x ¼ weak or x ¼ strong) of the mutant strain in co-infection kinetics
Our framework to evaluate the role of the mutant strain is divided into two main steps as outlined in the main text.
Step 1. We estimatedqðfa ij gÞ using the reference dataset n D r as described in §4.2. We used qðfa ij gÞ to generate models that belong to the weak or the strong category. The details regarding how these models were generated are given below. Weak models: A specific mutant strain is hypothesized to possess loss or gain of phenotype(s) at the design stage of the co-infection experiments (see table 1 for details). The hypothesized changes in the specific phenotypes for the mutant strain will result in changes in a subset of LV interaction parameters pertaining to the co-infection kinetics of the populations of the mutant strain and another bacterial species within the chinchilla host. This subset of LV interaction parameters is denoted by fa pq g , fa ij g, where p ¼ i and q ¼ j for each a ij in fa pq g. E.g. for the co-infection with NTHi (wt) þ Mcat (mcaB), the mcaB mutant strain increases the value of a 21 , thus the set fa pq g contains only one parameter a 21 , whereas for co-infection with NTHi (wt) þ Mcat (dtgt), the dtgt mutant strain increases both a 21 and a 22 , and fa pq g will contain two parameters, a 21 and a 22 . The models in the weak category are defined by the probability distribution function of the interactions parameters (q w (a;a)) in the models. The weak models are parametrized by fa k g which quantifies the extent by which each of the interaction parameters in fa pq g is perturbed inqðfa ij gÞ to generate q w (a;a), i.e. q w ða; aÞ ¼qðfa ij gÞ Y k[K Qða k , a pq Þ, ð4:7Þ where K ¼ f( p, q)g. Q(a k , a pq ) denotes the Heaviside theta function Q(a k 2a pq ) or Q(a pq 2 a k ). Q(a k 2 a pq ) is chosen when the mutation decreases the value of a pq such that the values a pq . a k are absent in the mutant, or Q(a pq 2 a k ) is chosen when the mutation increases the values of a pq such that the lower values a pq , a k are absent in the mutant. To illustrate, the weak model for the co-infection NTHi (wt) þ Mcat (mcaB) is parametrized by a 1 and is generated using q w (a; a 1 ) ¼qðfa ij gÞQ(a 21 À a 1 ). Strong models: The models in the strong category considered unanticipated changes in the LV interaction parameters for co-infection kinetics where a wild-type strain is replaced by a mutant strain. If a mutant strain in a co-infection is hypothesized to change a subset of LV interaction parameters fa pq g, then the strong models consider changes in fa pq g and additional LV interactions (or unanticipated changes) outside of fa pq g. Similar to the weak models, the strong models are defined by the probability distribution function of the interactions parameters (q s (a; a)) and are parametrized by fa k g which quantifies the extent by which each of the interaction parameters is perturbed in qðfa ij gÞ to generate q s (a; a), i.e. In the above equations, H , fa ij g\K, where H defines the set of LV interactions outside fa pq g. For example, the strong models for co-infection NTHi (wt) þ Mcat (mcaB) were generated from in total six types of models: three types that vary pairs of interactions simultaneously, namely, q s (a; a 1 , a 2 ) ¼qðfa ij gÞQ(a 21 À a 1 )Q(a 11 , a 2 ), q s (a; a 1 , a 3 ) ¼qðfa ij gÞQ(a 21 À a 1 )Q(a 12 , a 3 ), q s (a; a 1 , a 4 ) 1 qðfa ij gÞQ(a 21 À a 1 )Q(a 22 , a 4 ), and three types that vary a triplet of interactions simultaneously, namely, We first created a reference dataset analogous to a NTHi (wt) þ Mcat (wt) co-infection experiment. Using the simplified LV two-species model (equation (4.4)), we set the following ranges for each of the interaction parameters: a 11 e [2.74 Â 10 23 , 0.2], a 12 e [2200,5], a 21 e [25, 0.1] and a 22 e [1.89189, 140]. These ranges were guided by the values observed from the in vivo data so that the in silico reference set was roughly similar to in vivo data. To calculate a single data point of a (N 1 , N 2 ) pair of abundances, we randomly drew a value for each parameter (assuming a uniform distribution) and calculated the steady-state abundances using equation (4.4). For all the datasets, we generated n ¼ 10 5 points per dataset. The reference dataset n D r was calculated by collating all the (N 1 , N 2 ) pairs obtained by drawing the a parameters from random distributions as described above.

Mutation sets
All 'mutations' of the reference dataset were done by increasing the value of exactly one interaction parameter; specifically, by raising the minimum value in the parameter's range. The three levels of the mutation indicated the severity of the increase. A low level a 11 mutation has a 11 e [0.027, 0.2]; at moderate severity: a 11 e [0.1, 0.2] and at large severity: a 11 e [0.18, 0.2]. The ranges for all twelve possibilities are shown in electronic supplementary material, figures S9 -S12.
We also aimed to evaluate the framework with respect to different levels of a host's immune response. In cases with no host response, we use the above method. For a non-zero immune response, we used a simplified model wherein only N 2 induces the immune response, and only N 1 is affected (detrimentally). Specifically, we introduced a Michaelis -Menten term to the first ODE in the standard LV two-species model: In that term, K m is the amount of N 2 necessary for the immune response to be at half-max. The susceptibility of N 1 to an immune attack is governed by the k d parameter. So, with no immune response, we set k d ¼ 0. For a weak immune response k d e [0,1]. For a strong immune response, k d e [0,10]. As before, when choosing a value for k d , we assumed a uniform distribution.
Animal ethics. All chinchilla infections were performed according to the protocols approved by the Wake Forest Animal Care and Use Committee. This study was conducted according to the guidelines outlined by National Science Foundation Animal Welfare Requirements and the Public Health Service Policy on the Humane Care and Use of Laboratory Animals. The Institutional Animal Care and Use Committees (IACUC) at Wake Forest School of Medicine approved these animal studies. Wake Forest's IACUC oversees the welfare, well-being and proper care and use of all vertebrate animals used for research and educational purposes at Wake Forest School of Medicine. The approved protocol number for the project is A13-140 and the Wake Forest School of Medicine animal welfare Assurance Number is D16-00248.