Dynamic Models Of Biological Systems Research Articles

PyRates-A code-generation tool for modeling dynamical systems in biology and beyond.

The mathematical study of real-world dynamical systems relies on models composed of differential equations. Numerical methods for solving and analyzing differential equation systems are essential when complex biological problems have to be studied, such as the spreading of a virus, the evolution of competing species in an ecosystem, or the dynamics of neurons in the brain. Here we present PyRates, a Python-based software for modeling and analyzing differential equation systems via numerical methods. PyRates is specifically designed to account for the inherent complexity of biological systems. It provides a new language for defining models that mirrors the modular organization of real-world dynamical systems and thus simplifies the implementation of complex networks of interacting dynamic entities. Furthermore, PyRates provides extensive support for the various forms of interaction delays that can be observed in biological systems. The core of PyRates is a versatile code-generation system that translates user-defined models into "backend" implementations in various languages, including Python, Fortran, Matlab, and Julia. This allows users to apply a wide range of analysis methods for dynamical systems, eliminating the need for manual translation between code bases. PyRates may also be used as a model definition interface for the creation of custom dynamical systems tools. To demonstrate this, we developed two extensions of PyRates for common analyses of dynamic models of biological systems: PyCoBi for bifurcation analysis and RectiPy for parameter fitting. We demonstrate in a series of example models how PyRates can be used in combination with PyCoBi and RectiPy for model analysis and fitting. Together, these tools offer a versatile framework for applying computational modeling and numerical analysis methods to dynamical systems in biology and beyond.

Symmetries in Dynamic Models of Biological Systems: Mathematical Foundations and Implications

Symmetries are ubiquitous in nature. Almost all organisms have some kind of “symmetry”, meaning that their shape does not change under some geometric transformation. This geometrical concept of symmetry is intuitive and easy to recognize. On the other hand, the behavior of many biological systems over time can be described with ordinary differential equations. These dynamic models may also possess “symmetries”, meaning that the time courses of some variables remain invariant under certain transformations. Unlike the previously mentioned symmetries, the ones present in dynamic models are not geometric, but infinitesimal transformations. These mathematical symmetries can be related to certain features of the system’s dynamic behavior, such as robustness or adaptation capabilities. However, they can also arise from questionable modeling choices, which may lead to non-identifiability and non-observability. This paper provides an overview of the types of symmetries that appear in dynamic models, the mathematical tools available for their analyses, the ways in which they are related to system properties, and the implications for biological modeling.

AutoRepar: A method to obtain identifiable and observable reparameterizations of dynamic models with mechanistic insights

AbstractMechanistic dynamic models of biological systems allow for a quantitative and systematic interpretation of data and the generation of testable hypotheses. However, these models are often over‐parameterized, leading to nonidentifiability and nonobservability, that is, the impossibility of inferring their parameters and state variables. The lack of structural identifiability and observability (SIO) compromises a model's ability to make predictions and provide insight. Here we present a methodology, AutoRepar, that corrects SIO deficiencies of nonlinear ODE models automatically, yielding reparameterized models that are structurally identifiable and observable. The reparameterization preserves the mechanistic meaning of selected variables, and has the exact same dynamics and input‐output mapping as the original model. We implement AutoRepar as an extension of the STRIKE‐GOLDD software toolbox for SIO analysis, applying it to several models from the literature to demonstrate its ability to repair their structural deficiencies. AutoRepar increases the applicability of mechanistic models, enabling them to provide reliable information about their parameters and dynamics.

Scalable nonlinear programming framework for parameter estimation in dynamic biological system models.

We present a nonlinear programming (NLP) framework for the scalable solution of parameter estimation problems that arise in dynamic modeling of biological systems. Such problems are computationally challenging because they often involve highly nonlinear and stiff differential equations as well as many experimental data sets and parameters. The proposed framework uses cutting-edge modeling and solution tools which are computationally efficient, robust, and easy-to-use. Specifically, our framework uses a time discretization approach that: i) avoids repetitive simulations of the dynamic model, ii) enables fully algebraic model implementations and computation of derivatives, and iii) enables the use of computationally efficient nonlinear interior point solvers that exploit sparse and structured linear algebra techniques. We demonstrate these capabilities by solving estimation problems for synthetic human gut microbiome community models. We show that an instance with 156 parameters, 144 differential equations, and 1,704 experimental data points can be solved in less than 3 minutes using our proposed framework (while an off-the-shelf simulation-based solution framework requires over 7 hours). We also create large instances to show that the proposed framework is scalable and can solve problems with up to 2,352 parameters, 2,304 differential equations, and 20,352 data points in less than 15 minutes. The proposed framework is flexible and easy-to-use, can be broadly applied to dynamic models of biological systems, and enables the implementation of sophisticated estimation techniques to quantify parameter uncertainty, to diagnose observability/uniqueness issues, to perform model selection, and to handle outliers.

