Review of stochastic hybrid systems with applications in biological systems modeling and analysis

Stochastic hybrid systems (SHS) have attracted a lot of research interests in recent years. In this paper, we review some of the recent applications of SHS to biological systems modeling and analysis. Due to the nature of molecular interactions, many biological processes can be conveniently described as a mixture of continuous and discrete phenomena employing SHS models. With the advancement of SHS theory, it is expected that insights can be obtained about biological processes such as drug effects on gene regulation. Furthermore, combining with advanced experimental methods, in silico simulations using SHS modeling techniques can be carried out for massive and rapid verification or falsification of biological hypotheses. The hope is to substitute costly and time-consuming in vitro or in vivo experiments or provide guidance for those experiments and generate better hypotheses.     

Odefy - From discrete to continuous models

Background  

Phenomenological information about regulatory interactions is frequently available and can be readily converted to Boolean models. Fully quantitative models, on the other hand, provide detailed insights into the precise dynamics of the underlying system. In order to connect discrete and continuous modeling approaches, methods for the conversion of Boolean systems into systems of ordinary differential equations have been developed recently. As biological interaction networks have steadily grown in size and complexity, a fully automated framework for the conversion process is desirable.     

Information-Theoretic Analysis of the Dynamics of an Executable Biological Model

To facilitate analysis and understanding of biological systems, large-scale data are often integrated into models using a variety of mathematical and computational approaches. Such models describe the dynamics of the biological system and can be used to study the changes in the state of the system over time. For many model classes, such as discrete or continuous dynamical systems, there exist appropriate frameworks and tools for analyzing system dynamics. However, the heterogeneous information that encodes and bridges molecular and cellular dynamics, inherent to fine-grained molecular simulation models, presents significant challenges to the study of system dynamics. In this paper, we present an algorithmic information theory based approach for the analysis and interpretation of the dynamics of such executable models of biological systems. We apply a normalized compression distance (NCD) analysis to the state representations of a model that simulates the immune decision making and immune cell behavior. We show that this analysis successfully captures the essential information in the dynamics of the system, which results from a variety of events including proliferation, differentiation, or perturbations such as gene knock-outs. We demonstrate that this approach can be used for the analysis of executable models, regardless of the modeling framework, and for making experimentally quantifiable predictions.     

A continuous optimization approach for inferring parameters in mathematical models of regulatory networks

Background

The advances of systems biology have raised a large number of sophisticated mathematical models for describing the dynamic property of complex biological systems. One of the major steps in developing mathematical models is to estimate unknown parameters of the model based on experimentally measured quantities. However, experimental conditions limit the amount of data that is available for mathematical modelling. The number of unknown parameters in mathematical models may be larger than the number of observation data. The imbalance between the number of experimental data and number of unknown parameters makes reverse-engineering problems particularly challenging.

Results

To address the issue of inadequate experimental data, we propose a continuous optimization approach for making reliable inference of model parameters. This approach first uses a spline interpolation to generate continuous functions of system dynamics as well as the first and second order derivatives of continuous functions. The expanded dataset is the basis to infer unknown model parameters using various continuous optimization criteria, including the error of simulation only, error of both simulation and the first derivative, or error of simulation as well as the first and second derivatives. We use three case studies to demonstrate the accuracy and reliability of the proposed new approach. Compared with the corresponding discrete criteria using experimental data at the measurement time points only, numerical results of the ERK kinase activation module show that the continuous absolute-error criteria using both function and high order derivatives generate estimates with better accuracy. This result is also supported by the second and third case studies for the G1/S transition network and the MAP kinase pathway, respectively. This suggests that the continuous absolute-error criteria lead to more accurate estimates than the corresponding discrete criteria. We also study the robustness property of these three models to examine the reliability of estimates. Simulation results show that the models with estimated parameters using continuous fitness functions have better robustness properties than those using the corresponding discrete fitness functions.

Conclusions

The inference studies and robustness analysis suggest that the proposed continuous optimization criteria are effective and robust for estimating unknown parameters in mathematical models.

