Optimal control and Bayes inference applied to complex microbial communities

1.   Introduction

The rapid surge of resistant strains and the void in discovering of new antimicrobials has made antimicrobial resistance one of the most critical problems in public health nowadays. There has been a great effort from the scientific community to understand the basic mechanisms by which bacteria become resistant to several antibiotics; almost all efforts have focused on the study of a single strain, and it is just in recent times that we have been looking at how the interaction between members of complex communities are shaping the susceptibility profile.

In nature, microbes co-exist in multi-species communities with ecological dynamics driven by the complex interaction between individual strains and the environment. Syntrophic interactions, whereby bacterial cells exchange costly metabolites for the benefit of both interacting partners, are pervasive in bacterial communities as they enhance bacterial survival in hostile environments.

Availability of resources affects the nature of the interactions between members of a microbial consortium ^[1]. These interactions can be determinants of the stability of the consortia. Based on two-pair interactions, there are typically three outcomes; cooperation where both members of the community are beneficial to the other one; competition where the presence of one member is detrimental to the other one and neutral interaction; where the members of the community have no interaction.

Typically competitive interactions end up in competitive exclusion, decreasing the diversity of the consortia. For syntrophic interactions, if the competition for resources is very strong, one community member's collapse could collapse the total population.

The effect of changes in interactions have on the control of antibiotic resistance is unknown. There are some examples where cooperation can help a community to survive in stressful environments ^[2]. We model a synthetic community composed of two different strains; both are auxotrophic to an amino acid, an essential metabolite. This condition makes them incapable of growing on their own in an environment without amino acids. In contrast, they can grow in a minimal media when they are together due to the leak of all the other amino acids from their partner ^[3]. Auxotroph cells are essential for natural communities ^[4], they shape the composition and can be determinant of the stability of a microbial community. For instance, they have been found in aquatic communities ^[5], microbial soil communities^[6] and the human microbiome ^[7,8]. Although auxotrophies could increase the dependence of individuals in the communities, they could also increase fitness by reducing the nutrient requirements of an environment, which means that a community could support growth in a poorer context.

In the past years, auxotrophic strains have been used to study metabolic cross-feeding interactions in particular environmental context ^[1,9,10]. Pair-wise interactions are useful to determine changes in the interaction and to study division of labor inside the community, such as public goods. A relevant example of public goods is the production of enzymes able to degrade antibiotics.

Cross-feeding interactions in natural environments are important to increase diversity in the population by decreasing the fitness of the individual and increasing the fitness of the community^[11]. It has been shown that these interactions are dynamic and can change with respect to the environment. Auxotrophies are a particular example of this. It has also been shown that auxotrophic phenotypes appear very fast in a population growing in a rich nutrient environment ^[10]. This feature can be crucial in diverse communities, like the microbiota, particularly when an antimicrobial is present, which disrupts the community members ^[4]. If auxotrophic phenotypes appear, having functional redundancies could help preserve general function in a community. For this reason, understanding what is the effect of an antimicrobial in a syntrophic consortium could be of great value.

This paper explores a population-based model and control theory to understand the effect of antimicrobials in microbial communities with syntrophic interactions and how the environment minimizes resistant bacteria in the system.

We first propose an optimal control problem based on Hoek et al. in ^[1] of an auxotrophic consortium of two bacteria competing for the same resources in the presence of a gradient of antibiotic. Then we described our experimental data and performed numerical experiments. Finally, using Bayesian Inference, we estimate some interesting mathematical model parameters using our experimental data. Numerical experiments validate all theoretical results.

2.   Optimal control formulation

In order to better understand the different types of interaction between cross-feeding bacteria, we formulate a similar mathematical model to the proposed by Hoek et al. In ^[1], two strains of yeast, one auxotrophic to tryptophan (trp-) and the other auxotrophic to leucine (leu-) interact in an environment with different resource availability. In this formulation, we assume that strain trp- ($X$) and strain leu- ($Y$) are divided into susceptible and resistant bacteria ($X_s$, $X_r$ and $Y_s$, $Y_r$, respectively). We assume both strains interact with a bacterial infection in an individual's body. We are interested in modeling the dynamics of bacterial resistance acquisition through antimicrobial resistance induced by antibiotic-mediated selection. Thus, the following hypothesis is considered for the strain $X$ (similar assumptions can be considered for the strain $Y$): sensitive bacteria grows at a rate proportional to resistant bacteria (if $a$ is the growth rate of sensitive bacteria, then $a\bar \gamma_x$ is the growth rate of resistant bacteria) following the logistic equation with a carrying capacity $K$. $\bar F$ is the amount of supplemented amino acids, $\bar \kappa$ represents an effective Monod constant, and $\beta$ is a parameter that quantifies the asymmetry of benefit that each bacteria receives from its partner. The density of sensitive bacteria become resistant by antibiotic-mediated selection is given through the term $\epsilon q\lambda X_s$, where $\epsilon$ is the efficacy of the antimicrobial, $\lambda$ is the antimicrobial supply rate, $q$ is the proportion of mutation (therefore ($1-q$) represents the proportion of elimination). Both population turnover at a constant dilution rate $\bar \delta$. Thus, the per capita growth rate of both bacteria is adjusted by the mutualistic partner as well as the supplemented amino acids:

A complete description of the parameters involved in Model (2.1) can be found on Table 1.

Now, we define $\pmb \mu = (\mu_1(t), \mu_2(t))$ as a control variable associated with the mutations. Thus, both strains of susceptible bacteria mutate at a rate $(1-\mu_1(t))\epsilon q \lambda X_s$ and $(1-\mu_2(t))\epsilon q\lambda Y_s$, respectively, where for $i = 1, 2$, $\mu_i(t)\in [0, 1]$ ($\mu_i = 0$ represents no efficacy of the control, while $\mu_i = 1$ indicates that the use of the control is completely effective). Thus, we have that each control variable $\mu_i(t)$ provides information about the amount of bacteria that must not mutate at time $t$. To minimize the number of resistant bacteria, we define the cost function:

where the parameters $c_1$ and $c_2$ represent social cost and the parameters $d_1$ and $d_2$ represent relative weights associated to the controls. Additionally, we define the boundary conditions:

and we assume that the final time $T$ is a fixed implementation time of the control strategies, the final state ${\bf X}_f$ is variable, and the initial state ${\bf X}_0$ is a coexistence equilibrium of System (2.1). We suppose that each control is in the set of Lebesgue measurable functions with $0\leq \mu(t)\leq 1$, $t\in [0, T]$ ($\mathcal{U}$ called the set of admissible controls).

Let us define $\omega$ as a rescaling parameter with dimension $1/population\times time$. To obtain a dimensionless formulation of System (2.1), the state variables have been divided by the carrying capacity of the environment, the parameters and the time have also been adjusted to be dimensionless as follows

and

Thus, our optimal control problem can be written in the following dimensionless form

2.1.   Theoretical results

Let us set as $\textbf{X} = (S_x, R_x, S_y, R_y)$ the vector of states, $\textbf{Z} = (z_1, z_2, z_3, z_4)$ the vector of adjoint variables and $f_0(t, \textbf{X}, \pmb \mu)$ the integrand of the cost function. We will use the Pontryagin principle for bounded controls to compute the optimal control of Problem (2.2). First of all, following the classical results given in the references ^[12,13] it is easy to verify the existence of optimal controls (it is enough to check that the properties (i)–(iii) of Section 3.4 from reference ^[14] are fulfilled). Now, the Hamiltonian $H$ associated with the Problem (2.2) is defined by $H(t, \textbf{X}(t), \mu(t), \textbf{Z}(t)) = f_0(t, \textbf{X}, \mu)+\textbf{Z}(t)\cdot f(\textbf{X}, t, \mu)$, that is

The adjoint system and state equations define the optimal system. The following theorem summarizes the main result of this section.

Theorem 2.1. There is an optimal solution $\textbf{X}^*(t)$ that minimize $J[\pmb \mu]$ in $[0, T]$. Moreover, there exits a vector of adjoint variables $\textbf{Z}$ such that

with transversality condition $z_i(t) = 0$ for $i = 1, 2, 3, 4$ which satisfies:

To see a detailed proof of the previous theorem see Appendix A.

3.   Experimental data

We used Escherichia coli strains from the Keio collection containing either $\Delta$ilvA or $\Delta$tyrA deletions ^[15], each strain was transformed with the plasmid pBGT-1 ^[16] carrying blaTEM-1 gene that confers resistance to ampicillin and a GFP fluorescent marker inducible with arabinose, pBGT-1 is a non-conjugative plasmid with around 20 copies in average per cell. Strains were grown in LB medium (Lysogeny Broth) with 40 µg/mL of kanamycin and 100 µg/mL of ampicillin (for those with plasmid) at 30℃ for 16 hours after incubation strains were washed with M9 salts and re-suspended in M9 minimal medium supplemented with glucose.

To measure susceptibility to the antimicrobial, we performed dose-response experiments in 96-well plates. The ampicillin concentrations used were 0, 6, 15, 36, 89,220,540, 1326, 3257 and 8000 µg/mL, each well was inoculated with 20 µL of clean cells for a total of 200 µL per well and incubated in an ELx808 plate reader at 30℃ with continuous shaking, OD630 was measured every 20 minutes for 24 hours, data not shown.

