alexa Assessing Numerical Resolution Methods Performance for Kinetic Models of Receptors and Channels | Open Access Journals
ISSN: 0974-7230
Journal of Computer Science & Systems Biology
Like us on:
Make the best use of Scientific Research and information from our 700+ peer reviewed, Open Access Journals that operates with the help of 50,000+ Editorial Board Members and esteemed reviewers and 1000+ Scientific associations in Medical, Clinical, Pharmaceutical, Engineering, Technology and Management Fields.
Meet Inspiring Speakers and Experts at our 3000+ Global Conferenceseries Events with over 600+ Conferences, 1200+ Symposiums and 1200+ Workshops on
Medical, Pharma, Engineering, Science, Technology and Business

Assessing Numerical Resolution Methods Performance for Kinetic Models of Receptors and Channels

Merdan Sarmis1,2, Jean-Marie C Bouteiller1,3*, Nicolas Ambert1, Arnaud Legendre1,2, Serge Bischoff1, Olivier Haeberlé2 and Michel Baudry1,4*

1Rhenovia Pharma SA, Mulhouse, France

2Laboratoire MIPS EA2332, Université de Haute Alsace, Mulhouse, France

3Department of Biomedical Engineering, University of Southern California, Los Angeles, CA, USA

4Western University of Health Sciences, Pomona, CA, USA

*Corresponding Author:
Michel Baudry
Western University of Health Sciences
Pomona, CA, USA
Tel: (+33) 389 321 180
Fax: (+33) 389 555 145
E-mail:[email protected]
Jean-Marie C Bouteiller
Rhenovia Pharma SA, 20C Rue de Chemnitz
F-68200 Mulhouse, France
Tel: +33-389 321 180
Fax: +33-389 555 145
E-mail: jeanmarie [email protected]

Received date: June 17, 2013; Accepted date: July 19, 2013; Published date: July 22, 2013

Citation:Sarmis M, Bouteiller JMC, Ambert N, Legendre A, Bischoff S, et al. (2013) Assessing Numerical Resolution Methods Performance for Kinetic Models of Receptors and Channels. J Comput Sci Syst Biol 6:150-164. doi:10.4172/jcsb.1000112

Copyright: © 2013 Sarmis M, et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License,which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Visit for more related articles at Journal of Computer Science & Systems Biology


In systems biology, systems of kinetic reactions are generally used to model and simulate various biochemical pathways. These reactions are translated into ordinary differential equations, which are computationally resolved by numerical algorithms. Computation performance, defined by how fast the algorithm converges to a numerical solution of the system of ordinary differential equations, critically depends on the choice of the appropriate algorithm. In this paper, we compared several algorithms used to solve ordinary differential equations applied to several kinetic models that describe the dynamic behavior of receptors and ion channels found in chemical synapses of the Central Nervous System; we provide a simplified method to determine the performances of these ordinary differential equation solvers, in order to provide a benchmark for algorithm selection. This method will facilitate the choice of the most efficient algorithm for a given kinetic model with a minimum number of tests. Our results provide a tool for identifying optimal solvers for any biological bilinear kinetic models under various experimental conditions. This comparison also underscored the complexity of biological kinetic models and illustrates how their input dependency could interfere with performance. Despite these challenges, our simplified method helps to select the best solvers for any synaptic receptors kinetic models described, with a bilinear system with minimal a priori information on the solver structure and the model.


Algorithm selection; Synaptic receptor; Kinetic models; Ordinary differential equation; Bilinear systems; Java; Numerical resolution method.


The main goal of systems biology [1], consists in providing a quantitative and integrative description of living organisms by using simulation of complex and interconnected models of metabolic networks and signaling pathways. These complex processes are described by systems of biochemical reactions and quantitatively analyzed by corresponding sets of mathematical equations that relate reactant concentrations and elementary rate constants to changes in product concentrations with time. In order to analyze the temporal evolution of the systems of reactions, series of ordinary differential equations (ODEs) need to be computed. ODEs are computationally solved with ODE solver algorithms, providing numerical solutions for the temporal changes of the various variables of the biological systems. Over the years, mathematicians have developed numerous ODE algorithms. However, each ODE system or model is specific regarding to the number of states, input, stiffness and linearity, and therefore, each algorithm differs in terms of speed and accuracy. If different ODE solvers produce relatively similar results, their performances are highly variable, in term of speed to reach a stable numerical solution and accuracy of the final solution. It is therefore difficult to choose the best algorithm to solve particular sets of ODE without expertise analysis. In fact, algorithm selection is a complex problem, which depends on many criteria, as described by Rice [2].

According to Ewald [3] (Figure 1), the selection of the most appropriate algorithm to solve an ODE system requires knowing 1) the model and its inputs (problem space), 2) the model size and characteristics (feature space), 3) the ODE solver structure (algorithm space), and 4) user’s preferences (criteria space). Users usually have a specific expertise, either in biology, modeling, code development or mathematics, but rarely in all fields. There is, therefore, a real need for a tool guiding any users towards the choice of the best algorithm to solve ODE systems.


Figure 1: Algorithm selection problem. The algorithm selection problem was defined by Rice in 1976, and completed by Ewald [3]. As presented in this scheme ([3] page 22), users need to know model inputs, characteristics (size) and algorithms structure in order to select the right algorithm for a simulation.

