A Bayesian approach to analyzing phenotype microarray data enables estimation of microbial growth parameters

Biolog phenotype microarrays (PMs) enable simultaneous, high throughput analysis of cell cultures in di®erent environments. The output is high-density time-course data showing redox curves (approximating growth) for each experimental condition. The software provided with the Omnilog incubator/reader summarizes each time-course as a single datum, so most of the information is not used. However, the time courses can be extremely varied and often contain detailed qualitative (shape of curve) and quantitative (values of parameters) information. We present a novel, Bayesian approach to estimating parameters from Phenotype Microarray data, ¯tting growth models using Markov Chain Monte Carlo (MCMC) methods to enable high throughput estimation of important information, including length of lag phase, maximal \growth" rate and maximum output. We ¯nd that the Baranyi model for microbial growth is useful for ¯tting Biolog data. Moreover, we introduce a new growth model that allows for diauxic growth with a lag phase, which is particularly useful where Phenotype Microarrays have been applied to


Background
Biolog Phenotype Microarrays (PMs) are unique patented commercial products for assessment of cellular respiration of prokaryotic and eukaryotic cells in a wide range of conditions, including metabolism using di®erent carbon, nitrogen, phosphorous and sulphur sources, as well as osmotic, pH, antimicrobial and metal ion stresses. The PMs work by the reduction of a colorless tetrazolium dye in the growth media to a purple formazan by electrons from the NADH produced during cellular respiration. 1 For microbial PMs, 1920 di®erent phenotypes per organism (including controls) can be assessed simultaneously by using the full set of 20 di®erent 96 well plates. As the assays are performed for 24 h or longer, the output is high density time-course data for each well (growth condition), showing a measurement of the quantity of dye reduced. A typical experiment can contain as many as 450,000 data points, making the output especially suitable for mathematical and statistical modelling. The PM platform is°exible, allowing users to construct their own assays using plates, growth media and tetrazolium dyes. Thus, they have proven extremely°exible in terms of the experiments that can be carried out with them. 2 Among other uses, Biolog Omnilog PM technology can be used in research aimed at understanding and controlling the performance of biotechnological processes, through analysis of microbial metabolism in conditions relevant to industrial fermentations. 3,4 However, while PM technology can generate vast amounts of information about microbial growth in the form of time-series data, its usefulness is limited by a lack of robust, easy to use and°exible data analysis tools for such data. Often, the data (which in its raw form comprises up to several hundred data points) is being reduced to either a binary \growth/no-growth" distinction or, at best, a single datum, such as maximum signal reading, average signal height reading or area under curve (AUC). 5,6 Clearly, a lot of information is lost this way.
Biolog's own software allows rudimentary data analysis mainly focused at directly comparing two di®erent strains grown on PM plate types. As detailed in Bochner et al. 5 À À À cf. Box 1 therein À À À this consists of plotting curves from two strains (generally a mutant and a control strain) against each other and highlighting di®erences above a certain threshold. This results in e®ectively a ternary distinction between either of the strains showing higher respiration or there being no signi¯cant di®erence between them. Alternatively the software can give a numerical output of the di®erence in respiration rates between the two strains compared, as used for instance in Ref. 7.
Whether average readout, AUC or endpoint is used, the use of a single value to represent a time-series comprising up to several hundred data points entails losing valuable information about the shape of the underlying curve. For instance, the curves in Fig. 1 all have very similar average readout and AUC. Yet clearly they are qualitatively rather di®erent, with di®erent lengths of lag phases, maximal growth rates and carrying capacities. This is the problem at the heart of Biolog data analysis that we seek to address.
Some previous research has been carried out in this area, but has been limited in its scope. For instance, Ref. 8 focuses on visualizing the raw Biolog data without any parameter estimation. A more recent approach published in Ref. 30 provides both visualization as well as parameter estimation using the gro¯t R package. This package uses nonlinear least squares regression to¯t Gompertz and Richards models to growth curves, but also provides a model-free spline¯t. 10 The logistic growth model has also been found to be e®ective in¯tting Biolog data and has been used to facilitate normalization and data comparison. 11 In this paper, we describe a method to extract further information from these curves. We employ a Bayesian approach using Markov Chain Monte Carlo (MCMC) techniques to sample from the posterior distributions of the parameters of several di®erent growth models¯tted to given Biolog data. Such an approach not only o®ers robust estimates of both best-¯t model parameters but also model-independent characteristics such as lag time, maximal growth rate, and maximal carrying capacity. We have tested these methods using data from customized PM plates, developed to assess the potential of PM technology in the di®erentiation of 100 proprietary brewing yeast strains in di®erent worts. The focus of the paper is on the modeling and methodology and these data are brought as relevant examples.
We anticipate that these ideas have the potential to transform the capacity of research groups to obtain useful and meaningful information from Biolog Phenotype Array data.