Co-cultures were performed in 96-well plates with 180 µL of LB media with 6 µg/mL of ampicillin, and 20 µL of cells with $80\%$ of susceptible cells ($\Delta$ilvA or $\Delta$tyrA) and $20\%$ of resistant cells ($\Delta$tyrA-pBGT or $\Delta$ilvA-pBGT). Two equal plates were grown in a ELx808 plate reader at 30℃ with eight replicates for condition in each plate, OD630 was measured every 20 minutes for 24 hours, one plate was used for growth rate analysis, and the second one was used to measure colony forming units (CFU), relative abundance and cell viability, every 2 hours we collected one of the replicates of the second plate and divided the sample, as we had to open the well, this procedure was destructive to the sample, to have sampled every 2 hours, we completely used one of the plates, the other one was kept in the plate reader.

For relative abundance and cell viability, one replicate was collected every 2 hours and supplemented with 20 µL of arabinose ($0.5\%$), incubated at 30℃ for 4 hours, and stored at 4℃ for one day to induce GFP fluorescence. After storage, 100 µL of each sample was frozen, 5 µL from a 1:1000 dilution were incubated for 2 days at $30^\circ$C in selective agar with M9 minimal media, supplemented with 10.9 mg/L isoleucine and 7.15 mg/L tyrosine and $0.5\%$ of arabinose for CFU counting, the rest of the sample was used to measured abundance of resistant cells using a CytoFLEX S cytometer (20,000 events per sample).

For experimental data, we performed 16 replicas in total, 8 were kept at a plate reader to measure OD shown in Figure 1A, the other eight replicas were used to measured abundance of cells and viability, the experimental data measured in the flux cytometer every 2 hours splitting resistance and susceptible bacteria are shown in Figure 1B.

4.   Numerical simulations

In order to performed the optimal control numerical experiments, we used the Backward-Forward Sweep Method described on Lenhart et al. ^[13] in Matlab interface, and also we used our experimental data given in Section 3 and some data collected in the reference ^[1], from two non-mating budding yeast strains of Saccharomyces cerevisiae that were designed to be deficient in the biosynthesis of an essential amino acid and also to overproduce the amino acid required by the ir partner. In such laboratory experiments, the auxotrophic strain of leucine (leu-) overproduces tryptophan, while the auxotrophic strain of tryptophan (trp-) overproduces leucine. According to reference^[1], these strains have previously been shown to form a cross-feeding mutualism when grown on solid agar, with each strain losing the amino acid needed by its partner. We assume that the carrying capacity $K = 0.1\times 10^6$ and the temporal value is $\omega = 0.3\times 10^{-5}$. The values of the other parameters taken in this study are given in Table 2.

In these experiments we consider that the amount of supplementary amino acids and the proportion of mutations are variables. The mutation portion ($q$) values used were 0.1, 0.5 and 0.9, parameters for the amount of amino acids used are shown in Table 3.

In Figure 2 we show some numerical experiments of sensitive bacteria ($X$ and $Y$) of Model (2.2), when there are no antimicrobials, for different values of $F$. From this figure, we can see that when $F$ is small enough ($F$ = 0.01) there is the extinction of both strains. When $F$ grows ($F$ = 0.05) there is cooperation or mutualism between both strains, whereas, for large values of $F$ (1.5), there is competition between both strains, being evident in the increase in the strain $X$ and the extinction of the strain $Y$.

In Figure 3, we supply antimicrobials in our system, we keep fixed the value of $F$ = 0.01 (where the is the extinction of both bacteria), and vary the values of the mutation proportion $q$ (0.1, 0.5 and 0.9). In this case, we can see that the sensitive strains tend to zero regardless of the value of $q$, whereas for large values of $q$, the resistant bacteria of both types ($X$ and $Y$) increase in density at the beginning of the time.

A different scenario occurs for $F$ = 0.05 (where there is cooperation between both strains), which can be evident in Figure 4. Here we can see that both strains become extinct only for small values of $q$ ($q$ = 0.5 or less), while if $q$ increases (0.9), they coexist in density over time.

Finally, for $F$ = 1.5 (see Figure 5), we can see that the strain $Y$ decreases no matter what value of $q$ is taken. On the other hand, the strain $X$ stabilizes at increasing densities as the value of $q$ changes.

