CADLIVE toolbox for MATLAB: automatic dynamic modeling of biochemical networks with comprehensive system analysis

Mathematical modeling has become a standard technique to understand the dynamics of complex biochemical systems. To promote the modeling, we had developed the CADLIVE dynamic simulator that automatically converted a biochemical map into its associated mathematical model, simulated its dynamic behaviors and analyzed its robustness. To enhance the feasibility by CADLIVE and extend its functions, we propose the CADLIVE toolbox available for MATLAB, which implements not only the existing functions of the CADLIVE dynamic simulator, but also the latest tools including global parameter search methods with robustness analysis. The seamless, bottom-up processes consisting of biochemical network construction, automatic construction of its dynamic model, simulation, optimization, and S-system analysis greatly facilitate dynamic modeling, contributing to the research of systems biology and synthetic biology. This application can be freely downloaded from http://www.cadlive.jp/CADLIVE_MATLAB/ together with an instruction.     

Cyto-Sim: a formal language model and stochastic simulator of membrane-enclosed biochemical processes

MOTIVATION: Compartments and membranes are the basis of cell topology and more than 30% of the human genome codes for membrane proteins. While it is possible to represent compartments and membrane proteins in a nominal way with many mathematical formalisms used in systems biology, few, if any, explicitly model the topology of the membranes themselves. Discrete stochastic simulation potentially offers the most accurate representation of cell dynamics. Since the details of every molecular interaction in a pathway are often not known, the relationship between chemical species in not necessarily best described at the lowest level, i.e. by mass action. Simulation is a form of computer-aided analysis, relying on human interpretation to derive meaning. To improve efficiency and gain meaning in an automatic way, it is necessary to have a formalism based on a model which has decidable properties. RESULTS: We present Cyto-Sim, a stochastic simulator of membrane-enclosed hierarchies of biochemical processes, where the membranes comprise an inner, outer and integral layer. The underlying model is based on formal language theory and has been shown to have decidable properties (Cavaliere and Sedwards, 2006), allowing formal analysis in addition to simulation. The simulator provides variable levels of abstraction via arbitrary chemical kinetics which link to ordinary differential equations. In addition to its compact native syntax, Cyto-Sim currently supports models described as Petri nets, can import all versions of SBML and can export SBML and MATLAB m-files. AVAILABILITY: Cyto-Sim is available free, either as an applet or a stand-alone Java program via the web page (http://www.cosbi.eu/Rpty_Soft_CytoSim.php). Other versions can be made available upon request.     

COPASI--a COmplex PAthway SImulator

MOTIVATION: Simulation and modeling is becoming a standard approach to understand complex biochemical processes. Therefore, there is a big need for software tools that allow access to diverse simulation and modeling methods as well as support for the usage of these methods. RESULTS: Here, we present COPASI, a platform-independent and user-friendly biochemical simulator that offers several unique features. We discuss numerical issues with these features; in particular, the criteria to switch between stochastic and deterministic simulation methods, hybrid deterministic-stochastic methods, and the importance of random number generator numerical resolution in stochastic simulation. AVAILABILITY: The complete software is available in binary (executable) for MS Windows, OS X, Linux (Intel) and Sun Solaris (SPARC), as well as the full source code under an open source license from http://www.copasi.org.     

Design and implementation of a qualitative simulation model of lambda phage infection

MOTIVATION: Molecular biology databases hold a large number of empirical facts about many different aspects of biological entities. That data is static in the sense that one cannot ask a database 'What effect has protein A on gene B?' or 'Do gene A and gene B interact, and if so, how?'. Those questions require an explicit model of the target organism. Traditionally, biochemical systems are modelled using kinetics and differential equations in a quantitative simulator. For many biological processes however, detailed quantitative information is not available, only qualitative or fuzzy statements about the nature of interactions. RESULTS: We designed and implemented a qualitative simulation model of lambda phage growth control in Escherichia coli based on the existing simulation environment QSim. Qualitative reasoning can serve as the basis for automatic transformation of contents of genomic databases into interactive modelling systems that can reason about the relations and interactions of biological entities.      

Exact Hybrid Particle/Population Simulation of Rule-Based Models of Biochemical Systems

Detailed modeling and simulation of biochemical systems is complicated by the problem of combinatorial complexity, an explosion in the number of species and reactions due to myriad protein-protein interactions and post-translational modifications. Rule-based modeling overcomes this problem by representing molecules as structured objects and encoding their interactions as pattern-based rules. This greatly simplifies the process of model specification, avoiding the tedious and error prone task of manually enumerating all species and reactions that can potentially exist in a system. From a simulation perspective, rule-based models can be expanded algorithmically into fully-enumerated reaction networks and simulated using a variety of network-based simulation methods, such as ordinary differential equations or Gillespie''s algorithm, provided that the network is not exceedingly large. Alternatively, rule-based models can be simulated directly using particle-based kinetic Monte Carlo methods. This “network-free” approach produces exact stochastic trajectories with a computational cost that is independent of network size. However, memory and run time costs increase with the number of particles, limiting the size of system that can be feasibly simulated. Here, we present a hybrid particle/population simulation method that combines the best attributes of both the network-based and network-free approaches. The method takes as input a rule-based model and a user-specified subset of species to treat as population variables rather than as particles. The model is then transformed by a process of “partial network expansion” into a dynamically equivalent form that can be simulated using a population-adapted network-free simulator. The transformation method has been implemented within the open-source rule-based modeling platform BioNetGen, and resulting hybrid models can be simulated using the particle-based simulator NFsim. Performance tests show that significant memory savings can be achieved using the new approach and a monetary cost analysis provides a practical measure of its utility.     

DFBAlab: a fast and reliable MATLAB code for dynamic flux balance analysis

Background

Dynamic Flux Balance Analysis (DFBA) is a dynamic simulation framework for biochemical processes. DFBA can be performed using different approaches such as static optimization (SOA), dynamic optimization (DOA), and direct approaches (DA). Few existing simulators address the theoretical and practical challenges of nonunique exchange fluxes or infeasible linear programs (LPs). Both are common sources of failure and inefficiencies for these simulators.

Results

DFBAlab, a MATLAB-based simulator that uses the LP feasibility problem to obtain an extended system and lexicographic optimization to yield unique exchange fluxes, is presented. DFBAlab is able to simulate complex dynamic cultures with multiple species rapidly and reliably, including differential-algebraic equation (DAE) systems. In addition, DFBAlab’s running time scales linearly with the number of species models. Three examples are presented where the performance of COBRA, DyMMM and DFBAlab are compared.

Conclusions

Lexicographic optimization is used to determine unique exchange fluxes which are necessary for a well-defined dynamic system. DFBAlab does not fail during numerical integration due to infeasible LPs. The extended system obtained through the LP feasibility problem in DFBAlab provides a penalty function that can be used in optimization algorithms.

Electronic supplementary material

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

A stochastic simulator of a blood product donation environment with demand spikes and supply shocks

The availability of an adequate blood supply is a critical public health need. An influenza epidemic or another crisis affecting population mobility could create a critical donor shortage, which could profoundly impact blood availability. We developed a simulation model for the blood supply environment in the United States to assess the likely impact on blood availability of factors such as an epidemic. We developed a simulator of a multi-state model with transitions among states. Weekly numbers of blood units donated and needed were generated by negative binomial stochastic processes. The simulator allows exploration of the blood system under certain conditions of supply and demand rates, and can be used for planning purposes to prepare for sudden changes in the public's health. The simulator incorporates three donor groups (first-time, sporadic, and regular), immigration and emigration, deferral period, and adjustment factors for recruitment. We illustrate possible uses of the simulator by specifying input values for an 8-week flu epidemic, resulting in a moderate supply shock and demand spike (for example, from postponed elective surgeries), and different recruitment strategies. The input values are based in part on data from a regional blood center of the American Red Cross during 1996-2005. Our results from these scenarios suggest that the key to alleviating deficit effects of a system shock may be appropriate timing and duration of recruitment efforts, in turn depending critically on anticipating shocks and rapidly implementing recruitment efforts.     

WinBEST-KIT: Windows-based biochemical reaction simulator for metabolic pathways

We have implemented an efficient, user-friendly biochemical reaction simulator called Web-based BEST-KIT (Biochemical Engineering System analyzing Tool-KIT) for analyzing large-scale nonlinear networks such as metabolic pathways. Users can easily design and analyze an arbitrary reaction scheme through the Internet and an efficient graphical user interface without considering the mathematical equations. The reaction scheme can include several reaction types, which are represented by both the mass action law (mass balance) and approximated velocity functions of enzyme kinetics at steady state, such as Michaelis-Menten, Hill cooperative, Competitive inhibition. However, since all modules in Web-based BEST-KIT have been developed in Java applet style, users cannot optionally make use of original mathematical equations in addition to the prepared equations. In the present study, we have developed a new version of BEST-KIT (for Microsoft Windows called WinBEST-KIT) to allow users to define original mathematical equations and to customize these equations very easily as user-defined reaction symbols. The following powerful system-analytical methods are prepared for system analysis: time-course calculation, parameter scanning, estimation of the values of unknown kinetic parameters based on experimentally observed time-course data of reactants, dynamic response of reactants against virtual external perturbations, and real-time simulation (Virtual Dry Lab).     

Usefulness of a New Three-Dimensional Simulator System for Radiofrequency Ablation

Multipuncture radiofrequency ablation is expected to produce a large ablated area and reduce intrahepatic recurrence of hepatocellular carcinoma; however, it requires considerable skill. This study evaluated the utility of a new simulator system for multipuncture radiofrequency ablation. To understand positioning of multipuncture electrodes on three-dimensional images, we developed a new technology by expanding real-time virtual ultrasonography. We performed 21 experimental punctures in phantoms. Electrode insertion directions and positions were confirmed on computed tomography, and accuracy and utility of the simulator system were evaluated by measuring angles and intersections for each electrode. Moreover, to appropriately assess placement of the three electrodes, puncture procedures with or without the simulator were performed by experts and non-experts. Technical success was defined as maximum angle and distance ratio, as calculated by maximum and minimum distances between electrodes. In punctures using 2 electrodes, correlations between angles on each imaging modality were strong (ultrasound vs. simulator: r = 0.991, p<0.001, simulator vs. computed tomography: r = 0.991, p<0.001, ultrasound vs. computed tomography: r = 0.999, p<0.001). Correlations between distances in each imaging modality were also strong (ultrasound vs. simulator: r = 0.993, p<0.001; simulator vs. computed tomography: r = 0.994, p<0.001; ultrasound vs. computed tomography: r = 0.994, p<0.001). In cases with 3 electrodes, distances between each electrode correlated strongly (yellow-labeled vs. red-labeled: r = 0.980, p<0.001; red-labeled vs. blue-labeled: r = 0.953, p<0.001; yellow-labeled vs. blue-labeled: r = 0.953, p<0.001). Both angle and distance ratio (expert with simulator vs. without simulator; p = 0.03, p = 0.02) were significantly smaller in procedures performed by experts using the simulator system. The new simulator system appears to accurately guide electrode positioning. This simulator system could allow multipuncture radiofrequency ablation to be performed more effectively and comfortably.     

The high number of neurons contributes to the robustness of the locust flight-CPG against parameter variation

 Real pattern-generating networks often consist of more neurons than necessary for the production of a certain rhythm. We investigated the question of whether these neurons contribute to the robustness of a pattern-generating system of using the central pattern generator (CPG) for flight of the locust, generating the deafferented activity pattern of wing elevator and wing depressor motoneurons, as an example of a rhythm-generating system. The neuronal network was reconstructed, based on the known connectivity of the interneurons in the flight CPG, using a biologically orientated network simulator (BioSim 3.0). This simulator allows a physiologically realistic simulation of particular neurons as well as the synaptic connections between them. The flight CPG consists of at least five cyclic loops. The simulation shows that each of them is in principle able to produce a rhythm comparable to the rhythm produced by the whole network, i.e. the ‘deafferented’ flight pattern of elevator and depressor motoneurons. Varying the parameter ‘synaptic strength’ in each of these loops and in the complete system shows that this parameter can be changed within certain ranges without loosing the ability to produce oscillations. These ranges are much smaller in each of the subloops than in the whole network. This result demonstrates that the robustness of the system is increased by supranumerary neurons and connections. Changing the active properties of the simulated neurons so that they are able to produce plateau potentials has no effect on the robustness of the simulated network. Received: 13 April 1994/Accepted in revised form: 15 September 1994     

The SBML discrete stochastic models test suite

MOTIVATION: Stochastic simulation is a very important tool for mathematical modelling. However, it is difficult to check the correctness of a stochastic simulator, since any two realizations from a single model will typically be different. RESULTS: We have developed a test suite of stochastic models that have been solved either analytically or using numerical methods. This allows the accuracy of stochastic simulators to be tested against known results. The test suite is already being used by a number of stochastic simulator developers. AVAILABILITY: The latest version of the test suite can be obtained from http://www.calibayes.ncl.ac.uk/Resources/dsmts/ and is licensed under GNU Lesser General Public License.     

Accelerated maximum likelihood parameter estimation for stochastic biochemical systems

ABSTRACT: BACKGROUND: A prerequisite for the mechanistic simulation of a biochemical system is detailed knowledge of its kinetic parameters. Despite recent experimental advances, the estimation of unknown parameter values from observed data is still a bottleneck for obtaining accurate simulation results. Many methods exist for parameter estimation in deterministic biochemical systems; methods for discrete stochastic systems are less well developed. Given the probabilistic nature of stochastic biochemical models, a natural approach is to choose parameter values that maximize the probability of the observed data with respect to the unknown parameters, a.k.a. the maximum likelihood parameter estimates (MLEs). MLE computation for all but the simplest models requires the simulation of many system trajectories that are consistent with experimental data. For models with unknown parameters, this presents a computational challenge, as the generation of consistent trajectories can be an extremely rare occurrence. RESULTS: We have developed Monte Carlo Expectation-Maximization with Modified Cross-Entropy Method (MCEM2): an accelerated method for calculating MLEs that combines advances in rare event simulation with a computationally efficient version of the Monte Carlo expectation-maximization (MCEM) algorithm. Our method requires no prior knowledge regarding parameter values, and it automatically provides a multivariate parameter uncertainty estimate. We applied the method to five stochastic systems of increasing complexity, progressing from an analytically tractable pure-birth model to a computationally demanding model of yeast-polarization. Our results demonstrate that MCEM2 substantially accelerates MLE computation on all tested models when compared to a stand-alone version of MCEM. Additionally, we show how our method identifies parameter values for certain classes of models more accurately than two recently proposed computationally efficient methods. CONCLUSIONS: This work provides a novel, accelerated version of a likelihood-based parameter estimation method that can be readily applied to stochastic biochemical systems. In addition, our results suggest opportunities for added efficiency improvements that will further enhance our ability to mechanistically simulate biological processes.     

CADCS simulation of the closed-loop cardiovascular system

A pulsatile simulator of the closed-loop cardiovascular system, designed to solve simulation, identification and control problems in a research and education context, is presented. Its implimentation makes use of a command-driven interactive program for simulation of non-linear ordinary differential equations. The flexibility of the simulator is demonstrated by the results presented which refer to a basal steady-state circulatory condition as well as a transient induced by an abrupt change in peripheral resistance.     

Serial NetEvolve: a flexible utility for generating serially-sampled sequences along a tree or recombinant network

SUMMARY: Serial NetEvolve is a flexible simulation program that generates DNA sequences evolved along a tree or recombinant network. It offers a user-friendly Windows graphical interface and a Windows or Linux simulator with a diverse selection of parameters to control the evolutionary model. Serial NetEvolve is a modification of the Treevolve program with the following additional features: simulation of serially-sampled data, the choice of either a clock-like or a variable rate model of sequence evolution, sampling from the internal nodes and the output of the randomly generated tree or network in our newly proposed NeTwick format. AVAILABILITY: From website http://biorg.cis.fiu.edu/SNE Contacts: giri@cis.fiu.edu SUPPLEMENTARY INFORMATION: Manual and examples available from http://biorg.cis.fiu.edu/SNE.     

Parameter estimation using Simulated Annealing for S-system models of biochemical networks

MOTIVATION: High-throughput technologies now allow the acquisition of biological data, such as comprehensive biochemical time-courses at unprecedented rates. These temporal profiles carry topological and kinetic information regarding the biochemical network from which they were drawn. Retrieving this information will require systematic application of both experimental and computational methods. RESULTS: S-systems are non-linear mathematical approximative models based on the power-law formalism. They provide a general framework for the simulation of integrated biological systems exhibiting complex dynamics, such as genetic circuits, signal transduction and metabolic networks. We describe how the heuristic optimization technique simulated annealing (SA) can be effectively used for estimating the parameters of S-systems from time-course biochemical data. We demonstrate our methods using three artificial networks designed to simulate different network topologies and behavior. We then end with an application to a real biochemical network by creating a working model for the cadBA system in Escherichia coli. AVAILABILITY: The source code written in C++ is available at http://www.engg.upd.edu.ph/~naval/bioinformcode.html. All the necessary programs including the required compiler are described in a document archived with the source code. SUPPLEMENTARY INFORMATION: Supplementary material is available at Bioinformatics online.     

Predicting knee replacement damage in a simulator machine using a computational model with a consistent wear factor

Wear of ultrahigh molecular weight polyethylene remains a primary factor limiting the longevity of total knee replacements (TKRs). However, wear testing on a simulator machine is time consuming and expensive, making it impractical for iterative design purposes. The objectives of this paper were first, to evaluate whether a computational model using a wear factor consistent with the TKR material pair can predict accurate TKR damage measured in a simulator machine, and second, to investigate how choice of surface evolution method (fixed or variable step) and material model (linear or nonlinear) affect the prediction. An iterative computational damage model was constructed for a commercial knee implant in an AMTI simulator machine. The damage model combined a dynamic contact model with a surface evolution model to predict how wear plus creep progressively alter tibial insert geometry over multiple simulations. The computational framework was validated by predicting wear in a cylinder-on-plate system for which an analytical solution was derived. The implant damage model was evaluated for 5 million cycles of simulated gait using damage measurements made on the same implant in an AMTI machine. Using a pin-on-plate wear factor for the same material pair as the implant, the model predicted tibial insert wear volume to within 2% error and damage depths and areas to within 18% and 10% error, respectively. Choice of material model had little influence, while inclusion of surface evolution affected damage depth and area but not wear volume predictions. Surface evolution method was important only during the initial cycles, where variable step was needed to capture rapid geometry changes due to the creep. Overall, our results indicate that accurate TKR damage predictions can be made with a computational model using a constant wear factor obtained from pin-on-plate tests for the same material pair, and furthermore, that surface evolution method matters only during the initial "break in" period of the simulation.     

A simple unconstrained dynamic knee simulator

The design of a simple dynamic knee simulator is described. In the simulator the joint dynamics are reproduced in-vitro in a knee specimen by controlling the time-histories of the tensions in two flexible cables acting as lumped muscle group equivalents, without constraining the natural conjunct and passive motions of the specimen. The two cable tensions acting individually are used to control the active flexion/extension motion, while their simultaneous action is used to control joint compressive force. The characteristics of the electrohydraulic servo system acting under real-time microprocessor control are described. The system performance during simulation of an idealized level-walking function is evaluated.     

Simulation system of spinal cord motor nuclei and associated nerves and muscles,in a Web-based architecture

A Web-based simulation system of the spinal cord circuitry responsible for muscle control is described. The simulator employs two-compartment motoneuron models for S, FR and FF types, with synaptic inputs acting through conductance variations. Four motoneuron pools with their associated interneurons are represented in the simulator, with the possibility of inclusion of more than 2,000 neurons and 2,000,000 synapses. Each motoneuron action potential is followed, after a conduction delay, by a motor unit potential and a motor unit twitch. The sums of all motor unit potentials and twitches result in the electromyogram (EMG), and the muscle force, respectively. Inputs to the motoneuron pool come from populations of interneurons (Ia reciprocal inhibitory interneurons, Ib interneurons, and Renshaw cells) and from stochastic point processes associated with descending tracts. To simulate human electrophysiological experiments, the simulator incorporates external nerve stimulation with orthodromic and antidromic propagation. This provides the mechanisms for reflex generation and activation of spinal neuronal circuits that modulate the activity of another motoneuron pool (e.g., by reciprocal inhibition). The generation of the H-reflex by the Ia-motoneuron pool system and its modulation by spinal cord interneurons is included in the simulation system. Studies with the simulator may include the statistics of individual motoneuron or interneuron spike trains or the collective effect of a motor nucleus on the dynamics of muscle force control. Properties associated with motor-unit recruitment, motor-unit synchronization, recurrent inhibition and reciprocal inhibition may be investigated.