Symbolic kinetic models in python (SKiMpy): intuitive modeling of large-scale biological kinetic models.

MOTIVATION: Large-scale kinetic models are an invaluable tool to understand the dynamic and adaptive responses of biological systems. The development and application of these models have been limited by the availability of computational tools to build and analyze large-scale models efficiently. The toolbox presented here provides the means to implement, parameterize and analyze large-scale kinetic models intuitively and efficiently. RESULTS: We present a Python package (SKiMpy) bridging this gap by implementing an efficient kinetic modeling toolbox for the semiautomatic generation and analysis of large-scale kinetic models for various biological domains such as signaling, gene expression and metabolism. Furthermore, we demonstrate how this toolbox is used to parameterize kinetic models around a steady-state reference efficiently. Finally, we show how SKiMpy can implement multispecies bioreactor simulations to assess biotechnological processes. AVAILABILITY AND IMPLEMENTATION: The software is available as a Python 3 package on GitHub: https://github.com/EPFL-LCSB/SKiMpy, along with adequate documentation. SUPPLEMENTARY INFORMATION: Supplementary data are available at Bioinformatics online.

The effects of model complexity and size on metabolic flux distribution and control: case study in Escherichia coli.

BACKGROUND: Significant efforts have been made in building large-scale kinetic models of cellular metabolism in the past two decades. However, most kinetic models published to date, remain focused around central carbon pathways or are built around ad hoc reduced models without clear justification on their derivation and usage. Systematic algorithms exist for reducing genome-scale metabolic reconstructions to build thermodynamically feasible and consistently reduced stoichiometric models. However, it is important to study how network complexity affects conclusions derived from large-scale kinetic models built around consistently reduced models before we can apply them to study biological systems. RESULTS: We reduced the iJO1366 Escherichia Coli genome-scale metabolic reconstruction systematically to build three stoichiometric models of different size. Since the reduced models are expansions around the core subsystems for which the reduction was performed, the models are nested. We present a method for scaling up the flux profile and the concentration vector reference steady-states from the smallest model to the larger ones, whilst preserving maximum equivalency. Populations of kinetic models, preserving similarity in kinetic parameters, were built around the reference steady-states and their metabolic sensitivity coefficients (MSCs) were computed. The MSCs were sensitive to the model complexity. We proposed a metric for measuring the sensitivity of MSCs to these structural changes. CONCLUSIONS: We proposed for the first time a workflow for scaling up the size of kinetic models while preserving equivalency between the kinetic models. Using this workflow, we demonstrate that model complexity in terms of networks size has significant impact on sensitivity characteristics of kinetic models. Therefore, it is essential to account for the effects of network complexity when constructing kinetic models. The presented metric for measuring MSC sensitivity to structural changes can guide modelers and experimentalists in improving model quality and guide synthetic biology and metabolic engineering. Our proposed workflow enables the testing of the suitability of a kinetic model for answering certain study-specific questions. We argue that the model-based metabolic design targets that are common across models of different size are of higher confidence, while those that are different could be the objective of investigations for model improvement.

Publisher Correction: MEMOTE for standardized genome-scale metabolic model testing.

An amendment to this paper has been published and can be accessed via a link at the top of the paper.

Kinetic models of metabolism that consider alternative steady-state solutions of intracellular fluxes and concentrations.

Large-scale kinetic models are used for designing, predicting, and understanding the metabolic responses of living cells. Kinetic models are particularly attractive for the biosynthesis of target molecules in cells as they are typically better than other types of models at capturing the complex cellular biochemistry. Using simpler stoichiometric models as scaffolds, kinetic models are built around a steady-state flux profile and a metabolite concentration vector that are typically determined via optimization. However, as the underlying optimization problem is underdetermined, even after incorporating available experimental omics data, one cannot uniquely determine the operational configuration in terms of metabolic fluxes and metabolite concentrations. As a result, some reactions can operate in either the forward or reverse direction while still agreeing with the observed physiology. Here, we analyze how the underlying uncertainty in intracellular fluxes and concentrations affects predictions of constructed kinetic models and their design in metabolic engineering and systems biology studies. To this end, we integrated the omics data of optimally grown Escherichia coli into a stoichiometric model and constructed populations of non-linear large-scale kinetic models of alternative steady-state solutions consistent with the physiology of the E. coli aerobic metabolism. We performed metabolic control analysis (MCA) on these models, highlighting that MCA-based metabolic engineering decisions are strongly affected by the selected steady state and appear to be more sensitive to concentration values rather than flux values. To incorporate this into future studies, we propose a workflow for moving towards more reliable and robust predictions that are consistent with all alternative steady-state solutions. This workflow can be applied to all kinetic models to improve the consistency and accuracy of their predictions. Additionally, we show that, irrespective of the alternative steady-state solution, increased activity of phosphofructokinase and decreased ATP maintenance requirements would improve cellular growth of optimally grown E. coli.

