Efficient characterization of high-dimensional parameter spaces for systems biology

DyMMM-LEAPS: An ML-based framework for modulating evenness and stability in synthetic microbial communities

There have been a growing number of computational strategies to aid in the design of synthetic microbial consortia. A framework to identify regions in parametric space to maximize two essential properties, evenness and stability, is critical. In this study, we introduce DyMMM-LEAPS (dynamic multispecies metabolic modeling-locating evenness and stability in large parametric space), an extension of the DyMMM framework. Our method explores the large parametric space of genetic circuits in synthetic microbial communities to identify regions of evenness and stability. Due to the high computational costs of exhaustive sampling, we utilize adaptive sampling and surrogate modeling to reduce the number of simulations required to map the vast space. Our framework predicts engineering targets and computes their operating ranges to maximize the probability of the engineered community to have high evenness and stability. We demonstrate our approach by simulating five cocultures and one three-strain culture with different social interactions (cooperation, competition, and predation) employing quorum-sensing-based genetic circuits. In addition to guiding circuit tuning, our pipeline gives an opportunity for a detailed analysis of pockets of evenness and stability for the circuit under investigation, which can further help dissect the relationship between the two properties. DyMMM-LEAPS is easily customizable and can be expanded to a larger community with more complex interactions.

Deciphering Metabolic Pathways in High-Seeding-Density Fed-Batch Processes for Monoclonal Antibody Production: A Computational Modeling Perspective

Due to their high specificity, monoclonal antibodies (mAbs) have garnered significant attention in recent decades, with advancements in production processes, such as high-seeding-density (HSD) strategies, contributing to improved titers. This study provides a thorough investigation of high seeding processes for mAb production in Chinese hamster ovary (CHO) cells, focused on identifying significant metabolites and their interactions. We observed high glycolytic fluxes, the depletion of asparagine, and a shift from lactate production to consumption. Using a metabolic network and flux analysis, we compared the standard fed-batch (STD FB) with HSD cultivations, exploring supplementary lactate and cysteine, and a bolus medium enriched with amino acids. We reconstructed a metabolic network and kinetic models based on the observations and explored the effects of different feeding strategies on CHO cell metabolism. Our findings revealed that the addition of a bolus medium (BM) containing asparagine improved final titers. However, increasing the asparagine concentration in the feed further prevented the lactate shift, indicating a need to find a balance between increased asparagine to counteract limitations and lower asparagine to preserve the shift in lactate metabolism.

Deduction of signaling mechanisms from cellular responses to multiple cues

Cell signaling networks are complex and often incompletely characterized, making it difficult to obtain a comprehensive picture of the mechanisms they encode. Mathematical modeling of these networks provides important clues, but the models themselves are often complex, and it is not always clear how to extract falsifiable predictions. Here we take an inverse approach, using experimental data at the cell level to deduce the minimal signaling network. We focus on cells' response to multiple cues, specifically on the surprising case in which the response is antagonistic: the response to multiple cues is weaker than the response to the individual cues. We systematically build candidate signaling networks one node at a time, using the ubiquitous ingredients of (i) up- or down-regulation, (ii) molecular conversion, or (iii) reversible binding. In each case, our method reveals a minimal, interpretable signaling mechanism that explains the antagonistic response. Our work provides a systematic way to deduce molecular mechanisms from cell-level data.

Design centering enables robustness screening of pattern formation models

MOTIVATION

Access to unprecedented amounts of quantitative biological data allows us to build and test biochemically accurate reaction-diffusion models of intracellular processes. However, any increase in model complexity increases the number of unknown parameters and, thus, the computational cost of model analysis. To efficiently characterize the behavior and robustness of models with many unknown parameters remains, therefore, a key challenge in systems biology.

RESULTS

We propose a novel computational framework for efficient high-dimensional parameter space characterization of reaction-diffusion models in systems biology. The method leverages the Lp-Adaptation algorithm, an adaptive-proposal statistical method for approximate design centering and robustness estimation. Our approach is based on an oracle function, which predicts for any given point in parameter space whether the model fulfills given specifications. We propose specific oracles to efficiently predict four characteristics of Turing-type reaction-diffusion models: bistability, instability, capability of spontaneous pattern formation and capability of pattern maintenance. We benchmark the method and demonstrate that it enables global exploration of a model's ability to undergo pattern-forming instabilities and to quantify robustness for model selection in polynomial time with dimensionality. We present an application of the framework to pattern formation on the endosomal membrane by the small GTPase Rab5 and its effectors, and we propose molecular mechanisms underlying this system.

AVAILABILITY AND IMPLEMENTATION

Our code is implemented in MATLAB and is available as open source under https://git.mpi-cbg.de/mosaic/software/black-box-optimization/rd-parameter-space-screening.

SUPPLEMENTARY INFORMATION

Supplementary data are available at Bioinformatics online.

Analysis of cellular kinetic models suggest that physiologically based model parameters may be inherently, practically unidentifiable

Physiologically-based pharmacokinetic and cellular kinetic models are used extensively to predict concentration profiles of drugs or adoptively transferred cells in patients and laboratory animals. Models are fit to data by the numerical optimisation of appropriate parameter values. When quantities such as the area under the curve are all that is desired, only a close qualitative fit to data is required. When the biological interpretation of the model that produced the fit is important, an assessment of uncertainties is often also warranted. Often, a goal of fitting PBPK models to data is to estimate parameter values, which can then be used to assess characteristics of the fit system or applied to inform new modelling efforts and extrapolation, to inform a prediction under new conditions. However, the parameters that yield a particular model output may not necessarily be unique, in which case the parameters are said to be unidentifiable. We show that the parameters in three published physiologically-based pharmacokinetic models are practically (deterministically) unidentifiable and that it is challenging to assess the associated parameter uncertainty with simple curve fitting techniques. This result could affect many physiologically-based pharmacokinetic models, and we advocate more widespread use of thorough techniques and analyses to address these issues, such as established Markov Chain Monte Carlo and Bayesian methodologies. Greater handling and reporting of uncertainty and identifiability of fit parameters would directly and positively impact interpretation and translation for physiologically-based model applications, enhancing their capacity to inform new model development efforts and extrapolation in support of future clinical decision-making.

Review of applications and challenges of quantitative systems pharmacology modeling and machine learning for heart failure

Quantitative systems pharmacology (QSP) is an important approach in pharmaceutical research and development that facilitates in silico generation of quantitative mechanistic hypotheses and enables in silico trials. As demonstrated by applications from numerous industry groups and interest from regulatory authorities, QSP is becoming an increasingly critical component in clinical drug development. With rapidly evolving computational tools and methods, QSP modeling has achieved important progress in pharmaceutical research and development, including for heart failure (HF). However, various challenges exist in the QSP modeling and clinical characterization of HF. Machine/deep learning (ML/DL) methods have had success in a wide variety of fields and disciplines. They provide data-driven approaches in HF diagnosis and modeling, and offer a novel strategy to inform QSP model development and calibration. The combination of ML/DL and QSP modeling becomes an emergent direction in the understanding of HF and clinical development new therapies. In this work, we review the current status and achievement in QSP and ML/DL for HF, and discuss remaining challenges and future perspectives in the field.

A rationally engineered decoder of transient intracellular signals

