3.1.1 Automated cultivation platform.

All experiments were conducted on a robotic platform with an integrated small-scale cultivation system. In our setup, a BioLector Pro microbioreactor system (m2p-labs GmbH, Baesweiler, Germany), is integrated into a Tecan Freedom EVO liquid handling robot (Tecan, Männedorf, Switzerland). The BioLector Pro is a device to quasi-continuously observe biomass, pH and dissolved oxygen (DO) during cultivation of microorganisms in specialized microtiter plates (MTPs). These rectangular plates comprise multiple reaction cavities called “wells”, usually with volumes in microliter or milliliter scale. The BioLector allows to control temperature and humidity while shaking the plates at adjustable frequencies between 800 and 1500 rpm.

The liquid handler, which allows to take samples for at-line measurements during cultivation, is surrounded by a laminar flow hood to ensure sterile conditions for liquid transfer operations. Next to the BioLector Pro, various other devices are available on the platform, including an Infinite M Nano+ microplate photometer (Tecan, Männedorf, Switzerland), a cooling carrier and a Hettich Rotanta 460 robotic centrifuge (Andreas Hettich GmbH & Co. KG, Tuttlingen, Germany). The overall setup is similar to the one described by Unthan et al. 2015 [8]. The automated platform enables to perform growth experiments with different microorganisms, to autonomously take samples of the running process and to perform bioanalytical measurements, e.g. quantification of glucose. It is thus a device for miniaturised, automated bioprocess cultivation experiments.

In this work, we used specialized 48-well deepwell plates called FlowerPlates (MTP-48-B, m2p-labs GmbH, Baesweiler, Germany) that are commonly used for cultivations in the BioLector device. The plates are pre-sterilized and disposable and were covered with a gas-permeable sealing film with a pierced silicone layer for automation (m2p-labs GmbH, Baesweiler, Germany). The biomass was quasi-continuously detected via scattered light [24] at gain 3 with 4 minutes cycle time to obtain backscatter measurements. DO and pH were not measured since they are not relevant for the application examples. Both cultivation and biomass calibration experiments were conducted in the BioLector Pro at 30°C, 3 mm shaking diameter, 1400 rpm shaking frequency, 21% head-space oxygen concentration and ≥ 85% relative humidity.

3.1.2 Strain, media preparation and cell banking and cultivation.

The wild-type strain Corynebacterium glutamicum ATCC 13032 [25] was used in this study. If not stated otherwise, all chemicals were purchased from Sigma–Aldrich (Steinheim, Germany), Roche (Grenzach-Wyhlen, Germany) or Carl Roth (Karlsruhe, Germany) in analytical quality.

Cultivations were performed with CGXII defined medium with the following final amounts per liter of distilled water: 20 g D-glucose, 20 g (NH₄)₂SO₄, 1 g K₂HPO₄, 1 g KH₂PO₄, 5 g urea, 13.25 mg CaCl₂ · 2 H₂O, 0.25 g MgSO₄ · 7 H₂O, 10 mg FeSO₄ · 7 H₂O, 10 mg MnSO₄ · H₂O, 0.02 mg NiCl₂ · 6 H₂O, 0.313 mg CuSO₄ · 5 H₂O, 1 mg ZnSO₄ · 7 H₂O, 0.2 mg biotin, 30 mg protocatechuic acid. MOPS were used as buffering agent and the pH was adjusted to 7.0 using 4 M NaOH.

A working cell bank (WCB) was prepared from a shake flask culture containing 50 mL of the described CGXII medium and 10% (v/v) brain heart infusion (BHI) medium (). It was inoculated with 100 μl cryo culture from a master cell bank stored at -80°C. The culture was incubated for approximately 16 hours in an unbaffled shake flask with 500 ml nominal volume at 250 rpm, 25 mm shaking diameter and 30°C. The culture broth was then centrifuged at 4000 × g for 10 minutes at 4°C and washed once in 0.9% (w/v) sterile NaCl solution. After centrifugation, the pellets were resuspended in a suitable volume of NaCl solution to yield a suspension with an optical density at 600 nm (OD₆₀₀) of 60. The suspension was then mixed with an equal volume of 50% (w/v) sterile glycerol, resulting in cryo cultures of OD₆₀₀≈30. Aliquots of 1 mL were quickly transferred to low-temperature freezer vials, frozen in liquid nitrogen and stored at -80°C.

3.1.3 Algorithmic planning of dilution series.

All calibration experiments require a set of standards (reference samples) with known concentrations, spanning across sometimes multiple orders of magnitude. Traditionally such standards are prepared by manually pipetting a serial dilution with a 2x dilution factor in each step. This can result in a series of standards whose concentrations are evenly spaced on a logarithmic scale. While easily planned, a serial dilution generally introduces inaccuracies that accumulate with an increasing number of dilution steps. It is therefore desirable to plan a dilution series of reference standards such that the number of serial dilution steps is minimized.

To reduce the planning effort and allow for a swift and accurate preparation of the standards, we devised an algorithm that plans liquid handling instructions for preparation of standards. Our DilutionPlan algorithm considers constraints of a (R × C) grid geometry, well volume, minimum and maximum transfer volumes to generate pipetting instructions for human or robotic liquid handlers.

First, the algorithms reshapes a length R ⋅ C vector of sorted target concentrations into the user specified (R × C) grid typically corresponding to a microtiter plate. Next, it iteratively plans the transfer and filling volumes of grid columns which are subject to the volume constraints. This column-wise procedure improves the compatibility with multi-channel manual pipettes, or robotic liquid handlers. Diluting from a stock solution is prioritized over the (serial) dilution from already diluted columns. The result of the algorithm are (machine-readable) pipetting instructions to create R ⋅ C single replicates with concentrations very close to the targets. We open-sourced the implementation as part of the robotools library [26].

As the accuracy of the calibration model parameter estimate increases with the number of calibration points, we performed all calibrations with the maximum number of observations that the respective measurement system can make in parallel. The calibration with 96 glucose and 48 biomass concentrations is covered in the following chapters.

3.1.4 Glucose assay calibration.

For the quantification of D-glucose, the commercial enzymatic assay kit “Glucose Hexokinase FS” (DiaSys Diagnostic Systems, Holzheim, Germany) was used. For the master mix, four parts buffer and one part enzyme solution were mixed manually. The master mix was freshly prepared for each experiment and incubated at room temperature for at least 30 minutes prior to the application for temperature stabilization. All other pipetting operations were performed with the robotic liquid handler. For the assay, 280 μL master mix were added to 20 μL analyte in transparent 96-well flat bottom polystyrol plates (Greiner Bio-One GmbH, Frickenhausen, Germany) and incubated for 6 minutes, followed by absorbance measurement at 365 nm. To treat standards and cultivation samples equally, both were diluted by a factor of 10 (100 μL sample/standard + 900 μL diluent) as part of the assay procedure.

As standards for calibration, 96 solutions with concentrations between 0.075 and were prepared from fresh CGXII cultivation medium (Section 3.1.2) with a concentration of D-glucose. The DilutionPlan algorithm (Section 3.1.3) was used to plan the serial dilution procedure with glucose-free CGXII media as the diluent, resulting in 96 unique concentrations, evenly distributed on a logarithmic scale. Absorbance results from the photometer were parsed with a custom Python package (not published) and paired with the concentrations from the serial dilution series to form the calibration dataset used in Section 4.1.2. 83 of the 96 concentration/absorbance pairs lie below and were used to fit a linear model in Section 4.1.1.

3.1.5 Biomass calibration.

