A CellML simulation compiler and code generator using ODE solving schemes

  • Florencio Rusty Punzalan1Email author,

    Affiliated with

    • Yoshiharu Yamashita2,

      Affiliated with

      • Naoki Soejima1,

        Affiliated with

        • Masanari Kawabata1,

          Affiliated with

          • Takao Shimayoshi3,

            Affiliated with

            • Hiroaki Kuwabara2,

              Affiliated with

              • Yoshitoshi Kunieda2 and

                Affiliated with

                • Akira Amano1

                  Affiliated with

                  Source Code for Biology and Medicine20127:11

                  DOI: 10.1186/1751-0473-7-11

                  Received: 28 March 2012

                  Accepted: 31 July 2012

                  Published: 19 October 2012

                  Abstract

                  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.

                  Introduction

                  In recent years, the continued development in computer processing power paved the way for the increased use of biological function simulation. Computers have proven to be invaluable in analysing complex and nonintuitive biological models and biologists are turning to them to complement their experiments. Simulations enable the testing of experimentally unfeasible scenarios and can potentially reduce experimental costs. However, the number and complexity of physiological models has also grown with the increase in computing performance. This creates challenges in reproducing simulated behaviours of the published models and reuse of models by other researchers, hindering the dissemination of science and knowledge integration.

                  One way to address model complexity is to use markup language-based model descriptions. Some popular examples include CellML [1], SBML (Systems Biology Markup Language) [2] and insilicoML [3]. CellML is an open standard markup language capable of describing mathematical models of cellular functions. SBML is an open interchange machine-readable format for representing models of functions such as metabolism and cell signalling. Meanwhile, insilicoML describes mathematical models for biophysical objects and incorporates morphological information such as shape, angle and position. SED-ML (Simulation Experiment Description Language) [4] is another type of description language which can encode the information of simulation experiments. These markup languages allow researchers to take advantage of the vast amount of biological function models using a common set of easily readable and versatile description rules.

                  Biological and physiological function models are generally described by differential equations. A typical simulation of biological function models consists of three parts: a model equation, a boundary condition, and an ordinary differential equation (ODE) solver. Model equations and boundary conditions can be described using CellML, while ODE numerical solutions like Euler and Runge-Kutta methods are typically built into the simulation software. However, it is necessary to be flexible in using ODE schemes in order to strike a balance between computational stability and speed. In addition, those using special hardware environments such as massively parallel computer systems require dedicated proprietary software to support their numerical solution needs. Thus, description languages like CellML and dedicated simulation software are not suitable or practical for flexibly incorporating different ODE solving schemes.

                  To address the need for more flexibility in creating simulation software, we created Time Evolution Calculation Markup Language (TecML), a machine-readable format for encoding ODE numerical solutions. TecML is a description language based on the extensible markup language (XML). This description language is designed to specify and store the numerical methods that can be used for solving the ODEs in biological models. It also allows the assignment of boundary conditions into the simulation experiments. The following sections describe TecML and how it is integrated into the proposed code generation system, which automatically generates codes for biological simulations.

                  The target of this study is limited to the use of different ODE numerical solutions and their application to models described in CellML. We propose an algorithm that allows users to change the ODE solution and boundary conditions of the model according to the computational needs of their simulation. To verify the effectiveness of the proposed system, we generate executable codes for several CellML models using a number of ODE numerical solutions. The system can generate code in several programming languages and code that runs in both sequential and parallel computing environments. Simulations on GPU (Graphics Processing Units) were undertaken to show the effect of using parallel computing on processing time.

                  Biological simulation code generation system

                  Summary of simulation code generation system

                  The proposed method is composed of two stages (Figure 1). In the first stage, the system represents the biological model by incorporating an ODE numerical solution method into the model’s differential equations. This creates the equation sets that calculate the time evolution of the mathematical model. The second stage generates the simulation code for these sets of equations, allowing the user to run computer simulations of the model in machines with general-purpose compilers. This approach enables the flexible construction of code generation modules that can support complex sets of equations. By generating the code separately for each section of the mathematical model, parallel code execution can be easily integrated into the simulation.
                  http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Fig1_HTML.jpg
                  Figure 1

                  Code generation system inputs and outputs. Input and output of the Biological Simulation Code Generation System. The inputs are composed of a CellML, TecML, and RelML file. The system can generate simulation codes in C/C++, Java, and Cuda C programming language.

                  To illustrate the capabilities of the proposed system, we used the FitzHugh-Nagumo (FHN) excitable media model [5] as an example. The model, proposed by R. FitzHugh in 1961 and which J. Nagumo et al. created an equivalent circuit, is a simplification of the Hodgkin-Huxley model and can be described by the following equations:
                  r = x 3 , http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Equ1_HTML.gif
                  (1)
                  dx / dt = x r / 3 . 0 y + a , http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Equ2_HTML.gif
                  (2)
                  dy / dt = b ( x + c d · y ) , http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Equ3_HTML.gif
                  (3)

                  where r, x and y are variables, t indicates time and a, b, c, d are constants. Removable algebraic expressions such as equation (1) are often used in biological models to improve the model’s readability.

                  As for the ODE numerical solution, we used the Modified Euler method to solve the ODEs in (2) and (3). By combining a mathematical description of the model above with the ODE solution method and boundary conditions, we can derive the formula to calculate the time evolution of the independent variables x and y. This involves finding the solution to the variable vector ξ using its derivatives, temporal variable vector ι, and initial conditions at time t given by
                  d ξ / dt = f ( t , ξ , ι ) , http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Equ4_HTML.gif
                  (4)
                  ι = g ( t , ξ ) . http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Equ5_HTML.gif
                  (5)
                  The rest of this section shows how the proposed system incorporates the Modified Euler method into the model equations to calculate for ξ. The following set of equations represent the numerical solution with boundary conditions:
                  ξ 0 = ξ t , http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Equ6_HTML.gif
                  (6)
                  t 0 = t , http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Equ7_HTML.gif
                  (7)
                  ι 0 = g ( t 0 , ξ 0 ) , http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Equ8_HTML.gif
                  (8)
                  κ 1 = f ( t 0 , ξ 0 , ι 0 ) , http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Equ9_HTML.gif
                  (9)
                  ξ 1 = ξ 0 + κ 1 · δ , http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Equ10_HTML.gif
                  (10)
                  t 1 = t 0 + δ , http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Equ11_HTML.gif
                  (11)
                  ι 1 = g ( t 1 , ξ 1 ) , http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Equ12_HTML.gif
                  (12)
                  κ 2 = f ( t 1 , ξ 1 , ι 1 ) , http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Equ13_HTML.gif
                  (13)
                  ξ 2 = ξ 0 + 1 2 ( κ 1 + κ 2 ) · δ , http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Equ14_HTML.gif
                  (14)
                  ξ t + δ = ξ 2 , http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Equ15_HTML.gif
                  (15)

                  where ξ,κ and ι are the vectors representing the variables and differential equations in the biological model.

                  By applying the ODE solver described above, we can arrive at the set of equations detailing the numerical solution for the initial value problem of the FHN model. First, from equations (6) and (7), the current values of the differential variables at t are assigned as the initial values with
                  x 0 = x t , http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Equ16_HTML.gif
                  (16)
                  y 0 = y t , http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Equ17_HTML.gif
                  (17)
                  t 0 = t. http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Equ18_HTML.gif
                  (18)
                  Then with equations (8) and (9), the differential and nondifferential equations are expanded and evaluated using these initial values to obtain
                  r 0 = g 1 ( t 0 , [ x 0 y 0 ] T ) = x 0 3 , http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Equ19_HTML.gif
                  (19)
                  κ 1 , x = f 1 ( t 0 , [ x 0 y 0 ] T , r 0 ) = x 0 r 0 / 3 . 0 y 0 + a , http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Equ20_HTML.gif
                  (20)
                  κ 1 , y = f 2 ( t 0 , [ x 0 y 0 ] T , r 0 ) = b ( x 0 + c d · y 0 ) . http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Equ21_HTML.gif
                  (21)
                  Next, the new values of x and y are computed as functions of κ1 and δ (equations (10) and (11)) and given by
                  x 1 = x 0 + κ 1 , x · δ , http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Equ22_HTML.gif
                  (22)
                  y 1 = y 0 + κ 1 , y · δ , http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Equ23_HTML.gif
                  (23)
                  t 1 = t 0 + δ. http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Equ24_HTML.gif
                  (24)
                  This process is repeated for r1 and κ2 as shown in equations (12) and (13);
                  r 1 = g 1 ( t 1 , [ x 1 y 1 ] T ) = x 1 3 , http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Equ25_HTML.gif
                  (25)
                  κ 2 , x = f 1 ( t 1 , [ x 1 y 1 ] T , r 1 ) = x 1 r 1 / 3 . 0 y 1 + a , http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Equ26_HTML.gif
                  (26)
                  κ 2 , y = f 2 ( t 1 , [ x 1 y 1 ] t , r 1 ) = b ( x 1 + c d · y 1 ) . http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Equ27_HTML.gif
                  (27)
                  Finally, the value of the differential variables are obtained by advancing the solution from time t to t + δ (equations (14) and (15)) with
                  x 2 = x 0 + 1 2 ( κ 1 , x + κ 2 , x ) · δ , http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Equ28_HTML.gif
                  (28)
                  y 2 = y 0 + 1 2 ( κ 1 , y + κ 2 , y ) · δ , http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Equ29_HTML.gif
                  (29)
                  x t + δ = x 2 , http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Equ30_HTML.gif
                  (30)
                  y t + δ = y 2 . http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Equ31_HTML.gif
                  (31)

                  The variable r of the equation set corresponds to ι while x and y are the variables in the vector ξ. The derivatives dx/dt and dy/dt are expressed by κ in equations (20), (21), (26) and (27).

                  In order to automatically incorporate the ODE numerical solution into the CellML model, we introduced an ODE solving scheme description language called Time Evolution Calculation Markup Language (TecML). TecML is an XML-based description language that can be used to configure ODE solutions and implement the algorithm described in equations (6) – (15). Another input of the system is the Relation Markup Language (RelML) file. The basic role of a RelML file is to relate the variables and variable types of the CellML file into their equivalent in the TecML file.

                  The simulation code generation system uses these inputs to generate the set of equations describing the time evolution of the variables in biological models. This step creates equations (16) – (31) in the FHN example. Once the model equations are applied with the ODE numerical solution method, the next stage involves the generation of the executable code that will do the actual calculations.

                  Description language

                  The code generation system requires three inputs, namely, the CellML file, TecML file and RelML file. This section gives a short introduction of each and how it is used in the algorithm and consequently, in the FHN model example.

                  CellML model encoding standard

                  The CellML [1] language is an open standard based on XML for describing mathematical models. It was designed to allow scientists to share models even if they are using different model-building software. The majority of the models stored in its repository are cell representations and these include information about the cell structure, equations for underlying processes and in some cases, boundary conditions. Figure 2 shows how a CellML model describes the mathematical equations of the FHN cell model. The file lists the input and output variables as well as their initial values and underlying equations. Note that compact syntax representation is used to show its contents and is based on Alan Garny’s notation for CellML in the Cellular Open Resource (COR) software [6].
                  http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Fig2_HTML.jpg
                  Figure 2

                  Information in the CellML file of the FHN model (represented in Alan Garny’s COR notation).

                  TecML: an ODE solving scheme description language

                  TecML (Time Evolution Calculation Markup Language) is an XML-based description language designed to describe ODE numerical solutions that can be used in biosimulations. This proposed standard allows the mathematical description of solving schemes like the Euler and Runge-Kutta methods. It allows integration of numerical methods with other description languages like CellML. TecML categorizes variables into six different types (Table 1) namely, diffvar, derivativevar, arithvar, constvar, timevar, and deltatimevar. The variables determined by a rate of change with respect to time (ξ) are referred to as differential variables (diffvar) while their derivatives (κ) are called derivative variables (derivativevar). The removable variables (ι) are the arithmetic variables (arithvar) and variables that do not change in value (ζ) are the constants (constvar). In addition, the time (t) and time increment (δ) are referred to as timevar and deltatimevar, respectively. TecML also divides the mathematical equations into two types; namely, differential (f()) and non-differential (g()) equations. Equations of type (diffequ) are the derivatives of a function while (nondiffequ) are the arithmetic functions. The information and example of a TecML file for the Modified Euler method are shown in Figure 3 and Figure 4, respectively.
                  http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Fig3_HTML.jpg
                  Figure 3

                  Information in the TecML file of the Modified Euler method where ξ is the diffvar, κ is the derivativevar, ζ is the constvar, and variable type ι is the arithvar.

                  http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Fig4_HTML.jpg
                  Figure 4

                  TecML file example.

                  Table 1

                  Variable and function types in TecML

                  Type

                  Definition

                  diffvar

                  Differential variable (variable is a function of time)

                  arithvar

                  Temporal variable (can be substituted with math equations)

                  derivativevar

                  Derivative of diffvar

                  constvar

                  Constants

                  timevar

                  Time variable

                  deltatimevar

                  Variable denoting the change of time per step

                  diffequ

                  Differential equation

                  nondiffequ

                  Non-differential equation

                  Relation Markup Language (RelML)

                  RelML is a language for describing the correspondence between the variables in the CellML model file and the variable types in an ODE numerical solution scheme described in a TecML file. Figure 5 shows the correspondence of variables described in the FHN CellML file (Figure 2) and the variable types in the TecML file (Figure 3). For the FHN model in Figure 2, variables x and y are defined as diffvar and their respective derivatives dx/dt and dy/dt as derivativevar. In addition, a, b, c, and d are constvar type while r is considered a temporal variable or arithvar type. Equations (2) and (3) are the formulas for calculating the diffvar so the functions are categorized as diffequ while the arithmetic equation for r in equation (1) is a nondiffequ.
                  http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Fig5_HTML.jpg
                  Figure 5

                  Information contained in the RelML file for the FHN model and Modified Euler method.

                  An example RelML file where the Modified Euler method is applied to the FHN model is shown in Figure 6. The first two statements after the header indicate the filename and location of the corresponding CellML and TecML file. The succeeding lines enumerate all the variables used in the FHN model and their corresponding types. The complete RelML and TecML files used in this paper have been published on the Web [7].
                  http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Fig6_HTML.jpg
                  Figure 6

                  RelML file example for FHN model and Modified Euler method.

                  Executable code generation for equation sets

                  The system implementation of the code generator allows the creation of codes in a number of programming languages. Table 2 lists the tested inputs and outputs of the system. For a desktop CPU environment, the generator can produce source codes in C, Java, and Java with BigDecimal library, which can handle higher precision computing. Codes can be generated for both single cell and cell array simulation. The cell array code generator can produce both 1D and 2D excitation propagation codes. The generated codes for 2D simulation describes the excitation propagation in a rectangular array with N×M number of cells. Aside from generating single CPU codes, it can also create codes suited for parallel computing that runs on a GPU machine.
                  Table 2

                  Input and output code types in the system

                  Cell Physiology

                  FitzHugh-Nagumo model [5]

                  Model

                  LuoRudy 1991 model [8], Kyoto 2003 Model [9]

                  ODE Solution

                  Euler Method, Modified Euler Method

                  Scheme

                  1 st -, 2 nd -, 4 th -order Runge-Kutta Method

                  Generated codes

                  C code, Java, Java with BigDecimal

                   

                  Cuda code (GPU)

                  Methods: simulation code generation algorithm

                  Cell model description

                  The cell model is a collection of variables and equation and can be written as
                  k = d x / dt = f ( x , y , t , z ) y = g ( x , y , t , z ) , http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Equ32_HTML.gif
                  (32)
                  where x is a differential variables vector, k is a derivatives vector, y is a temporal variables vector, and z is a constants vector while t is the scalar variable for time. These variable vectors can be expressed as sets with
                  x = [ x 1 , x 2 , , x N x ] T , http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Equ33_HTML.gif
                  (33)
                  k = [ k 1 , k 2 , , k N x ] T , http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Equ34_HTML.gif
                  (34)
                  y = [ y 1 , y 2 , , y N y ] T , http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Equ35_HTML.gif
                  (35)
                  z = [ z 1 , z 2 , , z N z ] T , http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Equ36_HTML.gif
                  (36)
                  where N x ,N y ,N z are the variable count for k,y, and z, respectively. Furthermore, f(x,y,t,z) and g(x,y,t,z) are function vectors with k i  = f i (x,y,t,z) and y i  = g i (x,y,t,z). These function vectors are defined by
                  f ( x , y , t , z ) = [ f 1 ( x , y , t , z ) , , f N x ( x , y , t , z ) ] T , http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Equ37_HTML.gif
                  (37)
                  g ( x , y , t , z ) = [ g 1 ( x , y , t , z ) , , g N y ( x , y , t , z ) ] T . http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Equ38_HTML.gif
                  (38)

                  ODE solving scheme

                  Table 3 shows the elements inside a TecML file. Note that all the TecML variables are denoted with Greek letters. The differential equations for the cell model are represented by the terms d ξ/dt = Φ(ξ,ι,t,ζ) and ι = Γ(ξ,ι,t,ζ). Inside a TecML file, the dependence between the differential variables ξ 0 (time t) and ξ N ξ http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_IEq1_HTML.gif (time t + δ) is given by
                  ξ i = σ i ( Ξ , π , δ ) ( 0 i N ξ ) , http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Equ39_HTML.gif
                  (39)
                  κ i + 1 = Φ ( ξ i , ι i , τ i , ζ ) ( 0 i N ξ ) , http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Equ40_HTML.gif
                  (40)
                  ι i = Γ ( ξ i , ι i , τ i , ζ ) ( 0 i N ξ ) , http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Equ41_HTML.gif
                  (41)
                  τ i = T i ( t , δ ) ( 0 i N ξ ) , http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Equ42_HTML.gif
                  (42)
                  Table 3

                  Information written in a TecML file

                  Notation

                  Definition

                  ξ 0 = [ ξ 0 , 1 , ξ 0 , 2 , , ξ 0 , N x ] T http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_IEq2_HTML.gif

                  Input differential variable vector

                  ξ N ξ = [ ξ N ξ , 1 , ξ N ξ , 2 , , ξ N ξ , N x ] T http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_IEq3_HTML.gif

                  Output differential variable vector

                  ξ i (1 ≤ i ≤ N ξ )

                  Time differential variable vector

                  Ξ = [ ξ 0 , ξ 1 , , ξ N ξ ] T http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_IEq4_HTML.gif

                  Differential variable vector of ξ

                  t

                  Current time value

                  δ

                  Time step

                  ζ = [ ζ 0 , ζ 1 , , ζ N z ] T http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_IEq5_HTML.gif

                  Constants vector

                  τ = T i (t,δ)

                  Calculation of variable T i

                  Φ ( ξ , ι , τ , ζ ) = [ ϕ 1 ( ξ , ι , τ , ζ ) , , ϕ N x ( ξ , ι , t , ζ ) ] T http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_IEq6_HTML.gif

                  Differential equation vector

                  Γ ( γ , ι , τ , ζ ) = [ γ 1 ( ξ , ι , τ , ζ ) , , γ N y ( ξ , ι , t , ζ ) ] T http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_IEq7_HTML.gif

                  Temporal function vector

                  κi + 1 = Φ(ξ,ι,t,ζ)

                  Derivative variable vector

                  ι i  = Γ(ξ,ι,t,ζ)

                  Temporal variable derived from ξ

                  π = [ κ 1 , κ 2 , , κ N ξ ] T http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_IEq8_HTML.gif

                  Derivative vector of κ

                  ξ i  = σ ( Ξ , K , δ ) http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_IEq9_HTML.gif

                  Relation between ξ i and Ξ , K , δ http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_IEq10_HTML.gif

                  where the differential variable vector Ξ and derivative variable vector π are given by
                  Ξ = [ ξ 0 , ξ 1 , , ξ N ξ ] T , http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Equ43_HTML.gif
                  (43)
                  π = [ κ 1 , κ 2 , , κ N ξ ] T . http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Equ44_HTML.gif
                  (44)

                  CellML and TecML integration

                  The integration of the CellML model and TecML solving scheme involves the mapping of corresponding variables. The mapping shows how each physiological model variable written in CellML is replaced with its corresponding TecML variable. The differential variable vector ξ of TecML corresponds to the differential variable vector x of CellML and reads as
                  ξ k x k = [ x k , 1 , x k , 2 x k , N x ] T . http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Equ45_HTML.gif
                  (45)
                  TecML’s derivative variable κ equates to the CellML’s derivative variable d x/dt and the temporal variable ι corresponds to y and given by
                  κ k + 1 d x k / dt = d x k , 1 dt , d x k , 2 dt , , d x k , N x dt T , http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Equ46_HTML.gif
                  (46)
                  ι k y k = [ y k , 1 , y k , 2 , , y k , N y ] T . http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Equ47_HTML.gif
                  (47)
                  Furthermore, the TecML equations κ = Φ(ξ,ι,t,ζ) and ι = Γ(ξ,ι,t,ζ) correspond to the CellML equations f and g, respectively, as shown by
                  Φ ( ξ , ι , τ , ζ ) f ( x , y , t , z ) , http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Equ48_HTML.gif
                  (48)
                  Γ ( ξ , ι , τ , ζ ) g ( x , y , t , z ) . http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Equ49_HTML.gif
                  (49)

                  Replacement algorithm

                  Once the variables are mapped, the differential and arithmetic equations are transformed and expanded according to the numerical method described in the TecML file. First, the variables in the CellML model equations are replaced with their corresponding TecML variables. Each model equation is searched for any diffvar, derivativevar, arithvar, constvar, timevar, or deltatimevar and that variable is replaced with the corresponding TecML variable name. This replacement procedure is represented in Strachey brackets and expressed by the following function:
                  R v equ = R x R k R y R t R z R d equ , http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Equ50_HTML.gif
                  (50)
                  where ⟦equ⟧ is the TecML equation and, for all i = 1,2,…,N ξ , the replacement functions are
                  R x ξ i = x i , http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Equ51_HTML.gif
                  (51)
                  R k κ i = k i , http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Equ52_HTML.gif
                  (52)
                  R y ι i = y i , http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Equ53_HTML.gif
                  (53)
                  R t τ i = t i , http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Equ54_HTML.gif
                  (54)
                  R z ζ = z , http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Equ55_HTML.gif
                  (55)
                  R d ⟦δ⟧ = ⟦d⟧. http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Equ56_HTML.gif
                  (56)

                  Note that the Strachey brackets in R x ξ i = x i http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_IEq11_HTML.gifmeans that all the TecML diffvar (i.e. for i = 1…N ξ ) in the argument is replaced with the corresponding CellML diffvar x i . This function is true for all the other variable types in TecML.

                  The generated set of equations advances the solution of the differential variables from time t to t + δ. Figure 7 shows how the algorithm generates the simulation program from the input files. Note that the subroutine replace_v() in the algorithm replaces variables in the CellML model equation with the TecML variables as expressed in equation (50). The function replace_sj() generates one scalar equation from a vector equation and appends index j to the variables. In the equations which do not contain functions f or g, replace_d() generates multiple scalar equations from a vector equation. Finally, replace_f() and replace_g() unfold functions f and g, respectively. These subroutines can be expressed in the following transformations:
                  R s , j ⟦equ⟧ = makeEquation S y , j S k , j getLHS ⟦equ⟧⟧⟧ , S g , j S f , j getRHS ⟦equ⟧⟧⟧⟧ , http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Equ57_HTML.gif
                  (57)
                  http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Fig7_HTML.jpg
                  Figure 7

                  Algorithm for generating the set of ODE numerical solution equations.

                  R d , j equ = S k , j S x , j equ , http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Equ58_HTML.gif
                  (58)
                  R f ϕ j ( ξ , ι , τ , ζ ) = T y , l T x , k f j ( x , y , t , z ) ⟧⟧ , http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Equ59_HTML.gif
                  (59)
                  R g γ j ( ξ , ι , τ , ζ ) = T y , l T x , k g j ( x , y , t , z ) ⟧⟧ , http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Equ60_HTML.gif
                  (60)
                  where the index-appending functions S http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_IEq12_HTML.gifand unfolding functions T http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_IEq13_HTML.gifare defined as
                  S x , j x i = x i , j , http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Equ61_HTML.gif
                  (61)
                  S y , j y i = y i , j , http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Equ62_HTML.gif
                  (62)
                  S f , j Φ = ϕ j , http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Equ63_HTML.gif
                  (63)
                  S g , j Γ = γ j , http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Equ64_HTML.gif
                  (64)
                  T x , k x = x k , http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Equ65_HTML.gif
                  (65)
                  T y , k y = y k ⟧. http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Equ66_HTML.gif
                  (66)

                  Experiments and results

                  Code generation for the cell model simulation

                  The FHN cell model [5] was used to evaluate the proposed system. The system was also used to generate simulation codes for a number of cell models (LuoRudy1991 [8], LuoRudy1994 [10] and KyotoModel2003 [9]) and with other more accurate ODE numerical methods (e.g. 4th-order Runge-Kutta method). All the generated code used in the simulation experiments can be downloaded from the files section of the program generator site [7].

                  The traditional approach to creating simulation experiment codes also offers little flexibility once the software is created. Changing the numerical ODE solution method means making major revisions in the simulation code. The proposed method allows the flexibility of using different simulation models, boundary conditions, and numerical methods for a simulation experiment. Since the simulation codes are generated automatically, users can choose or change their desired ODE solver without making changes in the simulation codes themselves. Also, this can give clear information on what is calculated to generate the simulation results.

                  The different cell models and ODE numerical methods produced varying simulation code sizes. Table 4 lists the number of execution steps generated for the cell models using different ODE solutions. The table shows that the more complex the model becomes and the more equations it has, the larger the number of steps to compute the model equations in the code.
                  Table 4

                  The number of execution steps for the generated codes for different cell models and ODE numerical methods

                  Cell Model

                  Euler

                  Modified Euler

                  Runge-Kutta

                  FHN

                  11

                  16

                  26

                  LuoRudy1991

                  70

                  117

                  211

                  LuoRudy1994

                  123

                  211

                  387

                  Kyoto Model

                  335

                  560

                  1200

                  An issue in approximating ODE solutions is the accuracy of the numerical method used. A good approximation to the underlying differential equation needs to be achieved in order to arrive at accurate simulation results. We tested a number of commonly used ODE numerical methods to determine how the use of different solutions affect the accuracy of the calculations. The three methods used were Euler, Modified Euler and 4th-order Runge-Kutta. Each of these methods were used to generate an FHN model simulation code that runs in a single CPU. The generated code can run in different compilers and does not require third party software.

                  In order to test the accuracy of the ODE numerical methods, simulation codes were generated in Java using the BigDecimal class and numeric formatting. The Java BigDecimal can represent a large number of decimal places and help avoid rounding errors. It can offer higher precisions than the 16 decimal digits offered by floating point double. In the simulation codes, we used BigDecimal to represent the numbers in a 32-digit decimal point precision format for all the ODE solving methods. Different time steps were also used in testing the accuracy of these ODE methods, ranging from 10−1 ms to 10−5 ms. The simulation using a time step of 10−6 and Runge-Kutta as ODE solving scheme was used as the basis to compute for the root-mean-square error (RMSE) and evaluate the accuracy of the other calculations. The RMSE was computed for all the calculations using different ODE numerical methods and in varying time steps (Figure 8). Each level in the RMSE indicates a 1/10 less accuracy in the simulation results compared to the calculations using Runge-Kutta and 10−6 ms time step.
                  http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Fig8_HTML.jpg
                  Figure 8

                  Relationship between the time step and calculation error for each ODE numerical method used in generating the FHN model. The results and root-mean-square errors (RMSE) are computed using Java BigDecimal. The result of the Runge-Kutta method with 10-6 ms time step is used as reference for the RMSE computations.

                  The first-order Euler method gave the largest error (log error ≈10−3) while the Modified Euler resulted to a smaller error. The fourth-order Runge-Kutta method resulted with the best accuracy (log error ≈10−10). However, the Runge-Kutta error is almost constant from the time step 10−3 ms, . This can be attributed to the rounding error of digits over the used 32-digit precision.

                  Simulation codes using different ODE numerical methods were also generated for models more complicated than the FHN or Hodgkin-Huxley model. For this, we used the model introduced by C. Luo and Y. Rudy in 1991 [8]. It is a simple cell model of cardiac action potential that uses Hodgkin-Huxley type equations to calculate ionic currents. The RMSE computations undertaken for the FHN model were also applied for the Luo-Rudy 1991 simulations. Codes were generated for the three ODE solutions with time steps ranging from 10−1 to 10−4 ms. Meanwhile, the simulation using the Euler method with a time step of 10−5 ms was used as the reference for RMSE calculations (Figure 9).

                  Excitation propagation simulations

                  One possible application of 1- or 2-dimensional excitation propagation simulation is cardiac arrhythmia research. Computational models have provided new insights into the underlying mechanisms of re-entry like the role of ionic currents and ion channel mutations [11]. In addition, designs of defibrillation treatments can be optimized using simulations of excitation propagation to achieve good clinical results for patients.

                  Simulation codes for the Luo-Rudy 1991 model were generated to simulate cell excitation propagation on a 2D homogeneous sheet. The two-dimensional cell action potential propagation simulations were run for both CPU and GPU to compare the speed of calculations. The CPU simulation was performed with an intel Core i7 880 processor with 8 GB of memory running Windows 7. The GPU has a single Nvidia Tesla C2050 processor with 448 CUDA cores, and 3 GB memory in a CentOS 5.5 system. The programs are written in C for the CPU simulations and Cuda C for the GPU.

                  In the experiments, the size of the homogeneous sheet was varied from 10240 (10240 × 1) to 1433600 (10240 × 14) cells with 10240-cell increments. The computation time was measured for each increment. Figure 9 shows the CPU and GPU computational time at different cell array sizes. Results showed that by using the GPU, the computational time can be accelerated 50 times as compared to using a standard CPU.
                  http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Fig9_HTML.jpg
                  Figure 9

                  Relationship between the time step and calculation error for each ODE numerical method used in generating the LuoRudy 1991 model.

                  Discussion

                  Traditionally, biological function simulation programs have to be manually created from scratch. This is manageable with small models or simulation experiments but it becomes less practical once model complexity increases. The more complex the model becomes the larger the program becomes and the harder it is to create and maintain such programs. The Doty model [12] showed that the cost of software development is proportional to the square exponential of the number of computational steps. As seen from Table 4, the number of computational steps increases with the complexity of the model. By using this system, the creation of a simulation program from biological models becomes simpler. The system takes input from experts in biological functions through CellML models and mathematical experts through ODE numerical solutions and automatically generates the simulation code. The automatic code generation allows it to deal with large and complex models without the proportional increase in software development cost. This also keeps the software maintenance cost low since the programs are created with the same structure, regardless of the biological model under study.

                  The system’s support for multiple programming languages makes it easier for users to test the model in different programming environments and keep the cost of switching from one programming language to a low. Based on the codes generated for the purpose of this paper, science and engineering students typically create a Java simulation code from the C version within four days.

                  The system requires the formatted RelML file to do the first stage of the algorithm and the RelML information can also be used to automatically determine the boundary conditions. The RelML information can be extracted automatically from the CellML model file but boundary conditions are not always available. We need a mathematical representation of boundary conditions which are compatible with declarative description. This can be addressed by the inclusion of PEPML experimental protocol [13] in the system. PEPML is an XML-based language that can describe the initial conditions and procedures of an experimental protocol. It describes boundary conditions in a mathematical form and is purely represented in a declarative manner, making it easy to combine with the current implementation of our system.

                  The attainable order of accuracy for approximating the Luo-Rudy 1991 model equations was determined by running simulations of the model using different ODE numerical methods. The results of the comparisons between ODE solving schemes showed that the model can be accurately predicted by the first-order ODE solutions. The use of higher-order numerical solutions has a minimal effect on the values of the approximations (Figure 10). The Luo-Rudy model is a set of nonlinear differential equations that use higher order ODE solutions but it can be seen from the Taylor series expansion that equations from the second-order term has little or negligible effect. The structure of the system, where the ODE solver can be easily changed with the appropriate TecML file, allows one to confirm the accuracy of the model simulation such as the case with the Luo-Rudy 1991 model (see Additional file 1).
                  http://static-content.springer.com/image/art%3A10.1186%2F1751-0473-7-11/MediaObjects/13029_2012_Article_86_Fig10_HTML.jpg
                  Figure 10

                  CPU and GPU execution time ratio of the excitation propagation simulation.

                  Although the system is fully implemented, further improvements are still to be added. The generated excitation propagation code does not consider the optimization of the model variables for parallel processing. Therefore, code parallelization is low and may have contributed to slower execution in the GPU environment. Also, the system can only handle equations with a single term in the left-hand side. Other equation forms will be addressed in future implementations.

                  Furthermore, simultaneous calculations using differential algebraic equations (DAE) are still not available. Simultaneous equations appear in some models and implicit numerical methods. The handling of implicit methods is available in the inputs but not in the current implementation of our system. A simultaneous equation library, other ODE numerical solutions library like CVODE, or a linear algebra library like LINPACK could be incorporated into our system. This can be a possibility with our plan to design a library interface description language to handle numerical solution inputs.

                  Future implementations also need to address the need for a declarative design of describing procedural methods. Complex ODE numerical methods like the ones described in the Kinetic Simulation Algorithm Ontology (KiSAO) [14] are difficult or impossible to encode in the current version of TecML. We are in the process of redesigning TecML and the system to incorporate methods that are suitable for the types of problems scientists are addressing with manually-created code, including adaptive methods. This can greatly enhance the usability of the proposed method.

                  In addition, the system may be integrated with SED-ML. SED-ML allows description of the simulation environment, which includes the name of the numerical method to be used in the simulation. A starting point for integration would be to provide TecML information of a numerical method for the KiSAO entry in SED-ML to allow simulation software to refer to this description. Another point would be to create support for SED-ML files in the code generation system in order to create more complex experiments that use more than one model or experiments with different simulation methods applied. This will allow users the flexibility of choosing not only the ODE solving methods but also the experiment protocol when generating their simulation codes.

                  Conclusion

                  In this paper, we proposed a method to automatically generate executable simulation codes using CellML physiological models and ODE numerical solution methods. The generated code describes the time-evolution of the set of differential equations enumerated by the CellML model. The code generation system is composed of a two-stage approach that allows flexible generation of complex sets of equations.

                  To evaluate the effectiveness of the proposed system, several combinations of physiological cell models and ODE solving schemes were generated. The output of the numerical approximations were in accordance with the published results of the cell models. Results for the Luo-Rudy 1991 simulations also indicated that it only has first-order accuracy. The comparison of execution time for 2D excitation propagation also showed that the use of GPU can accelerate the processing time by 50 times as compared to a CPU.

                  The code generation system allows executable simulation codes to be easily generated from CellML and TecML files. This can be very useful in the field of biological model simulation since it provides the tools to quantitatively evaluate the mathematical equations in these models.

                  Author’s contributions

                  FRP has been involved in the drafting of the manuscript, programming and gathering of experimental data. YY was responsible for programming of the code generator software and gathering data. NS was the main programmer of the code generator. MK was involved in the simulation experiments and data gathering as well as programming. TS, HK and YK acted as consultants for the algorithm and system design. AA has been the main architect of the algorithms and system design and has given the final approval of the version to be published. All authors read and approved the final manuscript.

                  Declarations

                  Acknowledgements

                  The authors would like to thank the Japan Ministry of Education, Culture, Sports, Science and Technology (MEXT) for the Grant-in-Aid for Scientific Research (Kakenhi) funding.

                  Authors’ Affiliations

                  (1)
                  Graduate School of Life Sciences, Ritsumeikan University
                  (2)
                  Graduate School of Information Science and Engineering, Ritsumeikan University
                  (3)
                  ASTEM Research Institute of Kyoto

                  References

                  1. Cuellar AA, Lloyd CM, Nielsen PF, Bullivant DP, Nickerson DP, Hunter PJ: An overview of CellML 1.1, a biological model description language. Simulation. 2003, 79: 740-747. 10.1177/0037549703040939.View Article
                  2. Hucka M, Kitano H: The systems biology markup language (SBML): a medium for representation and exchange of biochemical network models. Bioinformatics. 2003, 19: 524-531. 10.1093/bioinformatics/btg015.View ArticlePubMed
                  3. Yoshiyuki A, Yasuyuki S, Yoshiyuki K: Specifications of insilicoML 1.0: a multilevel biophysical model description language. J Biol Sci. 2008, 58: 447-458.
                  4. Waltemath D, Adams R, Bergmann FT, Hucka M, Kolpakov F, Miller A, Moraru II, Nickerson DP, Sahle S, Snoep JL, Le Novere: Reproducible computational biology experiments with SED-ML - the simulation experiment description markup language. BMC Syst Biol. 2011, 5: 198-10.1186/1752-0509-5-198.PubMed CentralView ArticlePubMed
                  5. Fitzhugh R: Impulses and physiological states in theoretical models of nerve membrane. Biophys J. 1961, 1 (6): 445-466. 10.1016/S0006-3495(61)86902-6.PubMed CentralView ArticlePubMed
                  6. Garny A: Cellular Open Source (COR): a public cellml based environment for modeling biological function. Int J Bifurcation and Chaos. 2003, 13 (12): 3579-3590. 10.1142/S021812740300882X.View Article
                  7. The CellML Compiler and Code Generator. [http://​sourceforge.​net/​projects/​cellmlcompiler]
                  8. Luo CH, Rudy Y: A model of the ventricular cardiac action potential, depolarization, repolarization, and their interaction. Circ Res. 1991, 68: 1501-1526. 10.1161/01.RES.68.6.1501.View ArticlePubMed
                  9. Matsuoka S, Sarai N, Kuratomi S, Ono K, Noma A: Role of individual ionic current systems in ventricular cells hypothesized by a model study. Jpn J Physiol. 2003, 53: 105-123. 10.2170/jjphysiol.53.105.View ArticlePubMed
                  10. Luo CH, Rudy Y: A dynamic model of the cardiac ventricular action potential. I. Simulations of ionic currents and concentration changes. Circ Res. 1994, 74: 1071-1097. 10.1161/01.RES.74.6.1071.View ArticlePubMed
                  11. Ten Tusscher KH, Panfilov AV: Cell model for efficient simulation of wave propagation in human ventricular tissue under normal and pathological conditions. Phys Med Biol. 2006, 51: 6141-6156. 10.1088/0031-9155/51/23/014.View ArticlePubMed
                  12. Herd JR, Postak JN, Russell WE, Steward KR: Software cost estimation study: Study results, Final Technical Report, RADC-TR77-220, vol. I. Rockville. 1977, Inc., Doty Associates
                  13. Shimayoshi T, Amano A, Matsuda T: A generic representation format of physiological experimental protocols for computer simulation using ontology. Engineering in Medicine and Biology Society, 2007. EMBS 200. 29th Annual International Conference of the IEEE. 2007, 382-385.View Article
                  14. Kinetic Simulation Algorithm Ontology (KiSAO). [http://​www.​biomodels.​net/​kisao]

                  Copyright

                  © Punzalan et al.; licensee BioMed Central Ltd. 2012

                  This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://​creativecommons.​org/​licenses/​by/​2.​0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.