Rhenovia is repeatedly faced with this problem throughout the development of a simulation platform for hippocampal glutamatergic synapse [4], and its integration into complex neuronal network. To this purpose, kinetic models were built and implemented. Due to the nature of kinetic models, the ODE systems are bilinear. As this project was initiated several years ago using the Java programming language, the eight solvers compared are all Java-based. We decided to simplify the algorithm selection problem by comparing the eight ODE solver algorithms performances and determine the most appropriate for bilinear synaptic kinetic models, depending on the properties of the system under consideration. This study provides a priori knowledge on solvers performances to help users interested in launching a large number of simulations with a specific model, for example to run an analysis sensibility.

In the literature, bilinear systems are written in the mathematical form represented by equation (1), where x(t)∈ℜn is the state variable vector and u(t)∈ℜn is the system’s input vector, with ui(t) the ith command of the system and i B his associated constant matrix [5,6]. The bilinear systems are a special case of non linear affine systems.

       image                                         (1)

We compared eight ODE solvers (Table 1 and Table S-1), implemented in the Java language following good coding practice [7], in order to ensure adequate performance and high robustness. The ODE solvers are categorized into three groups: implicit, explicit and hybrid solvers (hybrid solver combines implicit and explicit methods at each integration step). Our kinetic models are usually bilinear systems and can be categorized as stiff or non-stiff according to their temporal evolution under a particular experimental protocol. The notion of stiffness was introduced by Curtiss and Hirschfelder [8], and has been refined several times. Since we do not intend to give a precise definition of stiffness here, we will consider a model to be non-stiff if an explicit solver is more efficient than an implicit one; otherwise we will consider it to be stiff [9]. The stiffness detection methods are commonly based on numerical resolution method stability bound [10-12]. This implies that the stiffness of a model is closely related to the selected numerical resolution method. According to Ekeland et al. [12], the precise definition of the stiffness of a system of equations is not crucial from a practical point of view. Therefore, in order to simplify the algorithm selection problem, we decide to follow Ekeland’s idea, and not to determine the stiffness degree of models used in this study.

Explicit Implicit Hybrid

Table 1: Classification of the 8 ODE solvers used in our study.

All evaluated solvers have variable step-sizes, except for the selected reference solver. The reference solver is a fourth-order explicit Runge-Kutta (RK) algorithm [13] (RK4), and assumed to be the most precise at the chosen (very small) step size. The other solvers were Runge-Kutta Fehlberg, TR-BDF2, Rosenbrock, JVODE ADAMS, JVODE BDF and JLSODE. The Runge-Kutta Fehlberg [14,15] (RKF) is a fourth/fifth-order explicit RK scheme. TR-BDF2 [16,17], an implicit solver, is composed of a trapezoidal rule and a second-order Backward Differential Formula (BDF scheme), and corresponds to the ode23tb solver [18] in the Matlab® software (see details in the supplemental document). The Rosenbrock solver [19] is a hybrid RK4, and IMEX (for IMplicit-EXplicit) [20,21] is another hybrid fourth-order Runge- Kutta. JVODE is a Java version of the CVODE solver [22,23]. JVODE has not been modified, as it was implemented in Java language by the BioUML platform developers [24]. JVODE includes two schemes: 1) a variable order (from 1 to 12), explicit Adams-Moulton method for non-stiff systems, and 2) a variable order (from 1 to 5), implicit Backward Differential Formula (BDF scheme) for stiff systems. We will differentiate these two schemes and call them JVODE ADAMS and JVODE BDF, respectively. The last ODE solver is JLSODE (Java version of LSODE [25]), which has the same Adams-Moulton and BDF schemes as JVODE. This last solver however adds stiffness detection, as written by Uteshev and Pennefather [26]. With this stiffness detection, JLSODE starts with the Adams-Moulton scheme and reversibly switches to the BDF scheme, if the system becomes stiff. With stiffness detection and the combination of two methods, one for non-stiff and one for stiff systems, we could assume that this solver would be the best performer and the most stable between all solvers. These eight solvers were implemented in the RHENOMS™ simulation platform [4]. Further details about these 8 ODE solvers are provided in the supplemental document.

Materials and Methods

As a common basis of benchmarking, we used the Rhenovia’s biosimulation platform RHENOMS and several models of synaptic receptors/channels. Figure 2 shows an example of a very simplified biological kinetic model of a generic glutamate receptor called elementary model. The glutamate neurotransmitter is represented by Glu, while R and RGlu represent the receptor in the free and liganded state, respectively. In this kinetic scheme, k1 and k2 represent the association and dissociation rate constants for the binding of the neurotransmitter on the receptor. The ordinary differential equations corresponding to this kinetic model consist in the system of equations (2), and the bilinearity is given by the term image


Figure 2: Example of an elementary kinetic model for ligand/receptor binding.

       image                               (2)