Calibration data for the biomass/backscatter calibration model (Section 4.1.2) was acquired by measuring 48 different biomass concentrations at cultivation conditions (Section 3.1.2) in the BioLector Pro. 100 mL C. glutamicum WT culture was grown overnight on glucose CGXII medium (Section 3.1.2) in two unbaffled 500 mL shake flasks with 50 mL culture volume each (N = 250 rpm, r = 25 mm). The cultures were harvested in the stationary phase, pooled, centrifuged and resuspended in 25 mL 0.9% (w/v) NaCl solution. The highly concentrated biomass suspension was transferred into a magnetically stirred 100 mL trough on the liquid handler for automated serial dilution with logarithmically evenly spaced dilution factors from 1× to 1000×. The serial dilution was prepared by the robotic liquid handler in a 6 × 8 (48-well square) deep well plate (Agilent Part number 201306–100) according to the DilutionPlan (Section 3.1.3). 6x 800 μL of biomass stock solution were transferred to previously dried and weighed 2 mL tubes, immediately after all transfers of stock solution to the 48 well plate had occurred. The 2 mL tubes were frozen at -80°C, lyophilized over night, dried again at room temperature in a desiccator over night and finally weighted again to determine the biomass concentration in the stock solution.

After a column in the 48 well plate was diluted with 0.9% (w/v) NaCl solution, the column was mixed twice by aspirating 950 μL at the bottom of the wells and dispensing above the liquid surface. The transfers for serial dilutions (columns 1 and 2) and to the 48 well FlowerPlate were executed right after mixing to minimize the effects of biomass sedimentation as much as possible. The FlowerPlate was sealed with a gas-permeable sealing foil (product number F-GP-10, m2p-labs GmbH, Baesweiler, Germany) and placed in the BioLector Pro device. The 1 h BioLector process for the acquisition of calibration data was programmed with shaking frequency profile of 1400, 1200, 1500, 1000 rpm while maintaining 30°C chamber temperature and measuring backscatter with gain 3 in cycles of 3 minutes.

The result file was parsed with the bletl Python package [27] to extract backscatter measurements made at 1400 rpm shaking frequency. A log(independent) asymmetric logistic calibration model was fitted as described in Section 4.1.2. The linear calibration model for comparison purposes (Section 4.2.3) was implemented with its intercept fixed to the background signal predicted by the asymmetric logistic model (). It was fitted to a subset of calibration points approximately linearly spaced at 30 different biomass concentrations from 0.01 to .

3.1.6 Microbial growth experiment.

Cultivations with C. glutamicum were performed in the automated cultivation platform (Section 3.1.1) under the described conditions. CGXII medium with glucose and without BHI was used as cultivation medium. To start the growth experiments, the liquid handler was used to transfer 20 μL of a WCB aliquot into the first column of FlowerPlate wells, which were pre-filled with 780 μL medium. These wells were run as a preculture. When precultures reached a backscatter readout of 15, which corresponds to a cell dry weight of approximately , the inoculation of the main culture wells was triggered. 780 μL medium were distributed into each main culture well (columns 2–8) and allowed to warm up for approximately 15 minutes. Preculture wells A01 and B01 were pooled and 20 μL culture broth was transferred to each main culture well, resulting in 800 μL final volume. The theoretical biomass concentration at the start of the main cultures is accordingly. This strategy was used to avoid a lag-phase with non-exponential growth.

Backscatter measurements of biomass concentration were acquired continuously, while main culture wells were harvested at predetermined time points to measure glucose concentrations in the supernatant. The time points were chosen between 0 and 15 hours after the inoculation of main cultures to cover all growth phases. For harvesting, the culture broth was transferred to a 1 mL deep-well plate by the liquid handler. The plate was centrifuged at 3190 × g at 4°C for 5 minutes and the supernatant was stored on a 1 mL deep well plate chilled to 4°C. The glucose assay was performed after all samples were taken.

3.2.1 Asymmetric logistic function.

The asymmetric, five-parameter logistic function (also known as Richard’s curve) was previously shown to be a good model for many applications [39], but it is often defined in a parameterization (Eq 10) that is non-intuitive. Some parametrizations even introduce a sixth parameter to make the specification more intuitive, but this comes at the cost of structural non-identifiability [40, 41]. Furthermore, in the most commonly found parametrization (Eq 10), one parameter is constrained to be strictly positive. We also found that structural non-identifiability between the parameters makes it difficult to define an initial guess and bounds to reliably optimize a model based on this parametrization. (10)

To make calibration model fitting more user-friendly, we reparameterized the commonly used form such that all five parameters are intuitively interpretable and structurally independent (Fig 5). With our reparameterization (Eq 11), the 5-parameter asymmetric logistic function is parameterized by lower limit , upper limit , inflection point x-coordinate , slope at inflection point and an asymmetry parameter . At c = 0, the y-coordinate of the inflection point lies centered between L_L and L_U. I_y moves closer to L_U when c > 0 and accordingly closer to L_L when c < 0 (Fig 5, black and gray curves). An interactive version of Fig 5 can be found in a Jupyter Notebook in the calibr8 GitHub repository [37].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Reparametrized asymmetric logistic function.

When parametrized as shown in Eq 11, each of the 5 parameters can be manipulated without influencing the others. Note that, for example, the symmetry parameter c can be changed without affecting the x-coordinate of the inflection point (I_x), or the slope S at the inflection point (gray vs. black).

https://doi.org/10.1371/journal.pcbi.1009223.g005

For readability and computational efficiency, we used SymPy [42] to apply common subexpression elimination to Eq 11 and our implementation respectively (S1 File). The step wise derivation from Eqs 10 to 11 is shown in S1 File and in a Jupyter Notebook in the calibr8 GitHub repository [37]. (11)

3.2.2 calibr8 package for calibration models and modular likelihoods.

With calibr8 we present a lightweight Python package that specializes on the definition and modular implementation of non-linear calibration models for calibration and modeling workflows.

The calibr8 application programming interface (API) was designed such that all calibration models are implemented as classes that inherit from calibr8.CalibrationModel, which implements properties and methods that are common to all calibration models (Fig 6). The common interface simplifies working with calibration models in a data analysis or modeling workflow. For a detailed explanation of classes and inheritance in calibr8 we refer the interested reader to the documentation [22]. For example, the CalibrationModel.objective can be used to create objective functions to optimize the model parameters. The objective relies on the loglikelihood method to compute the sum of log-likelihoods from independent and dependent variables. It uses the predict_dependent method internally to obtain the parameters of the probability distribution describing the dependent variables, conditioned on the independent variable.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. calibr8 class diagram.

All calibr8 models inherit from the same CalibrationModel and DistributionMixin classes that define attributes, properties and method signatures that are common to all calibration models. Some methods, like loglikelihood() or objective() are implemented by CalibrationModel directly, whereas others are implemented by the inheriting classes. Specifically the predict_* methods depend on the choice of the domain expert. With a suite of Base*T classes, calibr8 provides base classes for models based on Student-t distributed observations. A domain expert may start from any of these levels to implement a custom calibration model for a specific application.

https://doi.org/10.1371/journal.pcbi.1009223.g006

Through its inheritance-based design, the calibr8.CalibrationModel gives the domain expert full control over the choice of trend functions and probability distributions. Conveniently, calibr8 already implements functions such as polynomial, logistic and asymmetric_logistic, as well as base classes for commonly found noise models. By leveraging these base models, the implementation of a user-defined calibration model reduces to just a few lines of code (Box 1 Code 1 and Box 2 Code 2). In the current version we implemented base classes for continuous univariate and multivariate calibration models and anticipated models with discrete independent variables.

