The numerical solution of stiff differential equations

This paper first discusses the conditions in which a set of differential equations should give stable solutions, starting with linear systems assuming that these do not differ greatly in this respect from non-linear systems. Methods of investigating the stability of particular systems are briefly discussed. Most real biochemical systems are known from observation to be stable, but little is known of the regions over which stability persists; moreover, models of biochemical systems may not be stable, because of inaccurate choice of parameter values.The separate problem of stability and accuracy in numerical methods of approximating the solution of systems of non-linear equations is then treated. Stress is laid on the consistently unsatisfactory results given by explicit methods for systems containing "stiff" equations, and implicit multistep methods are particularly recommended for this class of problem, which is likely to include many biochemical model systems. Finally, an iteration procedure likely to give convergence both in multistep methods and in the steady-state approach is recommended, and areas in which improvement in methods is likely to occur are outlined.     

Parameter estimation for stiff equations of biosystems using radial basis function networks

Background  

The modeling of dynamic systems requires estimating kinetic parameters from experimentally measured time-courses. Conventional global optimization methods used for parameter estimation, e.g. genetic algorithms (GA), consume enormous computational time because they require iterative numerical integrations for differential equations. When the target model is stiff, the computational time for reaching a solution increases further.     

Practical aspects of the mathematical modeling of fermentation processes as a method of description,simulation, identification and optimization

Most mathematical models for describing the physiological state in fermentations lead to solutions of the so-called “stiff differential systems” during simulation on a digital computer There is no suitable conventional software for solving these systems As a result of a relatively extensive screening of suitable methods for the solution of “stiff” differential systems (about 200 methods) it may be concluded that the semiimplicit RUNGE -KUTTA -formulas of the ROSENBROCK type, which constitute a part of the collection of programmes STIFFSOLVER-80, are optimal for the simulation of fermentation processes For determining kinetic parameters from integral data the authors use the system of programmes BIOKIN Their practical application is discussed for 3 examples:

• 1 Growth of the yeast Saccharomyces cerevisiae and changes in the content of specific compounds (proteins and ergosterol)
• 2 Quantitative evaluation of “;direct oxygen transfer” in the submerged culture.
• 3 Biosynthesis of a new antibiotic substance mucidin

.     

Solution of large compartmental models using numerical transform inversion

Compartmental models of biological or physical systems are often described by a system of “stiff” differential equations. In this paper an algorithm for solving a system with linear coefficients is presented that employs numerical inversion of the Laplace transform of the model equations. The inversion algorithms and Gear's backward differentiation method are compared for two stiff test problems and a differential system governing a 27-compartment model of bile acid transport and metabolism. The inversion algorithm is reliable, requires modest computation time on a desktop computer and provides better accuracy than Gear's method, especially for the extremely stiff example.     

Function approximation approach to the inference of reduced NGnet models of genetic networks

Background  

The inference of a genetic network is a problem in which mutual interactions among genes are deduced using time-series of gene expression patterns. While a number of models have been proposed to describe genetic regulatory networks, this study focuses on a set of differential equations since it has the ability to model dynamic behavior of gene expression. When we use a set of differential equations to describe genetic networks, the inference problem can be defined as a function approximation problem. On the basis of this problem definition, we propose in this study a new method to infer reduced NGnet models of genetic networks.     

Parameter identification of the calvin photosynthesis cycle

Summary This article is concerned with the determination of kinetic parameters of the Calvin photosynthesis cycle which is described by seventeen nonlinear ordinary differential equations. It is shown that the task requires dynamic data for several sets of initial conditions. The numerical technique is based upon an algorithm for non-linear optimization and Gear's numerical integration scheme for stiff systems of differential equations. The sensitivity of the parameters to noise in the data is tested with a method adapted from Rosenbrook and Storey. A preliminary set of parameters has been obtained from a preliminary set of experimental data. The numerical methods are then tested with synthetic data derived from these parameters. The mathematical model and the results obtained in the simulation are used as an aid in designing new experiments.     