Cells can encode information about their environment by modulating signaling dynamics and responding accordingly. Yet, the mechanisms cells use to decode these dynamics remain unknown when cells respond exclusively to transient signals. Here, we approach design principles underlying such decoding by rationally engineering a synthetic short-pulse decoder in budding yeast. A computational method for rapid prototyping, TopoDesign, allowed us to explore 4122 possible circuit architectures, design targeted experiments, and then rationally select a single circuit for implementation. This circuit demonstrates short-pulse decoding through incoherent feedforward and positive feedback. We predict incoherent feedforward to be essential for decoding transient signals, thereby complementing proposed design principles of temporal filtering, the ability to respond to sustained signals, but not to transient signals. More generally, we anticipate TopoDesign to help designing other synthetic circuits with non-intuitive dynamics, simply by assembling available biological components.

Smart computational exploration of stochastic gene regulatory network models using human-in-the-loop semi-supervised learning

MOTIVATION

Discrete stochastic models of gene regulatory network models are indispensable tools for biological inquiry since they allow the modeler to predict how molecular interactions give rise to nonlinear system output. Model exploration with the objective of generating qualitative hypotheses about the workings of a pathway is usually the first step in the modeling process. It involves simulating the gene network model under a very large range of conditions, due to the large uncertainty in interactions and kinetic parameters. This makes model exploration highly computational demanding. Furthermore, with no prior information about the model behavior, labor-intensive manual inspection of very large amounts of simulation results becomes necessary. This limits systematic computational exploration to simplistic models.

RESULTS

We have developed an interactive, smart workflow for model exploration based on semi-supervised learning and human-in-the-loop labeling of data. The workflow lets a modeler rapidly discover ranges of interesting behaviors predicted by the model. Utilizing that similar simulation output is in proximity of each other in a feature space, the modeler can focus on informing the system about what behaviors are more interesting than others by labeling, rather than analyzing simulation results with custom scripts and workflows. This results in a large reduction in time-consuming manual work by the modeler early in a modeling project, which can substantially reduce the time needed to go from an initial model to testable predictions and downstream analysis.

AVAILABILITY AND IMPLEMENTATION

A python-package is available at https://github.com/Wrede/mio.git.

SUPPLEMENTARY INFORMATION

Supplementary data are available at Bioinformatics online.

Robustness and parameter geography in post-translational modification systems

Biological systems are acknowledged to be robust to perturbations but a rigorous understanding of this has been elusive. In a mathematical model, perturbations often exert their effect through parameters, so sizes and shapes of parametric regions offer an integrated global estimate of robustness. Here, we explore this “parameter geography” for bistability in post-translational modification (PTM) systems. We use the previously developed “linear framework” for timescale separation to describe the steady-states of a two-site PTM system as the solutions of two polynomial equations in two variables, with eight non-dimensional parameters. Importantly, this approach allows us to accommodate enzyme mechanisms of arbitrary complexity beyond the conventional Michaelis-Menten scheme, which unrealistically forbids product rebinding. We further use the numerical algebraic geometry tools Bertini, Paramotopy, and alphaCertified to statistically assess the solutions to these equations at ∼10⁹ parameter points in total. Subject to sampling limitations, we find no bistability when substrate amount is below a threshold relative to enzyme amounts. As substrate increases, the bistable region acquires 8-dimensional volume which increases in an apparently monotonic and sigmoidal manner towards saturation. The region remains connected but not convex, albeit with a high visibility ratio. Surprisingly, the saturating bistable region occupies a much smaller proportion of the sampling domain under mechanistic assumptions more realistic than the Michaelis-Menten scheme. We find that bistability is compromised by product rebinding and that unrealistic assumptions on enzyme mechanisms have obscured its parametric rarity. The apparent monotonic increase in volume of the bistable region remains perplexing because the region itself does not grow monotonically: parameter points can move back and forth between monostability and bistability. We suggest mathematical conjectures and questions arising from these findings. Advances in theory and software now permit insights into parameter geography to be uncovered by high-dimensional, data-centric analysis.

Biological organisms are often said to have robust properties but it is difficult to understand how such robustness arises from molecular interactions. Here, we use a mathematical model to study how the molecular mechanism of protein modification exhibits the property of multiple internal states, which has been suggested to underlie memory and decision making. The robustness of this property is revealed by the size and shape, or “geography,” of the parametric region in which the property holds. We use advances in reducing model complexity and in rapidly solving the underlying equations, to extensively sample parameter points in an 8-dimensional space. We find that under realistic molecular assumptions the size of the region is surprisingly small, suggesting that generating multiple internal states with such a mechanism is much harder than expected. While the shape of the region appears straightforward, we find surprising complexity in how the region grows with increasing amounts of the modified substrate. Our approach uses statistical analysis of data generated from a model, rather than from experiments, but leads to precise mathematical conjectures about parameter geography and biological robustness.

TopoFilter: a MATLAB package for mechanistic model identification in systems biology

Background

To develop mechanistic dynamic models in systems biology, one often needs to identify all (or minimal) representations of the biological processes that are consistent with experimental data, out of a potentially large set of hypothetical mechanisms. However, a simple enumeration of all alternatives becomes quickly intractable when the number of model parameters grows. Selecting appropriate dynamic models out of a large ensemble of models, taking the uncertainty in our biological knowledge and in the experimental data into account, is therefore a key current problem in systems biology.

Results

The TopoFilter package addresses this problem in a heuristic and automated fashion by implementing the previously described topological filtering method for Bayesian model selection. It includes a core heuristic for searching the space of submodels of a parametrized model, coupled with a sampling-based exploration of the parameter space. Recent developments of the method allow to balance exhaustiveness and speed of the model space search, to efficiently re-sample parameters, to parallelize the search, and to use custom scoring functions. We use a theoretical example to motivate these features and then demonstrate TopoFilter’s applicability for a yeast signaling network with more than 250’000 possible model structures.

Conclusions

TopoFilter is a flexible software framework that makes Bayesian model selection and reduction efficient and scalable to network models of a complexity that represents contemporary problems in, for example, cell signaling. TopoFilter is open-source, available under the GPL-3.0 license at https://gitlab.com/csb.ethz/TopoFilter. It includes installation instructions, a quickstart guide, a description of all package options, and multiple examples.

Uncertainty quantification, propagation and characterization by Bayesian analysis combined with global sensitivity analysis applied to dynamical intracellular pathway models

Motivation

Dynamical models describing intracellular phenomena are increasing in size and complexity as more information is obtained from experiments. These models are often over-parameterized with respect to the quantitative data used for parameter estimation, resulting in uncertainty in the individual parameter estimates as well as in the predictions made from the model. Here we combine Bayesian analysis with global sensitivity analysis (GSA) in order to give better informed predictions; to point out weaker parts of the model that are important targets for further experiments, as well as to give guidance on parameters that are essential in distinguishing different qualitative output behaviours.

Results

We used approximate Bayesian computation (ABC) to estimate the model parameters from experimental data, as well as to quantify the uncertainty in this estimation (inverse uncertainty quantification), resulting in a posterior distribution for the parameters. This parameter uncertainty was next propagated to a corresponding uncertainty in the predictions (forward uncertainty propagation), and a GSA was performed on the predictions using the posterior distribution as the possible values for the parameters. This methodology was applied on a relatively large model relevant for synaptic plasticity, using experimental data from several sources. We could hereby point out those parameters that by themselves have the largest contribution to the uncertainty of the prediction as well as identify parameters important to separate between qualitatively different predictions. This approach is useful both for experimental design as well as model building.

Availability and implementation

Source code is freely available at https://github.com/alexjau/uqsa.