Modeling aspects
A number of di®erent growth models have been proposed in the literature. 12 These have di®erent statistical properties 13 and it can be argued that there are limitations of the interpretations of these models. 14 These models are generally derived in the form of an ordinary di®erential equation or as a system of ODEs, though many of them also have a closed-form expression. We will focus on two models in particular, described in the sections below. As the Biolog Omnilog PM machines measure respiration rather than growth, it is not a priori clear that any of the traditional growth models will be able to accurately¯t the data produced in PM experiments. However, empirically we found that these models provide useful results.

The Baranyi model
We chose to focus mainly on the model developed by Baranyi and Roberts. This was rst introduced in Ref. 15, and discussed in more detail for instance in Refs. 16 and 17. The model is based on the Richards model, 18 but introduces another inhibition term to model the lag phase. The inclusion of a lag phase is important, as growth and metabolism of microorganisms in fresh media typically results in an initial period of delayed activity. Consequently, much of the data we analyzed includes a distinct lag phase; this can be seen in the example data shown in Fig. 1. Thus, models that do not include a lag phase (e.g. logistic growth), or models where there is insu±cient°exibility over the shape and time of the lag phase (e.g. Gompertz models 12 ) do not perform well. The form of the model we used is as in Ref. 16: where ðtÞ accounts for the inhibition at the beginning of growth. For an isothermal batch culture environment (e.g. as provided by Omnilog), the authors suggest setting With uðyÞ as in the Richard's model this gives a closed-form expression as follows: where In addition to the four parameters ðy 0 ; y max ; r; mÞ of the Richards model, this introduces another two parameters, and q 0 , that control the length and shape of the lag phase. This term is motivated biologically: q 0 is to be taken as the initial physiological state of the cells, while gives the rate at which they adapt to their new environment. This gives the full Baranyi model a high degree of freedom in accommodating a wide range of growth curve data; as others have noted 13 the Baranyi model performs very well in¯tting empirical growth curve data.
There are two things we would like to point out. First, the form of ðtÞ makes it entirely independent of the rest of the Baranyi model. That is, it is straightforward to incorporate this lag term into other models, and we have done so for instance with a simple diauxic model which we discuss below. Second, for any given lag phase length , there is an in¯nite number of combinations of q 0 and that satisfy the formula for given above. We have found it more convenient to parameterize the Baranyi model (and other models using ðtÞ) with and rather than q 0 and . We then derive q 0 by for use in the closed-form expression of the Baranyi model respectively ðtÞ. In this form, one parameter controls the duration of the lag phase itself, whereas controls its shape. Figure 2(b) shows AðtÞ for a¯xed and varying values of . Higher values of make the lag phase more pronounced, whereas lower values make the e®ect more subtle (in the¯gure, the topmost curve is the one with the lowest value for ). Baranyi suggested a slight simpli¯cation of the original model, setting m ¼ 1 and ¼ r 19 . The resulting model is simpler than the original one and still¯ts most of the data we have seen; On the other hand, the original model is more°exible in¯tting less typical growth curves. We thus adopted a two-fold strategy and used both the original and the simpli¯ed versions of the model, as well as two versions that incorporate either one of the proposed modi¯cations.

