Predicting microbial growth in a mixed culture from growth curve data

We present a model-based approach for prediction of microbial growth in a mixed culture and relative fitness using data solely from growth curve experiments, which are easier to perform than competition experiments. Our approach combines growth and competition models and utilizes the total densities of mixed cultures. We implemented our approach in an open-source software package, validated it using experiments with bacteria, and demonstrated its application for estimation of relative fitness. Our approach establishes that growth in a mixed culture can be predicted using growth and competition models. It provides a way to infer relative strain or species frequencies even when competition experiments are not feasible, and to determine how differences in growth affect differences in fitness.

Abstract

Determining the fitness of specific microbial genotypes has extensive application in microbial genetics, evolution, and biotechnology. While estimates from growth curves are simple and allow high throughput, they are inaccurate and do not account for interactions between costs and benefits accruing over different parts of a growth cycle. For this reason, pairwise competition experiments are the current “gold standard” for accurate estimation of fitness. However, competition experiments require distinct markers, making them difficult to perform between isolates derived from a common ancestor or between isolates of nonmodel organisms. In addition, competition experiments require that competing strains be grown in the same environment, so they cannot be used to infer the fitness consequence of different environmental perturbations on the same genotype. Finally, competition experiments typically consider only the end-points of a period of competition so that they do not readily provide information on the growth differences that underlie competitive ability. Here, we describe a computational approach for predicting density-dependent microbial growth in a mixed culture utilizing data from monoculture and mixed-culture growth curves. We validate this approach using 2 different experiments with Escherichia coli and demonstrate its application for estimating relative fitness. Our approach provides an effective way to predict growth and infer relative fitness in mixed cultures.

Microbial fitness is usually defined in terms of the relative growth of different microbial strains or species in a mixed culture (1). Pairwise competition experiments can provide accurate estimates of relative fitness (2), but can be laborious and expensive, especially when carried out with nonmodel organisms. Moreover, competition experiments cannot be used to estimate the effect of environmental perturbations on fitness, as the competing strains must be grown in a shared environment. Instead, comparisons of separate aspects of growth curves—for example, growth rates or lag times—are commonly used to estimate the fitness of individual microbial isolates, despite clear evidence that they provide an inadequate alternative (3, 4).

Growth curves describe the density of cell populations in liquid culture over time and are usually obtained by measuring the optical density (OD) of cell populations. The simplest way to infer fitness from growth curves is to estimate the growth rate (i.e., Malthusian parameter) during the exponential growth phase, using the slope of the log of the growth curve (5) (see example in Fig. 1). While relative growth rates are often used as a proxy for relative fitness (1, 6, 7), exponential growth rates do not capture the complete dynamics of typical growth curves, such as the duration of the lag phase and the cell density achieved at stationary phase (8) (Fig. 1A). Moreover, the maximal specific growth rate is not typical for the entire growth curve (Fig. 1B). Thus, growth rates are often poor estimators of relative fitness (3, 4).

Fig. 1.