In this study, we choose four elementary models with different kinetic structures and dynamic characteristics (fast or slow), and corresponding to ligand-gated receptor (activated by a ligand such as glutatmate or GABA) or voltage-dependent channels (activated by a change of the potential), two important elements for synaptic transmission. The first model tested in this study is represented in Figure 3, and could represent either an AMPA or NMDA synaptic receptor according to the set of parameters used [27,28]. To differentiate these two models, we will call them AMPA7 and NMDA7, respectively. Testing the algorithms with two sets of parameters applied to the same model structure allows differentiating the impact of the kinetic scheme and the impact of the parameters on the algorithms performances. Several other models developed by Rhenovia were then tested. The next tested model is a model of NMDA (NR1/NR2A) receptor [29], with 15 states variables (referred later as NMDA15), qualified as slow based on the response of the model (50-250 msec) in the referential of synaptic transmission. This NMDA receptor model gives more details and information on the receptor characteristics than the first NMDA7 model. The fourth model represents the GABAA synaptic receptor, which is considered as a fast model (20-50 msec). Lastly, a model of an N-type voltage-dependent calcium channel (VDCC) was tested. This VDCC model is relatively fast (0-5 msec for activation), and differs from the models previously described as it is based on the Hodgkin- Huxley formalism [30], instead of kinetic reactions. This VDCC model was parameterized and validated to fit Jaffe and Poirazi results [31,32]. The features of all tested models are summarized on Table 2.


Figure 3: Kinetic scheme of the AMPA7/NMDA7 model depending to the set of parameters. The rectangles represent the different state variables. The green ovals named Glu represent the input of the model and the location with where it interacts in the model.

Model State size Number of equations Number of parameters Number of inputs Dynamic rate (msec)
AMPA7 7 7 16 1 20-50
NMDA7 7 7 16 1 50-250
NMDA15 15 15 38 3 50-250
GABAA 8 8 18 1 20-50
N type VDCC 2 6 12 1 0-5
Glutamatergic synapse 144 326 432 1 unknown

Table 2: Features of the models used in this study.

We stimulated each model with two different protocols: 1) Protocol 1 (P1) is a single event (Figure 4A and 2) Protocol 2 (P2) is the same event repeated 4 times with a 10 msec interval (Figure 4B). These two protocols were selected because they correspond to the type of electrical activity that takes place in the brain under most conditions. In the resting state, neurons communicate at low frequency, with single action potential (the single event protocol). On the other hand, under a number of conditions (learning, exploration), neurons increased their firing frequency and often emit bursts of action potential at frequency between 5 and 100 Hz. As an example, we selected 10 Hz. The repeated events are protocols often used for testing synaptic receptor models [33]. In addition, this low frequency could show if input had an important impact on solver performances. For the ligand-gated receptors, the single event consists in a 1 msec pulse of 1 mM glutamate, while for the voltage-gated channel the single event is a 1 msec depolarization step from -70 mV to 0 mV. These two protocols were tested with two sets of tolerance, which are used to validate a computed step by all the numerical ODE resolution methods: 1) a common set of tolerance with the relative tolerance equal to 1e-3 (Rtol), and the absolute tolerance equal to 1e-6 (Atol), and 2) a set of restrictive tolerance with Rtol=1e-6 and Atol=1e-9. These two tolerance sets are respectively called Tol1 and Tol2. Although this choice may appear arbitrary, experience leads us to consider these values as relevant for numerical resolution methods [26,34].


Figure 4: Stimulation protocols. A: Single event protocol P1, B: Repeated event protocol P2. Inserts represent the stimulus embedded in the entire simulation run.

Explicit Runge-Kutta methods are stable algorithms, as long as the integration step-size remains small enough. Therefore, we selected the fourth-order explicit Runge-Kutta (RK4) algorithm solver with a constant step-size of 0.5 µsec, as our reference to ensure that the solution of the simulation would remain stable and accurate. To assess performance, we used the execution time (speed) and the number of points (which is proportional to memory consumption) necessary to complete the simulation (i.e. to reach a stable numerical solution), which satisfy the tolerance parameter sets. As the execution time differs a little for each simulation, we repeat it ten times and make an average for measuring the execution times. For each model and protocol, the calculated output of the simulation is the current generated by receptor/ channel activation as a function of time. We also evaluated the relative precision of each algorithm by determining the Mean Square Error (MSE), and the Normalized Root Mean Square Deviation (NRMSD) between the results produced by a solver and those generated by the reference solver (RK4). More details are available in the supplemental document for the measurement of these two accuracy criteria. These 4 performance criteria (execution time, memory consumption, MSE and NRMSD) were selected as indices for rapidly obtaining the most accurate result using the lowest computer resources.

To quantify each algorithm performance, we simplified the algorithm selection problem and used the ||p|| norm value, as illustrated in equation (3), where pi represents our selected criteria and wi the priority that user give for each criteria, respectively. Importantly, in this study, we decided to use the same priority level for all four criteria, and thus, the w vector is [1, 1, 1, 1]T. The norm value provides a rapid comparison of the overall performance of the various algorithms, as it generates a single value for the performance of each algorithm. Algorithms yielding a small norm value will thus be considered as more efficient than others yielding a larger norm value.

       image                               (3)

Furthermore, to verify the performance variation of algorithms, we computed the differences in norm values between the different sets of tolerance and protocols. The performance variations enable us to determine whether an algorithm is sensitive to input protocol, or to tolerance parameter sets. Simulations were performed on a workstation with a LINUX (Ubuntu 10.04) operating system and an Intel® Xeon® CPU at 2.67 GHz frequency equipped with 12 Gbytes of RAM and the version 1.6 of Java installed.

