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.     

Efficient solution of ordinary differential equations modeling electrical activity in cardiac cells

The contraction of the heart is preceded and caused by a cellular electro-chemical reaction, causing an electrical field to be generated. Performing realistic computer simulations of this process involves solving a set of partial differential equations, as well as a large number of ordinary differential equations (ODEs) characterizing the reactive behavior of the cardiac tissue. Experiments have shown that the solution of the ODEs contribute significantly to the total work of a simulation, and there is thus a strong need to utilize efficient solution methods for this part of the problem. This paper presents how an efficient implicit Runge-Kutta method may be adapted to solve a complicated cardiac cell model consisting of 31 ODEs, and how this solver may be coupled to a set of PDE solvers to provide complete simulations of the electrical activity.     

Modelling and convergence in arterial wall simulations using a parallel FETI solution strategy

Arterial walls are characterised by nearly incompressible, anisotropic, hyperelastic material behaviour. Several polyconvex material functions representing such materials are considered and adjusted to experimental data. For all of these functions and for different parameter sets numerical simulations using a three-dimensional model of a diseased artery are performed. A finite element tearing and interconnecting-dual primal domain decomposition algorithm is used to solve the linearised systems of equations. The numerical performance of the different models is discussed with respect to convergence of the linear and nonlinear solvers.     

Estimating a predator-prey dynamical model with the parameter cascades method

Summary .   Ordinary differential equations (ODEs) are widely used in ecology to describe the dynamical behavior of systems of interacting populations. However, systems of ODEs rarely provide quantitative solutions that are close to real field observations or experimental data because natural systems are subject to environmental and demographic noise and ecologists are often uncertain about the correct parameterization. In this article we introduce "parameter cascades" as an improved method to estimate ODE parameters such that the corresponding ODE solutions fit the real data well. This method is based on the modified penalized smoothing with the penalty defined by ODEs and a generalization of profiled estimation, which leads to fast estimation and good precision for ODE parameters from noisy data. This method is applied to a set of ODEs originally developed to describe an experimental predator–prey system that undergoes oscillatory dynamics. The new parameterization considerably improves the fit of the ODE model to the experimental data sets. At the same time, our method reveals that important structural assumptions that underlie the original ODE model are essentially correct. The mathematical formulations of the two nonlinear interaction terms (functional responses) that link the ODEs in the predator–prey model are validated by estimating the functional responses nonparametrically from the real data. We suggest two major applications of "parameter cascades" to ecological modeling: It can be used to estimate parameters when original data are noisy, missing, or when no reliable priori estimates are available; it can help to validate the structural soundness of the mathematical modeling approach.     

Moment Equations and Dynamics of a Household SIS Epidemiological Model

An SIS epidemiological model of individuals partitioned into households is studied, where infections take place either within or between households, the latter generally happening much less frequently. The model is explored using stochastic spatial simulations, as well as mathematical models which consist of an infinite system of ordinary differential equations for the moments of the distribution describing the proportions of individuals who are infectious among households. Various moment-closure approximations are used to truncate the system of ODEs to finite systems of equations. These approximations can sometimes lead to a system of ill-behaved ODEs which predict moments which become negative or unbounded. A reparametrization of the ODEs is then developed, which forces all moments to satisfy necessary constraints.Changing the proportion of contacts within and between households does not change the endemic equilibrium, but does affect the amount of time it takes to approach the fixed point; increasing the proportion of contacts within households slows the spread of the infection toward endemic equilibrium. The system of moment equations does describe this phenomenon, although less accurately in the limit as the proportion of between-household contacts approaches zero. The results indicate that although controlling the movement of individuals does not affect the long-term frequency of an infection with SIS dynamics, it can have a large effect on the time-scale of the dynamics, which may provide an opportunity for other controls such as immunizations to be applied.     

Parameter estimation of kinetic models from metabolic profiles: two-phase dynamic decoupling method

MOTIVATION: Time-series measurements of metabolite concentration have become increasingly more common, providing data for building kinetic models of metabolic networks using ordinary differential equations (ODEs). In practice, however, such time-course data are usually incomplete and noisy, and the estimation of kinetic parameters from these data is challenging. Practical limitations due to data and computational aspects, such as solving stiff ODEs and finding global optimal solution to the estimation problem, give motivations to develop a new estimation procedure that can circumvent some of these constraints. RESULTS: In this work, an incremental and iterative parameter estimation method is proposed that combines and iterates between two estimation phases. One phase involves a decoupling method, in which a subset of model parameters that are associated with measured metabolites, are estimated using the minimization of slope errors. Another phase follows, in which the ODE model is solved one equation at a time and the remaining model parameters are obtained by minimizing concentration errors. The performance of this two-phase method was tested on a generic branched metabolic pathway and the glycolytic pathway of Lactococcus lactis. The results showed that the method is efficient in getting accurate parameter estimates, even when some information is missing.     