General method to find the attractors of discrete dynamic models of biological systems.

Analyzing the long-term behaviors (attractors) of dynamic models of biological networks can provide valuable insight. We propose a general method that can find the attractors of multilevel discrete dynamical systems by extending a method that finds the attractors of a Boolean network model. The previous method is based on finding stable motifs, subgraphs whose nodes' states can stabilize on their own. We extend the framework from binary states to any finite discrete levels by creating a virtual node for each level of a multilevel node, and describing each virtual node with a quasi-Boolean function. We then create an expanded representation of the multilevel network, find multilevel stable motifs and oscillating motifs, and identify attractors by successive network reduction. In this way, we find both fixed point attractors and complex attractors. We implemented an algorithm, which we test and validate on representative synthetic networks and on published multilevel models of biological networks. Despite its primary motivation to analyze biological networks, our motif-based method is general and can be applied to any finite discrete dynamical system.

Using a system of differential equations that models cattle growth to uncover the genetic basis of complex traits.

The interplay between dynamic models of biological systems and genomics is based on the assumption that genetic variation of the complex trait (i.e., outcome of model behavior) arises from component traits (i.e., model parameters) in lower hierarchical levels. In order to provide a proof of concept of this statement for a cattle growth model, we ask whether model parameters map genomic regions that harbor quantitative trait loci (QTLs) already described for the complex trait. We conducted a genome-wide association study (GWAS) with a Bayesian hierarchical LASSO method in two parameters of the Davis Growth Model, a system of three ordinary differential equations describing DNA accretion, protein synthesis and degradation, and fat synthesis. Phenotypic and genotypic data were available for 893 Nellore (Bos indicus) cattle. Computed values for parameter k1 (DNA accretion rate) ranged from 0.005±0.003 and for α (constant for energy for maintenance requirement) 0.134±0.024. The expected biological interpretation of the parameters is confirmed by QTLs mapped for k1 and α. QTLs within genomic regions mapped for k1 are expected to be correlated with the DNA pool: body size and weight. Single nucleotide polymorphisms (SNPs) which were significant for α mapped QTLs that had already been associated with residual feed intake, feed conversion ratio, average daily gain (ADG), body weight, and also dry matter intake. SNPs identified for k1 were able to additionally explain 2.2% of the phenotypic variability of the complex ADG, even when SNPs for k1 did not match the genomic regions associated with ADG. Although improvements are needed, our findings suggest that genomic analysis on component traits may help to uncover the genetic basis of more complex traits, particularly when lower biological hierarchies are mechanistically described by mathematical simulation models.

DeBInfer: Bayesian inference for dynamical models of biological systems in R

Summary Understanding the mechanisms underlying biological systems, and ultimately, predicting their behaviours in a changing environment, requires overcoming the gap between mathematical models and experimental or observational data. Differential equations (DEs) are commonly used to model the temporal evolution of biological systems, but statistical methods for comparing DE models to data and for parameter inference are relatively poorly developed. This is especially problematic in the context of biological systems where observations are often noisy and only a small number of time points may be available. The Bayesian approach offers a coherent framework for parameter inference that can account for multiple sources of uncertainty, while making use of prior information. It offers a rigorous methodology for parameter inference, as well as modelling the link between unobservable model states and parameters, and observable quantities. We present deBInfer, a package for the statistical computing environment R, implementing a Bayesian framework for parameter inference in DEs. deBInfer provides templates for the DE model, the observation model and data likelihood, and the model parameters and their prior distributions. A Markov chain Monte Carlo (MCMC) procedure processes these inputs to estimate the posterior distributions of the parameters and any derived quantities, including the model trajectories. Further functionality is provided to facilitate MCMC diagnostics, the visualization of the posterior distributions of model parameters and trajectories, and the use of compiled DE models for improved computational performance. The templating approach makes deBInfer applicable to a wide range of DE models. We demonstrate its application to ordinary and delay DE models for population ecology.

Analysis of a dynamic model of guard cell signaling reveals the stability of signal propagation.