Results and Discussion

In computational neuroscience, most models are stimulusdependent and bilinear. The major goal of our study was to simplify the algorithm selection problem by benchmarking the solvers performances implemented in RHENOMS. This simplification could make a complex problem accessible to all users, and provide a recommendation for the possible use of biological bilinear kinetic models. Therefore, we will not discuss here the details of the solvers or of the models.

AMPA7/NMDA7 receptor model

This kinetic scheme, presented in Figure 3, comprises 7 states, 7 equations, 16 parameters and 1 input variable. It could model a fast receptor like AMPA or a slower one like NMDA, according to the set of parameters used. For this model, we choose the open state probability represented by the O4 rectangle in Figure 3 as readout. The open probability of the AMPA7 model is depicted on Figure 5, and the same state for the NMDA7 model is shown in Figure 6. As we can see, the two models (AMPA7 and NMDA7) do not give the same dynamics on the open probability state. The raw data of the AMPA7 model are depicted on Table S-2 and S-3 on the additional document for the P1 and P2 protocols. Similar data for the NMDA7 model are given in Table S-4 and S-5. The performance values are presented on Table 3 for the AMPA7 and on Table 4 for the NMDA7 model.


Figure 5: Normalized open state probability for AMPA7 model. Insert represents the results embedded in the entire simulation run.


Figure 6: Normalized open state probability for NMDA7 model.

Protocol 1 (single pulse)
Tol1 1.0 1.34e-3 1.46e-3 5.27e-4 5.23e-4 5.59e-4 1.02e-3 5.49e-4
Tol2 1.0 1.35e-3 4.31e-2 6.44e-4 5.98e-4 1.01e-3 1.97e-2 8.64e-4
ΔP1 Tol1/Tol2 - 1e-5 4.16e-2 1.17e-4 7.5e-5 4.5e-4 1.87e-2 3.15e-4
Protocol 2 (repeated pulses)
Tol1 1.0 2.39e-3 7.09e-3 6.14e-4 5.87e-4 7.85e-4 2.2e-3 7.19e-4
Tol2 1.0 2.41e-3 1.65e-2 1.12e-3 7.86e-4 2.61e-3 1.69e-2 1.79e-3
ΔP2 Tol1/Tol2 - 2.1e-5 9.41e-3 5.06e-4 1.99e-4 1.8e-3 1.47e-2 1.1e-3
ΔTol1 P1/P2 - 1.1e-3 5.6e-3 8.7e-5 7.5e-5 2.26e-4 1.2e-3 1.7e-4
ΔTol2 P1/P2 - 1.1e-3 2.66e-2 4.76e-3 7.11e-3 1.6e-3 2.8e-3 9.26e-4

Table 3: Quantification and comparison of normalized solver performances for the AMPA7 model.

Protocol 1 (single pulse)
Tol1 1.0 5.14e-4 1.83e-3 5.24e-4 5.21e-4 5.41e-4 8.36e-4 5.25e-4
Tol2 1.0 5.36e-4 3.93e-3 5.84e-4 5.63e-4 8.14e-4 9.54e-3 7.2e-4
ΔP1 Tol1/Tol2 - 2.2e-5 2.1e-3 6e-5 4.2e-5 2.73e-4 8.7e-3 1.95e-4
Protocol 2 (repeated pulses)
Tol1 1.0 5.29e-4 4.55e-3 5.64e-4 5.57e-4 6.08e-4 1.23e-3 6.28e-4
Tol2 1.0 5.79e-4 9.06e-3 7.1e-4 6.58e-4 1.13e-3 1.63e-2 1.03e-3
ΔP2 Tol1/Tol2 - 2.1e-5 9.41e-3 5.06e-4 1.99e-4 1.8e-3 1.47e-2 1.1e-3
ΔTol1 P1/P2 - 1.5e-5 2.72e-3 4e-5 3.6e-5 6.7e-5 3.9e-4 1.03e-4
ΔTol2 P1/P2 - 4.3e-5 5.13e-3 1.26e-4 9.5e-5 8.14e-4 6.76e-3 3.8e-4

Table 4: Quantification and comparison of normalized solver performances for the NMDA7 model.

In order to find the best algorithm for this model, we identified the solver(s) with the smallest norm values for the two stimulation protocols and the two sets of tolerance (Figure 7). Based on the first protocol (Figure 7B), the hybrid IMEX solver provides the best performances for the AMPA7 model whereas solvers with a BDF scheme are the worst performers. With the second protocol (Figure 7D), we observe the same tendancy. If we analyze the raw data (in the supplemental document), TR-BDF2 and JVODE BDF needed an important number of points to compute the results. Contrary, the IMEX algorithm is the one that had the lowest memory footprint and had the most accurate results. In addition, its execution time is very short.


Figure 7: Normalized norm values for the AMPA7 receptor model. A: Review of the protocol P1. B: The normalized norm values for the 8 solver algorithms with P1. C: Review of the protocol P2. D: The normalized norm values for the 8 solver algorithms with P2. All ||p|| norms are normalized to the reference RK4 ||pRK4|| norm.

According to performances obtained on this AMPA7 model, TRBDF2 and JVODE BDF appear to be the most sensitive to protocol and tolerance sets variations with this model. The RKF solver is the least sensitive to tolerance variations, whereas hybrid solvers were the least sensitive to protocol variations.