Electronic supplementary material

The online version of this article (doi:10.1186/1471-2105-15-256) contains supplementary material, which is available to authorized users.     

Preservation of dynamic properties in qualitative modeling frameworks for gene regulatory networks

Mathematical modeling often helps to provide a systems perspective on gene regulatory networks. In particular, qualitative approaches are useful when detailed kinetic information is lacking. Multiple methods have been developed that implement qualitative information in different ways, e.g., in purely discrete or hybrid discrete/continuous models. In this paper, we compare the discrete asynchronous logical modeling formalism for gene regulatory networks due to R. Thomas with piecewise affine differential equation models. We provide a local characterization of the qualitative dynamics of a piecewise affine differential equation model using the discrete dynamics of a corresponding Thomas model. Based on this result, we investigate the consistency of higher-level dynamical properties such as attractor characteristics and reachability. We show that although the two approaches are based on equivalent information, the resulting qualitative dynamics are different. In particular, the dynamics of the piecewise affine differential equation model is not a simple refinement of the dynamics of the Thomas model     

Simulation-based model checking approach to cell fate specification during Caenorhabditis elegans vulval development by hybrid functional Petri net with extension

Background  

Model checking approaches were applied to biological pathway validations around 2003. Recently, Fisher et al. have proved the importance of model checking approach by inferring new regulation of signaling crosstalk in C. elegans and confirming the regulation with biological experiments. They took a discrete and state-based approach to explore all possible states of the system underlying vulval precursor cell (VPC) fate specification for desired properties. However, since both discrete and continuous features appear to be an indispensable part of biological processes, it is more appropriate to use quantitative models to capture the dynamics of biological systems. Our key motivation of this paper is to establish a quantitative methodology to model and analyze in silico models incorporating the use of model checking approach.     

Evolutionary dynamics of continuous strategy games on graphs and social networks under weak selection

Understanding the emergence of cooperation among selfish individuals has been a long-standing puzzle, which has been studied by a variety of game models. Most previous studies presumed that interactions between individuals are discrete, but it seems unrealistic in real systems. Recently, there are increasing interests in studying game models with a continuous strategy space. Existing research work on continuous strategy games mainly focuses on well-mixed populations. Especially, little theoretical work has been conducted on their evolutionary dynamics in a structured population. In the previous work (Zhong et al., BioSystems, 2012), we showed that under strong selection, continuous and discrete strategies have significantly different equilibrium and game dynamics in spatially structured populations. In this paper, we further study evolutionary dynamics of continuous strategy games under weak selection in structured populations. By using the fixation probability based stochastic dynamics, we derive exact conditions of natural selection favoring cooperation for the death–birth updating scheme. We also present a network gain decomposition of the game equilibrium, which might provide a new view of the network reciprocity in a quantitative way. Finally, we make a detailed comparison between games using discrete and continuous strategies. As compared to the former, we find that for the latter (i) the same selection conditions are derived for the general 2 × 2 game; especially, the rule b/c > k in a simplified Prisoner's Dilemma is valid as well; however, (ii) for a coordination game, interestingly, the risk-dominant strategy is disfavored. Numerical simulations have also been conducted to validate our results.     

Discrete and continuous state population models in a noisy world

Simple ecological models operate mostly with population densities using continuous variables. However, in reality densities could not change continuously, since the population itself consists of integer numbers of individuals. At first sight this discrepancy appears to be irrelevant, nevertheless, it can cause large deviations between the actual statistical behaviour of biological populations and that predicted by the corresponding models. We investigate the conditions under which simple models, operating with continuous numbers of individuals can be used to approximate the dynamics of populations consisting of integer numbers of individuals. Based on our definition for the (statistical) distance between the two models we show that the continuous approach is acceptable as long as sufficiently high biological noise is present, or, the dynamical behaviour is regular (non-chaotic). The concepts are illustrated with the Ricker model and tested on the Tribolium castaneum data series. Further, we demonstrate with the help of T. castaneum's model that if time series are not much larger than the possible population states (as in this practical case) the noisy discrete and continuous models can behave temporarily differently, almost independently of the noise level. In this case the noisy, discrete model is more accurate [OR has to be applied].     