To observe the effect of controlling the antimicrobial resistance in the system where the strains $X$ and $Y$ interact in an environment with antimicrobials, we use the Forward-Backward Sweep Method proposed by Lenhart and Workman ^[13]. The values of the parameters associated with Problem (2.2) are shown in Table 4. We only introduce the control variables for the cases where there is a high proportion of mutation ($q$ = 0.9) and for the mutualism ($F$ = 0.05) and competition ($F$ = 1.5). In Figure 6 we show the results of control when $F = 0.05$ and $q = 0.9$ (where there is cooperation). In this case, with controls, both populations of bacteria are controlled immediately, maintaining the two controls at their maximum effort during the first 40 days and then rapidly decreasing to zero during the rest of the control campaign. In Figure 8 we show the results when $F = 1.5$ and $q = 0.9$ (when there is competition). Here, it can be seen that only the strain $X$ is controlled from the first day of the campaign. But the strain $Y$ increases with the control. The effort is only made for the first control at 100%, while the second control remains at zero throughout the campaign.

5.   Parameter estimation

For the parameter estimation of Model (2.1) we used the experimental data described in Section 3. To simplify notation, we will drop the bar symbol on the parameters from here. The experimental data is shown in Figure 9. Since the data size is very low, i.e., we have just six point for each time series of type of bacteria, we have implemented a re-sample with an interpolation for the values of odd hours in order to increase the data size. Specifically, we used a rich media called LB, and low ampicillin concentration. The values are described in Table 5.

The fitting curve or estimation of the parameters of a model is considered an inverse problem. Some references of works using Bayesian inference are available in references ^{[17,18,19,20,21]}. Let $F: \mathbb{R} ^{m} \rightarrow \mathbb{R}^{s\times k}$, denoted by $F(\theta)$ be the forward problem, where $\theta$ are the parameters of Model (2.1) to estimate, $m$ the number of parameter to estimate, and $k$ the number of state variables. The inverse problem is formulated as a standard optimization problem

with $y_{\rm{obs}}$. Given $y_{\rm{obs}} = (\tilde{I}, \tilde{A})$, which correspond to strain tyrA- ($X$) and strain ilvA- ($Y$) are divided into susceptible and resistant bacteria ($X_s$, $X_r$ and $Y_s$, $Y_r$, respectively), the conditional probability distribution $\pi (\theta | y_{\rm{obs}})$, called the posterior distribution of $\theta$ is given by the Bayes' theorem:

All the available information regarding the unknown parameter $\theta$ is codified into the prior distribution $\pi (\theta)$, which specifies our belief in a parameter before observing the data. All the available information regarding the way of how was obtained the measured data is codified into the likelihood distribution $\pi (y_{\rm{obs}} | \theta)$. This likelihood can be seen as an objective or cost function, as it punishes deviations of the model from the data. A Poisson distribution, $\mathcal{P}(y| \mu)$, with respect to the time, is typically used to account for the discrete nature of these counts, where µ is the mean of the random variable $y$, i.e., $\mathbb{E}[Y] = \mu$. We assume independent Poisson distributed noise $\eta$, i.e., all dependency in the data is codified into the model (2.1). In other words, the positive definite noise covariance matrix $\eta$ is assumed to be diagonal. The posterior distribution $\pi (\theta | y_{\rm{obs}})$ given by (5.2) does not have an analytical closed form since the likelihood function, which depends on the solution of the non-linear model (2.1), does not have an explicit solution. Then, we explore the posterior distribution using the $\texttt{Stan}$ Statistics package ^[22]. We have used the Automatic Differentiation Variational Inference method (ADVI), which is based on the automatic variational inference algorithm. Specifically, we have used the Full-Rank submethod of ADVI. We have used the interface in Python (PyStan) ^[22], for more details see this Github link. Table 5 shows the posterior mean, quantiles of all estimated parameters of model (2.1). Table 6 shows the prior distributions summary of the estimated parameters of Model (2.1) using the experiment data set. Gamma distributions were used for parameters $\kappa, \gamma_x, \gamma_y$ and $q$. Uniform distribution were used for the initial conditions $S_{x_0}, R_{x_0}, S_{y_0}, R_{y_0}$. Figure 11 shows the joint probability density distributions of the estimated parameters within $95\%$ (HPD) using the experimental data. The blue lines represent the medians. Figure 12 shows Model fit for sensitive and resistant bacteria of Model (2.1) using the experiment data. Blue and red dot points represent strain ilvA ($X$) data, i.e., the susceptible and resistant bacteria, called $S_{x}$ and $R_{x}$, respectively. Orange and purple solid lines represent the median of posterior distribution of the sensitive and resistant bacteria (ilvA-pBGT ($X$), respectively. Shaded area represent the $95\%$ probability bands for the expected value of sensitive (orange line), resistant bacteria (blue line). Figure 13 shows Model fit for sensitive and resistant bacteria of the model (2.1) using the experiment data. Blue and red dot points represent strain tyrA ($Y$) data, i.e., the susceptible and resistant bacteria, called $S_{y}$ and $R_{y}$, respectively. Orange and purple solid lines represent the median of the posterior distribution of the sensitive and resistant bacteria strain tyrA-pBGT ($Y$), respectively. Shaded area represent the $95\%$ probability bands for the expected value of sensitive (orange line), resistant bacteria (blue line).

6.   Discussion

In this work, we attempted to explain the ecological interactions between cross-feeding bacteria in strains that supply an essential amino acid for their mutualistic partner when they are exposed to antimicrobials. Although from the ecological point of view, the microbial interaction can occur between multiple bacteria, for illustration and ease of handling mathematical equations, we considered the interaction between two bacteria. We formulated a mathematical model using ODEs assuming that both strains interact in an environment with different resources availability. We assumed both strains are divided into susceptible and resistant. After that, we estimated the most important parameters of the ODEs model using Bayesian Inference.

To validate our theoretical results with numerical experiments, we used our experimental data and other data available in the literature. More specifically, we used Escherichia coli strains from the Keio collection containing either $\Delta$ilvA or $\Delta$tyrA deletions, where each strain was transformed with a plasmid that confers resistance to ampicillin. For some parameters, we also used data from Hoek et al. in ^[1] of auxotrophic consortia of two microbes competing for the same resources.

The theoretical and numerical results showed that when the strains are free of the antimicrobial, depending on the amount of amino acids freely available in the environment, the strains can exhibit extinction, mutualism, or competition, the availability of resources modulates the behavior of both species. In contrast, if the strains are exposed to antimicrobials, the population dynamics depends on the proportion of bacteria that presents resistance to the antimicrobial, finding that for low levels of the resource, the two species become extinct, whereas, for high levels of the resource, competition between both strains is given.

The optimal control results showed that both strains (sensitive and resistant) are immediately controlled under cooperation, while under competition, only the density of one of the strains decreases, whereas its mutualist partner with control increases. Finally, the growth rates of both strains ($\bar \gamma_x$ and $\bar \gamma_y$), the Monod constant ($\bar \kappa$) the mutation proportion ($q$) and the initial conditions of Model (2.1) were estimated using Bayesian Inference and the data set described in Section 3. From Figures 12 and 13, we could observe that the fitting is not entirely accurate. This could be due to the lack of data and/or also because Model (2.1) should be adjusted.

The results obtained with this study corroborated that the antimicrobial resistance phenomenon is a complex problem worldwide that the scientific community has extensively studied. This problem is not only related to biological aspects of microorganisms but also other aspects, including socioeconomic and governance factors of countries ^[23]. Even, some authors claim that a novel approach to antibiotic discovery would be based on the analysis of microbial consortia in their ecological context (see ^[24,25]). Others discuss the potential of microbial interactions to target and improve microbial dysbiosis as a strategy for the prevention or treatment of cancer ^[26]. Additionally, some of them affirm that exposure to sublethal concentrations of antimicrobials can indeed alter microbial metabolism and even change the behavior in beneficial ways, triggering reactions such as fleeing or hiding within the protective environment of a microbial aggregate ^[27].

Therefore, research questions are left open. A fundamental challenge will be to generate more laboratory experiments to obtain more data to allow better adjustments to the model, and probably adjust the model assuming additional hypotheses.

Acknowledgements

J. Romero thanks the support of Fundación Ceiba, Colombia. We thank A. San Millan for the plasmid p-BGT. D. Reyes-Gonzalez is a doctoral student in Programa de Doctorado en Ciencias Biomédicas, Universidad Nacional Autónoma de México, and received fellowship 572373 from CONACYT. A. Fuentes-Hernandez was funded by PAPIIT-UNAM (grant IN215920).

Conflict of interest

No potential conflict of interest was reported by the authors.

Appendix A. Proof of Theorem 2.1

By the Pontryagin principle, we can guarantee the existence of adjoint variables $z_i$, $i = 1, 2, 3, 4$ that satisfy:

From above, the adjoint system can be written as:

By doing the respective calculations in the previous equations, we obtain System (2.4). Now, the optimality condition for the Hamiltonian is $\partial H/\partial\pmb\mu^*$, or equivalently:

From the above, we obtain the characterization given on (2.5). In consequence, $\mu_1^*$ satisfies:

or equivalently:

Similar calculations can be done for $\mu_2$, and then we obtain the characterization given on Eq (2.5) which completes the proof.