In Figure 8, the same graphics are plotted for the parameters, which model the NMDA7 synaptic receptor to visualize the ||p|| performances values. With a slow dynamics on this kinetics scheme, the RKF solver became the best choice for all our tested protocols and tolerance sets. As for the AMPA7 model, solvers using BDF scheme are the worst performers. In fact, the raw data (in the supplemental data) show that these solvers needed a large number of points, and a long execution time. Regarding the performance variations on these models, RKF appears to be the least sensitive to protocols and tolerance variations, whereas JVODE BDF and TR-BDF2 are the most sensitive.

This kinetic scheme generates different dynamics with different parameters. With the parameters modeling the fast AMPA7 model, the IMEX solver appears to be the most suitable; and with the parameters modeling the slow NMDA7 model, the RKF solver appears to be the most suitable. With this observation, it appears that the dynamics of the system is more important than the structure of the kinetic scheme.

NMDA15 receptor model

The NMDA15 (NR1/NR2A) receptor is a relatively slow (50-250 msec) glutamate receptor. This model was previously calibrated and validated to fit a variety of experimental data [26]. The tested NMDA15 model had 3 inputs: glutamate concentration, glycine concentration and depolarization. For this study, we choose the glutamate as the protocol input. Glycine and depolarization are parameterized to be constant. The current generated by receptor activation using the single pulse protocol P1, and calculated with all the solvers is depicted on Figure 9. All algorithms provided similar results, and Figure 10 summarizes these performance values. The criteria of the NMDA15 model are quantified on Table S-6 and S-7 of the supplemental document. Algorithm performance values for the NMDA15 receptor model are quantified in Table 5.


Figure 8: Normalized norm values for the NMDA7 receptor model. A: Review of the protocol P1. B: The normalized norm values for the 8 solver algorithms with P1. C: Review of the protocol P2. D: The normalized norm values for the 8 solver algorithms with P2. All ||p|| norms are normalized to the reference RK4 ||pRK4|| norm.


Figure 9: Normalized NMDA15 receptor current. NMDA15 receptor’s current resulting from a single pulse protocol generated by all solver algorithms.

To find the best algorithm for this model, we identified the solver(s) with the smallest norm values for the two stimulation protocols and the two sets of tolerance (Figure 10B and 10D). Based on Figure 10B, for P1 protocol, the algorithm performances did not differ much between Tol1 and Tol2 tolerances except for the BDF scheme (TR-BDF2 and JVODE BDF solvers), for which a more than 10 times difference on ||p|| performance value is observed. For P2 protocol, the performance differences between Tol1 and Tol2 parameters were more pronounced with the last three algorithms. As for P1 protocol, BDF scheme performances deteriorated significantly with Tol2. For both stimulation protocols, the best solver was IMEX for both tolerance parameters, while the least performing algorithm was RKF for Tol1 and TR-BDF2 for Tol2 for both protocols. Considering the algorithm performance variations presented in Table 5, we concluded that TR-BDF2 solver was the most sensitive to protocol and tolerance parameters for the NMDA15 receptor model, whereas RKF solver has the overall most stable performance.


Figure 10: Normalized norm values for the NMDA15 receptor model. A: Review of the protocol P1. B: The normalized norm values for the 8 solver algorithms with P1. C: Review of the protocol P2. D: The normalized norm values for the 8 solver algorithms with P2. All ||p|| norms are normalized to the reference RK4 ||pRK4|| norm.

      BDF2 brock   ADAMS BDF  
Protocol 1 (single pulse)
Tol1 1.0 8.025e-3 3.587e-3 5.25e-4 5.22e-4 5.45e-4 8.83e-4 5.59e-4
Tol2 1.0 8.051e-3 1.8478 6.04e-4 5.76e-4 8.96e-4 1.163e-2 8.79e-4
ΔP1 Tol1/Tol2 - 2.6e-5 1.8444 7.9e-5 5.4e-5 3.51e-4 1.075e-2 3.29e-4
Protocol 2 (repeated pulses)
Tol1 1.0 8.034e-3 8.004e-3 5.74e-4 5.63e-4 6.76e-4 1.51e-3 6.6e-4
Tol2 1.0 8.067e-3 2.7826 8.04e-4 7.18e-4 1.554e-3 2.18e-2 1.35e-3
ΔP2 Tol1/Tol2 - 3.3e-5 2.7746 2.3e-4 1.55e-4 8.68e-4 2.029e-2 6.92e-4
ΔTol1 P1/P2 - 9.1e-6 4.417e-3 4.9e-5 4.1e-5 1.31e-4 6.27e-4 1.01e-4
Δ Tol2 P1/P2 - 1.6e-6 0.9347 2.4e-4 1.42e-4 6.48e-4 1.016e-2 4.73e-4

Table 5: Quantification and comparison of normalized solver performances for the NMDA15 model.

GABAA receptor model

The GABAA synaptic receptor is a fast model with approximately the same dynamics as the AMPA7 model. It is sensitive to the GABA neurotransmitter. Our GABAA model was developed using parameters is presented [35]. The tested GABAA model had 8 states, 8 equations, 18 parameters and 1 input. Figure 11 shows the GABAA model current generated with protocol P1. This current was our readout to determine performance values. Our results indicate that all solvers provide similar results. This model’s quantification criteria are depicted in details in the supplemental document on Table S-8 and S-9 (respectively for P1 and P2 protocols). The algorithm performance values for the tested GABAA receptor model are quantified in Table 6.