Reconstructing Causal Biological Networks through Active Learning

Reverse-engineering of biological networks is a central problem in systems biology. The use of intervention data, such as gene knockouts or knockdowns, is typically used for teasing apart causal relationships among genes. Under time or resource constraints, one needs to carefully choose which intervention experiments to carry out. Previous approaches for selecting most informative interventions have largely been focused on discrete Bayesian networks. However, continuous Bayesian networks are of great practical interest, especially in the study of complex biological systems and their quantitative properties. In this work, we present an efficient, information-theoretic active learning algorithm for Gaussian Bayesian networks (GBNs), which serve as important models for gene regulatory networks. In addition to providing linear-algebraic insights unique to GBNs, leading to significant runtime improvements, we demonstrate the effectiveness of our method on data simulated with GBNs and the DREAM4 network inference challenge data sets. Our method generally leads to faster recovery of underlying network structure and faster convergence to final distribution of confidence scores over candidate graph structures using the full data, in comparison to random selection of intervention experiments.     

Application of the stretched exponential function to fluorescence lifetime imaging

Conventional analyses of fluorescence lifetime measurements resolve the fluorescence decay profile in terms of discrete exponential components with distinct lifetimes. In complex, heterogeneous biological samples such as tissue, multi-exponential decay functions can appear to provide a better fit to fluorescence decay data than the assumption of a mono-exponential decay, but the assumption of multiple discrete components is essentially arbitrary and is often erroneous. Moreover, interactions, both between fluorophores and with their environment, can result in complex fluorescence decay profiles that represent a continuous distribution of lifetimes. Such continuous distributions have been reported for tryptophan, which is one of the main fluorophores in tissue. This situation is better represented by the stretched-exponential function (StrEF). In this work, we have applied, for the first time to our knowledge, the StrEF to time-domain whole-field fluorescence lifetime imaging (FLIM), yielding both excellent tissue contrast and goodness of fit using data from rat tissue. We note that for many biological samples for which there is no a priori knowledge of multiple discrete exponential fluorescence decay profiles, the StrEF is likely to provide a truer representation of the underlying fluorescence dynamics. Furthermore, fitting to a StrEF significantly decreases the required processing time, compared with a multi-exponential component fit and typically provides improved contrast and signal/noise in the resulting FLIM images. In addition, the stretched-exponential decay model can provide a direct measure of the heterogeneity of the sample, and the resulting heterogeneity map can reveal subtle tissue differences that other models fail to show.     

Qualitative theory of compartmental systems with lags

Dynamic models of many processes in the biological and physical sciences give systems of ordinary differential equations called compartmental systems. Often, these systems include time lags; in this context, continuous probability density functions (pdfs) of lags are far more important than discrete lags. There is a relatively complete theory of compartmental systems without lags, both linear and non-linear [SIAM Rev. 35 (1993) 43]. The authors extend their previous work on compartmental systems without lags to show that, for discrete lags and for a very large class of pdfs of continuous lags, compartmental systems with lags are equivalent to larger compartmental systems without lags. Consequently, the properties of compartmental systems with lags are the same as those of compartmental systems without lags. For a very large class of compartmental systems with time lags, one can show that the time lags themselves can be generated by compartmental systems without lags. Thus, such systems can be partitioned into a main system, which is the original system without the lags, plus compartmental subsystems without lags that generate the lags. The latter may be linear or non-linear and may be inserted into main systems that are linear or non-linear. The state variables of the compartmental lag subsystems are hidden variables in the formulation with explicit lags.     

A Multiple Binary Markov Chain Model for Explaining the Behaviour of Diverse Biological Phenomena