A Lagrange-based generalised formulation for the equations of motion of simple walking models

Simple 2D models of walking often approximate the human body to multi-link dynamic systems, where body segments are represented by rigid links connected by frictionless hinge joints. Performing forward dynamics on the equations of motion (EOM) of these systems can be used to simulate their movement. However, deriving these equations can be time consuming. Using Lagrangian mechanics, a generalised formulation for the EOM of n-link open-loop chains is derived. This can be used for single support walking models. This has an advantage over Newton–Euler mechanics in that it is independent of coordinate system and prior knowledge of the ground reaction force (GRF) is not required. Alternative strategies, such as optimisation algorithms, can be used to estimate joint activation and simulate motion. The application of Lagrange multipliers, to enforce motion constraints, is used to adapt this general formulation for application to closed-loop chains. This can be used for double support walking models. Finally, inverse dynamics are used to calculate the GRF for these general n-link chains. The necessary constraint forces to maintain a closed-loop chain, calculated from the Lagrange multipliers, are one solution to the indeterminate problem of GRF distribution in double support models. An example of this method’s application is given, whereby an optimiser estimates the joint moments by tracking kinematic data.     

Forcing Function Diagnostics for Nonlinear Dynamics

Summary .  This article investigates the problem of model diagnostics for systems described by nonlinear ordinary differential equations (ODEs). I propose modeling lack of fit as a time-varying correction to the right-hand side of a proposed differential equation. This correction can be described as being a set of additive forcing functions, estimated from data. Representing lack of fit in this manner allows us to graphically investigate model inadequacies and to suggest model improvements. I derive lack-of-fit tests based on estimated forcing functions. Model building in partially observed systems of ODEs is particularly difficult and I consider the problem of identification of forcing functions in these systems. The methods are illustrated with examples from computational neuroscience.     

Application of Petri net based analysis techniques to signal transduction pathways

Background  

Signal transduction pathways are usually modelled using classical quantitative methods, which are based on ordinary differential equations (ODEs). However, some difficulties are inherent in this approach. On the one hand, the kinetic parameters involved are often unknown and have to be estimated. With increasing size and complexity of signal transduction pathways, the estimation of missing kinetic data is not possible. On the other hand, ODEs based models do not support any explicit insights into possible (signal-) flows within the network. Moreover, a huge amount of qualitative data is available due to high-throughput techniques. In order to get information on the systems behaviour, qualitative analysis techniques have been developed. Applications of the known qualitative analysis methods concern mainly metabolic networks. Petri net theory provides a variety of established analysis techniques, which are also applicable to signal transduction models. In this context special properties have to be considered and new dedicated techniques have to be designed.     

On the use of qualitative reasoning to simulate and identify metabolic pathways

MOTIVATION: Perhaps the greatest challenge of modern biology is to develop accurate in silico models of cells. To do this we require computational formalisms for both simulation (how according to the model the state of the cell evolves over time) and identification (learning a model cell from observation of states). We propose the use of qualitative reasoning (QR) as a unified formalism for both tasks. The two most commonly used alternative methods of modelling biochemical pathways are ordinary differential equations (ODEs), and logical/graph-based (LG) models. RESULTS: The QR formalism we use is an abstraction of ODEs. It enables the behaviour of many ODEs, with different functional forms and parameters, to be captured in a single QR model. QR has the advantage over LG models of explicitly including dynamics. To simulate biochemical pathways we have developed 'enzyme' and 'metabolite' QR building blocks that fit together to form models. These models are finite, directly executable, easy to interpret and robust. To identify QR models we have developed heuristic chemoinformatics graph analysis and machine learning procedures. The graph analysis procedure is a series of constraints and heuristics that limit the number of ways metabolites can combine to form pathways. The machine learning procedure is generate-and-test inductive logic programming. We illustrate the use of QR for modelling and simulation using the example of glycolysis. AVAILABILITY: All data and programs used are available on request.     

Determining the three-dimensional relation between the skeletal elements of the human shoulder complex

In this paper, we present an inverse kinematics method to determining human shoulder joint motion coupling relationship based on experimental data in the literature. This work focuses on transferring Euler-angle-based coupling equations into a relationship based on the Denavit–Hartenberg (DH) method. We use analytical inverse kinematics to achieve the transferring. For a specific posture, we can choose points on clavicle, scapula, and humerus and represent the end-effector positions based on Euler angles or DH method. For both Euler and DH systems, the end-effectors have the same Cartesian positions. Solving these equations related to end-effector positions yields DH joint angles for that posture. The new joint motion coupling relationship is obtained by polynomial and cosine fitting of the DH joint angles for all different postures.     