Figure 11: Normalized GABAA receptor model current. GABAA receptor’s current resulting from the single pulse protocol generated by all solvers. Insert represents the results embedded in the entire simulation run.

Protocol 1 (single pulse)
Tol1 1.0 1.19e-3 5.66e-3 5.32e-4 5.27e-4 5.74e-4 8.62e-4 4.88e-4
Tol2 1.0 1.98e-3 2.11e-2 6.49e-4 6.06e-4 1.05e-3 1.47e-2 8.88e-4
ΔP1 Tol1/Tol2 - 7.9e-4 1.54e-2 1.17e-4 7.9e-5 4.76e-4 1.38e-2 4e-4
Protocol 2 (repeated pulses)
Tol1 1.0 1.95e-3 6.12e-3 6.07e-4 5.87e-4 7.48e-4 1.95e-3 7.53e-4
Tol2 1.0 2.09e-3 1.86e-2 1.03e-3 8.8e-4 2.46e-3 3.07e-2 1.7e-3
ΔP2 Tol1/Tol2 - 1.4e-4 1.25e-2 4.23e-4 2.93e-4 1.71e-3 1.12e-2 9.47e-3
ΔTol1 P1/P2 - 7.6e-4 4.6e-4 7.5e-5 6e-5 1.74e-4 1.08e-3 2.65e-4
ΔTol2 P1/P2 - 11e-4 2.5e-3 3.81e-4 2.74e-4 6.6e-4 1.6e-2 8.12e-4

Table 6: Quantification and comparison of normalized solver performances for the GABAA model.

The graphs presenting the ||p|| performance value in Figure 12 quantitatively helps to determine the most appropriate algorithm for the GABAA receptor model. A brief examination of Figure 12B and 12D suggests that TR-BDF2 and JVODE BDF are the worst solvers for this model, as they yield large ||p|| values. JLSODE appears as the best solver for this model with the first protocol (P1), and the first tolerance parameters (Tol1). For all others simulations, IMEX yields the best performances. Observing the performance variations on Table 6, the worst solvers (JVODE BDF and TR-BDF2) for this model are the most sensitive to input and tolerance variations. With the GABAA model and the tested situations, IMEX and RKF yield the most consistent performances.


Figure 12: Normalized norm values for the GABAA receptor model. A: review of the protocol P1. B: The normalized norm values for the 8 solver algorithms with P1. C: review of the protocol P2. D: The normalized norm values for the 8 solver algorithms with P2. All ||p|| norms are normalized to the reference RK4 ||pRK4|| norm.

In fact, from Table S-8 and S-9 in the supplemental file, we could observe that RKF is the faster solver. However, the IMEX was most accurate with the GABAA model with all tested simulation. The IMEX was clearly not competing with other solvers for the execution time on the second protocol, as it needs longer time to compute the results. But, its low memory consumption combined with its final accuracy background made it the best solver for this particular model.

N-Type VDCC model

In general, it has proven difficult to model voltage-dependent calcium channels (VDCC), due to the existence of a tail current at the end of the plateau current. For our comparative studies, we decided to simulate the N-Type VDCC, a high voltage-activated channel. This channel is very fast (0-5 msec activation), as shown in Figure 13. For this channel, we increased the depolarization duration to 10 msec. The interval between pulses was adjusted in protocol P2 to maintain the 10 Hz stimulation condition. All the solvers were able to accurately fit the references results. The various performance values parameters are summarized in Figure 14.


Figure 13: N-type VDCC current. N-type VDCC’s current resulting from a single event protocol generated by all ODE solvers Insert represents the results embedded in the entire simulation run.


Figure 14: Normalized norm values for the N type VDCC model. A: review of the protocol P1. B: The normalized norm values for the 8 solver algorithms with P1. C: review of the protocol P2. D: The normalized norm values for the 8 solver algorithms with P2. All ||p|| norms are normalized to the reference RK4 ||pRK4|| norm.

No significant differences were observed between the eight considered algorithms, except for TR-BDF2, which generated the highest ||p|| norm performance values for both protocols and tolerance sets. Comparing the overall performance values of the solvers (Table 7), JLSODE was the algorithm that provided the best performance for both stimulation protocols and tolerance sets. However, RKF performance values were very close to that of JLSODE (1e-6 difference).

Protocol 1 (single pulse)
Tol1 1.0 5.004e-2 5.154e-2 5.006e-2 5.007e-2 5.008e-2 5.008e-2 5.003e-2
Tol2 1.0 5.004e-2 0.683 5.006e-2 5.007e-2 5.008e-2 5.008e-2 5.003e-2  
ΔP1 Tol1/Tol2 - 2.5e-8 0.6314 9.84e-7 9.4e-8 1.056e-6 1e-6 1.2e-8
Protocol 2 (repeated pulses)
Tol1 1.0 5.001e-2 5.647e-2 5.018e-2 5.019e-2 5.037e-2 5.037e-2 5e-2
Tol2 1.0 5.001e-2 3.215 5.018e-2 5.019e-2 5.037e-2 5.037e-2 5e-2
ΔP2 Tol1/Tol2 - 2.7e-10 3.1587 2.8e-8 1.08e-7 6.4e-8 2e-9 1e-9
ΔTol1 P1/P2 - 5.976e-6 4.933e-3 1.198e-5 1.192e-5 2.988e-5 2.899e-5 5.99e-6
ΔTol2 P1/P2 - 6e-6 2.5322 1.102e-5 1.213e-5 3.006e-5 3e-5 6.01e-6