Supplementary information

Supplementary data are available at Bioinformatics online.

Global Sensitivity Analysis of the SCOPE Model in Sentinel-3 Bands: Thermal Domain Focus

Sentinel-3 satellite has provided simultaneous observations in the optical (visible, near infrared (NIR), shortwave infrared (SWIR)) and thermal infrared (TIR) domains since 2016, with a revisit time of 1–2 days. The high temporal resolution and spectral coverage make the data of this mission attractive for vegetation monitoring. This study explores the possibilities of using the Soil Canopy Observation, Photochemistry and Energy fluxes (SCOPE) model together with Sentinel-3 to exploit the two sensors onboard of Sentinel-3 (the ocean and land color instrument (OLCI) and sea and land surface temperature radiometer (SLSTR)) in synergy. Sobol’ variance based global sensitivity analysis (GSA) of top of atmosphere (TOA) radiance produced with a coupled SCOPE-6S model was conducted for optical bands of OLCI and SLSTR, while another GSA of SCOPE was conducted for the land surface temperature (LST) product of SLSTR. The results show that in addition to ESA level-2 Sentinel-3 products, SCOPE is able to retrieve leaf area index (LAI), leaf chlorophyll content (Cab), leaf water content (Cw), leaf senescent material (Cs), leaf inclination distribution (LAD). Leaf dry matter content (Cdm) and soil brightness, despite being important, were not confidently retrieved in some cases. GSA of LST in TIR domain showed that plant biochemical parameters—maximum carboxylation rate (Vcmax) and stomata conductance-photosynthesis slope (Ball-Berry m)—can be constrained if prior information on near-surface weather conditions is available. We conclude that the combination of optical and thermal domains facilitates the constraint of the land surface energy balance using SCOPE.

Computational analysis of viable parameter regions in models of synthetic biological systems

Background

Gene regulatory networks with different topological and/or dynamical properties might exhibit similar behavior. System that is less perceptive for the perturbations of its internal and external factors should be preferred. Methods for sensitivity and robustness assessment have already been developed and can be roughly divided into local and global approaches. Local methods focus only on the local area around nominal parameter values. This can be problematic when parameters exhibits the desired behavior over a large range of parameter perturbations or when parameter values are unknown. Global methods, on the other hand, investigate the whole space of parameter values and mostly rely on different sampling techniques. This can be computationally inefficient. To address these shortcomings ’glocal’ approaches were developed that apply global and local approaches in an effective and rigorous manner.

Results

Herein, we present a computational approach for ’glocal’ analysis of viable parameter regions in biological models. The methodology is based on the exploration of high-dimensional viable parameter spaces with global and local sampling, clustering and dimensionality reduction techniques. The proposed methodology allows us to efficiently investigate the viable parameter space regions, evaluate the regions which exhibit the largest robustness, and to gather new insights regarding the size and connectivity of the viable parameter regions. We evaluate the proposed methodology on three different synthetic gene regulatory network models, i.e. the repressilator model, the model of the AC-DC circuit and the model of the edge-triggered master-slave D flip-flop.

Conclusions

The proposed methodology provides a rigorous assessment of the shape and size of viable parameter regions based on (1) the mathematical description of the biological system of interest, (2) constraints that define feasible parameter regions and (3) cost function that defines the desired or observed behavior of the system. These insights can be used to assess the robustness of biological systems, even in the case when parameter values are unknown and more importantly, even when there are multiple poorly connected viable parameter regions in the solution space. Moreover, the methodology can be efficiently applied to the analysis of biological systems that exhibit multiple modes of the targeted behavior.

Transcriptional Regulation Contributes to Prioritized Detoxification of Hydrogen Peroxide over Nitric Oxide

Hydrogen peroxide (H2O2) and nitric oxide (NO·) are toxic metabolites that immune cells use to attack pathogens. These antimicrobials can be present at the same time in phagosomes, and it remains unclear how bacteria deal with these insults when simultaneously present. Here, using Escherichia coli, we observed that simultaneous exposure to H2O2 and NO· leads to prioritized detoxification, where enzymatic removal of NO· is impeded until H2O2 has been eliminated. This phenomenon is reminiscent of carbon catabolite repression (CCR), where preferred carbon sources are catabolized prior to less desirable substrates; however, H2O2 and NO· are toxic, growth-inhibitory compounds rather than growth-promoting nutrients. To understand how NO· detoxification is delayed by H2O2 whereas H2O2 detoxification proceeds unimpeded, we confirmed that the effect depended on Hmp, which is the main NO· detoxification enzyme, and used an approach that integrated computational modeling and experimentation to delineate and test potential mechanisms. Plausible interactions included H2O2-dependent inhibition of hmp transcription and translation, direct inhibition of Hmp catalysis, and competition for reducing equivalents between Hmp and H2O2-degrading enzymes. Experiments illustrated that Hmp catalysis and NAD(P)H supply were not impaired by H2O2, whereas hmp transcription and translation were diminished. A dependence of this phenomenon on transcriptional regulation parallels CCR, and we found it to involve the transcriptional repressor NsrR. Collectively, these data suggest that bacterial regulation of growth inhibitor detoxification has similarities to the regulation of growth substrate consumption, which could have ramifications for infectious disease, bioremediation, and biocatalysis from inhibitor-containing feedstocks.IMPORTANCE Bacteria can be exposed to H2O2 and NO· concurrently within phagosomes. In such multistress situations, bacteria could have evolved to simultaneously degrade both toxic metabolites or preferentially detoxify one over the other. Here, we found that simultaneous exposure to H2O2 and NO· leads to prioritized detoxification, where detoxification of NO· is hampered until H2O2 has been eliminated. This phenomenon resembles CCR, where bacteria consume one substrate over others in carbon source mixtures. Further experimentation revealed a central role for transcriptional regulation in the prioritization of H2O2 over NO·, which is also important to CCR. This study suggests that regulatory scenarios observed in bacterial consumption of growth-promoting compound mixtures can be conserved in bacterial detoxification of toxic metabolite mixtures.

Comprehensive Review of Models and Methods for Inferences in Bio-Chemical Reaction Networks

The key processes in biological and chemical systems are described by networks of chemical reactions. From molecular biology to biotechnology applications, computational models of reaction networks are used extensively to elucidate their non-linear dynamics. The model dynamics are crucially dependent on the parameter values which are often estimated from observations. Over the past decade, the interest in parameter and state estimation in models of (bio-) chemical reaction networks (BRNs) grew considerably. The related inference problems are also encountered in many other tasks including model calibration, discrimination, identifiability, and checking, and optimum experiment design, sensitivity analysis, and bifurcation analysis. The aim of this review paper is to examine the developments in literature to understand what BRN models are commonly used, and for what inference tasks and inference methods. The initial collection of about 700 documents concerning estimation problems in BRNs excluding books and textbooks in computational biology and chemistry were screened to select over 270 research papers and 20 graduate research theses. The paper selection was facilitated by text mining scripts to automate the search for relevant keywords and terms. The outcomes are presented in tables revealing the levels of interest in different inference tasks and methods for given models in the literature as well as the research trends are uncovered. Our findings indicate that many combinations of models, tasks and methods are still relatively unexplored, and there are many new research opportunities to explore combinations that have not been considered-perhaps for good reasons. The most common models of BRNs in literature involve differential equations, Markov processes, mass action kinetics, and state space representations whereas the most common tasks are the parameter inference and model identification. The most common methods in literature are Bayesian analysis, Monte Carlo sampling strategies, and model fitting to data using evolutionary algorithms. The new research problems which cannot be directly deduced from the text mining data are also discussed.

