A CellML simulation compiler and code generator using ODE solving schemes

Models written in description languages such as CellML are becoming a popular solution to the handling of complex cellular physiological models in biological function simulations. However, in order to fully simulate a model, boundary conditions and ordinary differential equation (ODE) solving schemes have to be combined with it. Though boundary conditions can be described in CellML, it is difficult to explicitly specify ODE solving schemes using existing tools. In this study, we define an ODE solving scheme description language-based on XML and propose a code generation system for biological function simulations. In the proposed system, biological simulation programs using various ODE solving schemes can be easily generated. We designed a two-stage approach where the system generates the equation set associating the physiological model variable values at a certain time t with values at t + Δt in the first stage. The second stage generates the simulation code for the model. This approach enables the flexible construction of code generation modules that can support complex sets of formulas. We evaluate the relationship between models and their calculation accuracies by simulating complex biological models using various ODE solving schemes. Using the FHN model simulation, results showed good qualitative and quantitative correspondence with the theoretical predictions. Results for the Luo-Rudy 1991 model showed that only first order precision was achieved. In addition, running the generated code in parallel on a GPU made it possible to speed up the calculation time by a factor of 50. The CellML Compiler source code is available for download at http://sourceforge.net/projects/cellmlcompiler

Time-dependent changes in membrane excitability during glucose-induced bursting activity in pancreatic β cells

In our companion paper, the physiological functions of pancreatic β cells were analyzed with a new β-cell model by time-based integration of a set of differential equations that describe individual reaction steps or functional components based on experimental studies. In this study, we calculate steady-state solutions of these differential equations to obtain the limit cycles (LCs) as well as the equilibrium points (EPs) to make all of the time derivatives equal to zero. The sequential transitions from quiescence to burst–interburst oscillations and then to continuous firing with an increasing glucose concentration were defined objectively by the EPs or LCs for the whole set of equations. We also demonstrated that membrane excitability changed between the extremes of a single action potential mode and a stable firing mode during one cycle of bursting rhythm. Membrane excitability was determined by the EPs or LCs of the membrane subsystem, with the slow variables fixed at each time point. Details of the mode changes were expressed as functions of slowly changing variables, such as intracellular [ATP], [Ca2+], and [Na+]. In conclusion, using our model, we could suggest quantitatively the mutual interactions among multiple membrane and cytosolic factors occurring in pancreatic β cells

Quantitative Decomposition of Dynamics of Mathematical Cell Models: Method and Application to Ventricular Myocyte Models.

Mathematical cell models are effective tools to understand cellular physiological functions precisely. For detailed analysis of model dynamics in order to investigate how much each component affects cellular behaviour, mathematical approaches are essential. This article presents a numerical analysis technique, which is applicable to any complicated cell model formulated as a system of ordinary differential equations, to quantitatively evaluate contributions of respective model components to the model dynamics in the intact situation. The present technique employs a novel mathematical index for decomposed dynamics with respect to each differential variable, along with a concept named instantaneous equilibrium point, which represents the trend of a model variable at some instant. This article also illustrates applications of the method to comprehensive myocardial cell models for analysing insights into the mechanisms of action potential generation and calcium transient. The analysis results exhibit quantitative contributions of individual channel gating mechanisms and ion exchanger activities to membrane repolarization and of calcium fluxes and buffers to raising and descending of the cytosolic calcium level. These analyses quantitatively explicate principle of the model, which leads to a better understanding of cellular dynamics

Each dynamic and total dynamics of intracellular calcium concentration, [Ca²⁺]_i of Priebe model after stimulus onset at two phases (A, B) with different scales, together with [Ca²⁺]_i (black) and the instantaneous equilibrium point (red) at the top of each figure.

<p>Each dynamic and total dynamics of intracellular calcium concentration, [Ca²⁺]_i of Priebe model after stimulus onset at two phases (A, B) with different scales, together with [Ca²⁺]_i (black) and the instantaneous equilibrium point (red) at the top of each figure.</p

Dynamics of the membrane potential <i>V</i>_<i>m</i> of Takeuchi Model after stimulus offset at different phases (A–D) with different scales, together with <i>V</i>_<i>m</i> (black) and the instantaneous equilibrium point (red) at the top of each figure on the same time base.

<p><i>I</i>_NaCa denotes the sum of <i>I</i>_NaCa p(E_1total), p(I₁) and p(I₂) dynamics of <i>V</i>_<i>m</i>. A notch of <i>V</i>_<i>m</i> dynamic around 185 ms and a tiny fluctuation of <i>I</i>_Na p(AP) dynamic around 235 ms in D are errors associated with division by nearly zero in numerical differentiation.</p

Time course of membrane potential (black) and its stable and unstable instantaneous equilibrium point (red and blue, respectively) of Priebe model.

<p>Time course of membrane potential (black) and its stable and unstable instantaneous equilibrium point (red and blue, respectively) of Priebe model.</p

Major currents and differential variables simulated using Takeuchi model.

<p>(A) <i>I</i>_Na, fast sodium current; <i>I</i>_Na p(AP), open probability of <i>I</i>_Na voltage gate. (B) <i>I</i>_K1, inward rectifier potassium current; <i>I</i>_K1<i>f</i>_O, open probability of <i>I</i>_K1 magnesium gate; <i>I</i>_K1<i>y</i>, open probability of <i>I</i>_K1 polyamine gate. (C) <i>I</i>_CaL, L-type calcium current; <i>I</i>_CaL p(AP), open probability of <i>I</i>_CaL voltage-dependent gate; <i>I</i>_CaL p(U), open probability of <i>I</i>_CaL calcium gate; <i>I</i>_CaL<i>y</i>, open probability of <i>I</i>_CaL ultra-slow gate. (D) <i>I</i>_CaT, T-type calcium current; <i>I</i>_CaT<i>y</i>₁, open probability of <i>I</i>_CaT activation gate. (E) <i>I</i>_NaCa, Na⁺/Ca²⁺ exchange current. (F) <i>I</i>_<i>l</i>(Ca), Ca²⁺-activated background cation current. (G) <i>I</i>_RyR, ryanodine receptor channel current; <i>I</i>_RyR p(Open), open probability of <i>I</i>_RyR gate. (H) <i>I</i>_SERCA, SR Ca²⁺ pump current; <i>I</i>_SERCA<i>y</i>, probability of the conformation state with the Ca²⁺-binding sites onto the SR side. (I) [Ca²⁺]_rel, calcium concentration in SR release site.</p

Derivative function of <i>V</i> of Priebe model at 257 ms and 340 ms (thick curve in A and B, respectively) assuming that all the other differential variables are constant.

<p>× is the instantaneous equilibrium point of <i>V</i>, which is the x-intercept of the tangent line (thin line) at the value of <i>V</i> at that time (+).</p