Table 7: Quantification and comparison of normalized solver performances for the N-type VDCC model.

When we analyzed the criteria quantification table (Table S-10 and S-11 in supplemental file), all algorithms required approximately the same number of points, except for TR-BDF2, which needed more points with Tol2 parameters (10 times more for P1 and 64 times more for P2). In addition, all solvers generated almost the same NRMSD, except for RKF (for both stimulation protocols) and TR-BDF2 (for Tol2 tolerance set).

As illustrated in the supplemental file (Figure S-1), when switching from Tol1 to Tol2, the TR-BDF2 solver yielded smaller error values (both MSE and NRMSD), although execution time and memory consumption (i.e. number of points) increased and exceeded RK4 values. In addition, a small relative tolerance value (Rtol<1e-6) was not required to obtain accurate results with TR-BDF2 solver. These results indicate that the JLSODE solver is the best choice for this VDCC model, as it provides the best performance. However, RKF results were very close. In addition, when observing performance variations, RKF performances were less sensitive to input protocols or tolerance sets.

Application to RHENOMS model

In order to complete our study, we then analyzed the performance of the different ODE solvers with a model of glutamatergic synapse that integrates a large number of synaptic elementary models, such as receptors, channels, transporters, enzymes and signaling pathway models. Overall, the glutamatergic synapse platform integrates more than 300 equations to be solved at each integration step, 144 variable states for 21 bilinear synaptic elementary kinetic models. A simplified representation of the glutamatergic model is illustrated in Figure 15, which specifically highlights the previously tested elementary models. In addition to the large number of equations, the neurotransmitter release model uses the Systems Biology Markup Language (SBML) events [36], which results in more non-linearity, stiffness and difficulties to solve equations during simulations.


Figure 15: Simplified schematic representation of the used glutamatergic synapse model.

Despite the power and available memory in the computer we used, computing the glutamatergic synapse model with TR-BDF2 algorithm with the second tolerance set (Tol2) was beyond the available capacity, and we, therefore, cannot present these results. Indeed such calculation would require a computer cluster with much more memory as the number of points computed for all elementary models becomes prohibitive. The second set of tolerance placed TR-BDF2 in its worst configuration, and this algorithm did not need a restrictive tolerance to provide accurate enough results. In addition, users do not often change the default tolerance parameters, which are commonly set close to our first set (Tol1). For these reasons, we present algorithm performances with the first tolerance set only (Tol1) with the two stimulation protocols (P1 and P2). Quantified criteria are presented in the supplemental file by Table S-12.

Considering the number of synaptic receptors integrated in this model, we selected the postsynaptic current as the main readout, as it corresponds to the sum of currents generated by activation of all receptors, and is a common readout in biological experiments. In addition, we replaced the number of points (used to evaluate the performance of previous tested models) by the memory consumption for the simulation, which seems more appropriate for this complex model. A two-second simulation is performed, with a 100 msec delay before the start of stimulation was applied, and the quantified performances values and variations for the glutamatergic synapse model are presented in Table 8.

      BDF2 brock   ADAMS BDF  
P1 1.0 0.3156 0.4989 0.1334 0.4327 0.1425 0.1456 0.1537
P2 1.0 0.3273 0.261 0.1663 0.4705 0.1884 0.1831 0.1926
Δ Tol1 P1/P2 - 1.16e-2 0.2379 3.291e-2 3.773e-2 4.585e-2 3.753e-2 3.886e-2

Table 8: Quantification and comparison of normalized solver performances for the glutamate synapse model.

As shown in Figure 16, the Rosenbrock solver provided the smallest norm values for both protocols. In contrast, IMEX, which is a hybrid solver too, with a fourth-order RK-scheme, generated a large norm value. Surprisingly, TR-BDF2 algorithm improved its performance between P1 and P2 protocol. This could be due to its step-size adaptation tstop (see supplemental file), with which TR-BDF2 readapts the integration step-size around an abrupt variation (i.e. a stimulation input). In theory, the more stimulation the model receives, the more efficient the algorithm becomes [16]. It was the sole algorithm with a better performance with P2 as compared to P1. TR-BDF2 was also the most sensitive to input protocol with the glutamatergic synapse model, whereas RKF was the least sensitive.


Figure 16: Normalized norm values for glutamatergic model. Normalized norm values for the 8 solver algorithms for the glutamate synapse model stimulated with P1 and P2 protocols. All ||p|| norms are normalized to the reference RK4 ||pRK4|| norm.