Ensemble optimization of microbial conversion of glycerol into 1, 3-propanediol by Klebsiella pneumoniae

Using mathematical model and computer simulation to predict biological processes and optimize the target production is an important strategy for optimizing fermentation process. However, the inherent uncertainty of the kinetic model severely limits the predictive capability. In this study, optimize target production, such as productivity and yield of 1, 3-propanediol produced by Klebsiella pneumoniae using glycerol as substrate, the ensemble modeling approach was used to reduce the model's uncertainty for fermentation process as much as possible, and effectively improve its prediction performance. Firstly, through sensitivity analysis, the parameters having significant influence on the model were determined as the adjustable parameters for the ensemble modeling. After comparison, the appropriate threshold coefficient of the model error was determined, and the sampling method was used to generate as many equivalent parameter sets as possible. Each set of parameters was separately applied for the simulation, and all the predicted values were integrated for the weighted average. Therefore, the expected value of the prediction was obtained. Compared with the traditional simulation using single parameter set, the ensemble modeling method achieved the lower relative error between the prediction and the experimental value and the greatly improved model prediction performance. Moreover, the optimal productivity and yield of 1, 3-propanediol and the corresponding operating conditions were obtained, respectively. The ensemble modeling approach effectively compensates for the uncertainties of the model, making its prediction performance more practical, which is important for computer simulations to predict and guide the actual production process.

Application of Parameter Optimization to Search for Oscillatory Mass-Action Networks Using Python

Biological systems can be described mathematically to model the dynamics of metabolic, protein, or gene-regulatory networks, but locating parameter regimes that induce a particular dynamic behavior can be challenging due to the vast parameter landscape, particularly in large models. In the current work, a Pythonic implementation of existing bifurcation objective functions, which reward systems that achieve a desired bifurcation behavior, is implemented to search for parameter regimes that permit oscillations or bistability. A differential evolution algorithm progressively approximates the specified bifurcation type while performing a global search of parameter space for a candidate with the best fitness. The user-friendly format facilitates integration with systems biology tools, as Python is a ubiquitous programming language. The bifurcation–evolution software is validated on published models from the BioModels Database and used to search populations of randomly-generated mass-action networks for oscillatory dynamics. Results of this search demonstrate the importance of reaction enrichment to provide flexibility and enable complex dynamic behaviors, and illustrate the role of negative feedback and time delays in generating oscillatory dynamics.

Parameter estimation in models of biological oscillators: an automated regularised estimation approach

BACKGROUND

Dynamic modelling is a core element in the systems biology approach to understanding complex biosystems. Here, we consider the problem of parameter estimation in models of biological oscillators described by deterministic nonlinear differential equations. These problems can be extremely challenging due to several common pitfalls: (i) a lack of prior knowledge about parameters (i.e. massive search spaces), (ii) convergence to local optima (due to multimodality of the cost function), (iii) overfitting (fitting the noise instead of the signal) and (iv) a lack of identifiability. As a consequence, the use of standard estimation methods (such as gradient-based local ones) will often result in wrong solutions. Overfitting can be particularly problematic, since it produces very good calibrations, giving the impression of an excellent result. However, overfitted models exhibit poor predictive power. Here, we present a novel automated approach to overcome these pitfalls. Its workflow makes use of two sequential optimisation steps incorporating three key algorithms: (1) sampling strategies to systematically tighten the parameter bounds reducing the search space, (2) efficient global optimisation to avoid convergence to local solutions, (3) an advanced regularisation technique to fight overfitting. In addition, this workflow incorporates tests for structural and practical identifiability.

RESULTS

We successfully evaluate this novel approach considering four difficult case studies regarding the calibration of well-known biological oscillators (Goodwin, FitzHugh-Nagumo, Repressilator and a metabolic oscillator). In contrast, we show how local gradient-based approaches, even if used in multi-start fashion, are unable to avoid the above-mentioned pitfalls.

CONCLUSIONS

Our approach results in more efficient estimations (thanks to the bounding strategy) which are able to escape convergence to local optima (thanks to the global optimisation approach). Further, the use of regularisation allows us to avoid overfitting, resulting in more generalisable calibrated models (i.e. models with greater predictive power).

Distilling Robust Design Principles of Biocircuits Using Mixed Integer Dynamic Optimization

A major challenge in model-based design of synthetic biochemical circuits is how to address uncertainty in the parameters. A circuit whose behavior is robust to variations in the parameters will have more chances to behave as predicted when implemented in practice, and also to function reliably in presence of fluctuations and noise. Here, we extend our recent work on automated-design based on mixed-integer multi-criteria dynamic optimization to take into account parametric uncertainty. We exploit the intensive sampling of the design space performed by a global optimization algorithm to obtain the robustness of the topologies without significant additional computational effort. Our procedure provides automatically topologies that best trade-off performance and robustness against parameter fluctuations. We illustrate our approach considering the automated design of gene circuits achieving adaptation.

Exploration of sexual dimorphism and inter-individual variability in multivariate parameter spaces for a pharmacokinetic compartment model

Pharmacokinetic models are particularly useful to study the underlying and complex physiological mechanisms contributing to clinical differences across patient subgroups or special populations. Unfortunately, the inherent variability of biological systems and knowledge gaps in physiological data limit confidence in model predictions for special populations. Sourcing data to reflect the desired physiologies can be resource intensive, particularly for a larger model. Thus, a critical step in model development for special populations involves an in-depth analysis of model inputs, which can be guided by Monte Carlo simulations. Such an approach enables the generation of parameter values by stochastic sampling that are subsequently restricted to the combinations that describe biologically plausible or target model output. Our approach utilized a published pharmacokinetic compartmental model to demonstrate how sampling in conjunction with global sensitivity analysis can be used to explore sexual dimorphism and inter-individual variability in multivariate parameter spaces for differentiation of model input and behavior across phenotypes. Despite limiting the model output to relatively narrow ranges, male and female phenotypes were associated with wide variability in both individual parameter values and combinations of parameters. Through an integrated approach using a support vector machine, principal component analysis and global sensitivity analysis, our approach revealed that specific combinations of parameters gave rise to a certain phenotype, while individual parameters influenced the shape of plasma concentration profile. Augmenting analysis of the model input with global sensitivity analysis enabled an understanding of both sexual dimorphism and inter-individual variability in pharmacokinetics. While the current study revealed how model input could be separated by sex for a simple compartment model, the approach could be extended to other patient factors, such as age or disease, and to a more complex physiologically-based model that describes absorption, distribution, metabolism, and elimination with more detail.

Bayesian uncertainty analysis for complex systems biology models: emulation, global parameter searches and evaluation of gene functions

Background

Many mathematical models have now been employed across every area of systems biology. These models increasingly involve large numbers of unknown parameters, have complex structure which can result in substantial evaluation time relative to the needs of the analysis, and need to be compared to observed data of various forms. The correct analysis of such models usually requires a global parameter search, over a high dimensional parameter space, that incorporates and respects the most important sources of uncertainty. This can be an extremely difficult task, but it is essential for any meaningful inference or prediction to be made about any biological system. It hence represents a fundamental challenge for the whole of systems biology.

Methods