BackgroundAnalyzing the long-term behaviors (attractors) of dynamic models of biological systems can provide valuable insight into biological phenotypes and their stability. In this paper we identify the allowed long-term behaviors of a multi-level, 70-node dynamic model of the stomatal opening process in plants.ResultsWe start by reducing the model’s huge state space. We first reduce unregulated nodes and simple mediator nodes, then simplify the regulatory functions of selected nodes while keeping the model consistent with experimental observations. We perform attractor analysis on the resulting 32-node reduced model by two methods: 1. converting it into a Boolean model, then applying two attractor-finding algorithms; 2. theoretical analysis of the regulatory functions. We further demonstrate the robustness of signal propagation by showing that a large percentage of single-node knockouts does not affect the stomatal opening level.ConclusionsCombining both methods with analysis of perturbation scenarios, we conclude that all nodes except two in the reduced model have a single attractor; and only two nodes can admit oscillations. The multistability or oscillations of these four nodes do not affect the stomatal opening level in any situation. This conclusion applies to the original model as well in all the biologically meaningful cases. In addition, the stomatal opening level is resilient against single-node knockouts. Thus, we conclude that the complex structure of this signal transduction network provides multiple information propagation pathways while not allowing extensive multistability or oscillations, resulting in robust signal propagation. Our innovative combination of methods offers a promising way to analyze multi-level models.Electronic supplementary materialThe online version of this article (doi:10.1186/s12918-016-0327-7) contains supplementary material, which is available to authorized users.

Enhanced parallel Differential Evolution algorithm for problems in computational systems biology

Many key problems in computational systems biology and bioinformatics can be formulated and solved using a global optimization framework. The complexity of the underlying mathematical models require the use of efficient solvers in order to obtain satisfactory results in reasonable computation times. Metaheuristics are gaining recognition in this context, with Differential Evolution (DE) as one of the most popular methods. However, for most realistic applications, like those considering parameter estimation in dynamic models, DE still requires excessive computation times.Here we consider this latter class of problems and present several enhancements to DE based on the introduction of additional algorithmic steps and the exploitation of parallelism. In particular, we propose an asynchronous parallel implementation of DE which has been extended with improved heuristics to exploit the specific structure of parameter estimation problems in computational systems biology. The proposed method is evaluated with different types of benchmarks problems: (i) black-box global optimization problems and (ii) calibration of non-linear dynamic models of biological systems, obtaining excellent results both in terms of quality of the solution and regarding speedup and scalability.

Hierarchical multivariate regression-based sensitivity analysis reveals complex parameter interaction patterns in dynamic models

Dynamic models of biological systems often possess complex and multivariate mappings between input parameters and output state variables, posing challenges for comprehensive sensitivity analysis across the biologically relevant parameter space. In particular, more efficient and robust ways to obtain a solid understanding of how the sensitivity to each parameter depends on the values of the other parameters are sorely needed.We report a new methodology for global sensitivity analysis based on Hierarchical Cluster-based Partial Least Squares Regression (HC-PLSR)-based approximations (metamodelling) of the input–output mappings of dynamic models, which we expect to be generic, efficient and robust, even for systems with highly nonlinear input–output relationships. The two-step HC-PLSR metamodelling automatically separates the observations (here corresponding to different combinations of input parameter values) into groups based on the dynamic model behaviour, then analyses each group separately with Partial Least Squares Regression (PLSR). This produces one global regression model comprising all observations, as well as regional regression models within each group, where the regression coefficients can be used as sensitivity measures. Thereby a more accurate description of complex interactions between inputs to the dynamic model can be revealed through analysis of how a certain level of one input parameter affects the model sensitivity to other inputs. We illustrate the usefulness of the HC-PLSR approach on a dynamic model of a mouse heart muscle cell, and demonstrate how it reveals interaction patterns of probable biological significance not easily identifiable by a global regression-based sensitivity analysis alone.Applied for sensitivity analysis of a complex, high-dimensional dynamic model of the mouse heart muscle cell, several interactions between input parameters were identified by the two-step HC-PLSR analysis that could not be detected in the single-step global analysis. Hence, our approach has the potential to reveal new biological insight through the identification of complex parameter interaction patterns. The HC-PLSR metamodel complexity can be adjusted according to the nonlinear complexity of the input–output mapping of the analysed dynamic model through adjustment of the number of regional regression models included. This facilitates sensitivity analysis of dynamic models of varying complexities.

Calibration of dynamic models of biological systems with KInfer

Methods for parameter estimation that are robust to experimental uncertainties and to stochastic and biological noise and that require a minimum of a priori input knowledge are of key importance in computational systems biology. The new method presented in this paper aims to ensure an inference model that deduces the rate constants of a system of biochemical reactions from experimentally measured time courses of reactants. This new method was applied to some challenging parameter estimation problems of nonlinear dynamic biological systems and was tested both on synthetic and real data. The synthetic case studies are the 12-state model of the SERCA pump and a model of a genetic network containing feedback loops of interaction between regulator and effector genes. The real case studies consist of a model of the reaction between the inhibitor kappaB kinase enzyme and its substrate in the signal transduction pathway of NF-kappaB, and a stiff model of the fermentation pathway of Lactococcus lactis.

