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

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.


Introduction
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.  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), is the state variable vector and m t u   ) ( is the s ste s i put e to ith ) (t u i the th i 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. We compared eight ODE solvers (Table 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 Curtis 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,11,12]. This implies that the stiffness of a model is closely related to the selected numerical resolution method. According to Ekeland, the precise definition of the stiffness of a system of equations is not crucial from a practical point of view [12]. 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.  [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 [24] platform developers. JVODE includes two schemes: 1) a variable order (from 1 to 12), explicit Adams-Moulton method for nonstiff 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 Linda R. Petzold [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 we e i ple e ted i the 'HENOMS™ simulation platform [4]. Further details about these 8 ODE solvers are provided in the supplemental document.

Material and Methods
As a common basis of benchmarking, e used the 'he o ia 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, k 1 and k 2 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 1 kR G l u  . Figure 2: Example of an elementary kinetic model for ligand/receptor binding. 12 21 dRGlu k R Glu k RGlu dt dR k RGlu k R Glu dt 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 GABA A 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.

Model
State size  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 Tol 1 and Tol 2 . Although this choice may appear arbitrary, experience leads us to consider these values as relevant for numerical resolution methods [26,34]. 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 satisfying the tolerance parameter sets.

Number of inputs
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 p i represents our selected criteria and w i 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.
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 stimulus-dependent 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 Table 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 Table S-5. The  performance values are presented on Table 3 for the AMPA7 and on Table 4 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. Protocol 2 (repeated pulse)  According to performances obtained on this AMPA7 model, TR-BDF2 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.

RK4 RKF TR-BDF2
In Figure 8, the same graphics are plotted for the parameters, which models 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 Tables S-6 and S-7 of the supplemental document. Algorithm performance values for the NMDA15 receptor model are quantified in Table 5. 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 ( Figures 10B and 10D). Based on Figure 10B, for P1 protocol, the algorithm performances did not differ much between Tol 1 and Tol 2 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 Tol 1 and Tol 2 parameters were more pronounced with the last three algorithms. As for P1 protocol, BDF scheme performances deteriorated significantly with Tol 2 . For both stimulation protocols, the best solver was IMEX for both tolerance parameters, while the least performing algorithm was RKF for Tol 1 and TR-BDF2 for Tol 2 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.

GABA A receptor model
The GABA A synaptic receptor is a fast model with approximately the same dynamics as the AMPA7 model. It is sensitive to the GABA neurotransmitter. Our GABA A model was developed using parameters presented in [35]. The tested GABA A model had 8 states, 8 equations, 18 parameters and 1 input. Figure 11 shows the GABA A 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 odel s ua tification criteria are depicted in details in the supplemental document on Tables S-8 and S-9 (respectively for P1 and P2 protocols). The algorithm performance values for the tested GABA A receptor model are quantified in Table 6.   The graphs presenting the ||p|| performance value in Figure 12 quantitatively helps to determine the most appropriate algorithm for the GABA A 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.

RK4 RKF TR-BDF2
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. In fact, from Table S-8 an S-9 in the supplemental file, we could observe that RKF is the faster solver. However, the IMEX was most accurate with the GABA A 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.  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).
As illustrated in the supplemental file ( Figure S-1), when switching from Tol 1 to Tol 2 , 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.
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 (Tol 2 ) 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 (Tol 1 ). For these reasons, we present algorithm performances with the first tolerance set only (Tol 1 ) 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.  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 t stop (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 stimulations 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.

Conclusion
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.
According to Table 9, 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 second-order Backward Differential Formula (BDF2). As shown in our results, a relative tolerance set below or equal to 1e -6 reduced this algo ith 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.  Table 9: Summary of the best and worst solvers. Summary of the best and worst solvers for the two types of models with two different stimulation protocols and tolerance values, and the applicative simulation with the glutamatergic synapse. P1 stands for the single event protocol and P2 stands for the repeated event protocol.
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 (ẋ=f(t,x,u)), 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 a 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.

The Numerical Resolution Methods
We implemented the fourth-order Runge-Kutta solver (RK4), as indicated in the book edited by Hairer, Norsett and Wanner (referenced in the main article). We did not modify the parameters of this algorithm, which is very easy to implement. The Runge-Kutta Fehlberg solver (RKF) is a fourth/fifth-order explicite RK-scheme. We implemented it at presented by Bogacki and Shampine (referenced in the main article). The c parameters for RKF solver at different stages are presented in Table S-1, and the stage equations are:  The TR-BDF2 algorithm corresponds to the ode23tb function in the Matlab software with the γ parameter equal to (2− √ 2), as explained by S. Dharmaraja who is referenced in the main article. In the original version of this algorithm, the authors describe a t stop parameter, which readjusts its calue for the next calculated step-size befor a large variation. Since in our case, we cannot forecast future variations, we transformed this equation into equation (S-7), which better fits our simulations.
The IMEX and Rosenbrock solvers are both fourth-order, RK-scheme-based hybrid solvers. For each computation step, hybrid solvers resolve the system twice, first with an explicit scheme, then with an implicit one. The two results are compared and if they are similar, the step is accepted. If not, the step-size is reduced and the computation step restarts. These two solvers are implemented as described in their respective references (see main article).
JVODE is a Java version of the CVODE solver. This is the only solver that we did not implement, as it is provided as open source files by the BioUML developers. This solver is in fact composed of two solvers: 1) an explicit Adams-Moulton method with variable order (from 1 to 12) for non-stiff systems and 2) an implicit Backward Differential Formula (BDF scheme) with variable order (from 1 to 5) for stiff systems. These two solvers are called JVODE ADAMS and JVODE BDF, respectively.
The last algorithm we consider is JLSODE, a Java version of the LSODE solver, which has the same Adams-Moulton and BDF schemes as JVODE. This last solver adds stiffness detection, as written by Linda R. Petzold (referenced in the main article). With this stiffness detection, JLSODE starts with the Adams-Moulton scheme and switches to the BDF scheme if the system becomes stiff, and vice-versa. In this case, we used the Adams-Moulton and BDF solvers integrated in JVODE and added stiffness detection, as described by Linda R Petzold (referenced in the main article).