Bayesian statistical methodology for the uncertainty analysis of complex models is introduced, which is designed to address the high dimensional global parameter search problem. Bayesian emulators that mimic the systems biology model but which are extremely fast to evaluate are embeded within an iterative history match: an efficient method to search high dimensional spaces within a more formal statistical setting, while incorporating major sources of uncertainty.

Results

The approach is demonstrated via application to a model of hormonal crosstalk in Arabidopsis root development, which has 32 rate parameters, for which we identify the sets of rate parameter values that lead to acceptable matches between model output and observed trend data. The multiple insights into the model’s structure that this analysis provides are discussed. The methodology is applied to a second related model, and the biological consequences of the resulting comparison, including the evaluation of gene functions, are described.

Conclusions

Bayesian uncertainty analysis for complex models using both emulators and history matching is shown to be a powerful technique that can greatly aid the study of a large class of systems biology models. It both provides insight into model behaviour and identifies the sets of rate parameters of interest.

Electronic supplementary material

The online version of this article (doi:10.1186/s12918-017-0484-3) contains supplementary material, which is available to authorized users.

Modeling the Sex Differences and Interindividual Variability in the Activity of the Hypothalamic-Pituitary-Adrenal Axis

Significant sex differences exist in the activity of the hypothalamic-pituitary-adrenal (HPA) axis. These differences are thought to contribute to the disparity in the prevalence of various autoimmune and infectious diseases between males and females. We used a mathematical model of the HPA axis to evaluate the hypothesis that differential sensitivity and negative feedback of the HPA axis network are causal factors for the observed sex differences in its activity. In doing so, we implicitly accounted for the differential influence of gonadal hormones on the HPA axis. Furthermore, we determined whether the putative mechanisms responsible for differences in basal HPA axis activity might also contribute to the observed differences in its stimulus-driven response. Model simulations predicted that the female HPA axis has greater adrenal sensitivity and weaker negative feedback. We identified two distinct sex-specific parameter spaces that generate corticosterone profiles in qualitative agreement with experimental results. We propose that these parameter subspaces indicate the interindividual variability in the regulatory mechanisms of the HPA axis. Furthermore, the model predicts that the maintenance of homeostatic rhythms in response to chronic stress requires specific regulatory adaptations resulting in a phenotype of allostatically driven chronic stress-sensitization. We propose that these adaptations indicate a physiological cost of adaptation to chronic stress. Model simulations suggest that individuals with high adrenal sensitivity are more vulnerable to chronic stress sensitization and might be more susceptible to the development of neuropsychiatric disorders. These results contribute to the study of sex differences in physiological feedback systems within a quantitative framework.

Formulation, construction and analysis of kinetic models of metabolism: A review of modelling frameworks

Kinetic models are critical to predict the dynamic behaviour of metabolic networks. Mechanistic kinetic models for large networks remain uncommon due to the difficulty of fitting their parameters. Recent modelling frameworks promise new ways to overcome this obstacle while retaining predictive capabilities. In this review, we present an overview of the relevant mathematical frameworks for kinetic formulation, construction and analysis. Starting with kinetic formalisms, we next review statistical methods for parameter inference, as well as recent computational frameworks applied to the construction and analysis of kinetic models. Finally, we discuss opportunities and limitations hindering the development of larger kinetic reconstructions.

The Design Space of the Embryonic Cell Cycle Oscillator

One of the main tasks in the analysis of models of biomolecular networks is to characterize the domain of the parameter space that corresponds to a specific behavior. Given the large number of parameters in most models, this is no trivial task. We use a model of the embryonic cell cycle to illustrate the approaches that can be used to characterize the domain of parameter space corresponding to limit cycle oscillations, a regime that coordinates periodic entry into and exit from mitosis. Our approach relies on geometric construction of bifurcation sets, numerical continuation, and random sampling of parameters. We delineate the multidimensional oscillatory domain and use it to quantify the robustness of periodic trajectories. Although some of our techniques explore the specific features of the chosen system, the general approach can be extended to other models of the cell cycle engine and other biomolecular networks.

How to deal with parameters for whole-cell modelling

Dynamical systems describing whole cells are on the verge of becoming a reality. But as models of reality, they are only useful if we have realistic parameters for the molecular reaction rates and cell physiological processes. There is currently no suitable framework to reliably estimate hundreds, let alone thousands, of reaction rate parameters. Here, we map out the relative weaknesses and promises of different approaches aimed at redressing this issue. While suitable procedures for estimation or inference of the whole (vast) set of parameters will, in all likelihood, remain elusive, some hope can be drawn from the fact that much of the cellular behaviour may be explained in terms of smaller sets of parameters. Identifying such parameter sets and assessing their behaviour is now becoming possible even for very large systems of equations, and we expect such methods to become central tools in the development and analysis of whole-cell models.

L p -Adaptation: Simultaneous Design Centering and Robustness Estimation of Electronic and Biological Systems

The design of systems or models that work robustly under uncertainty and environmental fluctuations is a key challenge in both engineering and science. This is formalized in the design-centering problem, which is defined as finding a design that fulfills given specifications and has a high probability of still doing so if the system parameters or the specifications fluctuate randomly. Design centering is often accompanied by the problem of quantifying the robustness of a system. Here we present a novel adaptive statistical method to simultaneously address both problems. Our method, L_p-Adaptation, is inspired by the evolution of robustness in biological systems and by randomized schemes for convex volume computation. It is able to address both problems in the general, non-convex case and at low computational cost. We describe the concept and the algorithm, test it on known benchmarks, and demonstrate its real-world applicability in electronic and biological systems. In all cases, the present method outperforms the previous state of the art. This enables re-formulating optimization problems in engineering and biology as design centering problems, taking global system robustness into account.

An integrated network analysis reveals that nitric oxide reductase prevents metabolic cycling of nitric oxide by Pseudomonas aeruginosa

Nitric oxide (NO) is a chemical weapon within the arsenal of immune cells, but is also generated endogenously by different bacteria. Pseudomonas aeruginosa are pathogens that contain an NO-generating nitrite (NO₂^-) reductase (NirS), and NO has been shown to influence their virulence. Interestingly, P. aeruginosa also contain NO dioxygenase (Fhp) and nitrate (NO₃^-) reductases, which together with NirS provide the potential for NO to be metabolically cycled (NO→NO₃^-→NO₂^-→NO). Deeper understanding of NO metabolism in P. aeruginosa will increase knowledge of its pathogenesis, and computational models have proven to be useful tools for the quantitative dissection of NO biochemical networks. Here we developed such a model for P. aeruginosa and confirmed its predictive accuracy with measurements of NO, O₂, NO₂^-, and NO₃^- in mutant cultures devoid of Fhp or NorCB (NO reductase) activity. Using the model, we assessed whether NO was metabolically cycled in aerobic P. aeruginosa cultures. Calculated fluxes indicated a bottleneck at NO₃^-, which was relieved upon O₂ depletion. As cell growth depleted dissolved O₂ levels, NO₃^- was converted to NO₂^- at near-stoichiometric levels, whereas NO₂^- consumption did not coincide with NO or NO₃^- accumulation. Assimilatory NO₂^- reductase (NirBD) or NorCB activity could have prevented NO cycling, and experiments with ΔnirB, ΔnirS, and ΔnorC showed that NorCB was responsible for loss of flux from the cycle. Collectively, this work provides a computational tool to analyze NO metabolism in P. aeruginosa, and establishes that P. aeruginosa use NorCB to prevent metabolic cycling of NO.

Bioprocess optimization under uncertainty using ensemble modeling