pyTFA and matTFA: a Python package and a Matlab toolbox for Thermodynamics-based Flux Analysis.

Summary: pyTFA and matTFA are the first published implementations of the original TFA paper. Specifically, they include explicit formulation of Gibbs energies and metabolite concentrations, which enables straightforward integration of metabolite concentration measurements. Motivation: High-throughput analytic technologies provide a wealth of omics data that can be used to perform thorough analyses for a multitude of studies in the areas of Systems Biology and Biotechnology. Nevertheless, most studies are still limited to constraint-based Flux Balance Analyses (FBA), neglecting an important physicochemical constraint: thermodynamics. Thermodynamics-based Flux Analysis (TFA) in metabolic models enables the integration of quantitative metabolomics data to study their effects on the net-flux directionality of reactions in the network. In addition, it allows us to estimate how far each reaction operates from thermodynamic equilibrium, which provides critical information for guiding metabolic engineering decisions. Results: We present a Python package (pyTFA) and a Matlab toolbox (matTFA) that implement TFA. We show an example of application on both a reduced and a genome-scale model of E. coli., and demonstrate TFA and data integration through TFA reduce the feasible flux space with respect to FBA. Availability and implementation: Documented implementation of TFA framework both in Python (pyTFA) and Matlab (matTFA) are available on www.github.com/EPFL-LCSB/. Supplementary information: Supplementary data are available at Bioinformatics online.

redGEM: Systematic reduction and analysis of genome-scale metabolic reconstructions for development of consistent core metabolic models.

Genome-scale metabolic reconstructions have proven to be valuable resources in enhancing our understanding of metabolic networks as they encapsulate all known metabolic capabilities of the organisms from genes to proteins to their functions. However the complexity of these large metabolic networks often hinders their utility in various practical applications. Although reduced models are commonly used for modeling and in integrating experimental data, they are often inconsistent across different studies and laboratories due to different criteria and detail, which can compromise transferability of the findings and also integration of experimental data from different groups. In this study, we have developed a systematic semi-automatic approach to reduce genome-scale models into core models in a consistent and logical manner focusing on the central metabolism or subsystems of interest. The method minimizes the loss of information using an approach that combines graph-based search and optimization methods. The resulting core models are shown to be able to capture key properties of the genome-scale models and preserve consistency in terms of biomass and by-product yields, flux and concentration variability and gene essentiality. The development of these "consistently-reduced" models will help to clarify and facilitate integration of different experimental data to draw new understanding that can be directly extendable to genome-scale models.

Rites of passage: requirements and standards for building kinetic models of metabolic phenotypes.

The overarching ambition of kinetic metabolic modeling is to capture the dynamic behavior of metabolism to such an extent that systems and synthetic biology strategies can reliably be tested in silico. The lack of kinetic data hampers the development of kinetic models, and most of the current models use ad hoc reduced stoichiometry or oversimplified kinetic rate expressions, which may limit their predictive strength. There is a need to introduce the community-level standards that will organize and accelerate the future developments in this area. We introduce here a set of requirements that will ensure the model quality, we examine the current kinetic models with respect to these requirements, and we propose a general workflow for constructing models that satisfy these requirements.

Spatial proteomic and phospho-proteomic organization in three prototypical cell migration modes.

BACKGROUND: Tight spatio-temporal signaling of cytoskeletal and adhesion dynamics is required for localized membrane protrusion that drives directed cell migration. Different ensembles of proteins are therefore likely to get recruited and phosphorylated in membrane protrusions in response to specific cues. RESULTS: HERE, WE USE AN ASSAY THAT ALLOWS TO BIOCHEMICALLY PURIFY EXTENDING PROTRUSIONS OF CELLS MIGRATING IN RESPONSE TO THREE PROTOTYPICAL RECEPTORS: integrins, recepor tyrosine kinases and G-coupled protein receptors. Using quantitative proteomics and phospho-proteomics approaches, we provide evidence for the existence of cue-specific, spatially distinct protein networks in the different cell migration modes. CONCLUSIONS: The integrated analysis of the large-scale experimental data with protein information from databases allows us to understand some emergent properties of spatial regulation of signaling during cell migration. This provides the cell migration community with a large-scale view of the distribution of proteins and phospho-proteins regulating directed cell migration.

Prediction stability in a data-based, mechanistic model of σF regulation during sporulation in Bacillus subtilis.

Mathematical modeling of biological networks can help to integrate a large body of information into a consistent framework, which can then be used to arrive at novel mechanistic insight and predictions. We have previously developed a detailed, mechanistic model for the regulation of σ(F) during sporulation in Bacillus subtilis. The model was based on a wide range of quantitative data, and once fitted to the data, the model made predictions that could be confirmed in experiments. However, the analysis was based on a single optimal parameter set. We wondered whether the predictions of the model would be stable for all optimal parameter sets. To that end we conducted a global parameter screen within the physiological parameter ranges. The screening approach allowed us to identify sensitive and sloppy parameters, and highlighted further required datasets during the optimization. Eventually, all parameter sets that reproduced all available data predicted the physiological situation correctly.