Diauxic growth
As will be seen below, the example data we analyze (from worts), typical of growth on complex mixtures of substrates, includes curves that show a clear diauxic e®ect. We modeled this using a simple diauxic growth model based on Monod-type substrate inhibition terms. As a starting point, we used the following model, in ODE form: While s 1 is large, this essentially gives logistic growth on substrate 1. As substrate 1 is being used up, the inhibition term k kþs 1 increases and logistic growth on substrate 2 starts. Smaller values of the inhibition constant k give a more pronounced diauxic e®ect due to the stronger in°uence of s 1 on the inhibition term. k kþs 1 is a hyperbolic inhibition term that is supposed to model the inhibitory e®ects happening within the cells as s 1 is being exhausted. We amended this model to include a Baranyi-like lag phase term. In its nonintegral form ðtÞ this transfers to our ODE model in a straightforward fashion, giving the following for our diauxic growth model: Arguably, this model is a rather simpli¯ed view of the biological processes happening but has empirically proven to be useful. A number of other models for diauxic growth have been proposed. 20,21 Our model is simpler (for example Kompala et al.'s 1986 model has¯ve parameters per substrate, ours has eight in total), but still su±ciently°exible to be able to describe virtually all of the diauxic growth curves we have seen.

Gompertz model
For comparison purposes, we also¯tted a Gompertz model 12 to the data. The Gompertz model does not include a lag phase, so it is useful to determine when the Baranyi model is helpful in identifying lag phase lengths.

Parameter estimation
We adopted a Bayesian approach 22 to infer model parameters from the time-series data recorded by Biolog machines. In particular we used a variant of the Metropolis-Hastings algorithm 23 to sample from the posterior distribution of these parameters. This algorithm starts with an initial state vector ð0Þ . At each iteration, a candidate state vector 0 is generated by drawing from a proposal distribution qð:j i Þ. With probability ð i ; 0 Þ that move is accepted, where If accepted, we set iþ1 ¼ 0 , otherwise set iþ1 ¼ i . The i form a Markov chain, and the stationary distribution of this chain is the desired posterior distribution pðÞ irrespective of the proposal distribution qð:j:Þ in Refs. 24 or 25. A number of modi¯cations and amendments to the original Metropolis-Hastings algorithm have been proposed to better explore the target distribution or to improve the algorithm's rate of convergence. In particular, signi¯cant attention has been drawn to the construction of adaptive algorithms, i.e. algorithms which do not require the speci¯cation of the tuning variance of the proposal distribution qð:j:Þ. One of the most commonly used algorithms is the Adaptive Metropolis (AM) algorithm 26 that has been discussed in detail in the literature [27][28][29][30] .
The use of an adaptive algorithm is important for the analysis of high throughput PM data because of the large number of data sets being analyzed. Each PM plate contains 96 wells, and with an experiment of 100 plates this leads to 9600 separate curves. It is not possible to manually tune the parameters for each curve so an automated, adaptive approach is necessary. We implemented an AM algorithm with global scaling as described in algorithm 4 in Andrieu and Thoms 2008 as follows 2 : (1) For an initial segment of the chain (i < i 0 , for some sensible choice of i 0 ) perform a \Random-Walk Metropolis with global scaling" with proposal covariance matrix i AE 0 , with AE 0 an initial \best guess" of the true covariance matrix. (2) For the remainder of the chain, do as above but use i AE i as the covariance matrix for the proposal distribution, where AE i is the sample covariance matrix of the previous history of the chain, and i is a varying scaling factor. We update AE i iteratively.
We used i 0 ¼ 2500 (chosen empirically). For i < i 0 we updated AE i and i only when we accept a move, as suggested in Haario et al. 2001. 26 We reset i to its original value of 2:4 2 = dimðÞ in iteration i 0 , where is the parameter vector. 2 We also added I to AE i in each step to keep it from becoming singular for some small 26 . If ever it still became singular due to rounding errors, we added I to it until it was not singular anymore. We did not use any form of thinning, always discarded the initial 50,000 iterations as burn-in, and in production runs used 500,000 iterations in total. We used uninformative or uniform priors on suitable regions, and utilized simple heuristics to give rough estimates of the initial parameter vector for the AM algorithm (see Appendix for details). These heuristics proved to be e®ective in that all Markov chains converged to a stationary distribution (see Sec. 2.6 on Performance below) after appropriate burn-in. Alternative approaches could be to use an iterative least-squares approach to obtain initial parameter values near to the stationary distribution. 11 To calculate the likelihood function, we assumed normally i.i.d. measurement errors, which we heuristically estimated from the raw data as detailed in the appendix.

Model choice
In addition to inferring parameters from them, we used the Deviance Information Criterion. 31,32 This is de¯ned as follows: where D is the deviance de¯ned by and " D is the mean of this deviance. The model complexity p V is given as Þ to account for model complexity. We have found however that this is highly dependent on the quality of the estimation of the posterior. In some cases where our posterior sample was not a good estimate, we saw the p D term dropping to arti¯cially low (often negative) values leading to a bad, arti¯cially low estimate of DIC. While it is possible to check that the posterior sample is reasonable before calculating DIC this way, we expect that in any high-throughput environment there will be individual data that slip through such checks. Thus, we have found that using p V was in practice more suitable for our purposes. For comparison purposes we have also calculated BIC for all models, 33 using the highest likelihood observed in the posterior sample for the BIC estimate.
We used the DIC to choose which model to use to derive estimated parameters for each well. Furthermore, we were looking to get an indication as to whether our models provide a reasonable¯t at all. To this end we compared them to a simple \dummy" model de¯ned by yðtÞ ¼ c for a constant parameter c. This is noninformative, and similar to a classical H 0 model. A comparison via DIC or BIC to the constant model is meant to show whether the model¯ts the data at all.

Parameter comparisons between models
One potential issue with¯tting a number of di®erent models to data is that parameters of one model do not necessarily relate easily to those of another. Therefore, in addition to the model-speci¯c parameters we used model-independent measurements of three key growth characteristics, which we derived numerically from the¯tted curves. This is to allow for simple quantitative comparisons between wells that show qualitatively di®erent behavior.
In the following paragraphs, we will take yðt; Þ to mean the signal level predicted by our model and the parameter vector at time t. First, for A, the maximum coloration change achieved, we simply used yðt last ; Þ À yð0; Þ, that is, the absolute increase in coloration our model predicts over the period of time that was recorded in the experiment. It may be argued that for all the models we are using we could also compute a similar A from the y 0 and y max (respectively y 0 and s 1 , s 2 ). However, we found that if the recorded respiration data stops before a maximum is attained, the maximum or substrate level parameters are e®ectively inestimable, and so would A be with such a de¯nition. Note that if a maximum is e®ectively reached within the experiment record, these two de¯nitions will for practical purposes be equivalent.
For max we took the steepest slope of yðt; Þ (again within the experiment period), as derived numerically from the¯tted curve. This is not a transformation of the model's rate parameter(s) alone, but is model-independent and coincides with a more natural de¯nition of the maximal respiration rate.
For the lag time L an essentially model-independent de¯nition is slightly trickier, and a number of possible de¯nitions could be used. 14 We de¯ne L to be the t-coordinate of the intersection of the tangent at the steepest point of yðt; Þ with the line y yð0; Þ. 13 For y 0 $ y max this will approximate the lag parameter in ðtÞ almost exactly. Figure 2(a) illustrates our de¯nition of the lag phase length. The tangent of the steepest point is shown in blue, and the derived lag length L in red. We note that for y 0 ( y max this estimate will di®er substantially from the lag parameter . However, so long as the relative di®erence between initial and¯nal cell concentration is constant across a dataset, our estimate of L will still be comparable between wells and, crucially, between di®erent models.

Identifying presence and absence of growth
In order to identify whether a given well exhibits signi¯cant levels of growth at all, we compared the maximum coloration attained A to a 95% quantile of the same parameter for a control well present on each plate.

Data preprocessing
In a small number of cases we observed anomalous behavior where coloration actually decreases signi¯cantly after attaining a maximum. The causes of this behavior are as yet unknown (in principle the reduction of the tetrazolium dye should be a oneway process), but we still want to be able to extract meaningful information from the data. We used a simple heuristic to remove the aberrant parts of the data. More precisely, we looked at the maximum y 0 max (attained at t max ) of a smoothed curve as above (using a 9-point window). We then removed the tail of the data series if in any interval ½t max ; t; t > t max at least 90% of data points (of the unsmoothed curve) were at least 0.5 standard deviations (using the numerically estimated measurement error) below y max . We would then remove the tail of the data series after the least such t. However, we always left at least the initial 40 data points.

Our models¯t the data
The model was applied to 40 Biolog arrays comprising a total of 3840 time courses. Of these, 2989 exhibited signi¯cant growth compared to a known control well. According to DIC 598 of these wells were¯tted best by the Baranyi model, 2382 by the diauxic model, 9 by the Gompertz model and none by the constant dummy model. The remaining 851 wells exhibited no signi¯cant growth, and model¯t between the Baranyi, diauxic and Gompertz models was almost arbitrary.
According to BIC the picture is similar, with a slight bias against the diauxic model (2160 diauxic cases with growth) and toward simpler models (793 Baranyi, 36 Gompertz). Again no wells were¯tted by the constant dummy model.
Furthermore, in all cases in which the Gompertz model was preferred over Baranyi or diauxic models, DIC scores of these models were very close to each other (within 1% of absolute values). The converse was not the case. Best-¯t models usually achieved a DIC score in the range of 200-400. For cases where the Gompertz model scored best, the mean di®erence between the Gompertz model and the second-best model was 1.7; Conversely, where the Gompertz model did not score best, the mean di®erence in BIC score between the best model and the Gompertz model was 658.
This demonstrates that the Baranyi model, and the diauxic model with Baranyi lag phase, are e®ective tools for the analysis and interpretation of Biolog PM data.
Example¯ts where the Baranyi model and the diauxic model with Baranyi lag¯t best are shown in Figs. 3(a) and 3(b), respectively. Plotted also is a¯tted Gompertz model, demonstrating the relevance of the models we have used. In these cases, the Baranyi and diauxic models receive much the best BIC scores (Table 1). Figure 3(c) shows one well with no growth, and a slight decrease in coloration, together with¯ts by the Baranyi and constant models.
Further evidence for the value of the Baranyi model and Baranyi lag term in the diauxic model in¯tting this type of data can be seen from the relative wide spread of the estimated mean for the curvature parameter of the lag phase, in both models, as well as the curvature parameter m in the Baranyi model (Fig. 4). Thus, we consider it likely that any model lacking such an extra parameter would fail to accommodate the range of curve shapes we have encountered.     (e)

Identi¯cation of curve features beyond AUC
There are two main features of our approach: the¯rst, as presented above, is that we are able to¯t curves that are qualitatively diverse, namely nondiauxic and diauxic growth. The second is that we are able to estimate key relevant parameters from Biolog data, in particular a length of lag phase, maximal growth rate and maximal output. In many biotechnological applications, the identi¯cation of strains or conditions with minimal lag phase or maximal growth could be particularly important, as such strains or conditions could speed production and/or cut costs. Figures 5(a Table 2 lists the mean derived parameters for the six di®erent strains shown in the beginning of the paper, clearly showing that AUC alone is not su±cient to allow comparisons between qualitatively di®erent data. The lag phases range from 12.3 h (P103-C01) to 38.9 h (P101-C01) and the maximal growth rates range from 11.7 per hour (P102-B12) to 35.0 per hour (P101-C01). Both of these model-derived parameters vary three-fold among the strains shown, demonstrating that our approach is considerably superior to using AUC.
We have found that a comparison to the 95% quantile of maximum coloration of a control well works reliably in identifying whether growth occurs.

The estimated parameters and model-independent measurements are reliable
The Bayesian approach has allowed us to use the posterior distributions to estimate standard deviations of the individual model parameters as well as the numerically derived model-independent growth measurements. These are generally low (typically < 1% of the parameter values) and are elevated only in cases where some parameters are not estimable from the data, e.g. estimation of y max when¯tting a Baranyi model to data that is cut o® before stationary phase is reached (standard deviation as high as 10% of the parameter), but also sometimes s 1 and s 2 when¯tting the diauxic model to data clearly depicting simple growth. We have found that A (the maximum coloration change) very closely matches the respective model parameter(s), and that max has given good comparability of results between wells even with di®erent best-¯t models. In cases of diauxic growth, particularly with small s 2 and the bulk of growth on s 1 (as is the case in the vast majority of diauxic cases observed) the Baranyi model would sometimes give an arti¯cially low estimate of max as it tries to compensate for the (in a sense slower) two-step approach to peak coloration by¯tting an altogether slower curve. Figure 3(e) shows one such case together with¯ts by the Baranyi and diauxic models.
The lag time L in the vast majority of cases similarly gives good comparability between wells and models. In one atypical case 3(f), a diauxic growth curve features a higher growth rate on s 2 than on s 1 . By our de¯nition of L as the time coordinate of the intersection of the steepest tangent with the°at line y y 0 , the lag time in such cases is determined by the growth on s 2 and could be signi¯cantly later than start of growth on s 1 . Arguably this is not what we want, as it would not be consistent with a physiological de¯nition of the lag phase, since the lag time identi¯ed includes the time grown on s 1 . On the other hand, in some circumstances this output may be preferential, e.g. in an industrial research application a slow initial growth on s 1 may be of considerably less interest than the time until peak growth rate. It would be possible in our framework to derive a specialized de¯nition of L tailored to diauxic models, e.g. using the steepest tangent before the point where s 1 is e®ectively depleted.
While our descriptive measurements L, max and A give us consistent de¯nitions and readily comparable results for the vast majority of cases, the various model parameters we are also inferring allow us a more explanatory analysis of subsets of (or even individual) cases. For instance, within a set of diauxic curves, our estimates of s 1 and s 2 allow us to identify cases in which the bulk of growth is on the second substrate. Again, in industrial research applications this may be of particular interest. Similar analysis could be carried out e.g. on the curvature parameters of the Baranyi model, to identify outliers in lag phase behavior.
As a further test of reliability, we have compared parameter estimates from wells measuring the same conditions: each condition on these particular phenotype arrays  Fig. 3, this provides three independent estimates for each of the three derived parameters for six models, a total of 18 comparisons. The parameter estimates were generally very consistent, with median percentage error of 5:8%. The worst case is for Fig. 3(d), where there is no growth, and the parameter estimates vary by approximately 30%; however, because there is no growth, this is not a problem. Full details of these comparisons are provided in Appendix D.

Performance
The MCMC methodology based on an Adaptive Metropolis algorithm performed well in combination with the models we used. We tested all the Markov chain outputs for the model¯ts that appear in all¯gures for convergence using the Heidelberger and Welch's convergence test as implemented in the CODA package in R (https://cran.r-project.org/web/packages/coda/index.html). The results from this test showed that all the chains have passed the test, i.e. have converged to the stationary distribution. The outputs are not especially interesting (many tables of nonsigni¯cant p-values) so we have placed one example in Appendix E and have not included the other outputs. The adaptation to di®erent target distributions in particular worked very well and our algorithm required no manual¯ne-tuning to explore di®erent data with high e±ciency. Figure 6 shows an initial segment of one AM chain we ran. It is clearly visible that the mixing of the chain improves rapidly after the initial discovery phase. Running our methodology implemented in C++ for a single plate (comprising 96 wells recorded every 15 min over 72 h) took around 50-60 min on a quad-core machine (Intel Xeon E3 1230v2).  6. Initial segment from the trace plot for for one AM chain. It is easy to see that the initial acceptance rate is far from optimal, but rapidly improves after the¯rst 2500 iterations.

Comparison to previous approaches
Vaas et al. 2012 discuss one particular case (their Figs. 4(a) and 4(b)) where their model-¯tting approach fares worse than their alternative spline-based method. It appears to us that the curve in question is simply diauxic in nature, and we expect that our approach would be able to¯t this time course and identify the behavior as diauxic. The second example they discuss (Fig. 4(c) same curve) could equally likely be¯tted by our methodology with an appropriate model. In either case Ref. 9 approach likely could not be extended easily to encompass such additional models, using a third-party software package to do the actual parameter estimation.

Conclusion
We have presented a Bayesian approach to estimating curve-parameter information from respiration data gathered in phenotype microarray experiments. We had aimed to extract meaningful information from complex respiration curves in a statistically sound way, and have shown that our approach succeeds in doing so.
Our solution has several key advantages. First, the Bayesian framework in which our approach operates a®ords us a great deal of°exibility in extracting information from respiration curves. As we are approximating the joint posterior distribution of all model parameters, we are free to perform further analysis on this and can for instance derive the distribution of arbitrary functions of our model parameters from them. Model choice criteria allow us to readily identify qualitatively di®erent behavior in the data on top of quantitative measurements. This gives us a \best of two worlds" solution encompassing both uni¯ed model-independent descriptive measurements as well as more explanatory model-speci¯c ones.
Second, the Baranyi model provides an excellent model to¯t a wide range of observed growth curves. The only exception we encountered were curves exhibiting diauxie. In these cases, our diauxic model provided us with good¯ts in these cases as well. In particular, these two models outperform traditional growth models in many cases. Our modular implementation allows us to adapt our code to new types of experimental data with minimal e®ort. As PM machinery is being used in a wide range of di®erent areas, this allows our solution to remain applicable as new types of data become available. Third, the model-independent measurements L, max and A we derive from thē tted curves allow an immediate comparison of key values between di®erent time series, even if qualitatively di®erent models were used to¯t the data. This allows us to match simpler (e.g. spline-based) approaches in their ability to perform quantitative comparisons between a wide range of growth behaviors. However in addition, our approach retains information from the richer model-dependent parameters. This gives us powerful tools to perform further analysis on speci¯c subsets of data. In particular, we can use these to compare wells with qualitatively similar behavior in more detail, for instance utilization of di®erent substrates in diauxic curves.
Last, our implementation is suited to high-throughput analysis of PM data. Compared to previous approaches to data analysis on PM experiments, we believe that our solution o®ers improvements in several areas. Even the relatively small number of models we have implemented so far allows our solution to successfully¯t a wide range of data.
We believe that our extensible approach is uniquely suited to \keep up" with PM data, as PM machinery is being used for an increasingly wide¯eld of di®erent types of experiments. As PM machinery continues to¯nd novel usage scenarios, the modular implementation we have chosen will allow our solution to remain°exible in accommodating data from these scenarios.

Program Code Availability
We have made the program code available on GitHub at URL https://github.com/ dovstekellab/mcmc-pma.git under license GNU GPLv3. We also include instructions on how to compile and run the code, as well as how to interpret the results¯les.
(2) For the lag parameter, we¯t lines to every 20-point interval of data points, and intersected the maximum-slope line with y ¼ y 0 , with y 0 derived as in the previous line. (3) For the growth rate r, we used a simple search heuristic to guess an initial value that we have found empirically to be e®ective. We start from a su±ciently large interval of possible values ½r min ; r max , and divided this into 10 equal parts r 0 ¼ r min ; r 1 ; . . . ; r 11 ¼ r max . We then compared the maximum slope of the modeled data with r ¼ r 1 ; . . . ; r 10 with the maximum slope as in the previous line. For r i giving the smallest di®erence in maximum slope, we recurse by setting r min ¼ r iÀ1 ; r max ¼ r iþ1 . We repeat this 100 times. (4) For the remaining parameters we used the same search algorithm, except we compared the sum of squares of di®erences between modeled and observed data instead of maximum slope.
Speci¯c to each of the models we did the following: We used ½0:01; 1, ½2; 10, ½0:3; 1 as initial intervals for r, m, in this search algorithm. We¯rst called this for r, setting both remaining parameters to 1, then similarly for the remaining ones. For the diauxic models, we proceed similarly, except we additionally needed to¯nd s 1 and s 2 , respectively the division of the total growth between the two. To do so, we found the maximum slope on a smoothed curve as in our estimation for the lag parameter, and then the minimum slope between that point and when the smoothed curve¯rst comes within 12 standard deviations of the maximum. (In other words, we looked for the characteristic intermittent slowing-down of growth.) We then subtracted y 0 from the value at this point, and took this as the initial guess for s 1 . For the remaining parameters we proceeded similarly as in the Baranyi model, using ½10 À8 ; 10 À4 ; ½0:0 5; 1; ½0:0001; 11 as initial intervals for r 1 , , k 1 in the search algorithm. We always use 0.0003 for r 2 , as performing our search heuristic for this parameter did not improve results.

Appendix C. Estimation of Measurement Error from Data
We took a smoothed curve (taking the average of every nine adjacent measurements; without taking logarithms) as reference to numerically estimate the variance. We however always assumed a minimum variance of 5 Biolog units. That is, we de¯ned ; :

Appendix D. Consistency Analysis
Consistency of parameter estimates was tested for the six best-¯t models shown in Fig. 3. These particular phenotype microarrays have triplicate wells for each condition, so the derived model parameters shown were compared with those for the two other triplicate wells. The full output is:

Appendix E. Convergence Tests
Convergence for all of the Markov chains for model¯ts used in the¯gures was carried out using the Heidelberger and Welch's convergence test as implemented in the CODA package in R (https://cran.r-project.org/web/packages/coda/index.html). All of the chains passed the test, i.e. demonstrating convergence to a stationary distribution. The output for each chain is in the form of a table with a p-value for each parameter. An example from one chain that is typical of all output is: Each row, labeled V1 to V8, represents one variable. The null hypothesis is that the chain is from a stationary distribution and so it can be seen that the null