A Bayesian non-parametric mixed-effects model of microbial growth curves.

Substantive changes in gene expression, metabolism, and the proteome are manifested in overall changes in microbial population growth. Quantifying how microbes grow is therefore fundamental to areas such as genetics, bioengineering, and food safety. Traditional parametric growth curve models capture the population growth behavior through a set of summarizing parameters. However, estimation of these parameters from data is confounded by random effects such as experimental variability, batch effects or differences in experimental material. A systematic statistical method to identify and correct for such confounding effects in population growth data is not currently available. Further, our previous work has demonstrated that parametric models are insufficient to explain and predict microbial response under non-standard growth conditions. Here we develop a hierarchical Bayesian non-parametric model of population growth that identifies the latent growth behavior and response to perturbation, while simultaneously correcting for random effects in the data. This model enables more accurate estimates of the biological effect of interest, while better accounting for the uncertainty due to technical variation. Additionally, modeling hierarchical variation provides estimates of the relative impact of various confounding effects on measured population growth.

This is a PLOS Computational Biology Methods paper.

Introduction

Microbial growth phenotypes inform studies in microbiology, including gene functional discovery, bioengineering process development, and food safety testing [1-3]. For example, recent advances in microbial functional genomics and phenotyping, or “phenomics”, have enabled transformative insights into gene functions, proving critical for mapping the genotype to phenotype relationship [4]. Methods such as genome-wide CRISPRi [5] and targeted genome-scale deletion libraries [6, 7] frequently rely upon accurate quantitation of microbial population growth as an assay to identify novel mutants with significant growth phenotypes. Population growth, as measured by the growth curve of a microbial culture, is an aggregate measure of all cellular processes and captures how microbial cells adapt and survive in their environmental niche [8]. Because microbial culturing is a necessary precursor to many experimental procedures in microbiology [9], reproducible results require accurate quantification of the variability in culture state measured through growth [9, 10].

Typical analyses of microbial population growth involve estimating parametric models under the assumptions of standard growth conditions comprised of three successive growth phases: (1) lag phase, in which the population adapts to a new environment, typically fresh growth medium at culture inoculation; (2) log phase, when the population grows exponentially at a rate dependent on nutrients in the environment; and (3) stationary phase, where measurable population growth terminates thereby reaching the culture carrying capacity [11]. Recent studies have shown that the estimates of parameters in these models are highly uncertain [12-14]. This uncertainty arises both from factors of biological interest, such as differences in genetic background and environment, as well as uncontrolled technical noise from experimental manipulation of microbial cultures. While such sources of variability can be modeled using fixed and random effects [15-19], parametric population growth models have additional limitations. Specifically, the parametric assumptions of these growth models require that growth measurements match the sigmoid shape expected for the growth curve under optimum conditions. When population growth deviates from the standard sigmoid shape assumed in these models, model extensions must be developed on a case by case basis for each new experimental perturbation [20, 21]. Additionally, we have shown in previous work that in cases such as extreme stress or strongly deleterious mutations, no parametric growth model accurately represents the growth curve, regardless of extension [19, 22, 23].

Factors affecting microbial growth measurements include both fixed and random effects [24]. Fixed effects are assumed to be drawn from a finite set of perturbations of interest, for example the effect of different concentrations of a chemical on growth that are entirely represented in the dataset. Random effects, conversely, can be viewed as a random sample from a larger set of interest. For example, repeating the same design over many experiments corresponds to sampling the random experimental effect from the theoretical set of all possible experiments that could be conducted with this design [3, 25]. Random effects arising from repeated experimental design are typically referred to as batch effects [26, 27]. Batch effects are often a significant component of measurement noise in high-throughput genomics experiments [28]. However, random effects are not always due to experimental noise, and may represent quantities of direct scientific interest; for example, assaying a set of genetic backgrounds may be viewed as sampling from the set of all possible genetic variants [29-33]. Models which include both fixed and random effects are referred to as mixed effects models.

In this study we present phenom, a general model for analysis of phenomic growth curve experiments based on a Bayesian non-parametric functional mixed effects model of microbial growth. We demonstrate the utility of phenom to analyze population growth measurements of two microorganisms: the hypersaline adapted archaeon, Halobacterium salinarum; and the opportunistic bacterial pathogen, Pseudomonas aeruginosa. H. salinarum is a model organism for transcriptional regulation of stress response in the domain of life Archaea [34-36]. H. salinarum is particularly well adapted to resisting oxidative stress (OS), which arises from the buildup of reactive oxygen species and causes damage to many critical cellular components, including DNA, protein, and lipids [37-43]. Population growth measurements of H. salinarum under OS have been used previously to quantify these harmful effects on physiology, as well as identify regulatory factors important for OS survival [22, 40–42]. The presence of batch effects in H. salinarum OS response was reported (and corrected for) previously [19], but these efforts did not include modeling of individual batch effects for each term in the model. This motivated the explicit deconstruction of batch effects between different factors (e.g. strain and stress), which we have implemented in phenom and reported here.

Pseudomonas aeruginosa is an opportunistic microbial pathogen and a growing problem in hospital-borne infections. Rising antimicrobial resistance of these organisms has necessitated the development of alternative treatment strategies. For example, topical treatment of infected burn wounds with acetic or organic acids (OAs) has been successful [44]. OA impact on growth depends on external pH levels—in acidic environments the OA does not dissociate, but rather freely traverses the cellular membrane as an uncharged particle. Within the neutral cytoplasm, the OA dissociates, and the protons released induce acid stress [45]. Here we apply phenom to the P. aeruginosa dataset, which is foundational for a larger study of P. aeruginosa strains responding to pH and OA perturbation as a potential novel treatment of pathogenic bacterial infections [23].

Stress occurs constantly in the environment: as conditions change, mild to severe cellular damage occurs, and cells must regulate their molecular components to survive [46-49]. Population growth measurements are particularly vital to the study of stress response by providing a quantitative measure of growth differences against a non-stressed control [1]. Our model recovers fixed effects due to high and low levels of OS in H. salinarum and interactions between organic acid concentration and pH in P. aeruginosa. Random effects from multiple sources are corrected, thus enabling more accurate estimates of the biological significance of the stress treatment effect. Notably, in cases where random effect and fixed effect sizes are comparable, we demonstrate that mixed modeling is critical for accurate quantification of model uncertainty. If random effects are not included in the model, the significance of the effect of stress treatments on population growth can be erroneously overestimated. We discuss the implications of these findings for multiple areas of microbiology research.

Materials and methods

Experimental growth data

H. salinarum growth was performed as described previously [22]. Briefly, starter cultures of H. salinarum NRC-1 Δura3 control strain [50] were recovered from frozen stock and streaked on solid medium for single colonies. Four individual colonies per strain were grown at 42°C with shaking at 225 r.p.m. to an optical density at 600 nm (OD600) ∼1.8–2.0 in 3 mL of Complete Medium (CM; 250 NaCl, 20 g/l MgSO4•7H2O, 3 g/l sodium citrate, 2 g/l KCl, 10 g/l peptone) supplemented with uracil (50 μg/ml). These starter cultures were diluted to OD600 ∼0.05 in 200 μl CM (“biological replicates”) then transferred in triplicate (“technical replicates”) into individual wells of a microplate. Cultures were grown in a high throughput microplate reader (Bioscreen C, Growth Curves USA, Piscataway, NJ), and culture density was monitored automatically by OD600 every 30 minutes for 48 hours at 42°C. High and low levels of OS were induced by adding 0.333 mM and 0.083 mM of paraquat to the media, respectively, at culture inoculation.