The implementations depicted in Fig 6 are fully compatible with aesara.Variable inputs, resulting in TensorVariable outputs. Aesara is a graph computation framework that auto-differentiates computation graphs written in Python and compiles functions that evaluate with high performance [43]. This way, the loglikelihood function of a CalibrationModel can be auto-differentiated and compiled, thus facilitating efficient computation with optimization or gradient-based MCMC sampling algorithms (Section 3.2.6). For more details about the implementation, please refer to the documentation and code of the calibr8 package [37].

Convenience features. To facilitate modeling workflows, calibr8 implements convenience functions for optimization (fit_scipy, fit_pygmo) and creation of diagnostic plots (calibr8.plot_model) as shown in Section 2.4 and Section 4.1.2. As explained in Section 2.4the residual plot on the right of the resulting figure is instrumental to judge the quality of the model fit.

Standard properties of the model, estimated parameters and calibration data can be saved to a JSON file via the CalibrationModel.save method. The saved file includes additional information about the type of calibration model and the calibr8 version number (e.g. S3 File) to support good versioning and data provenance practices. When the CalibrationModel.load method is called to instantiate a calibration model from a file, standard properties of the new instance are set and the model type and calibr8 version number are checked for compatibility.

3.2.3 Numerical inference.

To numerically infer the posterior distribution of the independent variable, given one or more observations, infer_independent implements a multi-step procedure. The three outputs of this procedure are a vector of posterior probability evaluations, densely resolved around the locations of high probability mass, and the bounds of the equal-tailed as well as the highest-density intervals (ETI, HDI) corresponding to a user-specified credible interval probability.

In the first step, the likelihood function is integrated in the user-specified interval [lower, upper] with scipy.integrate.quad. Second, we evaluate its cumulative density function (CDF) at 10 000 locations in [lower, upper] and determine locations closest to the ETI^99.999%. Next, we re-evaluate the CDF at 100 000 locations in the ETI^99.999% to obtain it with sufficiently high resolution in the region of highest probability. Both ETI and HDI with the (very close to) user-specified ci_prob are obtained from the high-resolution CDF. Whereas the ETI is easily obtained by finding the CDF evaluations closest to the corresponding lower and upper probability levels, the HDI must be determined through optimization (Eq 12). (12)

3.2.4 murefi package for building multi-replicate ODE models.

Models of biochemical processes are traditionally set up to describe the temporal progression of an individual system, such as a reaction vessel. Experimental data, however, is commonly obtained from multiple reaction vessels in parallel, often run under different conditions to maximize information gain. This discrepancy between the prototypical model of the biological system and the heterogeneous experimental data to be fitted is typically resolved by replicating the model for all realizations of the biological system in the dataset. Along with the replication of the model, some model parameters may be kept global, while others can be local to a subset of the replicates, for example due to batch effects or different start conditions.

With a Python package we call murefi (multi-replicate fitting), we implemented data structures, classes and auxiliary functions that simplify the implementation of models for such heterogeneous time series datasets. It seamlessly integrates with calibr8 to construct likelihood-based objective functions for optimization or Bayesian inference. To enable the application of efficient optimization or sampling algorithms, the use of automatic differentiation to obtain gradients of the likelihood w.r.t. input parameters is highly desirable. Various methods for automatic differentiation of ODE models are available, but their efficiency is closely connected to the implementation and size of the model [44]. In murefi we implemented support for sunode [45], a recently implemented Python wrapper around the SUNDIALS suite of nonlinear and differential/algebraic equation solvers [46]. When used in the context of a PyMC model, a process model created with calibr8 and murefi can therefore be auto-differentiated, solved, optimized and MCMC-sampled with particularly high computational efficiency.

Structure of time series data and models. To accommodate for the heterogeneous structure of time series experiments in biological applications, we implemented a data structure of three hierarchy levels. The murefi.Timeseries object represents the time and state vectors , of a single state variable or observation time series. To allow association of state and observed variables via calibr8 calibration models, the Timeseries is initialized with independent_key and dependent_key. Multiple Timeseries are bundled to a murefi.Replicate, which represents either the observations obtained from one reaction vessel, or the predictions made by a process model. Consequently, the murefi.Dataset aggregates replicates of multiple reaction vessels, or the model predictions made for them (Fig 7, center). To allow for a separation of data preprocessing and model fitting in both time and implementation, a murefi.Dataset can be saved as and loaded from a single HDF5 file [47, 48].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 7. Data structures and computation graph of murefi models.

Elements in a comprehensive parameter vector are mapped to replicate-wise model instances. In the depicted example, the model instances for both replicates “B1” and “B2” share θ_1,global as the first element in their parameter vectors. The second model parameter θ₂ is local to the replicates, hence the full parameter vector (left) is comprised of three elements. Model predictions are made such that they resemble the structure of the observed data, having the same number of time points for each predicted time series. An objective function calculating the sum of log-likelihoods is created by associating predicted and observed time series via their respective calibration models. By associating the calibration models based on the dependent variable name, a calibration model may be re-used across multiple replicates, or kept local if, for example, the observations were obtained by different methods.

https://doi.org/10.1371/journal.pcbi.1009223.g007

To describe a reaction system by a system of ODEs, a new class is implemented by subclassing the murefi.BaseODEModel convenience type. In the constructor of the class, the names and order of parameters and state variables are defined, whereas the differential equations are implemented in a dydt instance method. An example is shown in Box 3 Code 3 with the implementation of the Monod kinetics for microbial growth.

Parameter mapping and objective function.

In addition to a murefi model instance, a murefi.Dataset and calibration models, a murefi.ParameterMapping must be defined to facilitate the creation of an objective function. This mapping specifies whether parameters are local or global and the rules with which they are shared between replicates. The ParameterMapping may be represented as a table, assigning each element of replicate-wise parameter vectors to constants or names of parameters in a comprehensive parameter vector. In Fig 7, the parameter mapping is depicted by arrows mapping elements of a 3-element comprehensive parameter vector to 2-element parameter vectors of the replicate-wise models. A table representation of the parameter mapping used to fit the Monod model in Section 4.2 is shown in S4 File.

Model predictions are made such that the time points of the predicted time series match those of the observed data (Fig 7, center). Based on the (in)dependent_key, the predicted and observed Timeseries can be associated with each other and passed to the corresponding CalibrationModel.loglikelihood method to calculate . Note that this procedure conveniently allows for calibration models to be shared by multiple replicates, as well as making observations of one state variable with more than one analytical method.

An objective function performing the computation depicted in Fig 7 can be created with a single call to a convenience function. For compute-efficient optimization and specifically Bayesian parameter estimation, the elements in the parameter vector can also be Aesara tensors, resulting in the creation of a symbolic computation graph. The computation graph cannot only be statically compiled but also auto-differentiated, if all operations in the graph are also auto-differentiable. This is already the case for standard calibr8 calibration models and is also available for murefi -based process models when the sunode [45] package is installed. A corresponding example of obtaining gradients from an ODE model with calibr8, murefi and sunode is included in the murefi documentation [49].

3.2.5 Optimization.

In this work, optimization algorithms are involved at multiple steps of the workflow. Unless otherwise noted we used scipy.optimize.minimize with default settings to obtain the MLEs of calibration and process models. Our current implementation to compute HDIs (Section 3.2.3) uses scipy.optimize.fmin with an underlying Nelder-Mead simplex algorithm.

Initial guesses, as well as parameter bounds for maximum-likelihood optimization, were motivated from prior assumptions or exploratory plots of the data. Based on the intuitive parametrization of the asymmetric logistic (Section 3.2.1) we specified initial guesses for calibration models such that the model prediction from the guessed parameter vector was at least in the same order of magnitude as the data. For MLE of process model parameters, the guessed parameters were motivated from prior assumptions, e.g. the amount of used substrate or values obtained from literature research. Likewise, we specified bounds to be realistic both biologically and based on exploratory scatter plots of the data.