The performance of model-based bioprocess optimizations depends on the accuracy of the mathematical model. However, models of bioprocesses often have large uncertainty due to the lack of model identifiability. In the presence of such uncertainty, process optimizations that rely on the predictions of a single "best fit" model, e.g. the model resulting from a maximum likelihood parameter estimation using the available process data, may perform poorly in real life. In this study, we employed ensemble modeling to account for model uncertainty in bioprocess optimization. More specifically, we adopted a Bayesian approach to define the posterior distribution of the model parameters, based on which we generated an ensemble of model parameters using a uniformly distributed sampling of the parameter confidence region. The ensemble-based process optimization involved maximizing the lower confidence bound of the desired bioprocess objective (e.g. yield or product titer), using a mean-standard deviation utility function. We demonstrated the performance and robustness of the proposed strategy in an application to a monoclonal antibody batch production by mammalian hybridoma cell culture.

Biphasic regulation of the transcription factor ABORTED MICROSPORES (AMS) is essential for tapetum and pollen development in Arabidopsis

Viable pollen is essential for plant reproduction and crop yield. Its production requires coordinated expression at specific stages during anther development, involving early meiosis-associated events and late pollen wall formation. The ABORTED MICROSPORES (AMS) transcription factor is a master regulator of sporopollenin biosynthesis, secretion and pollen wall formation in Arabidopsis. Here we show that it has complex regulation and additional essential roles earlier in pollen formation. An inducible-AMS reporter was created for functional rescue, protein expression pattern analysis, and to distinguish between direct and indirect targets. Mathematical modelling was used to create regulatory networks based on wild-type RNA and protein expression. Dual activity of AMS was defined by biphasic protein expression in anther tapetal cells, with an initial peak around pollen meiosis and then later during pollen wall development. Direct AMS-regulated targets exhibit temporal regulation, indicating that additional factors are associated with their regulation. We demonstrate that AMS biphasic expression is essential for pollen development, and defines distinct functional activities during early and late pollen development. Mathematical modelling suggests that AMS may competitively form a protein complex with other tapetum-expressed transcription factors, and that biphasic regulation is due to repression of upstream regulators and promotion of AMS protein degradation.

Feedback, Mass Conservation and Reaction Kinetics Impact the Robustness of Cellular Oscillations

Oscillations occur in a wide variety of cellular processes, for example in calcium and p53 signaling responses, in metabolic pathways or within gene-regulatory networks, e.g. the circadian system. Since it is of central importance to understand the influence of perturbations on the dynamics of these systems a number of experimental and theoretical studies have examined their robustness. The period of circadian oscillations has been found to be very robust and to provide reliable timing. For intracellular calcium oscillations the period has been shown to be very sensitive and to allow for frequency-encoded signaling. We here apply a comprehensive computational approach to study the robustness of period and amplitude of oscillatory systems. We employ different prototype oscillator models and a large number of parameter sets obtained by random sampling. This framework is used to examine the effect of three design principles on the sensitivities towards perturbations of the kinetic parameters. We find that a prototype oscillator with negative feedback has lower period sensitivities than a prototype oscillator relying on positive feedback, but on average higher amplitude sensitivities. For both oscillator types, the use of Michaelis-Menten instead of mass action kinetics in all degradation and conversion reactions leads to an increase in period as well as amplitude sensitivities. We observe moderate changes in sensitivities if replacing mass conversion reactions by purely regulatory reactions. These insights are validated for a set of established models of various cellular rhythms. Overall, our work highlights the importance of reaction kinetics and feedback type for the variability of period and amplitude and therefore for the establishment of predictive models.

Rhythmic behavior is omnipresent in biology and has many crucial functions. In cells the activation levels and abundances of signaling molecules such as NF-κB, p53, EGFR or calcium repeatedly increase and decrease in response to stimuli. Such a dynamic behavior can also be observed monitoring the concentrations of mRNAs and proteins in the circadian clock and the cell cycle. Period and amplitude which are the time span between peaks and the peak height, respectively, as well as their variabilities are important features of oscillations. The circadian period is very stable allowing for proper time keeping, whereas in calcium signaling the period is very variable encoding different stimulation strengths. Our goal is to examine the origin of differences in sensitivities of periods and amplitudes using a computational approach. We use prototype oscillators and demonstrate that they can be used to derive general principles that explain the degree of robustness in period and amplitude for a set of commonly used models of cellular oscillators. Our findings imply that the robustness of oscillating systems can be influenced by feedback type and kinetic properties to which special attention should be paid when designing mathematical models of cellular rhythms.

An Orthogonal Permease-Inducer-Repressor Feedback Loop Shows Bistability

Feedback loops in biological networks, among others, enable differentiation and cell cycle progression, and increase robustness in signal transduction. In natural networks, feedback loops are often complex and intertwined, making it challenging to identify which loops are mainly responsible for an observed behavior. However, minimal synthetic replicas could allow for such identification. Here, we engineered a synthetic permease-inducer-repressor system in Saccharomyces cerevisiae to analyze if a transport-mediated positive feedback loop could be a core mechanism for the switch-like behavior in the regulation of metabolic gene networks such as the S. cerevisiae GAL system or the Escherichia coli lac operon. We characterized the synthetic circuit using deterministic and stochastic mathematical models. Similar to its natural counterparts, our synthetic system shows bistable and hysteretic behavior, and the inducer concentration range for bistability as well as the switching rates between the two stable states depend on the repressor concentration. Our results indicate that a generic permease-inducer-repressor circuit with a single feedback loop is sufficient to explain the experimentally observed bistable behavior of the natural systems. We anticipate that the approach of reimplementing natural systems with orthogonal parts to identify crucial network components is applicable to other natural systems such as signaling pathways.

A Statistical Approach Reveals Designs for the Most Robust Stochastic Gene Oscillators

The engineering of transcriptional networks presents many challenges due to the inherent uncertainty in the system structure, changing cellular context, and stochasticity in the governing dynamics. One approach to address these problems is to design and build systems that can function across a range of conditions; that is they are robust to uncertainty in their constituent components. Here we examine the parametric robustness landscape of transcriptional oscillators, which underlie many important processes such as circadian rhythms and the cell cycle, plus also serve as a model for the engineering of complex and emergent phenomena. The central questions that we address are: Can we build genetic oscillators that are more robust than those already constructed? Can we make genetic oscillators arbitrarily robust? These questions are technically challenging due to the large model and parameter spaces that must be efficiently explored. Here we use a measure of robustness that coincides with the Bayesian model evidence, combined with an efficient Monte Carlo method to traverse model space and concentrate on regions of high robustness, which enables the accurate evaluation of the relative robustness of gene network models governed by stochastic dynamics. We report the most robust two and three gene oscillator systems, plus examine how the number of interactions, the presence of autoregulation, and degradation of mRNA and protein affects the frequency, amplitude, and robustness of transcriptional oscillators. We also find that there is a limit to parametric robustness, beyond which there is nothing to be gained by adding additional feedback. Importantly, we provide predictions on new oscillator systems that can be constructed to verify the theory and advance design and modeling approaches to systems and synthetic biology.

Discovery and dissection of metabolic oscillations in the microaerobic nitric oxide response network of Escherichia coli