Global Versus Local Centrality in Evolution of Yeast Protein Network

It is shown here that in the yeast protein interaction network the global centrality measure (betweenness) depends on the protein evolutionary age (i.e., on historic contingency) more weakly than the local centrality measure (degree). This phenomenon is not observed in mutational duplication-and-divergence models. The network domains responsible for this difference deal with DNA/RNA information processing, regulation, and cell cycle. A selection vector can operate in these domains, which integrates the network activity and thus compensates for the process of mutational divergence. Electronic supplementary material  The online version of this article (doi:) contains supplementary material, which is available to authorized users.     

Mathematical simulation and analysis of cellular metabolism and regulation.

MOTIVATION: A better understanding of the biological phenomena observed in cells requires the creation and analysis of mathematical models of cellular metabolism and physiology. The formulation and study of such models must also be simplified as far as possible to cope with the increasing complexity demanded and exponential accumulation of the metabolic reconstructions computed from sequenced genomes. RESULTS: A mathematical simulation workbench, DBsolve, has been developed to simplify the derivation and analysis of mathematical models. It combines: (i) derivation of large-scale mathematical models from metabolic reconstructions and other data sources; (ii) solving and parameter continuation of non-linear algebraic equations (NAEs), including metabolic control analysis; (iii) solving the non-linear stiff systems of ordinary differential equations (ODEs); (iv) bifurcation analysis of ODEs; (v) parameter fitting to experimental data or functional criteria based on constrained optimization. The workbench has been successfully used for dynamic metabolic modeling of some typical biochemical networks (Dolgacheva et al., Biochemistry (Moscow), 6, 1063-1068, 1996; Goldstein and Goryanin, Mol. Biol. (Moscow), 30, 976-983, 1996), including microbial glycolytic pathways, signal transduction pathways and receptor-ligand interactions. AVAILABILITY: DBsolve 5. 00 is freely available from http://websites.ntl.com/ approximately igor.goryanin. CONTACT: gzz78923@ggr.co.uk     

A comparison of non-standard solvers for ODEs describing cellular reactions in the heart

Mathematical models for the electrical activity in cardiac cells are normally formulated as systems of ordinary differential equations (ODEs). The equations are nonlinear and describe processes occurring on a wide range of time scales. Under normal accuracy requirements, this makes the systems stiff and therefore challenging to solve numerically. As standard implicit solvers are difficult to implement, explicit solvers such as the forward Euler method are commonly used, despite their poor efficiency. Non-standard formulations of the forward Euler method, derived from the analytical solution of linear ODEs, can give significantly improved performance while maintaining simplicity of implementation. In this paper we study the performance of three non-standard methods on two different cell models with comparable complexity but very different stiffness characteristics.     

A comparison of non-standard solvers for ODEs describing cellular reactions in the heart

Mathematical models for the electrical activity in cardiac cells are normally formulated as systems of ordinary differential equations (ODEs). The equations are nonlinear and describe processes occurring on a wide range of time scales. Under normal accuracy requirements, this makes the systems stiff and therefore challenging to solve numerically. As standard implicit solvers are difficult to implement, explicit solvers such as the forward Euler method are commonly used, despite their poor efficiency. Non-standard formulations of the forward Euler method, derived from the analytical solution of linear ODEs, can give significantly improved performance while maintaining simplicity of implementation. In this paper we study the performance of three non-standard methods on two different cell models with comparable complexity but very different stiffness characteristics.     

Statistical model comparison applied to common network motifs

Background  

Network motifs are small modules that show interesting functional and dynamic properties, and are believed to be the building blocks of complex cellular processes. However, the mechanistic details of such modules are often unknown: there is uncertainty about the motif architecture as well as the functional form and parameter values when converted to ordinary differential equations (ODEs). This translates into a number of candidate models being compatible with the system under study. A variety of statistical methods exist for ranking models including maximum likelihood-based and Bayesian methods. Our objective is to show how such methods can be applied in a typical systems biology setting.     