3.2.6 MCMC sampling.

In contrast to optimization, MCMC sampling follows a very different paradigm. Whereas in MLE the likelihood function is iteratively evaluated to find its maximum, Bayesian inference aims to approximate the posterior probability distribution according to Eq 7.

Most sampling algorithms draw the posterior samples in the form of a Markov chain with a equilibrium distribution that matches the posterior probability distribution. While early MCMC algorithms, such as Random-walk Metropolis [50] are conceptually simple and easy to implement, they are computationally ineffective on problems with more than just a handful of dimensions [51, 52]. Instead of implementing inefficient algorithms by hand, practitioners can rely on state of the art libraries for Bayesian inference. These libraries apply automatic transformations, provide diagnostic tools and implement much more efficient sampling algorithms that often use gradients for higher computational efficiency.

Probabilistic programming languages / libraries (PPL), such as PyMC [53], Pyro [54], Stan [55] or Tensorflow Probability [56] use automatic differentiation and typically implement at least one of the gradient-based sampling algorithms Hamiltonian Monte Carlo (HMC) or No-U-Turn Sampling (NUTS) [52]. While PPLs typically require a model to be implemented using their API, other libraries such as emcee [57] provide, for example, gradient-free ensemble algorithms that can be applied to black-box problems. PyMC, the most popular Python-based PPL, implements both gradient-based (HMC, NUTS) as well as gradient-free algorithms, such as Differential Evolution MCMC (DE-MCMC) [58], DE-MCMC-Z [51] or elliptical slice sampling [59] in Python, allowing easy integration with custom data processing and modeling code. In this study, PyMC was used to sample posterior distributions with either DE-MCMC-Z (pymc.DEMetropolisZ) or NUTS.

MCMC sampling of the process model. Whereas in DE-MCMC, proposals are informed from a random pair of other chains in a “population”, the DE-MCMC-Z version selects a pair of states from its own history, the “Z”-dimension. Compared to DE-MCMC, DE-MCMC-Z yields good results with fewer chains that can run independently. The pymc.DEMetropolisZ sampler differs from the original DE-MCMC-Z in a tuning mechanism by which a tune_drop_fraction of by default 90% of the samples are discarded at the end of the tuning phase. This trick reliably cuts away unconverged “burn-in” history, leading to faster convergence.

pymc.DEMetropolisZ was applied to sample the process model in Section 4.2.3. MCMC chains were initialized at the MAP to accelerate convergence of the DE-MCMC-Z sampler in the tuning phase. 50 000 tuning iterations per chain were followed by 500 000 iterations to draw posterior samples for further analysis. The DEMetropolisZ settings remained fixed at ( (default), ϵ = 0.0001) for the entire procedure.

The diagnostic from ArviZ [29] was used to check for convergence (all , S1 Table).

3.2.7 Visualization techniques.

Plots were prepared from Python with a combination of Matplotlib [31], ArviZ and PyMC. We used POV-Ray to produce Figs 2 and 4 and https://diagrams.net for technical drawings. Probability densities were visualized with the pymc.gp.utils.plot_gp_dist helper function that overlays many polygons corresponding to percentiles of a distribution, creating the colorful bands of plots seen in Section 4.2.4 and others. Posterior predictive samples were obtained by randomly drawing observations from the calibration model, based on independent values sampled from the posterior distribution. If not stated otherwise, the densities plotted for MCMC prediction results were obtained from at least 1000 posterior samples. The pair plots of 2-dimensional kernel density estimates of posterior marginals (e.g. S2 Fig) were prepared with ArviZ.

4.1.1 Linear calibration model.

To acquire data for the establishment of a calibration model, 96 glucose standards between 0.001 and were subjected to the glucose assay. A frequent approach to calibration modeling in life sciences is to identify the linear range of an assay and to discard measurements outside this range. From a previous adaptation of the glucose assay for automation with liquid handling robotics, the linear range was expected to be up to (Holger Morschett, personal communication, 2019). Since samples are diluted by a factor of 10 before the assay, 83 glucose standards with concentrations below remain for a linear calibration model.

As described in Section 2.4, calibration models use a probability distribution to describe the relation between independent variable and measurement outcome, both subject to random noise. In this example, we chose a Student-t distribution, thus the change of location parameter μ over the independent variable determines the trend of the calibration model. calibr8 provides a convenience class BasePolynomialModelT that was used to implement a glucose calibration model with linear trend (Box 1 Code 1). For the spread parameter scale, we also chose a linear function dependent on μ to account for increasing random noise in dependency of the absorbance readout of the assay. Both can easily be adapted by changing the respective parameters mu_degree and scale_degree passed to the constructor of the convenience class. The degree of freedom ν in a BasePolynomialModelT is estimated from the data as a constant. The mathematical notation of this model can be found in S2 File.

The calibration model resulting from MLE of location and spread parameters was plotted with another calibr8 convenience function (Fig 8A–8C). The plot shows the calibration model and measurement data (Fig 8A), the same relation with a logarithmic x-axis (Fig 8B) and the relative residuals of data and model predictions (Fig 8C). As it is often recommended for biological assays, the glucose concentrations of the dilution series were evenly spaced on a logarithmic scale [60, 61], thus ensuring a higher accuracy of the model in the low-concentrated area (Fig 8B). To evaluate the quality of the linear calibration model, the residuals of data and model prediction were analyzed (Fig 8C). Overall, the residuals lie between ±5% of the observed values, demonstrating the high precision of the data. For concentrations higher than , an s-shaped trend is observed in the residuals, meaning that data first lies below and then above the linear model prediction. This indicates a lack-of-fit as described in Section 2.4. However, the discrepancy might also be caused by errors in the serial dilution that was pipetted with the robotic liquid handler, resulting in deviations from the expected linear relation. Moreover, it can be seen that the relative spread of residuals is quite constant, meaning that the absolute deviation increases with higher concentrations (Fig 8C). Although the linearly increasing scale parameter accounts for this rise of the underlying random noise, it can be seen that it is slightly overestimated by the model since all data points above lie within a 90% probability interval.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 8. Linear (top) and logistic (bottom) calibration model of glucose assay.

A calibration model comprising linear functions for both the location parameter μ_A365 and the scale parameter of a Student-t distribution was fitted to calibration data of glucose standard concentrations () and absorbance readouts by maximum likelihood estimation (A-C). The calibration data used to fit the linear model is the subset of standards that were spaced evenly on a log-scale up to (B, E). Likewise, a calibration model with a 5-parameter asymmetric logistic function for the μ parameter of the Student-t distribution was fitted to the full calibration dataset (D-E). In both models, the scale parameter was modeled as a 1st-order polynomial function of μ and the degree of freedom ν as a constant. The extended range of calibration standard concentrations up to reveals a saturation kinetic of the glucose assay (A, D) and depending on the glucose concentration, the residuals (C, F) with respect to the modeled location parameter are scattered by approximately 5%. Modeling the scale parameter of the distribution as a 1st-order polynomial function of μ describes the broadening of the distribution at higher concentrations (C).

https://doi.org/10.1371/journal.pcbi.1009223.g008

In comparison to simple linear regression, which is often evaluated by the coefficient of determination R² alone, the demonstrated diagnostics allow to judge whether the choice of model is appropriate. In this case, a more sophisticated model for the spread of the Student-t distribution could be chosen to reduce the lack-of-fit. Moreover, all data points lying above were not considered so far to allow for a linear model. In the following, we will therefore modify the existing calibration model to include a logistic function for the location parameter.

4.1.2 Logistic calibration model.