Computational neuroscience encompasses a wide range of kinetic models with very different characteristics. It is very difficult to select the appropriate algorithm to solve the ODE systems representing biological bilinear kinetic models. Indeed, in order to select the right algorithm, as suggested by Rice, users need to know the model input, model characteristics (size), and ODE solver structure. A benchmark comparing various ODE solver algorithms or a recommendation could help users to select the most appropriate algorithm for a given simulation. We simplified the algorithm selection problem and benchmarked eight ODE solvers performances using several types of kinetic models with two different stimulation protocols. Although this benchmark was done in Java programming language, it is important to note that the conclusions would not change with another programming language, as long as good coding practices are respected. Our overall conclusions are summarized in Table 9.

Model Protocol Tolerance Best ||p|| Worst ||p||
Fast model: P1 Tol1 IMEX TR-BDF2
  P2 Tol1 IMEX TR-BDF2
Slow model: P1 Tol1 RKF TR-BDF2
  P2 Tol1 RKF TR-BDF2
Slow model: P1 Tol1 IMEX RKF
  P2 Tol1 IMEX RKF
    Tol2 IMEX TR-BDF2
Fast model: P1 Tol1 IMEX TR-BDF2
Fast model: P1 Tol1 JLSODE TR-BDF2
Glutamatergic P1 Tol1 Rosenbrock TR-BDF2
synapse P2 Tol1 Rosenbrock IMEX

Table 9: Summary of the best and worst solvers.

According to Table 5, TR-BDF2 was the least appropriate algorithm for all the models used in this study, with the selected stimulation protocols and tolerance parameter sets. TR-BDF2 is an ODE solver combining a trapezoidal scheme, followed by a secondorder Backward Differential Formula (BDF2). As shown in our results, a relative tolerance set below or equal to 1e-6 reduced this algorithm’s performance, and produced the highest overall ||p|| norm values. This algorithm could produce better performances and accurate enough results with a larger tolerance, but a more economical (in terms of performance) solver could then give accurate enough results as well.

The IMEX algorithm was the overall best choice for the tested AMPA7, NMDA15 and GABAA models. RKF solver was appropriate for the NMDA7. JLSODE was the best choice for the fast N type Calcium Channel model and Rosenbrock solver was the most appropriate for the glutamatergic synapse model, which was considered as a complex model. However, algorithm performance variations indicated that the RKF algorithm was the most stable one, when comparing stimulation protocols and tolerance sets. In fact, RKF solver appeared to be more stable to input and tolerance sets variation compared to JLSODE, which possesses stiffness detection with our w preference criteria. This study clearly showed that for a kinetic model, changing the model dynamics will affect the solver performances. In fact, for our model, the dynamic appears to be more important than the underlying kinetic structure.

Synaptic elementary kinetic models are composed in most cases of time-varying continuous time control system with input image, also called affine in control or bilinear. Due to the dependency of ODEs with respect to inputs, it is important to note that model stiffness could differ if changing the stimulus protocol. Additionally, algorithm performance could change depending on the selected performance criteria or the criteria priority level (or weight), as described by Rice [2]. The presented method gives a priori information on solvers performance using an equal priority level for all four selected criteria. Users who need to run a large number of simulations with a model (for example to perform a sensitivity analysis or to test a drug concentration effect) could use this method to optimize the computational effort.

As a comparison with other simulation platforms or software, NEURON [37] uses the CVODE algorithm with a manual selection of either Adams-Moulton or BDF scheme. The Copasi [38] tool (biochemical network simulator) uses a fourth-order hybrid Runge- Kutta algorithm, in addition to the LSODE algorithm, which integrates Petzolds stiffness detection. BioUML [24] developers made the same choice as Copasi ones, using a fourth-order hybrid Runge-Kutta algorithm, in addition to the CVODE algorithm implemented in Java language. The Matlab software uses the same approach as ours, which consists in integrating several resolution algorithms, and allowing the users to choose their algorithm based on the model they intend on simulating.

PYTHIA [39] is software, which provides a strict selection algorithm and tries to select a pair of machine/algorithm, in order to optimize the computational effort with a given problem. However, to use this tool, users need some expertise due to the complexity of the algorithm selection problem. As synaptic receptor models are commonly designed using SBML [36], a possible extension of this work would be to generate an automatic algorithm selection, which uses the strength of SBML. The idea would be to use SBML standard for automatically extracting the model features, and make the process completely transparent to the user.


We would like to thank Lhassane Idoumghar at the University of Haute Alsace (France) and Tameem Albash at the University of Southern California (USA, CA), for constructive discussions; the operational team of Rhenovia Pharma, which helped for the elaboration of this study. Rhenovia obtained financing by the French National Agency of Research and Technology (ANRT), with a CIFRE scholarship (432/2010) for Merdan Sarmis. The authors would like thank both anonymous referees for their valuable comments and suggestions.


Select your language of interest to view the total content in your interested language
Post your comment

Share This Article

Relevant Topics

Recommended Conferences

Article Usage

  • Total views: 11539
  • [From(publication date):
    May-2013 - Aug 17, 2017]
  • Breakdown by view type
  • HTML page views : 7776
  • PDF downloads :3763

Post your comment

captcha   Reload  Can't read the image? click here to refresh

Peer Reviewed Journals
Make the best use of Scientific Research and information from our 700 + peer reviewed, Open Access Journals
International Conferences 2017-18
Meet Inspiring Speakers and Experts at our 3000+ Global Annual Meetings

Contact Us

© 2008-2017 OMICS International - Open Access Publisher. Best viewed in Mozilla Firefox | Google Chrome | Above IE 7.0 version