The virulence of many pathogens depends upon their ability to cope with immune-generated nitric oxide (NO·). In Escherichia coli, the major NO· detoxification systems are Hmp, an NO· dioxygenase (NOD), and NorV, an NO· reductase (NOR). It is well established that Hmp is the dominant system under aerobic conditions, whereas NorV dominates anaerobic conditions; however, the quantitative contributions of these systems under the physiologically relevant microaerobic regime remain ill defined. Here, we investigated NO· detoxification in environments ranging from 0 to 50 μM O2, and discovered a regime in which E. coli NO· defenses were severely compromised, as well as conditions that exhibited oscillations in the concentration of NO·. Using an integrated computational and experimental approach, E. coli NO· detoxification was found to be extremely impaired at low O2 due to a combination of its inhibitory effects on NorV, Hmp, and translational activities, whereas oscillations were found to result from a kinetic competition for O2 between Hmp and respiratory cytochromes. Because at least 777 different bacterial species contain the genetic requirements of this stress response oscillator, we hypothesize that such oscillatory behavior could be a widespread phenomenon. In support of this hypothesis,Pseudomonas aeruginosa, whose respiratory and NO· response networks differ considerably from those of E. coli, was found to exhibit analogous oscillations in low O2 environments. This work provides insight into how bacterial NO· defenses function under the low O2 conditions that are likely to be encountered within host environments.

Construction and Experimental Validation of a Quantitative Kinetic Model of Nitric Oxide Stress in Enterohemorrhagic Escherichia coli O157:H7

Enterohemorrhagic Escherichia coli (EHEC) are responsible for large outbreaks of hemorrhagic colitis, which can progress to life-threatening hemolytic uremic syndrome (HUS) due to the release of Shiga-like toxins (Stx). The presence of a functional nitric oxide (NO·) reductase (NorV), which protects EHEC from NO· produced by immune cells, was previously found to correlate with high HUS incidence, and it was shown that NorV activity enabled prolonged EHEC survival and increased Stx production within macrophages. To enable quantitative study of EHEC NO· defenses and facilitate the development of NO·-potentiating therapeutics, we translated an existing kinetic model of the E. coli K-12 NO· response to an EHEC O157:H7 strain. To do this, we trained uncertain model parameters on measurements of [NO·] and [O₂] in EHEC cultures, assessed parametric and prediction uncertainty with the use of a Markov chain Monte Carlo approach, and confirmed the predictive accuracy of the model with experimental data from genetic mutants lacking NorV or Hmp (NO· dioxygenase). Collectively, these results establish a methodology for the translation of quantitative models of NO· stress in model organisms to pathogenic sub-species, which is a critical step toward the application of these models for the study of infectious disease.

A Kinetic Platform to Determine the Fate of Hydrogen Peroxide in Escherichia coli

Hydrogen peroxide (H2O2) is used by phagocytic cells of the innate immune response to kill engulfed bacteria. H2O2 diffuses freely into bacteria, where it can wreak havoc on sensitive biomolecules if it is not rapidly detoxified. Accordingly, bacteria have evolved numerous systems to defend themselves against H2O2, and the importance of these systems to pathogenesis has been substantiated by the many bacteria that require them to establish or sustain infections. The kinetic competition for H2O2 within bacteria is complex, which suggests that quantitative models will improve interpretation and prediction of network behavior. To date, such models have been of limited scope, and this inspired us to construct a quantitative, systems-level model of H2O2 detoxification in Escherichia coli that includes detoxification enzymes, H2O2-dependent transcriptional regulation, enzyme degradation, the Fenton reaction and damage caused by •OH, oxidation of biomolecules by H2O2, and repair processes. After using an iterative computational and experimental procedure to train the model, we leveraged it to predict how H2O2 detoxification would change in response to an environmental perturbation that pathogens encounter within host phagosomes, carbon source deprivation, which leads to translational inhibition and limited availability of NADH. We found that the model accurately predicted that NADH depletion would delay clearance at low H2O2 concentrations and that detoxification at higher concentrations would resemble that of carbon-replete conditions. These results suggest that protein synthesis during bolus H2O2 stress does not affect clearance dynamics and that access to catabolites only matters at low H2O2 concentrations. We anticipate that this model will serve as a computational tool for the quantitative exploration and dissection of oxidative stress in bacteria, and that the model and methods used to develop it will provide important templates for the generation of comparable models for other bacterial species.

Efficient Characterization of Parametric Uncertainty of Complex (Bio)chemical Networks

Parametric uncertainty is a particularly challenging and relevant aspect of systems analysis in domains such as systems biology where, both for inference and for assessing prediction uncertainties, it is essential to characterize the system behavior globally in the parameter space. However, current methods based on local approximations or on Monte-Carlo sampling cope only insufficiently with high-dimensional parameter spaces associated with complex network models. Here, we propose an alternative deterministic methodology that relies on sparse polynomial approximations. We propose a deterministic computational interpolation scheme which identifies most significant expansion coefficients adaptively. We present its performance in kinetic model equations from computational systems biology with several hundred parameters and state variables, leading to numerical approximations of the parametric solution on the entire parameter space. The scheme is based on adaptive Smolyak interpolation of the parametric solution at judiciously and adaptively chosen points in parameter space. As Monte-Carlo sampling, it is “non-intrusive” and well-suited for massively parallel implementation, but affords higher convergence rates. This opens up new avenues for large-scale dynamic network analysis by enabling scaling for many applications, including parameter estimation, uncertainty quantification, and systems design.

In various scientific domains, in particular in systems biology, dynamic mathematical models of increasing complexity are being developed and analyzed to study biochemical reaction networks. A major challenge in dealing with such models is the uncertainty in parameters such as kinetic constants; how to efficiently and precisely quantify the effects of parametric uncertainties on systems behavior remains a question. Addressing this computational challenge for large systems, with good scaling up to hundreds of species and kinetic parameters, is important for many forward (e.g., uncertainty quantification) and inverse (e.g., system identification) problems. Here, we propose a sparse, deterministic adaptive interpolation method tailored to high-dimensional parametric problems that allows for fast, deterministic computational analysis of large biochemical reaction networks. The method is based on adaptive Smolyak interpolation of the parametric solution at judiciously chosen points in high-dimensional parameter space, combined with adaptive time-stepping for the actual numerical simulation of the network dynamics. It is “non-intrusive” and well-suited both for massively parallel implementation and for use in standard (systems biology) toolboxes.

Stochastic simulation of Boolean rxncon models: towards quantitative analysis of large signaling networks

Background

Cellular decision-making is governed by molecular networks that are highly complex. An integrative understanding of these networks on a genome wide level is essential to understand cellular health and disease. In most cases however, such an understanding is beyond human comprehension and requires computational modeling. Mathematical modeling of biological networks at the level of biochemical details has hitherto relied on state transition models. These are typically based on enumeration of all relevant model states, and hence become very complex unless severely – and often arbitrarily – reduced. Furthermore, the parameters required for genome wide networks will remain underdetermined for the conceivable future. Alternatively, networks can be simulated by Boolean models, although these typically sacrifice molecular detail as well as distinction between different levels or modes of activity. However, the modeling community still lacks methods that can simulate genome scale networks on the level of biochemical reaction detail in a quantitative or semi quantitative manner.

Results

Here, we present a probabilistic bipartite Boolean modeling method that addresses these issues. The method is based on the reaction-contingency formalism, and enables fast simulation of large networks. We demonstrate its scalability by applying it to the yeast mitogen-activated protein kinase (MAPK) network consisting of 140 proteins and 608 nodes.

Conclusion