DeepCME: A deep learning framework for computing solution statistics of the chemical master equation

Stochastic models of biomolecular reaction networks are commonly employed in systems and synthetic biology to study the effects of stochastic fluctuations emanating from reactions involving species with low copy-numbers. For such models, the Kolmogorov’s forward equation is called the chemical master equation (CME), and it is a fundamental system of linear ordinary differential equations (ODEs) that describes the evolution of the probability distribution of the random state-vector representing the copy-numbers of all the reacting species. The size of this system is given by the number of states that are accessible by the chemical system, and for most examples of interest this number is either very large or infinite. Moreover, approximations that reduce the size of the system by retaining only a finite number of important chemical states (e.g. those with non-negligible probability) result in high-dimensional ODE systems, even when the number of reacting species is small. Consequently, accurate numerical solution of the CME is very challenging, despite the linear nature of the underlying ODEs. One often resorts to estimating the solutions via computationally intensive stochastic simulations. The goal of the present paper is to develop a novel deep-learning approach for computing solution statistics of high-dimensional CMEs by reformulating the stochastic dynamics using Kolmogorov’s backward equation. The proposed method leverages superior approximation properties of Deep Neural Networks (DNNs) to reliably estimate expectations under the CME solution for several user-defined functions of the state-vector. This method is algorithmically based on reinforcement learning and it only requires a moderate number of stochastic simulations (in comparison to typical simulation-based approaches) to train the “policy function”. This allows not just the numerical approximation of various expectations for the CME solution but also of its sensitivities with respect to all the reaction network parameters (e.g. rate constants). We provide four examples to illustrate our methodology and provide several directions for future research.     

A gradual update method for simulating the steady-state solution of stiff differential equations in metabolic circuits

Numerical simulation of differential equation systems plays a major role in the understanding of how metabolic network models generate particular cellular functions. On the other hand, the classical and technical problems for stiff differential equations still remain to be solved, while many elegant algorithms have been presented. To relax the stiffness problem, we propose new practical methods: the gradual update of differential-algebraic equations based on gradual application of the steady-state approximation to stiff differential equations, and the gradual update of the initial values in differential-algebraic equations. These empirical methods show a high efficiency for simulating the steady-state solutions for the stiff differential equations that existing solvers alone cannot solve. They are effective in extending the applicability of dynamic simulation to biochemical network models. Electronic supplementary material  The online version of this article (doi:) contains supplementary material, which is available to authorized users.     

Approximate parameter inference in systems biology using gradient matching: a comparative evaluation

Background

A challenging problem in current systems biology is that of parameter inference in biological pathways expressed as coupled ordinary differential equations (ODEs). Conventional methods that repeatedly numerically solve the ODEs have large associated computational costs. Aimed at reducing this cost, new concepts using gradient matching have been proposed, which bypass the need for numerical integration. This paper presents a recently established adaptive gradient matching approach, using Gaussian processes (GPs), combined with a parallel tempering scheme, and conducts a comparative evaluation with current state-of-the-art methods used for parameter inference in ODEs. Among these contemporary methods is a technique based on reproducing kernel Hilbert spaces (RKHS). This has previously shown promising results for parameter estimation, but under lax experimental settings. We look at a range of scenarios to test the robustness of this method. We also change the approach of inferring the penalty parameter from AIC to cross validation to improve the stability of the method.

Methods

Methodology for the recently proposed adaptive gradient matching method using GPs, upon which we build our new method, is provided. Details of a competing method using RKHS are also described here.

Results

We conduct a comparative analysis for the methods described in this paper, using two benchmark ODE systems. The analyses are repeated under different experimental settings, to observe the sensitivity of the techniques.

Conclusions

Our study reveals that for known noise variance, our proposed method based on GPs and parallel tempering achieves overall the best performance. When the noise variance is unknown, the RKHS method proves to be more robust.

     

Computing flow pipe of embedded hybrid systems using deep group preserving scheme

In this paper, we propose a novel methodology of numerical approximation to analyze flow of a nonlinear embedded hybrid system. For proving that all trajectories of a hybrid system do not enter an unsafe region, many classic numerical approaches such as Euler, Runge–Kutta methods for ordinary differential equations (ODEs) are applied, whereas, there exist several defects, including so-called spurious solutions and ghost fixed points. Moreover, to approximate the proper solution as much as possible, step size selection becomes especially important. In comparison, integrating group preserving scheme (GPS) which calculates true circumstance getting rid of spurious solutions and ghost fixed points, with neural network model which reduces numerical errors, deep GPS (DGPS) eliminates aforementioned adverse factors and gains better numerical approximation using a large time step size. The experimental results show that the proposed method makes safety verification for an embedded hybrid system well.