Material and Methods
To simplify the study and its reproducibility, we used sequential computations, even if the platform works on parallel computing with MPI library. The two tolerances, Rtol and Atol, are used to validate a computed step for the numerical ODE resolution methods. All simulations are run ten times and the presented quantitative results are the averages of the ten simulation results. All solvers are variable step-size except for the reference solver RK4, which uses a 5 µs constant step-size. To compute the accuracy criteria, we did an interpolation of all solver's results for the reference times.
In addition, in order to determine the sensitivity to stimulation protocoland tolerance parameters, we computed the following performance differences: 1) the differences between Tol 1 and Tol 2 for the first and second protocol, respectively called ∆P 1 Tol1/T ol2 and ∆P 2 Tol1/T ol2 and 2) the differences between the first and second protocol for Tol 1 and Tol 2 parameters, respectively called ∆Tol 1 P 1/P 2 and ∆Tol 2 P 1/P 2 .

AMPA7/NMDA7 receptor model
The kinetic model presented in Figure 3 of the main document could model either an AMPA or a NMDA receptor, depending on the set of parameters. The quantified criteria of the AMPA7 model are shown in Table S-2 for the P 1 protocol and Table S-3 for the P 2 protocol. TR-BDF2 and JVODE BDF required a large number of points to achieve an accurate result. In contrast, the hybrid IMEX solver generated more accurate results with the minimal number of points.
For NMDA7 model, the quantified criteria are in Table S-4 for the P 1 protocol and Table S-5 for the P 2 protocol. As for AMPA7 model, TR-BDF2 and JVODE BDF required a large number of points to achieve accurate results. The RKF solver consumed less memory but the hybrid IMEX solver was more accurate, according to MSE and NRMSD columns. With our selected criteria priority, RKF takes the advantage on IMEX solver due to its speed.

NMDA15 Receptor Model
Based on MSE and NRMSD values (Table S-6 and Table S-7, respectively) all the algorithms are accurate enough for this particular type of model. The Rosenbrock and IMEX solvers required the smallest number of points, which means they used the least amount of memory for the simulations for both stimulation protocols and both tolerance sets.
Analyzing the two tables, TR-BDF2 and RKF solvers appear to be the worst solvers for NMDA15 model for both protocols and tolerance sets. In addition, the performance of JVODE BDF solver seems to be worst with the second tolerance parameter Tol 2 . Thus, BDF-scheme solvers perform poorly with a restrictive tolerance set.

GABA A Receptor Model
Table S-8 and S-9 present the raw data for the criteria for protocols P 1 and P 2, respectively. As for AMPA7 model, the more accurate solver is IMEX and the faster is RKF. Despite RKF's rapidity, IMEX is more efficient due to its lower memory usage.

N-type VDCC Model
As shown in Figure S-1.B, all solvers did not provide the same current peak. All solvers have specific step-size adaptation algorithms and criteria, which result in different computing time. The last three solvers (JVODE ADAMS, JVODE BDF and JLSODE) include a step-size adaptation, but also an algorithm to adapt the order of the ODE solver, in order to provide accurate results without decreasing step-size. This variable step-size, together with the order adapt is likely to account for the different results obtained for the peak current. However, since variations in MSE and NRMSD are small (see Table S-10 and Table S-11), we can neglect these differences and consider that all solvers provide similar accurate results.
For the N-type VDCC model, all solvers required approximately the same number of points, except for TR-BDF2, which needed more points with Tol 2 parameters. In addition, all solvers generated almost the same NRMSD, except for RKF (for both stimulation protocols) and TR-BDF2 (for Tol 2 set tolerance). Considering the MSE error criterion, hybrid solvers are the worst solvers for this model. Switching from Tol 1 to Tol 2 , TR-BDF2 solver decreased the error parameters (both MSE and NRMSD), although execution time and memory consumption (i.e. number of points) increased and exceeded RK4 values.

Application to glutamatergic synapse Model
The postsynaptic current generated by the simulations is presented in Figure S-2 and Table S-12 presents the quantified results for the glutamatergic synapse model.
As was the case for the N-type VDCC model, the applicative case output currents (postsynaptic current) observed on Figure S-2 are different, depending on the choice of ODE solver (step-size and order variation). However, the small variations in MSE and NRMSD shown in Table S-12 allow us to neglect these differences and to consider that all solvers were accurate and provided the same postsynaptic current.
Hybrid solvers generated the largest MSE and NRMSD, which appears to be their weak features for solving biological kinetic models. But, as indicated in the main article, IMEX is the best solver for the NMDA receptor model (slow model) and Rosenbrock for RHENOMS, which models a glutamatergic synapse.
From Table S-12, MSE and NRMSD values indicated that all solvers are accurate; however, they also indicated that choosing the right solver could save execution time and memory consumption. Without evaluating solver performance, we might have selected a solver that would require a long simulation time. In fact, for a 2-second simulation of the glutamatergic synapse model, the execution time varies approximately between 150 and 900 seconds on our workstation depending on the selected algorithm. This difference would obviously become more important with longer simulations, underlying the importance of the choice of the numerical solver and its underlying parameters to obtain accurate results in minimal time.