The probabilistic Boolean model can be generated and parameterized automatically from a rxncon network description, using only two global parameters, and its qualitative behavior is robust against order of magnitude variation in these parameters. Our method can hence be used to simulate the outcome of large signal transduction network reconstruction, with little or no overhead in model creation or parameterization.

Automated design of programmable enzyme-driven DNA circuits

Molecular programming allows for the bottom-up engineering of biochemical reaction networks in a controlled in vitro setting. These engineered biochemical reaction networks yield important insight in the design principles of biological systems and can potentially enrich molecular diagnostic systems. The DNA polymerase-nickase-exonuclease (PEN) toolbox has recently been used to program oscillatory and bistable biochemical networks using a minimal number of components. Previous work has reported the automatic construction of in silico descriptions of biochemical networks derived from the PEN toolbox, paving the way for generating networks of arbitrary size and complexity in vitro. Here, we report an automated approach that further bridges the gap between an in silico description and in vitro realization. A biochemical network of arbitrary complexity can be globally screened for parameter values that display the desired function and combining this approach with robustness analysis further increases the chance of successful in vitro implementation. Moreover, we present an automated design procedure for generating optimal DNA sequences, exhibiting key characteristics deduced from the in silico analysis. Our in silico method has been tested on a previously reported network, the Oligator, and has also been applied to the design of a reaction network capable of displaying adaptation in one of its components. Finally, we experimentally characterize unproductive sequestration of the exonuclease to phosphorothioate protected ssDNA strands. The strong nonlinearities in the degradation of active components caused by this unintended cross-coupling are shown computationally to have a positive effect on adaptation quality.

REDEMPTION: reduced dimension ensemble modeling and parameter estimation

Summary: Here, we present REDEMPTION (Reduced Dimension Ensemble Modeling and Parameter estimation), a toolbox for parameter estimation and ensemble modeling of ordinary differential equations (ODEs) using time-series data. For models with more reactions than measured species, a common scenario in biological modeling, the parameter estimation is formulated as a nested optimization problem based on incremental parameter estimation strategy. REDEMPTION also includes a tool for the identification of an ensemble of parameter combinations that provide satisfactory goodness-of-fit to the data. The functionalities of REDEMPTION are accessible through a MATLAB user interface (UI), as well as through programming script. For computational speed-up, REDEMPTION provides a numerical parallelization option using MATLAB Parallel Computing toolbox.

Availability and implementation: REDEMPTION can be downloaded from http://www.cabsel.ethz.ch/tools/redemption.

Contact:rudi.gunawan@chem.ethz.ch

Causal drift, robust signaling, and complex disease

The phenotype of many regulatory circuits in which mutations can cause complex, polygenic diseases is to some extent robust to DNA mutations that affect circuit components. Here I demonstrate how such mutational robustness can prevent the discovery of genetic disease determinants. To make my case, I use a mathematical model of the insulin signaling pathway implicated in type 2 diabetes, whose signaling output is governed by 15 genetically determined parameters. Using multiple complementary measures of a parameter's importance for this phenotype, I show that any one disease determinant that is crucial in one genetic background will be virtually irrelevant in other backgrounds. In an evolving population that drifts through the parameter space of this or other robust circuits through DNA mutations, the genetic changes that can cause disease will vary randomly over time. I call this phenomenon causal drift. It means that mutations causing disease in one (human or non-human) population may have no effect in another population, and vice versa. Causal drift casts doubt on our ability to infer the molecular mechanisms of complex diseases from non-human model organisms.

Quantitative Analysis of Robustness of Dynamic Response and Signal Transfer in Insulin mediated PI3K/AKT Pathway

Robustness is a critical feature of signaling pathways ensuring signal propagation with high fidelity in the event of perturbations. Here we present a detailed quantitative analysis of robustness in insulin mediated PI3K/AKT pathway, a critical signaling pathway maintaining self-renewal in human embryonic stem cells. Using global sensitivity analysis, we identified robustness promoting mechanisms that ensure (1) maintenance of a first order or overshoot dynamics of self-renewal molecule, p-AKT and (2) robust transfer of signals from oscillatory insulin stimulus to p-AKT in the presence of noise. Our results indicate that negative feedback controls the robustness to most perturbations. Faithful transfer of signal from the stimulating ligand to p-AKT occurs even in the presence of noise, albeit with signal attenuation and high frequency cut-off. Negative feedback contributes to signal attenuation, while positive regulators upstream of PIP₃ contribute to signal amplification. These results establish precise mechanisms to modulate self-renewal molecules like p-AKT.

A method for inverse bifurcation of biochemical switches: inferring parameters from dose response curves

Background

Within cells, stimuli are transduced into cell responses by complex networks of biochemical reactions. In many cell decision processes the underlying networks behave as bistable switches, converting graded stimuli or inputs into all or none cell responses. Observing how systems respond to different perturbations, insight can be gained into the underlying molecular mechanisms by developing mathematical models. Emergent properties of systems, like bistability, can be exploited to this purpose. One of the main challenges in modeling intracellular processes, from signaling pathways to gene regulatory networks, is to deal with high structural and parametric uncertainty, due to the complexity of the systems and the difficulty to obtain experimental measurements. Formal methods that exploit structural properties of networks for parameter estimation can help to overcome these problems.

Results

We here propose a novel method to infer the kinetic parameters of bistable biochemical network models. Bistable systems typically show hysteretic dose response curves, in which the so called bifurcation points can be located experimentally. We exploit the fact that, at the bifurcation points, a condition for multistationarity derived in the context of the Chemical Reaction Network Theory must be fulfilled. Chemical Reaction Network Theory has attracted attention from the (systems) biology community since it connects the structure of biochemical reaction networks to qualitative properties of the corresponding model of ordinary differential equations. The inverse bifurcation method developed here allows determining the parameters that produce the expected behavior of the dose response curves and, in particular, the observed location of the bifurcation points given by experimental data.

Conclusions

Our inverse bifurcation method exploits inherent structural properties of bistable switches in order to estimate kinetic parameters of bistable biochemical networks, opening a promising route for developments in Chemical Reaction Network Theory towards kinetic model identification.

Electronic supplementary material

The online version of this article (doi:10.1186/s12918-014-0114-2) contains supplementary material, which is available to authorized users.

Predictive sulfur metabolism - a field in flux

The key role of sulfur metabolites in response to biotic and abiotic stress in plants, as well as their importance in diet and health has led to a significant interest and effort in trying to understand and manipulate the production of relevant compounds. Metabolic engineering utilizes a set of theoretical tools to help rationally design modifications that enhance the production of a desired metabolite. Such approaches have proven their value in bacterial systems, however, the paucity of success stories to date in plants, suggests that challenges remain. Here, we review the most commonly used methods for understanding metabolic flux, focusing on the sulfur assimilatory pathway. We highlight known issues with both experimental and theoretical approaches, as well as presenting recent methods for integrating different modeling strategies, and progress toward an understanding of flux at the whole plant level.

Advancing metabolic models with kinetic information

Kinetic models are crucial to quantitatively understand and predict how functional behavior emerges from dynamic concentration changes of cellular components. The current challenge is on resolving uncertainties about parameter values of reaction kinetics. Additionally, there are also major structural uncertainties due to unknown molecular interactions and only putatively assigned regulatory functions. What if one or few key regulators of biochemical reactions are missing in a metabolic model? By reviewing current advances in building kinetic models of metabolism, we found that such models experience a paradigm shift away from fitting parameters towards identifying key regulatory interactions.