Predictive models for cellular signaling networks.

This chapter provides an introduction to the formulation and analysis of differential-equation-based models for biological regulatory networks. In the first part, we discuss basic reaction types and the use of mass action kinetics and of simplifying approximations in the development of models for biological signaling. In the second part we introduce phase plane and linear stability analysis to evaluate the time evolution and identify the long-term attractors of dynamic systems. We then discuss the use of bifurcation diagrams to evaluate the parameter dependency of qualitative network behaviors (i.e., the emergence of oscillations or switches), and we give measures for the sensitivity and robustness of the signaling output.

Analyzing and constraining signaling networks: parameter estimation for the user.

The behavior of most dynamical models not only depends on the wiring but also on the kind and strength of interactions which are reflected in the parameter values of the model. The predictive value of mathematical models therefore critically hinges on the quality of the parameter estimates. Constraining a dynamical model by an appropriate parameterization follows a 3-step process. In an initial step, it is important to evaluate the sensitivity of the parameters of the model with respect to the model output of interest. This analysis points at the identifiability of model parameters and can guide the design of experiments. In the second step, the actual fitting needs to be carried out. This step requires special care as, on the one hand, noisy as well as partial observations can corrupt the identification of system parameters. On the other hand, the solution of the dynamical system usually depends in a highly nonlinear fashion on its parameters and, as a consequence, parameter estimation procedures get easily trapped in local optima. Therefore any useful parameter estimation procedure has to be robust and efficient with respect to both challenges. In the final step, it is important to access the validity of the optimized model. A number of reviews have been published on the subject. A good, nontechnical overview is provided by Jaqaman and Danuser (Nat Rev Mol Cell Biol 7(11):813-819, 2006) and a classical introduction, focussing on the algorithmic side, is given in Press (Numerical recipes: The art of scientific computing, Cambridge University Press, 3rd edn., 2007, Chapters 10 and 15). We will focus on the practical issues related to parameter estimation and use a model of the TGFß-signaling pathway as an educative example. Corresponding parameter estimation software and models based on MATLAB code can be downloaded from the authors's web page ( http://www.bsse.ethz.ch/cobi ).

Plasticity of TGF-ß signaling.

BACKGROUND: The family of TGF-ß ligands is large and its members are involved in many different signaling processes. These signaling processes strongly differ in type with TGF-ß ligands eliciting both sustained or transient responses. Members of the TGF-ß family can also act as morphogen and cellular responses would then be expected to provide a direct read-out of the extracellular ligand concentration. A number of different models have been proposed to reconcile these different behaviours. We were interested to define the set of minimal modifications that are required to change the type of signal processing in the TGF-ß signaling network. RESULTS: To define the key aspects for signaling plasticity we focused on the core of the TGF-ß signaling network. With the help of a parameter screen we identified ranges of kinetic parameters and protein concentrations that give rise to transient, sustained, or oscillatory responses to constant stimuli, as well as those parameter ranges that enable a proportional response to time-varying ligand concentrations (as expected in the read-out of morphogens). A combination of a strong negative feedback and fast shuttling to the nucleus biases signaling to a transient rather than a sustained response, while oscillations were obtained if ligand binding to the receptor is weak and the turn-over of the I-Smad is fast. A proportional read-out required inefficient receptor activation in addition to a low affinity of receptor-ligand binding. We find that targeted modification of single parameters suffices to alter the response type. The intensity of a constant signal (i.e. the ligand concentration), on the other hand, affected only the strength but not the type of the response. CONCLUSIONS: The architecture of the TGF-ß pathway enables the observed signaling plasticity. The observed range of signaling outputs to TGF-ß ligand in different cell types and under different conditions can be explained with differences in cellular protein concentrations and with changes in effective rate constants due to cross-talk with other signaling pathways. It will be interesting to uncover the exact cellular differences as well as the details of the cross-talks in future work.

A computational analysis of the dynamic roles of talin, Dok1, and PIPKI for integrin activation.

Integrin signaling regulates cell migration and plays a pivotal role in developmental processes and cancer metastasis. Integrin signaling has been studied extensively and much data is available on pathway components and interactions. Yet the data is fragmented and an integrated model is missing. We use a rule-based modeling approach to integrate available data and test biological hypotheses regarding the role of talin, Dok1 and PIPKI in integrin activation. The detailed biochemical characterization of integrin signaling provides us with measured values for most of the kinetics parameters. However, measurements are not fully accurate and the cellular concentrations of signaling proteins are largely unknown and expected to vary substantially across different cellular conditions. By sampling model behaviors over the physiologically realistic parameter range we find that the model exhibits only two different qualitative behaviors and these depend mainly on the relative protein concentrations, which offers a powerful point of control to the cell. Our study highlights the necessity to characterize model behavior not for a single parameter optimum, but to identify parameter sets that characterize different signaling modes.