The behaviour of many biological systems can be attributed to that of a large number of units, with each unit swinging between two competing states. During the past few years efforts have been made (e.g., Chung and Kennedy , 1996) to describe such discrete systems using a multiple binary Markov chain model. Here we explore the gamut of these models and classify their behaviour into five qualitatively distinct types, corresponding to subregions of the parameter space. It is suggested that these model behaviours may correspond to behaviours observed in nature. A simple method for fitting the model to data is presented.     

Dynamic partitioning for hybrid simulation of the bistable HIV-1 transactivation network

MOTIVATION: The stochastic kinetics of a well-mixed chemical system, governed by the chemical Master equation, can be simulated using the exact methods of Gillespie. However, these methods do not scale well as systems become more complex and larger models are built to include reactions with widely varying rates, since the computational burden of simulation increases with the number of reaction events. Continuous models may provide an approximate solution and are computationally less costly, but they fail to capture the stochastic behavior of small populations of macromolecules. RESULTS: In this article we present a hybrid simulation algorithm that dynamically partitions the system into subsets of continuous and discrete reactions, approximates the continuous reactions deterministically as a system of ordinary differential equations (ODE) and uses a Monte Carlo method for generating discrete reaction events according to a time-dependent propensity. Our approach to partitioning is improved such that we dynamically partition the system of reactions, based on a threshold relative to the distribution of propensities in the discrete subset. We have implemented the hybrid algorithm in an extensible framework, utilizing two rigorous ODE solvers to approximate the continuous reactions, and use an example model to illustrate the accuracy and potential speedup of the algorithm when compared with exact stochastic simulation. AVAILABILITY: Software and benchmark models used for this publication can be made available upon request from the authors.     

Temporal constraints of a gene regulatory network: Refining a qualitative simulation

The modelling of gene regulatory networks (GRNs) has classically been addressed through very different approaches. Among others, extensions of Thomas’s asynchronous Boolean approach have been proposed, to better fit the dynamics of biological systems: genes may reach different discrete expression levels, depending on the states of other genes, called the regulators: thus, activations and inhibitions are triggered conditionally on the proper expression levels of these regulators. In contrast, some fine-grained propositions have focused on the molecular level as modelling the evolution of biological compound concentrations through differential equation systems. Both approaches are limited. The first one leads to an oversimplification of the system, whereas the second is incapable to tackle large GRNs. In this context, hybrid paradigms, that mix discrete and continuous features underlying distinct biological properties, achieve significant advances for investigating biological properties. One of these hybrid formalisms proposes to focus, within a GRN abstraction, on the time delay to pass from a gene expression level to the next. Until now, no research work has been carried out, which attempts to benefit from the modelling of a GRN by differential equations, converting it into a multi-valued logical formalism of Thomas, with the aim of performing biological applications.     

A Computer Model of Engineered Cardiac Monolayers

Engineered monolayers created using microabrasion and micropatterning methods have provided a simplified in vitro system to study the effects of anisotropy and fiber direction on electrical propagation. Interpreting the behavior in these culture systems has often been performed using classical computer models with continuous properties. However, such models do not account for the effects of random cell shapes, cell orientations, and cleft spaces inherent in these monolayers on the resulting wavefront conduction. This work presents a novel methodology for modeling a monolayer of cardiac tissue in which the factors governing cell shape, cell-to-cell coupling, and degree of cleft space are not constant but rather are treated as spatially random with assigned distributions. This modeling approach makes it possible to simulate wavefront propagation in a manner analogous to performing experiments on engineered monolayer tissues. Simulated results are compared to previously published measured data from monolayers used to investigate the role of cellular architecture on conduction velocities and anisotropy ratios. We also present an estimate for obtaining the electrical properties from these networks and demonstrate how variations in the discrete cellular architecture affect the macroscopic conductivities. The simulations support the common assumption that under normal ranges of coupling strength, tissues with relatively uniform distributions of cell shapes and connectivity can be represented using continuous models with conductivities derived from random discrete cellular architecture using either global or local estimates. The results also reveal that in the presence of abrupt changes in cell orientation, local estimates of tissue properties predict smoother changes in conductivity that may not adequately predict the discrete nature of propagation at the transition sites.     