Although linear calibration models are useful in many cases, some relations in datasets are non-linear in their nature. Moreover, restricting analytical measurements to an approximately linear range instead of calibrating all concentrations of interest can be limiting. If the order of magnitude of sample concentrations is unknown, this leads to laborious dilution series or repetition of measurements to ensure that the linear range is met. In contrast, non-linear calibration models allow to describe complex relationships and, in case of biological assays, to reduce these time- and material-consuming workflows.

Many recommendations for experimental design in calibration can be found in literature (e.g. [60]). Having determined the range of interest for the calibration model, it should be exceeded in both directions if possible, thus ensuring that the relevant concentrations are well-covered. This way, all model parameters, including limits where applicable, can be identified from the observed data. Afterwards, the expected relationship between dependent and independent variable is to be considered. Since the glucose assay readout is based on absorbance in a plate reader (Section 3.1.4), which has a lower and upper detection limit, a saturation effect at high glucose concentrations is expected. In our demonstration example, glucose concentrations of up to were targeted to cover relevant concentration for cultivation (Section 4.2) and at the same time to exceed the linear range towards the upper detection limit.

Sigmoidal shapes in calibration data, e.g. often observed for immunoassays, can be well-described by logistic functions [39]. In the calibr8 package, a generalized logistic function with five parameters is used in an interpretable form (Section 3.2.1). It was used to implement a calibration model where the location parameter μ is described by a logistic function dependent on the glucose concentration. A respective base class BaseAsymmetricLogisticT is provided by calibr8 (S1 File). The mathematical notation of the resulting model is given in S2 File. Using the whole glucose dataset up to , parameters of the new calibration model were estimated (Fig 8D–8F).

Overall, the logistic trend of the location parameter matches the extended calibration data well (Fig 8D and 8E). Since the scale parameter of the Student-t distribution is modeled as a linear function dependent on μ, the width of the likelihood bands approaches a limit at high glucose concentrations (Fig 8F). For concentrations greater than , no residuals lie outside of the 90% probability interval, indicating that the distribution spread is overestimated as it was before. Importantly, a direct comparison between the two calibration models (Fig 8C and 8F) reveals a high similarity in the reduced range (). This demonstrates how a non-linear model extends the range of concentrations available for measurement and modeling while improving the quality of the fit. For the glucose assay, truncating to a linear range thus becomes obsolete.

While non-linear models were so far shown to be useful to extend the usable concentration range of an assay, other applications do not allow to linearly approximate a subrange of measurements. Such an example is the online detection of biomass in growth experiments, where the non-invasive backscatter measurement (Section 3.1.1) does not allow for dilution of the cell suspension during incubation. To model the distribution of backscatter observations as a function of the underlying biomass concentration, a structure similar to the glucose calibration model was chosen (S2 File). In contrast, the location parameter μ was modeled by a polynomial function of the logarithmic cell dry weight (CDW). The final CDW/backscatter calibration model was implemented using the calibr8.BaseLogIndependentAsymmetricLogisticT convenience class (Box 2 Code 2).

Two independent experiments were conducted to obtain calibration data as described in Section 3.1.5. The model was fitted to pooled data using the calibr8.fit_pygmo convenience function. As shown in Fig 9, the model accurately describes the nonlinear correlation between biomass concentration and observed backscatter measurements in the cultivation device (Fig 9A and 9B). Non-linearity is particularly observed for biomass concentrations below (Fig 9A). Moreover, the residual plot (Fig 9C) mainly shows a random distribution; solely residuals between 1 and indicate a lack-of-fit. To assess the potential influence, the resulting uncertainty in estimated biomass concentrations has to be considered, which will be further discussed in Section 4.1.3. Overall, the chosen logistic calibration model describes the calibration data well and is thus useful to transform backscatter measurements from the BioLector device into interpretable quantitative biomass curves.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 9. Calibration model of biomass-dependent backscatter measurement.

Backscatter observations from two independent calibration experiments (1400 rpm, gain = 3) on the same BioLector Pro cultivation device were pooled. A non-linearity of the backscatter/CDW relationship is apparent already from the data itself (A). The evenly spaced calibration data (B) are well-described with little lack-of-fit error (C). At low biomass concentrations the relative spread of the measurement responses starts at ca. 20% and reduces to approximately 2% at concentrations above .

https://doi.org/10.1371/journal.pcbi.1009223.g009

In summary, this section illustrated how calibration models can be built conveniently with calibr8 and showed that the asymmetric logistic function is suitable to describe many relationships in biotechnological measurement systems. Note that the only requirement to estimate CalibrationModel parameters with calibr8 is to provide a numpy.ndarray with data for the independent and dependent variable. The code is agnostic of bioprocess knowledge and can thus be transferred to other research fields without further adaptations.

Having demonstrated how concentration/readout relations can be described by different calibration models, a remaining question is how to apply those calibration models. An important use-case is to obtain an estimate of concentrations in unknown samples, where uncertainty quantification is a crucial step.

4.1.3 Uncertainty quantification on independent variables.

After establishing a calibration model, the practitioner can in most cases consider the parameters of the model as fixed. Characterization of measurement reproducibility is thereby externalized into the calibration procedure, where random noise is inherently described by the spread of a probability distribution. The calibration model can then be put into application for the quantification of the independent variable from new observations. As introduced before, not only a single value of the independent variable is desired, but also a measure of uncertainty about it.

Quantifying the uncertainty in the independent variable as a continuous probability density is not only intuitive to visually interpret (Section 2.4), but also flexible with respect to the question of interest. To quantify the uncertainty numerically, various kinds of credible intervals (Section 2.3) can be obtained. For example, one might estimate the equal-tailed interval in which the independent variable lies with 90% probability, or alternatively the probability that it lies above a certain threshold.

In calibr8, the CalibrationModel.infer_independent method is used to perform the uncertainty quantification from one or more observations (Section 3.2.3). Internally, it uses the loglikelihood method of the calibration model and numerically integrates the sum of log-likelihoods over a user-specified range of plausible independent variables (Section 3.2.3). The resulting calibr8.ContinuousUnivariateInference is equivalent to Eq 13, where the prior p(x) is specifying the plausible range. (13)

For convenience, the CalibrationModel.infer_independent method automatically determines median and credible interval (ETI and HDI) bounds. It determines vectors for the independent variable and the conditional probability density that can be plotted without further processing.

In Fig 10, various inferences obtained with infer_independent are illustrated with a biomass calibration model. For illustration purposes, the calibration model from Section 4.1.2 was slightly modified and the degree of freedom ν was modified to an extreme value of 1. Different sets of simulated observations were subjected to the infer_independent method and the resulting credible intervals are shown in subplots A and B.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 10. Independent variable PDFs in various observation scenarios.

Posterior densities inferred from various numbers of observations corresponding to different biomass concentrations are shown (A). The ends of the drawn lines in A indicate the 95% equal-tailed interval. Near biomass concentrations of 0, the posterior density is asymmetric (A, blue), indicating that very low concentrations cannot be distinguished. As the number of observations grows, the probability mass is concentrated and the ETIs shrink (A, oranges). The choice of a Student-t distribution model can lead to a multi-modality of the inferred posterior density when observations lie far apart (B). For asymmetric distributions, the median (dashed line) does not necessarily coincide with a mode and equal-tailed and highest-density intervals (ETI, HDI) can be different. Maximum likelihood estimates from individual observations, as obtained via predict_independent are shown as arrows. Note: and the model’s ν parameter were chosen at extreme values for illustrative purposes.

https://doi.org/10.1371/journal.pcbi.1009223.g010

In Fig 10A, the orange curves show the 95% ETIs for one, two or three observations of backscatter signals corresponding to . As expected, the probability mass is concentrated around with higher number of observations and the ETIs get narrower. Furthermore, observations in the lower or upper saturation of the measurement system typically result in one-sided probability densities (blue).