For P. aeruginosa, laboratory strain PAO1 (ATCC 15692) was grown as described in reference [23]. Briefly, cultures were grown in M9 minimal media supplemented with 0.4% (w/v) glucose and 0.2% (w/v) casamino acids and buffered with 100 mM each of MES and MOPS buffers. Initial cultures were diluted to a starting OD of 0.05 before growth in a microplate reader at a total volume of 200 μl per well. Population growth was measured with a CLARIOstar automated microplate reader (BMG Labtech) at 37°C with 300 rpm continuous shaking. The OD600 was recorded automatically every 15 minutes for a total of 24 hours. A full factorial design of pH and OA concentration was performed for benzoate, citric acid, and malic acid. An experimental batch corresponded to two repetitions of the experiment on separate days with a minimum of three biological replicates of each condition on each day. Two batches for each OA were performed.

All data generated or analyzed during this study are included in this published article (see github repository associated with this study, https://github.com/ptonner/phenom).

Parametric growth curve estimation

For comparison with our non-parametric methods, parametric growth curve models were estimated using the grofit package in R with default parameters [51]. The logistic model was used to fit each curve. Kernel density estimates of parameter distributions were calculated with the Python scipy package with default kernel bandwidth parameters [52].

phenom: A hierarchical Gaussian process model of microbial growth

Gaussian processes

A Gaussian process (GP) defines a non-parametric distribution over functions f(t), defined by the property that any finite set of observations of f follow a multivariate normal distribution [53]. A GP is fully defined by a mean function and a covariance function κ(t, t′) (Eq (1)):

GPs are commonly used for non-parametric curve fitting [53] where is typically set to 0, which we do here. Similarly, we use a common choice for covariance function defined by a radial basis function (RBF) kernel (Eq (2)). Where σ2 is the variance and ℓ is the length-scale. The parameter σ2 controls the overall magnitude of fluctuation in the population of functions described in the GP distribution, while ℓ controls the expected smoothness, with larger ℓ making smoother, slower varying functions more likely. In the process of non-parametric modeling of growth curves, these parameters are adaptively estimated from the dataset.

Fixed effects

We first define the fixed effects models used in this study; these will be augmented with random effects in the next section. We consider fixed effects models of increasing complexity: a mean growth phenotype, a single treatment phenotype, and a combinatorial phenotype with interactions between treatments. All of these models fall under the functional analysis of variance (ANOVA) framework [22, 54]. To estimate a mean growth profile, as in the case of measuring a single condition, a mean function m(t) is estimated from the data by modeling each replicate y(t) for 1 ≤ r ≤ R as consisting of an unknown mean function observed with additive noise (Eq (3)). Where m(t) ∼ GP(0, κ(t, t′)) provides a prior distribution over m, and κ is an RBF kernel with hyperparameters . Here is Gaussian white noise.

When estimating the effect of a perturbation on growth, as in the case of OS, we add a second function δ(t) that represents the effect of the stress being considered. The model then becomes (Eq (4)): where δ(t) ∼ GP(0, κ(t, t′)) also follows a GP prior independently of m, and κ has hyperparameters .

When incorporating possible interaction effects such as those between pH and organic acids in the P. aeruginosa dataset, the model becomes (Eq (5)): for pH p and molar acid concentration c, with α(t) representing the main effect of pH, β(t) the main effect of acid concentration, and (αβ)(t) the interaction between them. Each effect is drawn from a treatment specific GP prior (Eq (6)).

Again, each covariance function is specified by a RBF kernel with corresponding variance and lengthscale hyperparameters that adapt to the observed data. All models in this section correspond to Mnull for their respective analyses, as they do not include any random effects.

Random effects

The first random effects added to the model were those used to account for batch effects, in the model Mbatch. Under this model, each fixed functional effect is modified by a GP describing the population of possible batch-specific curves. For example, under the model of interaction effects on growth (Eq 5), replicate r from batch k is modeled as (Eq 7)): where the functions shared with Mnull are highlighted, the functions m((t), , , and are the corresponding random batch effects, and . Similar to the fixed effects, the batch effect functions are drawn from shared GP priors (Eq (8)): with kernel hyperparameters , , and that are distinct from those for the corresponding fixed effects, allowing for different variance and lengthscales between fixed and random effects. Other Mnull models are converted to Mbatch similarly, with each fixed effect becoming a mean of a GP prior for each batch effect.

Mfull develops the hierarchy one step deeper by adding replicate effects to Mbatch (Eq 7). Specifically, the error model ϵ is now described by a GP: ϵ ∼ GP(0, κ(t, t′)) with corresponding hyperparameters, accounting for replicate-specific variability rather than simply white noise.

Bayesian inference

The unknown functions (m(t), δ(t), α(t), β(t), and (αβ)(t)), kernel hyperparameters ( and ℓ) for each group of latent functions, and observation noise parameters () are all estimated by Bayesian statistical inference. In Bayesian inference, prior distributions on unknown quantities (e.g. p(m(t), θ) = p(m(t)|θ) × p(θ)) are combined with the likelihood, p(y(t)|m(t)) to obtain the posterior distribution p(m(t), θ|y(t)).

Latent functions are grouped by shared kernel hyperparameters into related sets (e.g. treatment effects, interaction effects, batch effects), which then provide the GP prior for the latent function. For each group, is assigned a Gamma(α, β) prior, with fixed effects assigned as a Gamma(10, 10) prior and random effects assigned a Gamma(7, 10) prior. ℓ was assigned an inverse- Gamma(α, β) prior, with parameter α = 6 and β = 1 for H. salinarum, and α = 2, β = 3 for P. aeruginosa fixed effects and α = 10, β = 1 for P. aeruginosa random effects. Noise variance was also assigned a gamma prior.

Bayesian inference was then performed, with the posterior distribution obtained by sampling using Markov chain Monte Carlo (MCMC) implemented with the Stan library, which uses a Hamilitonian Monte-Carlo procedure with No-U-turn sampling [55]. Multiple chains were run to diagnose convergence, with all parameter posterior means confirmed to have converged within as recommended [56].

Results

Hierarchical batch effects typical in phenomics datasets render parametric models ineffective

In the dataset used here, population growth for each of P. aeruginosa and H. salinarum cultures was monitored under standard (non-stressed) conditions vs. stress conditions (see Materials and methods and references [22, 23] for precise definition of “standard conditions” for each organism). Specifically, cultures were grown in liquid medium in a high throughput growth plate reader that measured population density at 30 minute intervals over the course of 24 hours (P. aeruginosa) or 48 hours (H. salinarum); the resulting data are shown in Fig 1. Experimental designs for each organism included biological replicates (growth curves from different colonies on a plate), technical replicates (multiple growth curves from the same colony), varying conditions (stress vs standard), and batches. Throughout, we define “batch” as a single run of the high throughput growth plate reader. In each run, this plate reader measures the growth of 200 individual cultures across a range of perturbations, including varying stress conditions and genetic mutations (see Methods). H. salinarum was grown under high (0.333 mM paraquat (PQ)) and low (0.083 mM PQ) levels of oxidative stress (OS); the data are combined from published [19, 22, 41] and unpublished studies (Fig 1A). The OS responses of H. salinarum were compared to a control of standard growth in rich medium, representing optimal conditions for the population. The experimental design was replicated in biological quadruplicate and technical triplicate, across nine batches (Fig 1A, individual curves and axes). P. aeruginosa was grown in the presence of increasing concentrations of three different organic acid (OA) chemicals (0–20mM; benzoate, citric acid, and malic acid), each combined with a gradient of pH (5.0–7.0) [23]. Each P. aeruginosa growth condition was repeated across 3 biological replicates and two batches (Fig 1B). The different P. aeruginosa and H. salinarum experimental designs with varying numbers of replicates at each level provides a rich test bed for modeling the impact of random effects with phenom (Fig 1B, S1 and S2 Figs).

Fig 1

Batch variation in high throughput phenomics studies.

A: Population growth measurements of H. salinarum under standard conditions (blue), and low (orange) and high (green) levels of OS. Individual measurement curves are replicates and each graph panel is a different batch. B: Growth of P. aeruginosa strain PA01 under gradient of pH (5–7) and citric acid (0–20 mM). Colors represent different batches.

Figs 1 and 2 demonstrate the two key issues described above and addressed in this paper. First, batch effects are present in both H. salinarum and the P. aeruginosa datasets. For H. salinarum, clear differences in growth under both standard and stress conditions are observed in the raw data across experimental batches (i.e. separate runs of the growth plate reader instrument; Fig 1). Some batches show a different phenotype, with either a complete cessation of growth or an intermediate effect with decreased growth relative to standard conditions. For example, in some batches, populations stressed with low OS grow at the same rate and reach the same carrying capacity as populations grown under standard conditions. For P. aeruginosa, a clear difference between batches grown under 10 mM citric acid at pH = 5.5 is observed [Fig 1B (graph in fourth column, third row) and Fig 2D]. Like with citric acid, batch effects were also found in some of the other conditions considered (e.g. growth under malic acid, S1 and S2 Figs).

Fig 2

Batch effects are prevalent in microbial phenomic datasets.

(A) Parametric fits to H. salinarum growth curves. (B) Residuals of parametric growth curve fit. (C) Growth of H. salinarum under standard conditions (blue), low (orange) and high (green) OS across three batches. (D) Measurement of P. aeruginosa growth under 10mM citric acid at 5.5 pH. Measurements for each condition vary significantly with batch.

Second, standard parametric growth curve models fail to describe experimental measurements adequately (Fig 2A and 2B), as we have shown previously with both datasets [19, 22, 23]. In Fig 2, we examined the impact of batch and replicate effects on our data by considering how they change parameters estimated from a mixed effects parametric model of population growth [32]. We focused on calculating μmax, the maximum instantaneous growth rate attained by the population, as this is a commonly used parameter for comparisons between conditions [19, 57]. Variation in μmax estimates were observed both on the replicate and batch level, as shown by the kernel density estimates (KDE) of μmax for each stress level (S3 Fig). The variance in μmax is remarkably high: the 95% confidence interval for μmax under standard growth is 0.050–0.141, a nearly 3-fold change between the lower and upper interval limits. Thus, while the t-test conducted on μmax estimates between standard conditions and each stress level is statistically significant (S3 Fig), it is difficult to conclude: (a) what the true magnitude of the stress effects may be; and (b) to what degree the variation due to replicate and batch should inform biological conclusions. The error of the logistic growth model under each PQ condition was also examined. Error increased under high OS (S4 Fig). High OS induces a growth phenotype that deviates heavily from the sigmoidal growth curve assumed in the logistic model as well as in other commonly used growth models. This leads to a poor fit under the high OS condition as has been shown previously (S4 Fig, [19]). The residuals under standard, low, and high OS conditions also appear to be dependent. Our previous work also demonstrated poor fits to the P. aeruginosa data using parametric models [23]. Taken together, the initial assessment of these two datasets indicates that: (a) technical variation due to batch and replicate in growth curve data can be high; and (b) commonly used standard parametric models are not able to adequately capture or correct for these sources of variability. These sources of error need to be corrected in order to model true growth behavior and inform biological conclusions from the data.

A hierarchical Bayesian model of functional random effects in microbial growth

We previously established the ability of non-parametric Bayesian methods to improve the modeling of growth phenotypes [19, 22, 23]. Here, we describe phenom, a fully hierarchical Bayesian non-parametric functional mixed effects model for population growth data. We highlight the utility of phenom to correct for confounding, random effects in growth phenotypes.

In order to model both biological and technical variation in microbial growth (Fig 3), we first assume that a set of population growth measurements are driven by an (unobserved) population curve m(t) (Fig 3A, blue curve) of unknown shape. For example, m(t) might represent the average growth behavior of an organism under standard conditions. This mean growth behavior may be altered by a treatment effect, represented by an additional unknown curve δ(t) (Fig 3A, orange curve). For example δ(t) may represent the effects on growth induced by low or high levels of OS (Fig 2A). The average growth behavior of a population under stress conditions would then be described by the curve f(t) = m(t) + δ(t).

Fig 3

Hierarchical model of functional data.

Representative diagram of hierarchical variation present in microbial growth data. Each tier of graphs represents a different variation source, and lines indicate relationship between them: experimental condition is the true growth behavior of interest, with the condition repeated across batches, and replicates repeated within each batch. (A) Functional phenotypes m(t) (blue), m(t) + δ(t) (orange), and δ(t) (green curve in inset). (B) Batch effects on m(t) and m(t) + δ(t). Each plot is a different batch, solid lines are the true functions as in (A), and the dashed lines are the observed batch effect of m(t) and m(t) + δ(t) for the corresponding batch. (C) Replicate effect within batches. Each axis is a different replicate, solid and dashed lines as in (B), dotted-dashed line is the observed replicate function. (D) Observations from the model described in (A-C). Each curve is sampled with a mean drawn from the global mean, with added batch and replicate effects (dotted-dashed lines in C) and iid observation noise. Each axis is a different batch. The smooth solid lines are the true functions m(t) and m(t) + δ(t) in (A).

When considering a combinatorial experimental design, such as that described for P. aeruginosa growth (Fig 1B), we model independent effects of different treatments as well as their interaction via the form (Eq (9)):

Here, y(t, i, j) denotes the observed population size at time t with treatments i and j of two independent stress conditions. Additionally, α(t) and β(t) are the independent effects of each stress condition, and (αβ)(t) is their interaction. This model corresponds to a functional analysis of variance [58], which we have previously used to estimate independent and interaction effects of microbial genetics and stress [22]. For the analysis of P. aeruginosa, we model the effect of pH (α for pH = p), organic acid combination (β for concentration = c) and their interaction ((αβ)), as well as their random functional effect equivalents (see Section “phenom: A hierarchical Gaussian process model of microbial growth”).

Variability around these fixed effect growth models is described by additional, random curves associated with two major sources of variation: batch and replicate (Fig 3B and 3C). Batches correspond to a single high-throughput growth experiment and replicates are the individual curve observations within a batch. Using phenom throughout this study, we only compare replicates that are contained within the same batch. This is due to the nested structure between batch and replicates (Fig 3). Noise due to both replicate and batch do not appear to be independent identically distributed (iid), as observed in the correlated residuals around the mean for each experimental variate (S5A and S5B Fig). Each observed growth curve is therefore described by a combination of the fixed effects and the corresponding batch and replicate effects (Fig 3D). Both replicate and batch variation are modeled as random effects because the variation due to both sources cannot be replicated, i.e. a specific batch effect cannot be purposefully re-introduced in subsequent experiments. Instead, these variates are assumed to be sampled from a latent distribution [59]. Combining the fixed and random effects, we arrive at a mixed-effects model of microbial phenotypes.

We adopted a hierarchical Bayesian framework to model these mixed effects. In this framework, batch effects are described by a shared generative distribution, allowing them to take on distinct values while still pooling across replicates for accurately estimating the generating distribution [60]. We use Gaussian process (GP) distributions for all groups in the model. GPs are flexible, non-parametric distributions suitable for smooth functions [53]. To assess the impact of incorporating random effects on estimation of the treatment effect of interest, we analyze three models of increasing complexity: Mnull excludes all hierarchical random effects, Mbatch incorporates batch variation only, and Mfull incorporates both batch and replicate variation. These models, collectively called phenom, were implemented using the probabilistic programming language Stan [55] to perform Bayesian statistical inference for all unknown functions and model parameters through Hamiltonian Monte Carlo sampling (see Materials and methods).

In previous work [19] we identified and corrected for batch effects in a single transcription factor mutant’s stress response, but this model did not provide an explicit deconstruction of batch effects between different factors (e.g. strain and stress) and could therefore not determine which factors were most strongly impacted by batch effects. Moreover, this approach utilized a standard GP regression framework, which has well-established limitations on dataset size, limiting its applicability to the large datasets we consider here. In reference [22] we described a functional ANOVA model for microbial growth phenotypes, which corresponds to the Mnull model in the phenom case. Again, a global batch effects term was included but individual batch effects were not modeled, and the computational approach utilized (Gibbs sampling) was prohibitively slow for the complete phenom model. phenom represents a significant advance on these previous modeling approaches and computational methodologies.

In order to demonstrate the impact of batch effects on the conclusions drawn from the analysis of microbial growth data, we estimated the latent functions driving both H. salinarum and P. aeruginosa growth using the Mnull model of phenom, with each batch analyzed separately (Fig 4). This corresponds to the analysis that would be conducted after generating any single set of experiments from a batch, without considering or controlling for batch effects, and therefore provides a test of the impact of ignoring batch effects.

Fig 4

Mnull model estimates are confounded by batch effects.

Posterior intervals of functions are shown for different analyses where phenom Mnull was fit using data from each batch separately. In all plots, solid line represents posterior mean, shaded region indicates 95% credible region, and each color corresponds to a different posterior conditioned on data from a single batch. (A) Posterior intervals of m(t) under standard growth, and δ(t) under low and high OS response for H. salinarum. (B) Posterior interval of interaction function (αβ)(t) for P. aeruginosa growth in indicated pH and acid concentration.

For H. salinarum, growth data under standard conditions was used to estimate a single mean function, m(t), and fixed effects were estimated for differential growth under low and high OS as δ(t) (Fig 4A). For the P. aeruginosa dataset, batch effects on the interaction between pH and organic acid concentration was represented by a function (αβ)(t), again estimated non-parametrically (Fig 4B). However, rather than reporting (αβ)(t) directly, we report its time derivative, which has the interpretation of instantaneous growth rate rather than absolute amount of growth, and provides an alternative metric for assessing the significance of a treatment effect on growth [61]. Specifically, assessing growth curve models can benefit from the estimates of derivatives as they may more accurately represent the differences between growth curves [58].

Fitting the Mnull model to each separate batch reveals that the posterior distributions obtained for each function of interest (m(t), δ(t), and (αβ)(t)) are highly variable across batches (Fig 4). This is observed in both the H. salinarum and P. aeruginosa datasets, where the experimental conditions, and therefore the underlying true mean functions, remain constant across batches in each case. Such variability can impact conclusions. We specifically assess the changes in statistically significant treatment effects, i.e. at time points where the effect (δ(t) or (αβ)(t)) has a 95% posterior credible interval excluding zero, indicating high confidence that the treatment effect at that time-point differs from the control. For example, in the low OS condition in the H. salinarum dataset, both the statistical significance of δ(t) and the sign (improved vs. impaired growth) differs between batches (Fig 4A, center). Additionally, the effect of low oxidative stress at time zero is estimated to be non-zero for many of the batches. This is due to technical variation that introduces an artificial offset in OD measurements at the beginning of the growth experiment. Such variation can arise from various factors, including variation between growth state in starter cultures and technical variation in plate reader measurements at low OD (Fig 1A). A similar batch variability was observed under high OS, but due to the stronger effect of the stress perturbation, estimates of δ(t) are less affected by batch and replicate variation (Fig 4A, right). Similarly, the batch variability observed in the raw P. aeruginosa growth data (Fig 1B) results in significantly different posterior estimates of the interaction effect (αβ)(t) across batches, as seen by the lack of overlap between 95% credible intervals (Fig 4B). Differences observed include the timing and length of negative growth impact (benzoate and citric acid), and completely opposite effects with either strong or no interaction (malic acid). In addition, the posterior variance of each function, which indicates the level of uncertainty remaining, is low for each batch modeled separately. This indicates high confidence in the estimated function despite observed differences across batches. These analyses suggest that use of a single experimental batch leads to overconfidence in explaining the true underlying growth behavior.

Hierarchical models correct for batch effects in growth data

To demonstrate the use of phenom to combat the impact of batch effects on growth curve analysis, we combined data across all batches and performed the analysis using each of the Mnull, Mbatch, and Mfull models (Fig 5). Estimates of m(t) between each model were largely similar, likely due to the abundance of data present to estimate this variable (S6 Fig). Instead, we focus on the estimates of δ(t) for low and high OS response of H. salinarum (Fig 5A) and the interaction (αβ) between pH and OA concentration effects on P. aeruginosa growth (Fig 5B). In cases where Mnull differs from Mbatch and Mfull, this indicates an inability for this model to correctly represent the uncertainty due to random effects in the data, which have been shown to be prevalent across different batches of experiments (Fig 4).

Fig 5

Hierarchical models of growth control for batch effects.

Posterior intervals of functions estimated by models of increasing hierarchical complexity: Mnull (blue), Mbatch (orange), and Mfull (green). Solid line indicates posterior mean and shaded regions indicate 95% credible regions. We also plot Mnull and Mfull in separate axes to highlight the posterior intervals estimated by these models. (A) Posterior interval of δ(t) for low (left) and high (right) OS response by H. salinarum. (B) Posterior interval of interaction function (αβ)(t) for P. aeruginosa growth in indicated pH and organic acid concentration.

Growth impairment in the presence of low OS relative to standard conditions (i.e. δ(t)) is estimated to be significant during the time points of ∼0–13 and ∼19–40 hours under Mnull. In contrast, only time points ∼23–31 are significantly non-zero under Mbatch, and no significant effect is identified under Mfull(Fig 5A, left). Conversely, due to the stronger stress effect in the high OS condition (Fig 5A, right), estimates of δ(t) were significantly non-zero under all three models, with only minor differences between the three model estimates. This highlights the importance of controlling for batch and replicate variability as in Mfull: even when estimating the low OS treatment effect under Mnull with all available data, without accounting for batch and replicate random effects the posterior estimates of δ(t) are overconfident and do not accurately represent the uncertainty with respect to the true treatment effect. The lack of significance under Mfull suggests that additional data are needed to confidently identify the true treatment effect in the presence of batch and replicate variation. Modeling the batch effects has also corrected for variability in treatment effects due to technical variation at the inoculation of the growth plate (time zero of the experiments) (S7 Fig).

The impact of modeling hierarchical variation on estimating interaction effects in P. aeruginosa growth was condition dependent (Fig 5B). Across all conditions, a decrease in posterior certainty on the true shape of the underlying function was again observed under Mbatch and Mfull. For benzoate and malic acid, the interaction between pH and acid concentration no longer appears to be a significant effect after accounting for batch and replicate variation, while the larger interaction under citric acid remains significant. As in the comparison of oxidative stress treatments in H. salinarum, stronger effect sizes are required to be confidently distinguished in the face of batch and replicate variability. Finally, the relative conclusions made for the absolute function scale are comparable to those of the derivative estimates for P. aeruginosa, highlighting the flexibility with which treatment effects can be analyzed as most relevant to the researcher (S8 Fig).

For both H. salinarum response to OS and P. aeruginosa growth under pH and OA exposure, an increase in posterior variance was observed under Mbatch and Mfull compared to Mnull (S9 Fig). However, posterior variance of δ(t) in the H. salinarum OS response was higher under Mbatch compared to Mfull. In this case, controlling for replicate effects appears to increase the signal needed to identify δ(t). In contrast, these variances are equal in the P. aeruginosa data, indicating that the relative improvement in variance afforded by modeling batch vs. replicate effects may be dataset dependent.

Variance components demonstrate the importance of controlling for batch effects

Variance components, which correspond to the estimated variance of each effect in the model, can be used to compare the impact each group has on the process of interest [24]. To better understand sources of variability in growth curve studies, we used phenom to estimate the variance components for each dataset above. In our hierarchical non-parametric setup, these variance components are the variance hyperparameters (e.g. σ2) of the Gaussian process kernels for each fixed and random effect group. These parameters control the magnitude of function fluctuations modeled by the GP distribution. Larger variance implies higher effect sizes and therefore a larger impact on the observations.

We show the value of variance components by considering the effects identified by Mfull for H. salinarum under low OS (Fig 6). The variance of the data is partitioned between the mean growth (m(t)), the OS (δ(t)), batch effects (batch curves of m(t) and δ(t)), biological noise (e.g. replicate variability) and instrument noise (). This analysis confirms that batch effects, compared to the other sources of experimental variability in the dataset (replicate noise and measurement error), are between 2 to 10 times more impactful on the phenotype measurements. Additionally, variance components enable comparisons between the experimental and treatment factors in the data. Of particular note is that the variance of the treatment of interest, δ(t), and the batch effects are similar in magnitude, at least in the case of a low-magnitude stress such as 0.083 PQ for H. salinarum. This suggests that proper modeling of this treatment requires both sufficient batch replication and accurate modeling of batch effects in those data. Future studies of similar phenotypes can be guided by these estimates in experimental design, choosing an appropriate batch replication for the degree of noise expected [62]. However, the extent of replication required may depend upon the dataset (factorial design, treatment severity, etc). Taken together, variance components provide an aggregated view of the contribution by various factors and guide future experimentation. Iterative rounds of phenom analysis and growth experiments can then determine the best suited designs, for example by leveraging the estimated batch effect variance from pilot experiments to determine the number of batch replications necessary to reliably estimate a treatment effect of given magnitude. The use of phenom for such formal statistical experimental design calculations represents an exciting direction for future work.

Fig 6

Posterior variance components in the phenom hierarchical phenotype model.

Posterior intervals are shown for the kernel variance hyperparameter for different groups of effects from phenom estimated on H. salinarum growth under low OS. Groups correspond to m(t) (mean), δ(t) (stress), batch effects (batch), replicate noise (biological), and measurement error (noise).

Discussion

We have provided a framework to test and control for random effects in microbial growth data using the hierarchical non-parametric Bayesian model, phenom (Fig 3). Analysis with phenom indicates that random effects (both batch and replicate) appear in the two microbial population growth datasets studied here, and constitute significant portions of the variability (Fig 1). Failure to correct for these effects confounds the interpretation of growth phenotypes for factors of interest in a large scale phenotyping analysis (Fig 4). phenom controls for these random effects and provides accurate estimates of the growth behavior of interest (Fig 5). Additionally, phenom can be used to estimate variance components, providing information about the relative impact of various sources of noise in the data (Fig 6). Controlling for batch effects in these datasets was therefore key to making accurate biological conclusions.

Related fields of functional genomics, such as transcriptomics, have seen considerable interest in controlling for different experimental sources of variation, broadly labeled as batch effects [28, 62–67]. These studies have shown that differences between batches first need to be corrected to avoid erroneous conclusions [68]. Here we have shown that, like in transcriptomics data, controlling for sources of variation in phenomics data—particularly due to batch—are an important step in making accurate biological conclusions regarding population growth. Additionally, the use of random batch effects in phenom highlights cases where additional information may be gained by further experimentation. Specifically, in cases where treatment effects differ strongly across batches, there may be underlying biological differences driving the variation. Follow-up experimental designs can then aim to delineate these effects directly in a way not confounded by batch. phenom establishes a complete and general method of controlling batch effects in microbial growth phenotypes, overcoming significant weaknesses of previously developed techniques.

Although we focus here on replicate and batch variation, the phenom model is easily extended to incorporate alternative or additional random and fixed effects appropriate for settings with other sources of variation. For example, depending on the experimental design, phenom could control for variation among labs, experimental material, culture history, or genetic background [25, 69–75]. phenom flexibly incorporates additional sources of variation and/or interaction between design variables, as demonstrated with the two different designs analyzed for H. salinarum and P. aeruginosa here. This flexibility allows phenom to be applied to control for many sources of technical variation within microbial population growth data, thereby improving the analysis and resulting conclusions regarding quantitative microbial phenotypes. We therefore expect our model to find broad applications in fields such as bioprocess control, microbial bioengineering, and microbial physiology.

P. aeruginosa growth under benzoate and pH gradient.

Growth of P. aeruginosa strain PAO1 under gradient of pH (7–5) and benzoate (0–20). Colors represent different batches.

(EPS)

Click here for additional data file.

P. aeruginosa growth under malic acid and pH gradient.

Growth of P. aeruginosa strain PAO1 under gradient of pH (7–5) and malic acid (0–20). Colors represent different batches.

(EPS)

Click here for additional data file.

KDE of μmax for H. salinarum growth across batches.

Crosses indicate significant difference between μmax standard conditions and each OS level (one-sided t-test, p < 0.05).

(EPS)

Click here for additional data file.

Error in parametric growth models.

Distribution of error (MSE) for each condition when fit with a logistic growth curve. The box show shows the inter-quartile range, red line is the median, whiskers show the 1.5 inter-quartile range, and the individual points are outliers.

(EPS)

Click here for additional data file.

Residual structure of microbial growth data across batches.

(A) Individual replicate curve residuals around the mean of the respective batch. Only standard conditions are shown. (B) Residual of the mean behavior for each batch around the global mean (standard condition only).

(EPS)

Click here for additional data file.

Posterior comparison of m(t) for H. salinarum growth across batches.

Posterior interval of m(t) for H. salinarum standard growth.

(EPS)

Click here for additional data file.

Posterior intervals of low oxidative stress batch effects estimated from Mfull.

Full estimates of the δ(t) batch effect under Mfull are shown, with solid lines representing posterior mean and shaded region representing 95% credible intervals (left). 95% posterior credible interval of batch effects for δ(t) at time zero are shown (right), with crosses marking posterior means. Many of the batch effects for δ(t) are estimated to be non-zero at the start of the experiment, reflecting the impact of technical variation in the high-throughput readings at the start of the growth curves.

(EPS)

Click here for additional data file.

Posterior intervals of interactions for P. aeruginosa on an absolute growth scale.

Posterior intervals of interactions in Fig 5B, but reported here on an absolute (log OD) scale. The same data is reported on the derivative (d log OD / dt) scale in the main text Fig 5B.

(EPS)

Click here for additional data file.

Posterior variance of function estimates under different models.

Each plot shows the posterior variance of a function at each time point under each of Mbatch and Mfull versus Mnull. (A) δ(t) estimated for H. salinarum growth under low (left) and high (right) OS. (B) (αβ)(t) at pH = 5, mM malic acid = 10.

(EPS)

Click here for additional data file.

23 Jun 2020

Dear Dr. Schmid,

Thank you very much for submitting your manuscript "A Bayesian Non-parametric Mixed-Effects Model of Microbial Phenotypes" for consideration at PLOS Computational Biology.

As with all papers reviewed by the journal, your manuscript was reviewed by members of the editorial board and by several independent reviewers. In light of the reviews (below this email), we would like to invite the resubmission of a significantly-revised version that takes into account the reviewers' comments. The reviewers bring up important points, but they all seem quite addressable.

We cannot make any decision about publication until we have seen the revised manuscript and your response to the reviewers' comments. Your revised manuscript is also likely to be sent to reviewers for further evaluation.

When you are ready to resubmit, please upload the following:

[1] A letter containing a detailed list of your responses to the review comments and a description of the changes you have made in the manuscript. Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out.

[2] Two versions of the revised manuscript: one with either highlights or tracked changes denoting where the text has been changed; the other a clean version (uploaded as the manuscript file).

Important additional instructions are given below your reviewer comments.

Please prepare and submit your revised manuscript within 60 days. If you anticipate any delay, please let us know the expected resubmission date by replying to this email. Please note that revised manuscripts received after the 60-day due date may require evaluation and peer review similar to newly submitted manuscripts.

Thank you again for your submission. We hope that our editorial process has been constructive so far, and we welcome your feedback at any time. Please don't hesitate to contact us if you have any questions or comments.

Sincerely,

Jason Papin

Editor-in-Chief

PLOS Computational Biology

***********************

Reviewer's Responses to Questions

Comments to the Authors:

Please note here if the review is uploaded as an attachment.

Reviewer #1: Summary:

The authors present a new Bayesian modeling approach for cell growth that is adaptable, allows for incorporation of several effectors, and includes batch-to-batch variation. While there is novelty to the overall approach, the manuscript is difficult to read and the implementation seems convoluted. If these issues can be addressed, the approach described could be adopted widely.

Major Points:

Overall, I found the modeling approach intriguing but the manuscript somewhat difficult to read. This is unfortunate given that growth models of microbes should be a topic accessible by a wide audience. The modeling approach described here has the potential to be adopted somewhat widely, as many engineers and biologists are looking for better models to describe microbial growth. The ability to incorporate batch-to-batch effects is certainly a benefit, and the ability to add in additional effects, inhibitors, etc. into the framework may make this approach universally adaptable. With that, there are problems that should be addressed. Currently, it is not clear how one would apply this modeling framework to build a model for monitoring bioprocesses. What phenomics experiments should be run a priori? How are model parameters found? How many batches and replicates are needed? Is batch-to-batch variability reduced by including more variables in the model? The code is available through GitHub, but the implementation in Python is not for the inexperienced user. From this, my main suggestion is to improve accessibility and provide user-friendly instructions and guidelines. This could lead to a large number of citations.

The introduction does not mention the contribution of “genome-scale” kinetic and dynamic flux-based modeling strategies in determining growth rates, effects of inhibitors, multiple substrates, genetic perturbations, etc. There is a wealth of knowledge in the literature that should be referenced. How does the Bayesian modeling approach proposed in this paper compare to the other dynamic genome-scale approaches that try to incorporate mechanistic effects?

Eq. 6 contains several conditional statements, and these can be expanded to incorporate other conditions, I’m assuming. How is this different from other parametric Monod-type models that include terms for inhibition, multiple substrates, etc.? Don’t they all contain constants that must be found through experimental data? Can’t a Monod-type model contain a batch-to-batch variability parameter? An effective approach would be to compare several different model types given the same experimental data generated in this study. There is some reference to previously published studies in this regard, but the findings remain unclear. It seems some direct comparisons can be made in this manuscript.

It’s not clear why a single batch can contain multiple growth curves. Are these all of the same medium lot? All of the curves seem to be related even though they were grown separately. A batch is usually from a single fermenter. I'm not sure how to address this better, but it is another somewhat confusing aspect in the manuscript.

Reviewer #2: This manuscript by Tonner et al. describes a much-needed contribution to the rigorous modeling of bacterial growth. They extend their previous work on Bayesian modeling of microbial growth curves by creating a fully non-parametric framework that accounts for experimental and technical variability. They demonstrate how incorporation of batch modeling terms can change interpretation of results from growth experiments performed in varying media conditions for P. aeruginosa and H. salinarum. Overall, the manuscript is very well written and logically builds on previous work. This approach will advance studies of microbial growth curves, most of which currently apply heuristic methods for summarizing growth curves into a single metric (growth rate, AUC, etc.). I have essentially no concerns related to the utility or implementation of the phenom method, but highlight one issue with the data source for P. aeruginosa that warrants reprocessing and reanalysis. I will use the methods the authors have developed in my own projects and I am very enthusiastic about the line of research being developed by the manuscript authors.

Major concerns

1. Unrealistic values of log(OD) are plotted in most figures describing growth of P. aeruginosa in the manuscript. One concrete demonstration of the problem is in the P. aeruginosa plots in Figure 1B; log(OD) = 5 is reported in many panels. Even assuming a base of 2, that is a raw OD of ~30, an order of magnitude higher than the maximum values that might be expected in very favorable growth conditions. What is the base of the logarithm? In digging through reference 23, it was reported that raw OD values were log2 transformed; however, the raw values in the github repository for that paper are between around 0-1.5 yet the processed data are in the range of around 0-6. In the paper itself, all plotted OD values seem to be from the raw data rather than processed data (or a different version of the processed data). In the publication’s github repo, it looks like there may be a problem with the scaling process (https://github.com/amyschmid/pseudomonas-organic-acids/blob/master/preprocess.ipynb).

Specifically, the problem seems to be in this loop:

for k,index in group.groups.items():

temp = data.loc[:,index]

od = temp.values[:5,:].ravel()

coeff = np.polyfit(time.tolist()*temp.shape[1],od,2)

temp = temp - np.polyval(coeff,data.index.values[0])

A second order polynomial is fit to the first 5 time points, and then the polynomial is evaluated at the first time point. This results in subtraction of a large negative number (-6.17 for the example growth curve I used), inflating the OD for the entire growth curve. The motivation for the polynomial fitting procedure on only the first five time points, then the scaling by the polynomial fit value for only the first time point, is unclear to me. Why not remove the first five time points, and scale the entire curve by subtracting by the real value of the first time point? The choice to scale based on a within-sample polynomial fit may also introduce artificial technical variation between samples.

Please resolve this issue with the data, rerun all analyses for P. aeruginosa, and include the base of the logarithm in the axis titles or the figure legend. I do not think that this issue points to any underlying problems with the phenom method. The H. salinarum data do not raise any red flags for me, but I would encourage you to inspect the data with the same rigor that you should for P. aeruginosa, given the potential errors.

2. Starting around the text describing Figure 4, “significant” starts to be used in a way that isn’t defined in the text (line 207, “significantly different posterior estimates”, line 220, “... is estimated to be significant during the time points...”; again on lines 221, 222, 224). Please define the procedure used to establish significance in these cases, or use a different term if a formal procedure was not used (which seems completely fair to me given the use of fully probabilistic models).

3. This is not a concern but a strong recommendation: the points about shrinkage of the credible interval with Mfull in Figure 5C are perhaps the most important in the manuscript, but are hard for the reader to grasp because of the overlapping intervals. The point is clear in the text but could be drawn out more clearly in the figure. There are many ways you could emphasize this, but I think the easiest would be having subpanels beneath each panel in Figure 5C which individual plots the estimates of Mnull, Mbatch, and Mfull (e.g., 3 plots with individual models and one plot with estimates from all models shown together for context).

Minor concerns

1. Lines 106-112, “technical”, “biological”, and “batch” should be defined by describing the experiments performed. Based on my own terminology usage, I assume for H. salinarum that 9 experiments were performed in which a frozen stock was struck on plates, a colony/lawn was gathered and pre-cultured, and the preculture was standardized via OD and used to inoculate each well (“9 batches”). Within each of these 9 separate experiments, each environmental condition/media condition was replicated in quadruplicate (i.e., “4 biological replicates”), and at each time point OD was measured from in each of the biological replicate wells three times (i.e., “3 technical replicates”). Many other researchers consider the “biological replicates” I describe here to be “technical replicates”. We also recommend indicating whether media and/or supplements (e.g., organic acids) were made separately for each “batch”, or whether the same solutions were added to media ingredients form the same lot, etc.

2. Line 292; “Cultures were then diluted to OD600~0.05 in a high throughput microplate reader” - was the culture diluted to an OD600 of ~0.05 as measured in a 1cm cuvette and then transferred to the microplate, or was culture diluted such that within each well of the microplate the OD600 was 0.05 at the beginning of the growth experiments as measured by the platereader? Please clarify, and if the latter is the case, please indicate the format of the microplates (e.g., 96-well flat-bottom) and the total volume used in each well. Please provide a similar level of detail in the description of culture setup for P. aeruginosa as well.

3. All figure captions and supplemental figure captions: mu and delta are occasionally expressed as a function of x (instead of t, as in the rest of the manuscript).

4. Labelling some panels of Figure 4 with species information would make the figure more interpretable; e.g., having only two batches in panel C makes it tough to figure out that the two colors are still referring to different batches and not some new concept.

5. In Figure 4C, why is the posterior estimate of the interaction term expressed as the derivative of the growth curve rather than using log(OD) as for all the other estimate plots?

6. In Figure 5A, label the left/right panels as low/high OS, respectively, on the figure itself

7. Lines 71-73, “The presence…” is incomplete, seems like the middle of the sentence should have the following inserted: “..., but [these efforts did not include modeling of] individual batch effects for each term in the model.”

Below are recommendations that we believe would improve the readability of the manuscript but do not think necessary to revise as part of the peer review process.

The broad readership of PLOS Computational Biology would benefit from explicit definition of several statistical modeling terms used throughout the paper. We recommend defining the following terms at first mention in the introduction. Some of these terms are defined in the text but could use additional clarification. The following are just a few of these terms; it may help to have an application-oriented computational biologist read your manuscript and identify terms that they don’t understand which limit their ability to understand the purpose and findings of the paper.

1. fixed effects. This definition seems suitable

2. random effects. This definition seems suitable but would benefit from using a statistical synonym for “population” since population implies something different in “population growth model”. Conversely, “population growth model” could use a biological synonym for “population”.

3. parametric model. Clearer indication that the parameters that make standard population growth models “parametric” are those that bake in the assumptions of lag-log-stationary, rather than another assumption about a distribution.

4. secondary model. The average reader will not be able to define “secondary model” based on the context provided alone.

Reviewer #3: This paper presents phenom, a method and software tool for fitting non-parametric mixed-effects models to microbial growth curves. The method is a generalized version of previous work by the same authors with a focus on handling batch and replicate effects without relying on a parametric growth curve. My main concerns are twofold: 1.) are interesting parts of the stress response swept into batch effects, and 2.) do the perturbation functions extracted by phenom correspond to the underlying biology?

Major Comments

My first concern is largely technical. The phenom method assumes that batch effects are simply random effects that perturb the true underlying \\mu(t). As a microbiologist, I'm not sure this is true. Take, for example, the PA growth curves if Fig. 1B with pH 5.5 and 10 mM acid. I have a hard time believing that the \\mu(t) underlying these curves are the same across the two batches. Bacteria do not have a single response to every combination of stressors. Small (random) perturbations early on could lead to vastly different transcriptional states and therefore two different functions \\mu(t). My interpretation of these data is the conditions (pH 5.5 and 10 mM acid) are an unstable region of the stress response network. There are actually two \\mu(t), or at least two functions \\delta(t), that the organism is pushed into randomly. Said another way, this is not a batch effect but rather a phase transition.

One way to analyze this is to look at how the function \\delta(t) behaves over the experimental region. Does phenom correctly place a steep transition around pH=5.5 and 10 mM acid? If not, then I would be concerned that the batch effects are hiding large changes in the stress response.

My second, and most significant, point is that the presentation of phenom feels disconnected from the underlying biology. As the authors mention, growth curves are rich sources of data regarding phenotypes and stress responses. The main claim of phenom is that it can dissect different sources of variability to uncover the biologically important signal. The paper shows many examples where batch effects supposedly hide perturbations, but there is no ground truth by which we can assess if the uncovered perturbation functions are correct. Many of these data were collected from other papers that study the underlying biology. It is important to show that the differences between the improved models more accurately reflect the stress response.

For example, the authors note that M_null finds significant perturbations from 10-40 hours, but M_batch reports the significant interval is 20-40 hours (line 220). Which is correct? Without comparing these intervals to the biological stress response, we cannot say which method is better; we can only say that they are different.

Ultimately, phenom will only be useful if it can provide insight into the underlying biology. I suggest the authors place a greater emphasis on how well their method translates OD measurements into latent biology. It would be fantastic if the authors can show how phenom could reduce the need for detailed intracellular characterization of the stress response, or at least identify interesting timepoints for further study.

Minor Comments

I thought the paper was very well written. It was clear and concise.

I suggest changing the title from "microbial phenotypes" to "microbial growth curves". "Phenotypes" is a broad term that includes almost any measurable feature of a microbe, e.g. toxin secretion or sporulation. This paper focuses on only one phenotype, i.e. growth. Growth is a very important phenotype, so I don't believe the more precise title diminishes the work.

The notation can be confusing, although I believe the confusion stems from notational conventions in two different fields rather than poor choices by the authors. The variable \\mu is used in microbiology as the instantaneous growth rate of a microbe. The authors use this convention when describing their parametric models (\\mu_{max}, line 126). Later (line 148) the authors define \\mu(t) as the "average growth behavior of an organism". From Figure 3 and eq (1), we see that this definition of \\mu(t) implies the units of \\mu are population size or [OD]. My suggestion is to change \\mu(t) to some other name. This may be strange for statisticians, but I believe research focused solely on microbial growth curves should not repurpose the standard names for the parameters in a sigmoid growth curve.

**********

Have all data underlying the figures and results presented in the manuscript been provided?

Large-scale datasets should be made available via a public repository as described in the PLOS Computational Biology data availability policy, and numerical data that underlies graphs or summary statistics should be provided in spreadsheet form as supporting information.

Reviewer #1: Yes

Reviewer #2: Yes

Reviewer #3: Yes

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: Yes: Greg Medlock

Reviewer #3: No

Figure Files:

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, . PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at .

Data Requirements:

Please note that, as a condition of publication, PLOS' data policy requires that you make available all data used to draw the conclusions outlined in your manuscript. Data must be deposited in an appropriate repository, included within the body of the manuscript, or uploaded as supporting information. This includes all numerical values that were used to generate graphs, histograms etc.. For an example in PLOS Biology see here: http://www.plosbiology.org/article/info%3Adoi%2F10.1371%2Fjournal.pbio.1001908#s5.

Reproducibility:

To enhance the reproducibility of your results, PLOS recommends that you deposit laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. For instructions, please see

10 Aug 2020

Submitted filename: 2020-08-10-phenom-review-response.pdf

Click here for additional data file.

30 Aug 2020

Dear Dr. Schmid,

We are pleased to inform you that your manuscript 'A Bayesian Non-parametric Mixed-Effects Model of Microbial Growth Curves' has been provisionally accepted for publication in PLOS Computational Biology.

Before your manuscript can be formally accepted you will need to complete some formatting changes, which you will receive in a follow up email. A member of our team will be in touch with a set of requests.

Please note that your manuscript will not be scheduled for publication until you have made the required changes, so a swift response is appreciated.

IMPORTANT: The editorial review process is now complete. PLOS will only permit corrections to spelling, formatting or significant scientific errors from this point onwards. Requests for major changes, or any which affect the scientific understanding of your work, will cause delays to the publication date of your manuscript.

Should you, your institution's press office or the journal office choose to press release your paper, you will automatically be opted out of early publication. We ask that you notify us now if you or your institution is planning to press release the article. All press must be co-ordinated with PLOS.

Thank you again for supporting Open Access publishing; we are looking forward to publishing your work in PLOS Computational Biology.

Best regards,

Jason A. Papin

Editor-in-Chief

PLOS Computational Biology

Jason Papin

Editor-in-Chief

PLOS Computational Biology

***********************************************************

Reviewer's Responses to Questions

Comments to the Authors:

Please note here if the review is uploaded as an attachment.

Reviewer #1: The authors have addressed my concerns regarding their manuscript. I encourage the authors to continue improving the accessibility and implementation of their method through shared computational resources in GitHub. The method must be easy to implement for microbiologists with minimal programming expertise. This will result in considerably more citations for this article, in my opinion.

Reviewer #2: The authors have thoroughly addressed all of my comments.

Reviewer #3: The authors have sufficiently addressed my concerns. The revised paper more accurately describes phenom as a statistical tool, and less-so as a method for understanding biology. I think it will be a useful contribution to the field.

**********

Have all data underlying the figures and results presented in the manuscript been provided?

Large-scale datasets should be made available via a public repository as described in the PLOS Computational Biology data availability policy, and numerical data that underlies graphs or summary statistics should be provided in spreadsheet form as supporting information.

Reviewer #1: Yes

Reviewer #2: Yes

Reviewer #3: Yes

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: Yes: Greg Medlock

Reviewer #3: No

19 Oct 2020

PCOMPBIOL-D-20-00079R1

A Bayesian Non-parametric Mixed-Effects Model of Microbial Growth Curves

Dear Dr Schmid,

I am pleased to inform you that your manuscript has been formally accepted for publication in PLOS Computational Biology. Your manuscript is now with our production department and you will be notified of the publication date in due course.

The corresponding author will soon be receiving a typeset proof for review, to ensure errors have not been introduced during production. Please review the PDF proof of your manuscript carefully, as this is the last chance to correct any errors. Please note that major changes, or those which affect the scientific understanding of the work, will likely cause delays to the publication date of your manuscript.

Soon after your final files are uploaded, unless you have opted out, the early version of your manuscript will be published online. The date of the early version will be your article's publication date. The final article will be published to the same URL, and all versions of the paper will be accessible to readers.

Thank you again for supporting PLOS Computational Biology and open-access publishing. We are looking forward to publishing your work!

With kind regards,

Matt Lyles

PLOS Computational Biology | Carlyle House, Carlyle Road, Cambridge CB4 3DN | United Kingdom ploscompbiol@plos.org | Phone +44 (0) 1223-442824 | ploscompbiol.org | @PLOSCompBiol