A model for one-dimensional morphoelasticity and its application to fibroblast-populated collagen lattices

The mechanical behaviour of solid biological tissues has long been described using models based on classical continuum mechanics. However, the classical continuum theories of elasticity and viscoelasticity cannot easily capture the continual remodelling and associated structural changes in biological tissues. Furthermore, models drawn from plasticity theory are difficult to apply and interpret in this context, where there is no equivalent of a yield stress or flow rule. In this work, we describe a novel one-dimensional mathematical model of tissue remodelling based on the multiplicative decomposition of the deformation gradient. We express the mechanical effects of remodelling as an evolution equation for the effective strain, a measure of the difference between the current state and a hypothetical mechanically relaxed state of the tissue. This morphoelastic model combines the simplicity and interpretability of classical viscoelastic models with the versatility of plasticity theory. A novel feature of our model is that while most models describe growth as a continuous quantity, here we begin with discrete cells and develop a continuum representation of lattice remodelling based on an appropriate limit of the behaviour of discrete cells. To demonstrate the utility of our approach, we use this framework to capture qualitative aspects of the continual remodelling observed in fibroblast-populated collagen lattices, in particular its contraction and its subsequent sudden re-expansion when remodelling is interrupted.     

Hybrid simulation of cellular behavior

MOTIVATION: To be valuable to biological or biomedical research, in silico methods must be scaled to complex pathways and large numbers of interacting molecular species. The correct method for performing such simulations, discrete event simulation by Monte Carlo generation, is computationally costly for large complex systems. Approximation of molecular behavior by continuous models fails to capture stochastic behavior that is essential to many biological phenomena. RESULTS: We present a novel approach to building hybrid simulations in which some processes are simulated discretely, while other processes are handled in a continuous simulation by differential equations. This approach preserves the stochastic behavior of cellular pathways, yet enables scaling to large populations of molecules. We present an algorithm for synchronizing data in a hybrid simulation and discuss the trade-offs in such simulation. We have implemented the hybrid simulation algorithm and have validated it by simulating the statistical behavior of the well-known lambda phage switch. Hybrid simulation provides a new method for exploring the sources and nature of stochastic behavior in cells.     

Sensitivity analysis of Markov models for communities of competing sessile organisms

1. Communities of competing sessile organisms are often modelled using Markov chains. Sensitivity analysis of the stationary distribution of these models tells us how we expect the abundance of each organism to respond to changes in interactions between species. This is important for conservation and management. 2. Markov models for such communities have usually been formulated in discrete time. Each column of the discrete-time transition matrix must sum to 1 (column stochasticity). Sensitivity analysis therefore involves defining a pattern of compensation that maintains column stochasticity as a single transition probability changes. There is little biological theory about the appropriate compensation pattern, but the usual choices involve changing only the elements of a single column of the transition matrix. 3. I argue that if the underlying dynamics occur in continuous time, each transition probability is the net outcome of direct and many indirect interactions. 4. Determining the consequences of changing a single direct interaction will often be of interest. I show how this can be achieved using a continuous-time model. The resulting discrete-time compensation pattern is quite different from those that have been considered elsewhere, with changes occurring in many columns. 5. I also show how to determine which direct interactions are being changed under any discrete-time compensation pattern.     

Neural system modeling and simulation using Hybrid Functional Petri Net

The Petri net formalism has been proved to be powerful in biological modeling. It not only boasts of a most intuitive graphical presentation but also combines the methods of classical systems biology with the discrete modeling technique. Hybrid Functional Petri Net (HFPN) was proposed specially for biological system modeling. An array of well-constructed biological models using HFPN yielded very interesting results. In this paper, we propose a method to represent neural system behavior, where biochemistry and electrical chemistry are both included using the Petri net formalism. We built a model for the adrenergic system using HFPN and employed quantitative analysis. Our simulation results match the biological data well, showing that the model is very effective. Predictions made on our model further manifest the modeling power of HFPN and improve the understanding of the adrenergic system. The file of our model and more results with their analysis are available in our supplementary material.