When the calibration model assumes the possibility of outliers (Student-t distributed measurement responses), the observation of drastically different measurement responses can translate into a multi-modal posterior belief in the independent variable. The intuition behind this multi-modality is that a subset of observations are “outliers” from the perspective of the remaining observations and vice versa. In the example shown in Fig 10B, the three observations around 0.5 could be “outliers”, or the ones around 1.3, but from the data alone both are equally likely. Hence the posterior belief in the biomass concentration is bimodal. The green arrows indicate the MLE estimates obtained via predict_independent with individual observations. The plot also reveals slight differences between ETI and HDI in this scenario. Note that both the degree of freedom and the simulated observations have been set to extreme values to illustrate the properties of the Student-t distribution.

The Bayesian or likelihood-based perspective on uncertainty in the independent variable (Eq 13) allows for quantification of uncertainty even with single observations, close to zero, or close to saturation limits of the measurement system. Calibration models built with calibr8 are easy to set up, visualize and diagnose and can thus be flexibly integrated into existing data analysis workflows of various domains. Moreover, the setup in a versatile, object-oriented programming language such as Python allows to use calibr8 in high-throughput, automated experiments where hundreds of calibration models must be fitted. Next, we will build upon the presented biomass and glucose calibration models and demonstrate how combining them with a bioprocess model enables to gain insight into the growth phenotype of a biotechnological model organism, Corynebacterium glutamicum.

4.2.1 Building an ODE process model for bacterial growth experiments.

The simplest model for microbial growth is the Monod kinetics differential equation model of substrate-limited exponential growth [64]. Similar to how the famous Michaelis-Menten kinetics describe enzymatic reaction rates, the Monod kinetics model the specific growth rate as a function of substrate concentration. Under the assumptions of homogeneous mixing, unlimited nutrient supply and constant ambient conditions, the Monod model can be applied to batch cultivations of bacterial, fungal, plant or cell cultures that grow with a maximum growth rate μ_max until a substrate, typically a carbon source, is depleted.

The Monod model (Eq 14) has five parameters including the initial conditions for substrate concentration S₀ and biomass concentration X₀. The maximum growth rate μ_max specifies the specific exponential growth rate that the organism can achieve under the modeled conditions. The actual specific growth rate μ(t) is modeled as a function of μ_max, the current substrate concentration S and a parameter K_S that corresponds to the substrate concentration at which . The last parameter Y_XS, called biomass yield, describes the amount of substrate consumed per unit of formed biomass. (14)

The experiment to be modeled in this application example was devised such that Monod-like growth behavior of C. glutamicum wild-type could be expected (Section 3.1.2). We grew 28 parallel batch cultures that were sampled to measure glucose concentrations in addition to the high-resolution backscatter time series. The resulting dataset comprises 28 replicates, each with backscatter time series of varying length and a time series of length 1 for the glucose absorbance readout. Building upon our Python package murefi for flexible multi-replicate fitting, we loaded the raw observations into a murefi.Dataset object (Section 3.2.4). The package was designed to simplify the definition and parameter estimation of process models that describe all replicates in a dataset simultaneously.

To build such elaborate process models with murefi, the user must specify the process model corresponding to a single replicate, as well as a set of rules that describe how parameters of this model are shared across replicates. The Monod kinetics in this application example were implemented in just a few lines of code by subclassing from murefi.BaseODEModel (Box 3 Code 3).

For heterogeneous datasets, the rules for sharing process model parameters across replicates can be complex and hard to implement and most modeling workflows require the practitioner to often change the parametrization. In murefi, the ParameterMapping class supports the modeler by specializing in the tedious translation of parameter sharing rules into a function (.repmap(…)) that takes a single parameter vector and transforms it into replicate-specific parameter vectors. At the same time, it provides mechanisms for specifying fixed parameters, initial guesses and bounds on the parameters. Reading a spread sheet with parameters into Python is an easy way of initializing the ParameterMapping (Table 1).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Tabular representation of a parameter mapping.

With columns corresponding to the parameter names of a naive Monod process model, the parametrization of each replicate, identified by a replicate ID (rid) is specified in a tabular format. Parameter identifiers that appear multiple times (e.g. S0) correspond to a parameter shared across replicates. Accordingly, replicate-local parameters names simply do not appear multiple times (e.g. X0_A06). Numeric entries are interpreted as fixed values and will be left out of parameter estimation. Columns do not need to be homogeneously fixed/shared/local, but parameters can only be shared within the same column. The parameter mapping can be provided as a DataFrame object.

https://doi.org/10.1371/journal.pcbi.1009223.t001

Unique names specify that a parameter is only estimated from the indicated replicate (e.g. X0_A02) while shared names correspond to global parameters (e.g. S0). For the application example at hand, a parameter mapping was defined such that the parameter X₀ is local to each replicate while S₀, μ_max and Y_XS are shared across all replicates. For the Monod substrate affinity constant K_S, literature reports values of approximately for Escherichia coli [65]), while no data is available for C. glutamicum. Because it is practically non-identifiable at the resolution of our dataset, K_S was fixed to an arbitrary, but numerically harmless value of . In Table 1, this is expressed by the numerical column entries.

A likelihood function for parameter estimation was created using the murefi.objectives.for_dataset convenience function (Box 4 Code 4). The objective is independent of the parameter estimation paradigm and was applied for optimization via MLE (Section 4.2.2) and sampling by MCMC (Section 4.2.3) in the scope of this work.

4.2.2 Estimating ODE process model parameters by maximum likelihood.

First, we determined maximum likelihood estimates of the process model parameters through optimization. In few lines of code, the calibration models from Section 4.1 and dataset are loaded (Box 4 Code 4, ll. 2–4), the process model is instantiated (Box 4 Code 4, l. 1) and the ParameterMapping is specified with bounds and guesses (Box 4 Code 4, ll. 7–21). The objective (Box 4 Code 4, ll. 22–27) can directly be used for an optimization algorithm (Box 4 Code 4, ll. 28–32), in this case one from the popular Python library scipy. A table with MLE parameters can be found in S1 Table.

Fig 11 shows the observations alone (A) and combined with MLE results (B) for glucose (absorbance) and biomass (backscatter). The replicates were sampled at different times to measure glucose concentrations; the end of a time series is indicated by an arrow and the replicate name (Fig 11A). Overall, the backscatter time series show a very high reproducibility, which demonstrates the effect of pooling precultures before inoculation (Section 3.1.2). The model describes the observations so accurately that they can only be distinguished in the inset plot (Fig 11B). Here, a small difference between replicates can be observed, which is caused by varying initial biomass concentrations due to inevitable pipetting errors in the automated inoculation of the main cultures. It becomes evident that replicate-wise X₀ parameters were necessary to account for this effect. The different initial biomasses are also visible from the spread of data points at the beginning of the growth curve (Fig 11B). For the biomass, the only period of systematic deviation between model prediction and observations is at the time of entry into the stationary phase, the phase where substrate is depleted and growth stops. Here, the biomass signal overshoots while the Monod kinetics predict a rapid change to a constant signal. This effect in the growth curve of C. glutamicum is also known from other experiments with the BioLector [66] and cannot be accounted for by the otherwise useful textbook process model.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 11. Measurements and maximum likelihood estimate of C. glutamicum growth Monod model.