Improved discretisation and linearisation of active tension in strongly coupled cardiac electro-mechanics simulations

Mathematical models of cardiac electro-mechanics typically consist of three tightly coupled parts: systems of ordinary differential equations describing electro-chemical reactions and cross-bridge dynamics in the muscle cells, a system of partial differential equations modelling the propagation of the electrical activation through the tissue and a nonlinear elasticity problem describing the mechanical deformations of the heart muscle. The complexity of the mathematical model motivates numerical methods based on operator splitting, but simple explicit splitting schemes have been shown to give severe stability problems for realistic models of cardiac electro-mechanical coupling. The stability may be improved by adopting semi-implicit schemes, but these give rise to challenges in updating and linearising the active tension. In this paper we present an operator splitting framework for strongly coupled electro-mechanical simulations and discuss alternative strategies for updating and linearising the active stress component. Numerical experiments demonstrate considerable performance increases from an update method based on a generalised Rush–Larsen scheme and a consistent linearisation of active stress based on the first elasticity tensor.     

Energetics and dynamics of global integrals modeling interaction between stiff filaments

The attractive and spacing interaction between pairs of filaments via cross-linkers, e.g. myosin oligomers connecting actin filaments, is modeled by global integral kernels for negative binding energies between two intersecting stiff and long rods in a (projected) two-dimensional situation, for simplicity. Whereas maxima of the global energy functional represent intersection angles of ‘minimal contact’ between the filaments, minima are approached for energy values tending to −∞, representing the two degenerate states of parallel and anti-parallel filament alignment. Standard differential equations of negative gradient flow for such energy functionals show convergence of solutions to one of these degenerate equilibria in finite time, thus called ‘super-stable’ states. By considering energy variations under virtual rotation or translation of one filament with respect to the other, integral kernels for the resulting local forces parallel and orthogonal to the filament are obtained. For the special modeling situation that these variations only activate ‘spring forces’ in direction of the cross-links, explicit formulas for total torque and translational forces are given and calculated for typical examples. Again, the two degenerate alignment states are locally ‘super-stable’ equilibria of the assumed over-damped dynamics, but also other stable states of orthogonal arrangement and different asymptotic behavior can occur. These phenomena become apparent if stochastic perturbations of the local force kernels are implemented as additive Gaussian noise induced by the cross-link binding process with appropriate scaling. Then global filament dynamics is described by a certain type of degenerate stochastic differential equations yielding asymptotic stationary processes around the alignment states, which have generalized, namely bimodal Gaussian distributions. Moreover, stochastic simulations reveal characteristic sliding behavior as it is observed for myosin-mediated interaction between actin filaments. Finally, the forgoing explicit and asymptotic analysis as well as numerical simulations are extended to the more realistic modeling situation with filaments of finite length.

---

Electronic supplementary material  The online version of this article (doi:) contains supplementary material, which is available to authorized users.      

Bifurcation analysis of continuous biochemical reactor models

The validity of a biochemical reactor model often is evaluated by comparing transient responses to experimental data. Dynamic simulation can be a rather inefficient and ineffective tool for analyzing bioreactor models that exhibit complex nonlinear behavior. Bifurcation analysis is a powerful tool for obtaining a more efficient and complete characterization of the model behavior. To illustrate the power of bifurcation analysis, the steady-state and transient behavior of three continuous bioreactor models consisting of a small number of ordinary differential equations are investigated. Several important features, as well as potential limitations, that are difficult to ascertain via dynamic simulation are disclosed through the bifurcation analysis. The results motivate the use of dynamic simulation and bifurcation analysis as complementary tools for analyzing the nonlinear behavior of bioreactor models.     

CompuCell, a multi-model framework for simulation of morphogenesis