Common approach to analyzing growth curve data using an exponential model. Growth rates are commonly estimated from growth curve data by taking the log of the growth curve and performing linear regression around the time of maximum growth, t_max (see *Materials and Methods* for specific details). Implicitly, this is equivalent to fitting an exponential growth model N(t) = N₀e^rt to the growth curve. (A) The circles represent N(t), the mean cell density of 22 growth curves of strain A1 growing in monoculture (see details of experiment A). The solid line represents a fit of a polynomial of degree 3 through the points. The dashed line represents the exponential model N₀e^rt fitted to the data, with r = 0.35 and N₀ = 0.088. The dotted vertical line denotes t_max. (B) The solid curve shows dN/*dt,* the derivative of the mean density (calculated as the derivative of the solid line in A). The dotted vertical line denotes t_max. Data in this figure correspond to the growth of strain A1 (red markers in Fig. 3A1).

By contrast, competition experiments can produce estimates of relative fitness that account for all growth phases (9). In pairwise competition experiments, 2 strains—a reference strain and a strain of interest—are grown in a mixed culture. The density or frequency of each strain in the mixed culture is measured during the course of the experiment using specific markers, either by counting colonies formed by competitors that differ in drug resistance, resource utilization, or auxotrophic phenotypes (9); by monitoring fluorescent markers with flow cytometry (2); or by counting DNA barcode reads using deep sequencing (10, 11). The relative fitness of the strain of interest is then estimated from changes in the densities or frequencies of the strains during the competition experiment. Such competition experiments allow relative fitness to be inferred with high precision (2), as they directly estimate fitness from changes in strain frequencies over time.

Competition experiments are often more demanding and expensive than simple growth curve experiments, especially in laboratories where they are not routinely performed. They require the strains of interest to be engineered with genetic or phenotypic markers (see ref. 3 and references therein), which is difficult or impossible in some nonmodel organisms or when measuring the fitness effect of environmental change. Moreover, many markers incur costs that can affect the outcome of competitions (2). Therefore, many investigators prefer to use proxies for fitness, such as growth rates estimated from growth curves of monocultures. However, it is difficult to infer how differences in growth during the growth phases affect relative fitness in competition (12, 13), even when competition experiments are a plausible approach [e.g., in microbial lineages with established markers (9)].

Here, we present a computational approach that provides a framework for estimating growth parameters from growth dynamics and for predicting relative growth in mixed cultures. We provide 2 different experimental validations of this approach and demonstrate its application to estimating the effect of protein expression on relative fitness.

Results

Our approach consists of 3 stages: (i) fitting growth models to monoculture growth curve data (Fig. 2A), (ii) fitting competition models to mixed-culture growth curve data (Fig. 2B), and (iii) predicting relative growth in a mixed culture using the estimated growth and competition parameters. Independent experimental validations of this approach include the use of fluorescent Escherichia coli strains, and the use of E. coli strains that previously evolved under metabolic challenges. In both of these experimental approaches, we measured growth of 2 strains in monoculture and mixed culture, predicted growth in the mixed culture, and compared these predictions to the empirical results. Finally, we describe an application of our method to estimating the effect of lac operon expression on relative fitness.

Fig. 2.

Illustration of the computational approach. (A) Two strains (here marked with red and green circles) are grown separately in monocultures while the cell culture density (e.g., OD) is measured. Then, growth models are fitted to the measured growth curves from each experiment, providing estimates for the growth parameters of each strain: initial cell density N₀, maximum cell density K, specific growth rate r, deceleration parameter ν, initial adjustment q₀, and adjustment rate m. (B) The same 2 strains are grown in a mixed culture, and the total cell density of the mixed culture is measured. The resulting growth curve is fitted to a competition model, providing estimates for the competition parameters c₁ and c₂. Importantly, at this stage, the cells are not identified by their strain; hence, the gray circles.

Experimental Validation Design.

Fluorescence experiments.

Three fluorescence experiments (denoted A, B, and C) were performed with 2 pairs of E. coli strains marked with green and red fluorescent proteins (GFP and RFP, respectively). The same pair of strains was used in experiment A and B, and a different pair was used in experiment C. Experiment A started by diluting stationary-phase bacteria from strains 1 and 2 into fresh media, yielding cultures in which lag phase was observably longer for strain 2. In experiment B, strains 1 and 2 were pregrown in fresh media for 4 h, allowing them to reach early exponential growth phase, and then diluted into fresh media, so that there was no observable lag phase. Experiment C was similar to experiment A but with different strains, denoted 3 and 4. Each experiment consisted of 3 subexperiments: 32 replicate monocultures of the GFP strain; 30 replicate monocultures of the RFP strain; and 32 replicate mixed cultures containing the GFP and RFP strains together. These subexperiments were performed under the same experimental conditions in a single 96-well plate format. The OD of every well (i.e., in all subexperiments, both in monoculture and mixed-culture wells) was measured using an automatic plate reader (Figs. 3 and 4). In addition, samples were collected from the mixed-culture subexperiment wells and the relative frequencies of the 2 strains were measured by flow cytometry (Fig. 5). See Materials and Methods for additional details (14).

Fig. 3.

Fitting growth models to monoculture growth curves. The panels show monoculture growth curve data (markers) and best-fit growth models (lines; Eqs. 2a and 2b). In panel labels, letters denote the experiment (A, B, and C) and numbers denote the strain (1–4; red strains in *Top* row, green strains in *Bottom* row); see *Results*, *Experimental Validation Design*. See Table 1 for estimated growth parameters. Thirty to thirty-two replicates per strain. The figure omits the stationary phase.

Fig. 4.

Fitting competition models to mixed-culture growth curves. The panels show mixed-culture growth curves (blue symbols show measured total densities) and best-fit competition models (solid blue lines; Eqs. 3a and 3b) from experiments A–C (i.e., from the mixed growth subexperiments). Thirty-two replicates per experiment. The dashed black lines show the prediction of an exponential model with $N_{1} (t) + N_{2} (t) = N_{0} (e^{r_{1} t} + e^{r_{2} t})$ , where $r_{i}$ are estimated by fitting an exponential model to monoculture growth curves (Fig. 1).

Fig. 5.

Predicting competition results. Comparison of experimental data (markers) and model prediction (dashed lines; see *SI Appendix*, Fig. S3 for confidence intervals) of relative strain frequencies in a mixed culture from experiments *A–C*. The red and green dashed lines show our model predictions for the red and green strains (red strains eventually win in all experiments). The dashed black lines show exponential model predictions ( $f_{1} (t) = 1 / (1 + (N_{2} (0) / N_{1} (0)) e^{- (r_{1} - r_{2}) t})$ and $f_{2} (t) = 1 - f_{1} (t)$ ; Fig. 1). Error bars show SD (hardly seen in A and C).

LacI experiments.

Eight E. coli strains were isolated from populations that previously evolved in lactose-containing environments (15). These strains maintained the ancestral allele (rather than having fixed a loss-of-function mutation) at the lacI gene, which represses the lac operon. These strains were then mutated at the lacI gene. For each pair of lacI+ and lacI− strains, growth curves were measured in a monoculture, and in competition experiments conducted in mixed culture (Fig. 6). See Materials and Methods for additional details.

Fig. 6.

Predicting relative fitness. Comparison of experimental and model estimates of the fitness (Eq. 4) of *lacI*− mutants relative to *lacI+* wild types for 8 different *E. coli* strains (strain identifiers on x axis). The open circles (○) show experimental estimates. The filled circles (●) show model predictions based only on monoculture growth curve data. Errors bars show 95% confidence intervals. Due to noise in the growth curve data, model fitting was unstable, and some predicted competitions resulted in no growth, or even a decrease in OD, for one strain. Such outliers were excluded from the analysis.

Estimating Growth Parameters.

Growth model.

The Baranyi–Roberts model (16) can be used to model growth that comprises several phases: lag phase, exponential phase, deceleration phase, and stationary phase (5). The model assumes that cell growth accelerates as cells adjust to new growth conditions, then decelerates as resources become scarce, and finally halts when resources are depleted (17). The model is described by the following ordinary differential equation (see equations 1c, 3a, and 5a in ref. 16; for a derivation of Eq. 1 and for further details, see Appendix A):

\frac{d N}{d t} = r \cdot α (t) \cdot N (1 - {(\frac{N}{K})}^{v}),

[1]

where

t

is time,

N = N (t)

is the cell density at time

t

,

r

is the specific growth rate at low density,

K

is the maximum cell density,

v

is a deceleration parameter, and

α (t)

is the adjustment function

(0 \leq α (t) \leq 1)

, which describes the fraction of the population that has adjusted to the new growth conditions by time

t

. In microbial experiments, an overnight liquid culture of microorganisms that has reached stationary phase is typically diluted into fresh media. Following dilution, cells enter lag phase until they adjust to the new growth conditions. We chose the specific adjustment function suggested by Baranyi and Roberts (16), which is both computationally convenient and biologically interpretable:

α (t) = q_{0} / q_{0} + e^{- m t}

, where

q_{0}

characterizes the physiological state of the initial population, and

m

is the rate at which the physiological state adjusts to the new growth conditions.

The Baranyi–Roberts differential equation (Eq. 1) has a closed-form solution:

N (t) = {K / [1 - (1 - {(\frac{K}{N_{0}})}^{v}) e^{- r \cdot v \cdot A (t)}]}^{\frac{1}{v}},

[2a]

where

A (t) = \int_{0}^{t} α (w) d w = t + \frac{1}{m} \log (\frac{e^{- m \cdot t} + q_{0}}{1 + q_{0}}),

[2b]

with parameters as defined above and with

N_{0} = N (0)

, the initial population density. For a derivation of Eqs. 2a and 2b from Eq. 1, see Appendix A. We used this growth model (Eqs. 2a and 2b) to estimate growth parameters, which we then used in a competition model (Eqs. 3a and 3b below) to infer relative growth in a mixed culture. Note that alternative models could be used with our approach, for example when analyzing biphasic growth curves (18).

Model fitting.

Growth model parameters were estimated by fitting the growth model (Eqs. 2a and 2b) to the monoculture growth curve data of each strain (Fig. 1A). The best-fit models (lines) and experimental data (markers) are shown in Fig. 3; see SI Appendix, Table S1 for the estimated growth parameters. From these best-fit models, we also estimated the maximum specific growth rate

(\max ((1 / N) \cdot (d N / d t)))

, the minimal specific doubling time (minimal time required for cell density to double), and the lag duration (time required to adjust to new environment); see SI Appendix, Table S1 and Materials and Methods. Different strains differ in their growth parameters: for example, strain A1 (red strain in experiment A) grows 41% faster than the strain A2 (green), has 23% higher maximum density, and has a 60% shorter lag phase (Fig. 3).

Estimating Competition Coefficients.

Competition model.

To model growth in a mixed culture, we assume that interactions between strains are density dependent, for example due to resource competition. This excludes frequency-dependent interactions, which may arise due to production of toxins (19) or public goods (20) (see Fig. 8 for a deviation from this assumption). Therefore, all interactions are described by the deceleration of the growth rate of each strain in response to the increased density of both strains. We have derived a 2-strain Lotka–Volterra competition model (21) based on resource consumption (Appendix B):

{\begin{matrix} \frac{d N_{1}}{d t} = r_{1} \frac{q_{0, 1}}{q_{0, 1} + e^{- m_{1} t}} N_{1} (1 - \frac{N_{1}^{v_{1}}}{K_{1}^{v_{1}}} - c_{2} \cdot \frac{N_{2}^{v_{2}}}{K_{1}^{v_{1}}}) \\ \frac{d N_{2}}{d t} = r_{2} \frac{q_{0,2}}{q_{0,2} + e^{- m_{2} t}} N_{2} (1 - c_{1} \cdot \frac{N_{1}^{v_{1}}}{K_{2}^{v_{2}}} - \frac{N_{2}^{v_{2}}}{K_{2}^{v_{2}}}), \end{matrix}

[3a, 3b]

where

N_{i}

is the density of strain

i = 1, 2

,

r_{i}, K_{i}, v_{i}, q_{0, i},

and

m_{i}

are the values of the corresponding parameters for strain

i

, obtained from fitting the growth model (Eqs. 2a and 2b) to monoculture growth curve data, and

c_{i}

are competition coefficients, the ratios between interstrain and intrastrain competitive effects. Note that in resource competition, each strain can be limited by a different resource, and strains may vary in their efficiency of resource utilization (i.e., uptake and conversion rates; Appendix A).

Model fitting.

The competition model (Eqs. 3a and 3b) was fitted to growth curve data from the mixed culture, in which the total OD of both strains in mixed culture was recorded over time (i.e., the bulk density, not the frequency or density of individual strains; Fig. 2B). The growth parameters

(r_{i}, K_{i}, v_{i}, q_{0, i})

were fixed to the values estimated from the monoculture growth curves at the previous stage, and the fitting at this stage only provided estimates for the competition coefficients

c_{i}

. Fitting was performed by minimizing the squared differences between

N_{1} + N_{2}

(the sum of the solutions of Eqs. 3a and 3b; integrals solved numerically using LSODA solver) and the total OD from the mixed culture (Fig. 4). Part of the strength of this approach stems from its use of measurements of the total density of mixed cultures, which is usually ignored when estimating fitness from growth curves (5). However, when such measurements are not available, competition coefficients can be set to

c_{i} = 1

. This is the case for the lacI experiments. See Materials and Methods for additional details.

Prediction and Validation.

Model prediction.

We solved the competition model (Eqs. 3a and 3b) using estimates of all of the competition model parameters, and numerical integration (LSODA solver), thereby providing a prediction for the cell densities

N_{1} (t)

and

N_{2} (t)

of the 2 strains growing in a mixed culture. From these predicted densities, the relative frequencies of each strain over time were estimated as

f_{i} (t) = (N_{i} (t)) / (N_{1} (t) + N_{2} (t))

.

Experimental validation: Relative growth.

We compared the model predictions

f_{i} (t)

to experimental relative frequencies obtained using flow cytometry from mixed-culture samples in the fluorescent experiments. Experimental results and model predictions are shown in Fig. 5, together with the exponential model predictions (which effectively use the growth rate as a proxy for fitness). Our model performed well and clearly improved upon the exponential model for predicting competition dynamics in a mixed culture.

Experimental validation: Relative fitness.

We validated the use of this approach for estimation of relative fitness using 8 pairs of lacI+ and lacI− strains. In each pair, the lacI+ strain had previously evolved in a lactose-containing environment (15) and the lacI− strain was produced by introducing a mutation that causes lac genes to be constitutively expressed at a high level. For each lacI+/lacI− strain pair, growth curves were measured in a monoculture and used to predict growth in mixed culture (Eqs. 3a and 3b; competition parameters were set to

c_{i} = 1

rather than estimated from growth in a mixed culture). In addition, competition experiments were conducted in mixed culture. The relative fitness W of lacI− strain relative to lacI+ strain was estimated from the experimental and predicted densities N_lacI− and N _lacI+ following Lenski et al. (7), where N_i(t) is the density of strain i after t hours, and

W = \log (\frac{N_{l a c I -} (24)}{N_{l a c I -} (0)}) / \log (\frac{N_{l a c I +} (24)}{N_{l a c I +} (0)}) .

[4]

Fig. 6 compares the relative fitness W of lacI− mutants from competition experiments and from model predictions in 8 strain pairs. Clearly, estimates from experiments and from model predictions are very similar. This suggests it was reasonable to assume competitions parameters can be fixed at

c_{i} = 1

.

Application: Predicting the Effect of lac Operon Expression on Relative Fitness.

We next tested an application of our computational approach by estimating the effect of lac operon expression on relative fitness in a strain of E. coli that evolved in a distinct lactose-containing environment (15). Quantitative manipulation of lac expression level is done by changing either the genotype or the environmental concentration of an inducer. With the latter, it is not possible to perform direct fitness competitions between strains expressing the lac operon at different levels.

To estimate the effect of lac operon expression on relative fitness, growth curves of a lacI+ strain, namely strain GL2, were measured in monoculture at a range of concentrations of isopropyl-β-d-thiogalactoside (IPTG), a molecular analog of allolactose that induces the lac operon (experimental conditions were similar to the lacI experiments used for experimental validation; see above). The growth medium contained glycerol as a sole carbon source so that expression of the lac operon was expected to confer a fitness cost. The effect of each IPTG concentration on lac expression was determined directly using Miller assays. The monoculture growth curves were used to predict growth in a mixed culture and then to estimate the relative fitness W of cells growing with each level of IPTG, and therefore lac expression, relative to cells growing without IPTG that did not express the lac operon (Fig. 7).

Fig. 7.

Predicting relative fitness of *lac* operon expression. Expression of the *lac* operon is induced by IPTG concentration. Thus, competing strains that grow in a shared environment cannot have different expression levels. Our approach was therefore applied to growth curves from monocultures with and without IPTG to predict the fitness W of cells expressing *lac* relative to cells that do not express *lac* (Eq. 4). Symbols denote the estimated relative fitness (y axis) for different expression levels (x axis, relative to a reference ancestor strain). Error bars show SEs.

Discussion

We developed a computational approach to predict relative growth in a mixed culture from growth curves of mono- and mixed cultures (Fig. 2). This approach removes the need to measure the frequencies of single isolates within a mixed culture. The approach performed well in 2 different experimental setups (Figs. 5 and 6), with results far superior to the current approach commonly used (3, 5). The 2 experimental validations provide strong support for the idea that our computational approach provides a way to simplify and reduce the cost of analyzing relative fitness. Indeed, this approach has already been used to estimate relative fitness of an E. coli strain in which the arginine codons CGU and CGC were edited to CGG in 60 highly expressed genes (22).

Our approach assumes that the assayed strains will grow in accordance with the density-dependent growth and competition models, which are appropriate when growth depends on the availability of a limiting resource (Appendix A and Appendix B). Therefore, this approach can be applied to data from a variety of organisms, experiments, and conditions. However, our approach is not appropriate if growth is frequency dependent, for example due to production of public goods (23–25) and toxins (19) or due to cross-feeding (26). Fig. 8 demonstrates the applicability of our model to simulated experimental results from 4 different frequency-dependent dynamics (1). When density- and frequency-dependent interactions work in the same direction, e.g., due to exploitation of the slow-growing strain (green) by the fast-growing strain (red), our approach is consistent with the simulated experiments: The competition model fits the total density in mixed culture quite well (Fig. 8A), and its mixed growth prediction is consistent with the final outcome after 10 h, but not with the full frequency trajectories (Fig. 8E). However, this is not the case when density- and frequency-dependent interactions do not agree so that the slow-growing strain benefits from the presence of the fast-growing strain, e.g., due to mutualism, competition, or exploitation by the slow-growing strain. In these cases, the fit of our competition model to total density in a mixed culture is poor (Fig. 8 B–D), and the model can fail to predict even the final outcome of pairwise competition (Fig. 8H). Future work will determine whether such divergences between experimental results and model predictions could be used to detect frequency-dependent interactions.

Fig. 8.

Frequency-dependent growth. Results from simulated experiments, which include frequency-dependent interactions, analyzed using our competition model, which assumes density-dependent growth. *Top* row shows the fit of competition models (solid lines) to total densities from mixed growth (markers), as in Fig. 4. *Bottom* row shows actual (markers) and predicted (dashed lines) relative strain frequencies, as in Fig. 5. Each column corresponds to a different type of frequency dependence. Growth of strain i in the presence of strain j follows $d N_{i} / dt = N_{i} r_{i} (1 - ((N_{i}^{ν_{i}} + c_{j} N_{j}^{ν_{j}}) / K_{i}^{ν_{i}})) + γ_{i} (N_{i} N_{j} / (N_{i} + N_{j}))$ , with $γ_{1} = γ_{2} > 0$ for mutualism (A and E); $- γ_{1} = γ_{2} > 0$ for exploitation of the green strain by the red strain (B and F); $γ_{1} = γ_{2} < 0$ for competition (C and G); and $γ_{1} = - γ_{2} > 0$ for exploitation of the red strain by the green strain (D and H); growth parameters as estimated from experiment A; competition coefficients $c_{1} = c_{2} = 1$ . Ten replications per pair of $γ_{1}$ and $γ_{2}$ .

Growth curve experiments, in which only OD is measured, require less effort and fewer resources than pairwise competition experiments, in which the cell frequency or count of each strain must be determined (2, 3, 9, 27). Current approaches to estimate fitness from growth curves only incorporate measurements from monoculture experiments. In contrast, our approach infers actual competition by directly incorporating measurements from mixed-culture experiments. Moreover, current approaches mostly use the growth rate and/or the maximum population density as a proxy for fitness (5), but proxies for fitness based on a single growth parameter cannot capture the full scope of effects that contribute to differences in overall fitness (13, 28). Most obviously, they fail to account for the lag and deceleration phases of growth. In contrast, our approach integrates several growth phases, allowing more accurate estimation of relative growth and fitness from growth curve data. Different growth phases also can be integrated into a single parameter by measuring or calculating the area under the curve (AUC) for the monoculture growth curves (29) (Fig. 3). This approach is easy to understand and to implement, and the AUC seems to correlate with both the growth rate and the maximum density (29). However, the biological interpretation of the AUC, how it is affected by the different growth parameters, and how it affects relative fitness and competition results, is unclear.

Our approach is useful even for laboratories that have considerable experience performing competition experiments. First, it can predict the results of hypothetical competition experiments. We demonstrated this by measuring growth of E. coli strains at different concentrations of IPTG, an inducer of the lac operon. We used our computational approach to predict how 2 populations of this strain would grow, if it were possible for them to compete in a mixed culture while keeping their IPTG exposures different. We then used this prediction to estimate the effect of protein expression on relative fitness (Fig. 7). We suggest that our approach can be similarly applied to predict the relative growth of strains experiencing different drug or nutrient concentrations. Second, it can be used to predict mixed growth, even if it is very hard or impossible to insert phenotypic or genetic markers into the strains in question, e.g., with some nonmodel organisms. Third, even when competition experiments can be performed, they are rarely designed in a way that gives insight into how differences in growth underlie differences in fitness (12, 13): Our approach can highlight whether strain 1 is more fit than strain 2 due to faster growth rate, or due to a shorter lag phase, for example. By inferring relative fitness from growth parameters, this approach sheds light on the source of differences in fitness. Furthermore, one can change specific growth parameters and simulate competition, thereby predicting the adaptive potential of such changes.

Another interesting approach to relating differences in growth during different growth phases to fitness has recently been described by Li et al. (11), who assumed that if a strain grows faster in a specific growth phase, prolonging that phase while keeping other phases fixed will increase the strain’s relative fitness. “Fitness profiles”—measurements of relative fitness with systematically varied growth phase durations—were characterized and used to find the underlying cause of fitness gain in strains that previously evolved in a glucose-limited environment. While the fitness profiles approach is very promising, it is also very labor intensive and expensive compared with ours.

We have released Curveball, an open-source software package that implements our approach (http://curveball.yoavram.com). This software is written in Python (30), an open-source and free programming language, and includes a user interface that does not require prior knowledge of programming. Curveball has already been used successfully to estimate relative fitness in E. coli (22). It is free and open (i.e., libre and gratis), so that additional data formats, growth and competition models, and other analyses can be added by the community to extend its utility.

Conclusions

We developed and tested an approach to analyzing growth curve data and applied it to predict the relative growth and fitness of individual strains within a mixed culture. Competitive fitness is defined as the relative change in frequency during growth in mixed culture. Therefore, any process that affects relative growth in a mixed culture might affect competitive fitness. Current approaches use growth curve experiments because they are easy to obtain, despite their clear deficiencies. Our approach allows the use of such growth curve data, incorporating growth curves measured in a mixed culture, and thus incorporates various processes that occur in a mixed culture, including actual competition dynamics. By predicting growth in mixed culture and estimating competitive fitness, our approach can improve the understanding of competitive fitness in microbes.

Materials and Methods

Strains and Plasmids.

Fluorescent experiments.

E. coli strains 1 and 2, used in both experiment A and B, were DH5α-GFP (J.B. Laboratory, Tel Aviv University, Tel Aviv, Israel) and TG1-RFP (E. Ron Laboratory, Tel Aviv University, Tel Aviv, Israel), respectively; E. coli strains 3 and 4, used in experiment C, were JM109-GFP (N. Ohad Laboratory, Tel Aviv University, Tel Aviv, Israel) and K12 MG1655-Δfnr-RFP (E. Ron Laboratory, Tel Aviv University, Tel Aviv, Israel), respectively. GFP or RFP genes were introduced using plasmids that also included genes conferring resistance to kanamycin (Kan^R) and chloramphenicol (Cap^R) [R. Milo Laboratory, Weizmann Institute of Science, Rehovot, Israel (31)]. lacI experiment: E. coli strains were selected from populations previously evolved by Cooper and Lenski (15).

Media.

Fluorescent experiment.

Experiments were performed in LB media [5 g/L Bacto yeast extract (BD, 212750), 10 g/L Bacto tryptone (BD; 211705), 10 g/L NaCl (Bio-Lab; 190305), and 1 L DDW] with 30 μg/mL kanamycin (Caisson Labs; K003) and 34 μg/mL chloramphenicol (Duchefa Biochemie; C0113). Green or red fluorescence of each strain was confirmed by fluorescence microscopy (Nikon Eclipe Ti; SI Appendix, Fig. S1). LacI experiments: Experiments with lacI strains were performed in DM (Davis–Mingioli minimal broth) with 0.021% lactose (Fig. 6) or 0.2% glycerol (Fig. 7).

Growth and Competition Experiments.

Fluorescent experiments.

Strains were inoculated into 3 mL of LB+Cap+Kan and grown overnight with shaking. Saturated overnight cultures were diluted into fresh media so that the initial OD was detectable above the OD of media alone (1:1–1:20 dilution rate). In experiment B, to avoid a lag phase, cultures were pregrown until the exponential growth phase was reached as determined by OD measurements (3–6 h). Cells were then inoculated into 100 μL LB+Cap+Kan in a 96-well flat-bottom microplate (Costar) in 3 subexperiments: 32 wells contained a monoculture of the GFP-labeled strain; 30 wells contained a monoculture of the RFP-labeled strain; 32 wells containing a mixed culture of both GFP and RFP-labeled strains. Two wells contained only growth medium.

The cultures were grown at 30 °C in an automatic microplate reader (Tecan Infinite F200 Pro), shaking at 886.9 rpm, until they reached stationary phase. OD₅₉₅ readings were taken every 15 min with continuous shaking between readings.

Samples were collected from the incubated microplate at the beginning of the experiment and once an hour for 6–8 h: 1–10 µL were removed from 4 wells (different wells for each sample), and diluted into cold PBS buffer (DPBS with calcium and magnesium; Biological Industries; 02-020-1). These samples were analyzed with a fluorescent cell sorter (Miltenyi Biotec; MACSQuant VYB). GFP was detected using a 488-nm/520(50)-nm FITC laser. RFP was detected with a 561-nm/615(20)-nm dsRed laser. Samples were diluted further to eliminate “double” event (events detected as both “green” and “red” due to high cell density) and noise in the cell sorter (2). LacI experiments: Strains were inoculated into 1 mL of LB media and grown overnight. Saturated overnight cultures were diluted and preconditioned to the DM media supplemented with lactose or glycerol by transferring 1 μL into 1 mL of said growth media and incubating for 24 h. The next day, 2 μL of the preconditioned culture was transferred into 89 μL of the same media, with variable IPTG concentrations, in a 96-well microplate. The microplate then incubated in a microplate reader (VersaMax) at 37 °C until cells reached stationary phase. OD₄₅₀ readings were taken every 5 min.

Data Analysis.

Fluorescent experiments.

Fluorescent cell sorter output data were analyzed using R (32) with the flowPeaks package that implements an unsupervised flow cytometry clustering algorithm (33). Growth curve data were analyzed using Curveball, an open-source software written in Python (30) that implements the approach presented in this manuscript. The software includes both a programmatic interface (API) and a command line interface (CLI), and therefore does not require programming skills. The source code makes use of several Python packages: NumPy (34), SciPy (35), Matplotlib (36), Pandas (37), Seaborn (38), LMFIT (39), Scikit-learn (40), and SymPy (41). LacI experiments: Growth curves of the lacI strains were analyzed using the same models but different software implementation. We note specific differences in the analysis wherever these apply.

Fitting Growth Models.

To fit growth models (Eqs. 2a and 2b) to monoculture density data, we used the least-squares nonlinear curve fitting procedure in SciPy’s least_squares function (35). We then calculate the Bayesian information criteria (BIC) of several nested models, defined by fixing the value of specific growth parameters (Appendix A and SI Appendix, Table S1 and Fig. S2). BIC is given by the following:

BIC = n \cdot \log (\frac{\sum_{i = 1}^{n} {(N (t_{i}) - \hat{N} (t_{i}))}^{2}}{n}) + k \cdot \log n,

where

k

is the number of model parameters,

n

is the number of data points,

t_{i}

are the time points,

N (t_{i})

is the OD at time point

t_{i}

, and

\hat{N} (t_{i})

is the expected density at time point

t_{i}

according to the model. We selected the model with the lowest BIC (42, 43). Other metrics for model selection can be used, but BIC was chosen for its simplicity and flexibility. Lag duration, max specific growth rate, and min doubling time (Table 1) were estimated from the fitted growth model,

N (t)

(Eqs. 2a and 2b). Lag duration is the time at which the line tangent to

N (t)

at the point of maximum derivative (i.e., the inflection point) intersects with N₀, the initial population size (44). The maximum specific (i.e., per-capita) growth rate is

\max_{t} ((1 / N) \cdot (d N / d t))

; specific growth rate in useful as a metric to compare different strains or treatments as it does not depend on the population density. Min doubling time is the minimal time required to double the population density,

N (t)

. LacI experiments: Model selection was not performed. Rather, we fitted the growth models (Eqs. 2a and 2b) but assumed that the rate at which the physiological state adjusts to the new growth conditions is equal to the specific growth rate at low density ((

m = r

)), to achieve more stable model fitting, as suggested by Baranyi (45). This was not necessary for the fluorescence experiment in which model fitting was stable.

Table 1.

Estimated growth parameters

Strain parameter	A1 (red)	A2 (green)	B1 (red)	B2 (green)	C1 (red)	C2 (green)
	Experiment A		Experiment B		Experiment C
Initial density $(N_{0})$	0.124	0.125	0.23	0.286	0.204	0.188
Max density $(K)$	0.650 (0.643, 0.658)	0.528 (0.525, 0.532)	0.628 (0.624, 0.632)	0.619 (0.612, 0.625)	0.741 (0.735, 0.746)	0.633 (0.627, 0.638)
Max specific growth rate	0.376 (0.371, 0.382)	0.268 (0.262, 0.275)	0.369 (0.355, 0.384)	0.256 (0.251, 0.261)	0.420 (0.391, 0.426)	0.228 (0.226, 0.231)
Min doubling time	1.844 (1.809, 1.88)	2.695 (2.636, 2.77)	2.451 (2.397, 2.506)	4.372 (4.269, 4.481)	2.075 (2.035, 2.124)	3.117 (3.087, 3.147)
Lag duration	1.578 (1.513, 1.64)	3.930 (3.82, 4.028)	0.014 (0.002, 0.029)	0.004 (0.002, 0.013)	0.039 (0.033, 0.081)	0.711 (0.684, 0.749)

Parentheses provide 95% confidence intervals (bootstrap, 1,000 samples). Min doubling time is the minimal time required to double the population density. Densities are in OD₅₉₅; growth rate is in hours⁻¹; doubling time and lag duration are in hours. See SI Appendix, Table S2, for additional parameter estimates.

Fitting Exponential Models.

The following represents a common approach to estimating growth rates from growth curve data and was used as a benchmark for our approach (Fig. 2 and black dashed lines in Fig. 4 and Fig. 5). A polynomial p(t) is fitted to the mean of the growth curve data N(t). The time of maximum growth rate t_max is found by differentiating the fitted polynomial and finding the time at which the maximum of the derivative max(p(t)) occurs. Then, a linear function at+b is fitted to the log of the growth curve logN(t) in the neighborhood of t_max (e.g., at 5 surrounding time points). The parameters a and b are then interpreted as the growth rate r = a and the log of the initial population density b = logN₀.

Fitting Competition Models.

To fit competition models (Eqs. 3a and 3b), we used the Nelder–Mead simplex method (also called downhill simplex method) from SciPy’s minimize function (35) to find the competition parameters

c_{i}

that minimize the difference between

N_{1} + N_{2}

(Eqs. 3a and 3b) and the total OD of mixed cultures. Other model parameters were fixed to the values estimated from monoculture growth curves.

N_{1}

and

N_{2}

were calculated using numerical integration of Eqs. 3a and 3b with SciPy’s odeint function (35). LacI experiments: To estimate the effect of lac operon expression on relative fitness, strains must grow in the presence of different IPTG concentrations, which is impossible in a mixed culture. Therefore, we did not perform mixed-culture experiments, and competition parameters were set to

c_{i} = 1

rather than estimated from growth in a mixed culture.

Data Availability

Data have been deposited on Figshare (DOI: 10.6084/m9.figshare.3485984.v1).

Code Availability.

Source code is available at https://github.com/yoavram/curveball; an installation guide, tutorial, and documentation are available at http://curveball.yoavram.com.

Figure Reproduction.

Data were analyzed and Figs. 1–6 were produced using a Jupyter notebook (46) that is available as a supporting file and at https://github.com/yoavram/curveball_ms.

Appendices

Appendix A: Mono-Culture Model.

We derive our growth models from a resource consumption perspective (21, 47). We denote by R the density of a limiting resource, and by N the density of the cell population, both in total mass per unit of volume.

We assume that the culture is well-mixed and homogeneous and that the resource is depleted by the growing cell population without being replenished. Therefore, the intake of resources occurs when cells meet resource via a mass action law with resource uptake rate h. Once inside the cell, resources are converted to cell mass at a conversion rate of ε

ε

. Cell growth is assumed to be proportional to

R \cdot N

, whereas resource intake is proportional to a power of cell density,

R \cdot N^{v}

. We set

Y ≔ N^{v}

.

We can describe this process with differential equations for R and N:

{\begin{cases} \frac{d R}{d t} = - h R N^{v} \\ \frac{d N}{d t} = ɛ h R N . \end{cases}

[A1a, A1b]

These equations can be converted to equations in R and Y:

\begin{array}{l} Y = N^{v} \Rightarrow \\ \frac{d Y}{d t} = v N^{v - 1} \frac{d N}{d t} = \\ v N^{v - 1} \cdot ɛ h RN = v ɛ h R N^{v}, \end{array}

which yields

{\begin{cases} \frac{d R}{d t} = - h R Y \\ \frac{d Y}{d t} = μ h R Y, \end{cases}

[A2a, A2b]

with

μ = v

To solve this system, we use a conservation law approach by setting

M = μ R + Y

(48). We find that M is constant

\frac{d M}{d t} = μ \frac{d R}{d t} + \frac{d Y}{d t} \equiv 0,

and we can substitute

μ R = M - Y

in Eq. A2b to get

\frac{d Y}{d t} = h Y (M - Y) = h M Y (1 - \frac{Y}{M})

[A3]

Substituting again

N^{v} = Y, \frac{d Y}{d t} = v N^{v - 1} \frac{d N}{d t},

and defining

K = M^{\frac{1}{v}}, r = \frac{h}{v} K^{v}

, we get

\frac{d N}{d t} = r \cdot N \cdot (1 - {(\frac{N}{K})}^{v}),

[A4]

which is the Richards differential equation (49), with the maximum population density K and the specific growth rate at low density r.

We solve Eq. A4 via Eq. A3, which is a logistic equation and therefore has a known solution. Setting the initial cell density

N (0) = N_{0}

we have

N (t) = \frac{K}{{(1 - {(1 - (\frac{K}{N_{0}}))}^{v} e^{- r v t})}^{\frac{1}{v}}} .

Eq. A4 is an autonomous differential equation (

d N / d t

doesn't explicitly depend on t). To include a lag phase, Baranyi and Roberts (16) suggested adding an adjustment function

α (t)

, which makes the equation nonautonomous (explicitly dependent on t):

\frac{d N}{d t} = r \cdot α (t) \cdot N \cdot (1 - {(\frac{N}{K})}^{v})

[A5]

Baranyi and Roberts suggested a Michaelis–Menten type of function (45)

α (t) = \frac{q_{0}}{q_{0} + e^{- m t}},

which has 2 parameters: q₀ is the initial physiological state of the population, and m is the rate at which the physiological state adjusts to growth conditions. Integrating

α (t)

gives

A (t) ≔ \int_{0}^{t} α (s) d s = \int_{0}^{t} \frac{q_{0}}{q_{0} + e^{- m s}} d s = t + \frac{1}{m} \log (\frac{e^{-mt} + q_{0}}{1 + q_{0}}) .

Therefore, integrating Eq. A5 produces Eqs. 2.

The term

1 - {(N / K)}^{v}

in Eq. A5 is used to describe the deceleration in the growth of the population as it approaches the maximum density K. When

v = 1

, the deceleration is the same as in the standard logistic model

(\frac{d N}{d t} = r \cdot N \cdot (1 - \frac{N}{K}))

and the density at the time of the maximum population growth

(\frac{d^{2} N}{d t^{2}} (t) = 0)

is half the maximum density,

\frac{K}{2}

. When

v > 1

or

1 > v

, the deceleration is slower or faster, respectively, and the density at the time of the maximum growth rate is

{K / (1 + v)}^{1 / v}

(see ref. 49, who uses different variables:

W = N, A = K, m = v + 1, k = r \cdot v

).

We use 6 forms of the Baranyi–Roberts model (see SI Appendix, Fig. S2 and Table S1). The full model is described by Eqs. 2 and has 6 parameters. A five-parameter form of the model assumes

v = 1

, such that the curve is symmetric as in the standard logistic model, but still incorporates the adjustment function

α (t)

and therefore includes a lag phase. Another five-parameter form has both rate parameters set to the same value

(m = r)

, which was suggested in order to make the fitting procedure more stable (45, 50). A four-parameter form has both of the previous constraints, setting

m = r

and

v = 1

(45). Another four-parameter form of the model has no lag phase, with

1 / m = 0 \Rightarrow α (t) \equiv 1

, which yields the Richards model (49), also called the θ-logistic model (51), or the generalized logistic model. This form of the model is useful in cases where there is no observed lag phase: Either because the population adjusts very rapidly or because it was already adjusted prior to the growth experiment, possibly by pregrowing it in fresh media before the beginning of the experiment. The last form is the standard logistic model (52), in which

v = 1

and

1 / m = 0

.

Appendix B: Mixed Culture Model.

We consider the case in which 2 species or strains grow in the same culture, competing for a single limiting resource, similarly to Eq. A1:

{\begin{cases} \frac{d R}{d t} = - h_{1} R N_{1}^{v_{1}} - h_{2} R N_{2}^{v_{2}} \\ \frac{d N_{1}}{d t} = ɛ_{1} h_{1} R N_{1} \\ \frac{d N_{2}}{d t} = ɛ_{2} h_{2} R N_{2} \end{cases}

[B1a–B1c]

We define

Y_{i} = N_{i}^{v_{i}}

, and

M_{i} = ɛ_{i} v_{i} R + Y_{i} + \frac{ɛ_{i} v_{i}}{ɛ_{j} v_{j}} Y_{j}

(where j is 1 when i is 2 and vice versa) and find that

\frac{d M_{i}}{d t} \equiv 0

and

M_{i}

is constant. We then substitute

ɛ_{i} v_{i} R = M_{i} - Y_{i} - \frac{ɛ_{i} v_{i}}{ɛ_{j} v_{j}} Y_{j}

into the differential equations for

\frac{d Y_{i}}{d t}

. Denoting

K_{i} = M_{i}^{\frac{1}{v_{i}}}

and

r_{i} = \frac{h_{i}}{v_{i}} K_{i}^{v_{i}}

, we get

{\begin{cases} \frac{d N_{1}}{d t} = r_{1} N_{1} (1 - \frac{N_{1}^{v_{1}}}{K_{1}^{v_{1}}} - c_{2} \cdot \frac{N_{2}^{v_{2}}}{K_{1}^{v_{1}}}) \\ \frac{d N_{2}}{d t} = r_{2} N_{2} (1 - c_{1} \frac{N_{1}^{v_{1}}}{K_{2}^{v_{2}}} - \frac{N_{2}^{v_{2}}}{K_{2}^{v_{2}}}), \end{cases}

[B2a, B2b]

where

C_{j} = \frac{ɛ_{i} v_{i}}{ɛ_{j} v_{j}}

. To get Eqs. 3 from Eqs. B2, we include a lag phase by adding the adjustment function

α_{i} (t) = \frac{q_{0, i}}{q_{0, i} + e^{- m_{i} t}}

; see details in A5 in Appendix A.

We get a similar result if the strains are limited by 2 resources R₁ and R₂ that both strains consume:

{\begin{cases} \frac{d R_{1}}{d t} = - h_{1} R_{1} N_{1}^{v_{1}} - h_{2} R_{1} N_{2}^{v_{2}} \\ \frac{d R_{2}}{d t} = - h_{1} R_{2} N_{1}^{v_{1}} - h_{2} R_{2} N_{2}^{v_{2}} \\ \frac{d N_{1}}{d t} = ɛ_{1} h_{1} R_{1} N_{1} \\ \frac{d N_{2}}{d t} = ɛ_{2} h_{2} R_{2} N_{2} \end{cases}

[B3a–B3d]

Here, we notice first that

\frac{d}{d t} \log (R_{1}) = \frac{d}{d t} \log (R_{2})

and therefore

ρ = \frac{R_{1}}{R_{2}}

is a constant. We then substitute

R_{1} = R, R_{2} = ρ R

in Eqs. B3 and continue as above. This changes the definition of the competition coefficients to

c_{j} = \frac{ɛ_{i} v_{i} R_{i}}{ɛ_{j} v_{j} R_{j}}

.

If the uptake rates h_i depend on the resource R_i rather than the strain N_i then

{\begin{cases} \frac{d R_{1}}{d t} = - h_{1} R_{1} N_{1}^{v_{1}} - h_{1} R_{1} N_{2}^{v_{2}} \\ \frac{d R_{2}}{d t} = - h_{2} R_{2} N_{1}^{v_{1}} - h_{2} R_{2} N_{2}^{v_{2}} \end{cases}

[B4a, B4b]

Then we define

H = h_{1} / h_{2}

and

ρ = \frac{R_{1}}{R_{2}^{H}}

and again continue as above.

Data Availability

Data deposition: The data reported in this paper have been deposited on Figshare (DOI: https://doi.org/10.6084/m9.figshare.3485984.v1). Source code is available at https://github.com/yoavram/curveball; an installation guide, tutorial, and documentation are available at http://curveball.yoavram.com.

Acknowledgments

We thank Y. Pilpel, D. Hizi, I. Françoise, I. Frumkin, O. Dahan, A. Yona, T. Pupko, A. Eldar, I. Ben-Zion, E. Even-Tov, H. Acar, J. Friedman, J. Masel, and E. Rosenberg for helpful discussions and comments, and L. Zelcbuch, N. Wertheimer, A. Rosenberg, A. Zisman, F. Yang, E. Shtifman Segal, I. Melamed-Havin, and R. Yaari for sharing materials and experimental advice. This research has been supported in part by Israel Science Foundation Grants 1568/13 (L.H.) and 340/13 (J.B.); the Minerva Center for Lab Evolution (L.H.); Manna Center Program for Food Safety and Security, the Israeli Ministry of Science and Technology, and Stanford Center for Computational, Evolutionary, and Human Genomics (Y.R.); Tel Aviv University Global Research and Training Fellowship in Medical and Life Science and the Naomi Foundation (M.B.); European Research Council FP7/2007-2013/ERC Grant 340087 (J.B.); National Science Foundation Grant DEB-1253650 (T.F.C.); and John Templeton Foundation/St. Andrews University Grant 13337 (Y.R. and M.W.F.)

Supporting Information

Appendix (PDF)

Download
539.08 KB

References

1

L.-M. Chevin, On measuring selection in experimental evolution. Biol. Lett. 7, 210–213 (2011).

Featured Topics

Articles By Topic

Featured Topics

Articles By Topic

Featured Topic

Articles By Topic

Significance

Abstract

Sign up for PNAS alerts.

Results

Experimental Validation Design.

Fluorescence experiments.

LacI experiments.

Estimating Growth Parameters.

Growth model.

Model fitting.

Estimating Competition Coefficients.

Competition model.

Model fitting.

Prediction and Validation.

Model prediction.

Experimental validation: Relative growth.

Experimental validation: Relative fitness.

Application: Predicting the Effect of lac Operon Expression on Relative Fitness.

Discussion

Conclusions

Materials and Methods

Strains and Plasmids.

Fluorescent experiments.

Media.

Fluorescent experiment.

Growth and Competition Experiments.

Fluorescent experiments.

Data Analysis.

Fluorescent experiments.

Fitting Growth Models.

Fitting Exponential Models.

Fitting Competition Models.

Data Availability

Code Availability.

Figure Reproduction.

Appendices

Appendix A: Mono-Culture Model.

Appendix B: Mixed Culture Model.

Data Availability

Acknowledgments

Supporting Information

References

Information

Published in

Classifications

Copyright

Data Availability

Submission history

Keywords

Acknowledgments

Authors

Affiliations

Notes

Competing Interests

Metrics

Altmetrics

Citations

Cite this article

Cited by

View options

PDF format

Login options

Recommend to a librarian

Purchase options

Restore content access

Share

Share article link

Share on social media

Further reading in this issue

Long-term expansion and differentiation of adult murine epidermal stem cells in 3D organoid cultures

Gas phase formation of c-SiC3 molecules in the circumstellar envelope of carbon stars

A modern maximum-likelihood theory for high-dimensional logistic regression

Molecular design of a therapeutic LSD analogue with reduced hallucinogenic potential

Bodily maps of emotions

Gas phase formation of c-SiC₃ molecules in the circumstellar envelope of carbon stars