Original measurement responses of online biomass (backscatter) and at-line endpoint glucose assay measurements (absorbance) are shown in (A). Glucose measurements were obtained by sacrificing culture wells, hence each backscatter time series terminates at the time of glucose assay observations. The time and well ID of sacrifices are marked by arrows, colored by row in the cultivation MTP. The inset plot shows a typical layout of the cultivation plate (FlowerPlate). The preculture wells are highlighted in green, main cultures in black. In B, the observations and MLE predictions of the ODE process model are shown in SI units. Observations were transformed from original units using the predict_independent method of the respective calibration model. Whereas all curves start at the same global initial substrate concentration S₀, each well has individual initial biomass concentrations, resulting in the time shifts visible in the zoomed-in inset plot. Biomass observations in the inset plot (●) correspond to the median posterior inferred from each backscatter observation individually.

https://doi.org/10.1371/journal.pcbi.1009223.g011

The glucose data shows more deviation, but follows the expected textbook behaviour of exponential decay (Fig 11B). Interestingly, the predictions for glucose concentrations at the end of cultivation lie slightly above , showing that the corresponding calibration model is not describing this range of concentrations well. The deviation could be caused by other components in the used cultivation medium that distort the measurement compared to calibration with fresh medium as diluent. However, this was not further investigated since the substrate data has little influence on the parameter estimation compared to the high-resolution backscatter measurements.

From a first inspection of MLE results, we can see that the simple Monod process model describes the high-resolution data very well. For more insight, we will take a look at the parameter estimation, correlations and systematic deviations using a Bayesian approach.

4.2.3 Hierarchical Bayesian ODE models with calibr8 and murefi.

The results presented in the previous chapter show that the Monod model, when combined with non-linear calibration models for the observations, can describe the observed biological process with high accuracy. However, the precision of the parameter set obtained by the maximum likelihood method is still unknown. Particularly, when decisions are made from model-based inferences and predictions, the uncertainty about these variables is a key factor.

The combination of (forward) sensitivity analysis with Gaussian error propagation could be applied to learn about the precision of the maximum likelihood estimate. Instead of maximum likelihood optimization of a parameter set, Bayes’ rule can be used to infer a posterior probability distribution of parameters. In comparison to the maximum likelihood method, the Bayesian approach allows to incorporate prior knowledge and inherently quantifies uncertainty and parameter correlations. Bayesian posteriors can in some (rare) cases be obtained analytically, or numerically as shown in Section 4.1.3. However, in most practical applications Markov chain Monte Carlo (MCMC) algorithms are applied. MCMC offers convergence guarantees as the number of iterations approaches infinity and can give satisfactory results with competitive computational performance when modern algorithms are used.

To build a Bayesian process model, one must explicitly state prior beliefs in the model parameters in the form of probability distributions. For our hierarchical Monod model application example, we must specify prior beliefs in the ODE parameters μ_max, Y_XS and initial conditions S₀ and X_0,well. Prior distributions for these parameters were specified to reflect biologically reasonable, but uninformative assumptions about the experiment (Eq 15). The initial substrate concentration S₀ was expected at approximately with a conservative 10% relative error. For C. glutamicum wild-type, our priors for biomass yields with and for maximum growth rates with are uninformative and based on literature [67]. Our process model describes initial biomass concentrations on a per-well basis (Section 4.2.1), but can still infer the mean initial biomass concentration X_0,μ as a hyperprior by modeling well-specific offsets w.r.t. the group mean as X_0,well = X_0,μ ⋅ F_offset,well. Through X_0,μ the priors for all initial biomass concentrations are parametrized by a common parameter, allowing each individual X_0,well to vary while concentrating around their common group mean. For more intuition and details about Bayesian hierarchical modeling in particular, we refer to [68]. While the experiment performed here is modeled with a hierarchical prior on a process model parameter, one may also use hierarchical priors for calibration model parameters. A corresponding example is shown in the calibr8 documentation [22].

The experiment was programmed to inoculate main cultures to approximately (Section 3.1.2), therefore the prior for X_0,μ was centered at with a 10% relative error. Our prior belief in the well-specific relative offset F_offset,well was also modeled by a Lognormal distribution with mean 0, corresponding to the expectation that the offset is centered around 1 in the transformed space. A standard deviation of 20% was chosen to account for random and systematic inaccuracy of the automated liquid handler at the low pipetting volume of 20 μL [66]. (15)

When modeling with calibr8 and murefi, this specification of prior beliefs is the only overhead compared to the MLE method. The API of both packages was designed to be fully compatible with the probabilistic programming library PyMC, such that calibr8 and murefi models can become fully Bayesian with little programming effort.

Concretely, the objective function created by murefi accepts Aesara tensors (e.g. PyMC random variables) as inputs, resulting in a symbolic TensorVariable likelihood instead of a numeric one. The PyMC model for the hierarchical ODE process model in our application example builds upon the previously established objective function (Box 4 Code 4, l. 22). The model code (Box 5 Code 5) resembles the mathematical notation of the same model shown in Eq 15.

After the PyMC process model was defined, its parameters were estimated by MCMC as described in Section 3.2.6. Two-dimensional marginals of the posterior samples obtained from MCMC sampling are shown in Fig 12 for two replicates and in S2 Fig for the whole dataset.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 12. Parameter correlations, data and posterior distributions of growth curves.

Each kernel density estimate (KDE) in the top half shows a 2-dimensional cross section of the full posterior, visualizing correlations between some of the model parameters. For example, the topmost KDE shows that the posterior samples of F_offset,D04 are correlated with X_0,μ. Axis labels correspond to the lower and upper bound of 90% HDIs. The large pair plot shows just the marginals that are relevant for the replicates D04 and D06, whereas the small pair plot shows the dimensions for all parameters (high resolution in S2 Fig). In the bottom half of the figure, the kinetics of replicates D04 and D06 are drawn. The red (substrate) and green (biomass) densities correspond to the distribution of predictions obtained from posterior samples, as described in Section 3.2.7. The red violins visualize the posterior inferred from single glucose measurement responses without the use of the process model. Likewise, the green vertical bars on the biomass concentrations show the 90% HDI.

https://doi.org/10.1371/journal.pcbi.1009223.g012

The pair plot visualization of the posterior reveals that some model parameters are strongly correlated with each other. Among those strong correlations are the pair of initial substrate concentration S₀ and biomass yield Y_XS. Interestingly, correlations were found even in the very narrow HDIs of X_0,well and μ_max, which is particularly clear for replicate D06. An interpretation is that when the initial biomass concentration is estimated at a smaller value, the maximum growth rate of cells must be higher to reach the same biomass level. The correlation is thus a natural consequence of the underlying process. Similarly, a lower initial substrate concentration results in a higher yield.

From a modeling point of view, the plot reveals how identifiable the model parameters are from the data. Furthermore, strong correlations, as observed for Y_XS and S₀, can be problematic for some optimization or MCMC sampling algorithms. In this case, the applied algorithm DE-Metropolis-Z [51] proved beneficial to sample the 32-dimensional parameter space with highly correlated parameters (Fig 12, top right). Interestingly, the strength of the correlation depends on the amount of data that was available for a particular replicate (S1 File). The more data available, the stronger the correlation between X₀ and μ_max; this can also be observed for wells D04 and D06. The parameter estimates by MCMC are also tabulated in S1 Table.

In the lower part of Fig 12, the observations as well as model predictions with parameters sampled from the posterior are shown. Each line in the density plot corresponds to one set of parameters sampled with the MCMC algorithm. The small width of the density bands express how well the parameters could be estimated from the data, which is in accordance to the pair plot above. The violins around the substrate data visualize the uncertainty of glucose concentration inferred with the calibration model alone, instead of using the process model with all evidence. The violin is wider than the posterior band from the process model accordingly. Similar to the the MLE results, it becomes obvious that the Monod model estimate is well-suited to describe the biological dataset. With calibr8 and murefi, building and sampling the Bayesian model needs a similar effort as MLE and the user can focus on structural requirements rather than cumbersome implementation. To assess the benefits of the Bayesian model in more detail, the role of different calibration models, the residuals and the hierarchical parameter X₀ are investigated in more detail in the next section.