MOTIVATION: CompuCell is a multi-model software framework for simulation of the development of multicellular organisms known as morphogenesis. It models the interaction of the gene regulatory network with generic cellular mechanisms, such as cell adhesion, division, haptotaxis and chemotaxis. A combination of a state automaton with stochastic local rules and a set of differential equations, including subcellular ordinary differential equations and extracellular reaction-diffusion partial differential equations, model gene regulation. This automaton in turn controls the differentiation of the cells, and cell-cell and cell-extracellular matrix interactions that give rise to cell rearrangements and pattern formation, e.g. mesenchymal condensation. The cellular Potts model, a stochastic model that accurately reproduces cell movement and rearrangement, models cell dynamics. All these models couple in a controllable way, resulting in a powerful and flexible computational environment for morphogenesis, which allows for simultaneous incorporation of growth and spatial patterning. RESULTS: We use CompuCell to simulate the formation of the skeletal architecture in the avian limb bud. AVAILABILITY: Binaries and source code for Microsoft Windows, Linux and Solaris are available for download from http://sourceforge.net/projects/compucell/     

Multiscale dynamic modeling and simulation of a biorefinery

A biorefinery comprises a variety of process steps to synthesize products from sustainable natural resources. Dynamic plant-wide simulation enhances the process understanding, leads to improved cost efficiency and enables model-based operation and control. It is thereby important for an increased competitiveness to conventional processes. To this end, we developed a Modelica library with replaceable building blocks that allow dynamic modeling of an entire biorefinery. For the microbial conversion step, we built on the dynamic flux balance analysis (DFBA) approach to formulate process models for the simulation of cellular metabolism under changing environmental conditions. The resulting system of differential-algebraic equations with embedded optimization criteria (DAEO) is solved by a tailor-made toolbox. In summary, our modeling framework comprises three major pillars: A Modelica library of dynamic unit operations, an easy-to-use interface to formulate DFBA process models and a DAEO toolbox that allows simulation with standard environments based on the Modelica modeling language. A biorefinery model for dynamic simulation of the OrganoCat pretreatment process and microbial conversion of the resulting feedstock by Corynebacterium glutamicum serves as case study to demonstrate its practical relevance.     

RuleMonkey: software for stochastic simulation of rule-based models

Background  

The system-level dynamics of many molecular interactions, particularly protein-protein interactions, can be conveniently represented using reaction rules, which can be specified using model-specification languages, such as the BioNetGen language (BNGL). A set of rules implicitly defines a (bio)chemical reaction network. The reaction network implied by a set of rules is often very large, and as a result, generation of the network implied by rules tends to be computationally expensive. Moreover, the cost of many commonly used methods for simulating network dynamics is a function of network size. Together these factors have limited application of the rule-based modeling approach. Recently, several methods for simulating rule-based models have been developed that avoid the expensive step of network generation. The cost of these "network-free" simulation methods is independent of the number of reactions implied by rules. Software implementing such methods is now needed for the simulation and analysis of rule-based models of biochemical systems.     

Steady-State Kinetic Modeling Constrains Cellular Resting States and Dynamic Behavior

A defining characteristic of living cells is the ability to respond dynamically to external stimuli while maintaining homeostasis under resting conditions. Capturing both of these features in a single kinetic model is difficult because the model must be able to reproduce both behaviors using the same set of molecular components. Here, we show how combining small, well-defined steady-state networks provides an efficient means of constructing large-scale kinetic models that exhibit realistic resting and dynamic behaviors. By requiring each kinetic module to be homeostatic (at steady state under resting conditions), the method proceeds by (i) computing steady-state solutions to a system of ordinary differential equations for each module, (ii) applying principal component analysis to each set of solutions to capture the steady-state solution space of each module network, and (iii) combining optimal search directions from all modules to form a global steady-state space that is searched for accurate simulation of the time-dependent behavior of the whole system upon perturbation. Importantly, this stepwise approach retains the nonlinear rate expressions that govern each reaction in the system and enforces constraints on the range of allowable concentration states for the full-scale model. These constraints not only reduce the computational cost of fitting experimental time-series data but can also provide insight into limitations on system concentrations and architecture. To demonstrate application of the method, we show how small kinetic perturbations in a modular model of platelet P2Y₁ signaling can cause widespread compensatory effects on cellular resting states.