Hierarchical analysis of piecewise affine models of gene regulatory networks

A key point in the analysis of dynamical models of biological systems is to handle systems of relatively high dimensions. In the present paper we propose a method to hierarchically organize a certain type of piecewise affine (PWA) differential systems. This specific class of systems has been extensively studied for the past few years, as it provides a good framework to model gene regulatory networks. The method, shown on several examples, allows a qualitative analysis of the asymptotic behavior of a PWA system, decomposing it into several smaller subsystems. This technique, based on the well-known strongly connected components decomposition, is not new. However, its adaptation to the non-smooth PWA differential equations turns out to be quite relevant because of the strong discrete structure underlying these equations. Its biological relevance is shown on a 7-dimensional PWA system modeling the gene network responsible for the carbon starvation response in Escherichia coli.

Novel metaheuristic for parameter estimation in nonlinear dynamic biological systems.

BackgroundWe consider the problem of parameter estimation (model calibration) in nonlinear dynamic models of biological systems. Due to the frequent ill-conditioning and multi-modality of many of these problems, traditional local methods usually fail (unless initialized with very good guesses of the parameter vector). In order to surmount these difficulties, global optimization (GO) methods have been suggested as robust alternatives. Currently, deterministic GO methods can not solve problems of realistic size within this class in reasonable computation times. In contrast, certain types of stochastic GO methods have shown promising results, although the computational cost remains large. Rodriguez-Fernandez and coworkers have presented hybrid stochastic-deterministic GO methods which could reduce computation time by one order of magnitude while guaranteeing robustness. Our goal here was to further reduce the computational effort without loosing robustness.ResultsWe have developed a new procedure based on the scatter search methodology for nonlinear optimization of dynamic models of arbitrary (or even unknown) structure (i.e. black-box models). In this contribution, we describe and apply this novel metaheuristic, inspired by recent developments in the field of operations research, to a set of complex identification problems and we make a critical comparison with respect to the previous (above mentioned) successful methods.ConclusionRobust and efficient methods for parameter estimation are of key importance in systems biology and related areas. The new metaheuristic presented in this paper aims to ensure the proper solution of these problems by adopting a global optimization approach, while keeping the computational effort under reasonable values. This new metaheuristic was applied to a set of three challenging parameter estimation problems of nonlinear dynamic biological systems, outperforming very significantly all the methods previously used for these benchmark problems.

A C++ implementation of an individual/landscape model

Dynamic models of biological systems are diverse, but have components in common. The description of a system may involve multiple databases which change and interact (individuals, populations, the landscape). A clock is required to advance the system description through time. In addition, the generation of random variates from different probability distributions, interpolation of values from tabular data, and the collection of statistics on system components may be required. All of these tasks can be effectively achieved through object oriented programming (OOP). The C++ programming language is the object oriented language utilized in the implementation described. Object oriented programming can be utilized to facilitate the organization and development of dynamic models of biological systems. A C++ hierarchy of classes describing objects in the system can start at the level of machine processing. A binary object can use individual bits of memory to represent binary data. Computer memory requirements can be minimized by using a linked list class to place objects in a list whose composition and size vary as the program executes. The C++ feature of runtime linking can connect a general base class with the data required for a specific application. These machine oriented classes can then be inherited into an individual animal class with binary data (male/female, diseased/nondiseased) and a population class which is a linked list of individuals. The binary arithmetic and memory pointers involved in these base classes can be completely transparent to the user of the individual and population classes which are described in biological rather than machine terms. OOP can facilitate the development of a description of the individual, population, landscape and their interactions. A landscape object can be described by a grid of rectangular cells linked on four sides to adjoining cells by memory pointers. These memory pointers to adjoining cells can be utilized to describe flows and movement of objects over the grid. This base class grid description can be inherited into a derived habitat class with additional data. The location or coordinates of objects located on the grid can be utilized to access the data for the cell in which it is located. Since grids describing landscapes can be large and heterogeneous, GIS can be utilized for grid initialization and for analysis of spatial data in the simulation output. A simulation clock is required to advance the description of the objects in the system and their interactions through time. A base clock class advances time in specified increments, or between events, or in specified increments until an event occurs. Events can be either scheduled long term or triggered by a change in the system and stored in an event queue by using the same base linked list class used to describe populations. The type of events which occur in a system are specified in a derived clock class which is accessed by the base clock class when the program executes through run-time linking. An example application is developed to describe winter kill of deer. Analysis of the simulation output in a GIS showed simple rules for individual animal behavior affected the pattern of landscape usage through time. Source code and documentation are available.