4.2.4 Process and model insight through Bayesian uncertainty quantification.

From the process model fit and the uncertainty estimates in particular, one can draw conclusions about the choice of model and the underlying biological process. First, to emphasize that the elaborate non-linear calibration model was required, we compare the process model fits obtained with a non-linear versus a linear calibration model. The more traditional linear biomass/backscatter correlation was fitted to calibration data as described in Section 3.1.5 and used to fit the D06 replicate from our dataset. For comparison, the asymmetric logistic calibration model from Section 4.1.1 was used to estimate parameters of the same process model and data.

On a first glance, the fit of the Monod process model using the linear biomass calibration model looks like a good description of the data (Fig 13A), but does not hold up to closer inspection. The residual plots (Fig 13B and 13C) reveal that using the linear calibration model results in systematically larger residuals of the process model compared to using the logistic calibration model. A thorough inspection of the linear calibration model itself (Fig 13D) also reveals that it already has a lack-of-fit of the location parameter (green line), similar to the depiction in Fig 3. We would like to point out that also the maximum growth rate estimated from a process model with linear biomass/backscatter calibration () is systematically overestimated compared to the one with the logistic model (). Regarding the choice of calibration model for the biomass/backscatter relationship, we conclude that the linear model should no longer be used, as it results in biased parameter estimates.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 13. Comparison of Monod model fit with linear error model.

Two Monod kinetic process models were fitted to the same observations from culture well D06 utilizing either a linear calibration model for the biomass/backscatter relationship (orange in A, calibration in D) or the previously established logistic model (blue in A). In A the posterior distribution of backscatter observations (density bands) is overlaid with actual backscatter observations. A linear calibration (D) model with fixed intercept (Section 3.1.5) was fitted to the subset of calibration data points up to such that it covers the range of biomass concentrations expected in the experiment. Residual plots of the observations compared to the posterior predictive distribution of backscatter observations (B, C) show that the fit obtained with the logistic calibration model (blue) has much less lack-of-fit compared to the one with the linear model (orange). Note that the backscatter residuals of ±1% are small compared to the amplitude of the absolute values going from close to 0 to approximately 20. The discrepancy between the two models is also evident from the 90% HDI of the maximum growth rate μ_max of [0.414, 0.423] h⁻¹ in the logistic and [0.480, 0.530] h⁻¹ in the linear case.

https://doi.org/10.1371/journal.pcbi.1009223.g013

Having chosen a suitable calibration model for the variables, the choice of the Monod model itself can be investigated. Fig 14 shows the high-resolution biomass data and predictions from MCMC on a logarithmic y-scale (Fig 14A) as well as the residuals in backscatter units (Fig 14B). In the left subplot, the data was transformed to biomass concentrations with the logistic biomass calibration model. The orange intervals represent the inferred from a single observation using only the calibration model. In contrast, the blue density represents the posterior of the process model, which contains all observations. Naturally, the posterior from all evidence, combined through the process model, is much narrower than the posterior from any single observation. The plot reveals that the exponential growth assumed by the Monod model is generally suitable for the growth on glucose, since the blue density is describing the trend of observations well.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 14. Predictions, observations and residuals of Monod model fitted to backscatter data.

A: Through a logarithmic y-axis, the plot A shows that both process model (blue density) and the obtained from the biomass calibation model with individual observations (orange) describe an exponentially increasing biomass concentration up to approximately 9 hours. B: The residuals between prediction and observed backscatter (black) and the posterior predictive backscatter distribution (green density) show that the lack-of-fit is consistently less than ±0.25 backscatter units with the exception of a fluctuation at the time of substrate depletion.

https://doi.org/10.1371/journal.pcbi.1009223.g014

To evaluate a lack-of-fit, the residual plot (Fig 14B) should be considered. Here, the residuals between the process model posterior and the observed backscatter are shown in black, the respective posterior predictive distribution of measurement responses (Section 3.2.3) is shown in green. The posterior predictive is the distribution of measurement responses that the model predicts. First, biomass concentrations are drawn from the posterior distribution. At each biomass concentration, another sample is taken from the Student-t distribution predicted by the biomass calibration model.

A large deviation that cannot be explained with the uncertainty of the estimate can be observed after 8 hours. Looking at the data, e.g. in Section 4.2.4, it can be seen that it accounts for the previously described overshoot of the backscatter signal at the beginning of the stationary phase (Section 4.2.2). This phenomenon cannot be explained by the Monod model, which assumes a constant biomass concentration after substrate depletion. Further investigations are needed to identify whether the change is morphological, e.g. a shrinking of cells due to carbon source depletion, or a decrease of biomass, e.g. by cell lysis.

Before 8 hours, an s-shaped systematic deviation can be observed, meaning that the observations first lie above and then below the prediction. Apart from the influence of the overshoot, which distorts the fit, this might be explained by a different growth rate. It was previously shown that C. glutamicum exhibits a higher specific growth rate on protocatechuic acid (PCA), which is a component of the cultivation medium CGXII [67]. Upon depletion of PCA after the first hours of cultivation, the growth rate decreases accordingly. This is not accounted for in the Monod kinetics, which describe an effectively constant growth rate at substrate concentrations much higher than the K_S value. To cover this effect, PCA must be measured, e.g. by sampling and liquid chromatography, and a more elaborate process models with several substrates must be utilized. Nevertheless, the very simple Monod kinetics describe the overall growth behaviour well and residuals are low.

In Fig 11, we have seen that the time differences in the exponential phases between replicates are well explained by the well-wise initial biomass concentrations . The choice of a hierarchical process models is further evaluated in Fig 15, which shows the estimated with uncertainties for all replicates. For replicates with more evidence (longer time series), the posterior probability for their initial biomass concentration is concentrated in a narrow interval, whereas X₀ in wells with little evidence was estimated with considerably more uncertainty. The posterior for the group mean X_0,μ is concentrated at , close to the theoretical concentration () expected from the experimental design.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 15. Posterior group mean and well-specific initial biomass concentrations X₀.

Variability between the growth curves in separate wells is described by well-specific initial biomass concentrations X_0,well. Their posterior probability distribution is wide if the well was sacrificed early (left) and narrows down with the number of observed time points (right). Their common hyperprior (a.k.a group mean prior) X_0,μ for the mean of each X_0,well was updated to a posterior with .

https://doi.org/10.1371/journal.pcbi.1009223.g015

Overall, the well-wise modeling of initial biomass concentrations as well as the separate modeling of replicates allowed us to account for inevitable differences between wells, while inferring the key process model parameters from all data. The combination of calibr8 and murefi made it possible to construct a process models of our application example with little code and apply both optimization (MLE) and Bayesian inference (MCMC) without needing to change any implementation details (Box 4 Code 4, Box 5 Code 5). Our application example showed that Bayesian inference with ODE-based process models to 28 parallel cultures with hundreds of observations is not only technically feasible, but also accessible without deep understanding of probabilistic programming frameworks.

As implied in the famous quote by George E.P. Box—“All models are wrong, but some are useful.”—also our Monod kinetics process model does not describe every aspect of the data, but is a powerful tool to quantify key process parameters under uncertainty. From its (in)accuracies, we can gain insight into the bioprocess and generate new hypotheses about the biological process or measurement system that are yet to be understood. In our case, the uncertainty quantification of process model parameters can become the cornerstone of bioprocess development by facilitating robust and intuitive statistical analysis or Bayesian decision-making.