Efficient Implicit Runge-Kutta Methods for Fast-Responding Ligand-Gated Neuroreceptor Kinetic Models

Neurophysiological models of the brain typically utilize systems of ordinary differential equations to simulate single-cell electrodynamics. To accurately emulate neurological treatments and their physiological effects on neurodegenerative disease, models that incorporate biologically-inspired mechanisms, such as neurotransmitter signalling, are necessary. Additionally, applications that examine populations of neurons, such as multiscale models, can demand solving hundreds of millions of these systems at each simulation time step. Therefore, robust numerical solvers for biologically-inspired neuron models are vital. To address this requirement, we evaluate the numerical accuracy and computational efficiency of three L-stable implicit Runge-Kutta methods when solving kinetic models of the ligandgated glutamate and γ-aminobutyric acid (GABA) neurotransmitter receptors. Efficient implementations of each numerical method are discussed, and numerous performance metrics including accuracy, simulation time steps, execution speeds, Jacobian calculations, and LU factorizations are evaluated to identify appropriate strategies for solving these models. Comparisons to popular explicit methods are presented and highlight the advantages of the implicit methods. In addition, we show a machinecode compiled implicit Runge-Kutta method implementation that possesses exceptional accuracy and superior computational efficiency. Keywords-implicit Runge-Kutta; neuroreceptor model; numerical stiffness; ODE simulation


I. INTRODUCTION
Mathematical modeling and computational simulation provide an in silico environment for investigating cerebral electrophysiology and neurological therapies including neurostimulation.Traditionally, volume-conduction models have been used to emulate electrical potentials and currents within the head cavity.In particular, these models can reproduced electroencephalograph (EEG) surface potentials [1]- [3], and have been successful in predicting cerebral current density distributions from neurostimulation administrations [1], [4]- [7].As these models become more refined, their utility in diagnosing, treating, and comprehending neurological disorders greatly increases.
Progress in field of computational neurology has motivated a migration towards models that incorporate cellular-level bioelectromagnetics.For example, bidomain based models have been used to simulate the effects of extracellular electrical current on cellular transmembrane voltage(s) [8]- [13].In addition, multiscale models have reproduced EEG measurements originating from action potentials [14], [15], and have also demonstrated an ability to simulate the influence of transcranial electrical stimulation on neuronal depolarization [16].
These models typically utilize a system of ordinary differential equations (ODEs) to emulate cellular-level electrophysiology.While the computational expense of simulating a single cell is essentially negligible, this is not the case with large-scale applications that may include hundreds of millions of cells; in multiscale applications, solving this set of ODEs is the computational bottleneck [17].In these applications, choosing an appropriate numerical solver and using efficient implementation approaches become paramount.
Alterations in neurotransmitter signalling is a hallmark of many neurodegenerative conditions and treatments.Parkinson's disease (PD), for example, which affects approximately one million individuals in the United States alone [18], culminates with pathological glutamate and γ-aminobutyric acid (GABA) binding activity throughout the basal ganglia-thalamocortical network [19], [20].As a treatment for PD, deep brain stimulation (DBS) electrically stimulates areas of the basal ganglia, such as the subthalamatic nucleus (STN) [21], to restore normal glutamate and GABA synaptic concentrations [22]- [24].Therefore, models that incorporate fundamental neurotransmitter-based signalling provide utility to the neurological research community.
Models of metabotropic and slow-responding ligand-gated receptors, such as the GABA B and N-methyl-D-aspartate (NMDA) glutamate receptors, can be efficiently solved with explicit Runge-Kutta (ERK) methods [25].On the contrary, fast-responding ionotropic receptors, such as the α-amino-3-hydroxy-5-methyl-4isoxazolepropionic acid receptor (AMPAR) and the GABA A receptor (GABA A R) result in models that are classified as stiff [26], which is an attribute of an ODE system that demands relatively small step sizes in portions of the numerical solution [27].For these ODE systems, L-stable implicit Runge-Kutta (IRK) solvers with adaptive timestepping are ideal given their exceptional stability properties [28].
In this paper, we examine L-stable IRK methods when solving models that represent the AMPA and GABA A neuroreceptors.Three L-stable IRK methods that are highly effective at solving stiff ODE systems were selected and implemented with custom Matlab [29] programming.Features including adaptive step-sizing, embedded error estimation, error-based step size selection, and simplified Newton iterations are incorporated [30].Numerical experiments were then used to identify the optimal maximum number of inner Newton iterations for each method.Then, for both the AMPAR and GABA A R models, simulation time step results of each IRK method are compared to commonly used ERK methods.In addition, the numerical accuracy and computational efficiency of each IRK method is compared to one other, as well as the highly-popular fifth order, variable step size Dormand-Prince method.Finally, a C++ based IRK implementation demonstrates exceptionally accurate and expedient performances, showcasing its potential to support large-scale multi-cellular brain simulations.

II. MATERIALS AND METHODS
A. Neuroreceptor models 1) AMPA: Glutamate is the single most abundant neurotransmitter in the human brain [31].It is produced by glutamatergic neurons, and is classified as excitatory in the sense that it predominately depolarizes post-synaptic neurons towards generating action potentials [32].Given the large concentration of glutamate in the nervous system, alterations in its production are associated with many neurodegenerative diseases and treatments.In PD patients, for example, stimulating the STN with DBS causes a cascade of cellular effects within the basal-ganglia thalamocortical pathway through its afferent and efferent projections, including increased glutamate secretion to the globus pallidus external (GPe), globus pallidus internal (GPi), and substantia nigra pars reticulata (SNr) [23].
Ligand-gated AMPA receptors for glutamate are permeable to sodium and potassium, have a reversal potential of 0 mV, and possess fast channel opening rates.Therefore, these receptors produce fast excitatory post-synaptic currents [33].Figure 1a displays the Markov kinetic binding model for the ligand-gated AMPAR that was utilized in this paper [34].In this network, there is the unbound AMPAR form C 0 , singly and doubly bound receptor forms C 1 and C 2 , which can lead to desensitized states D 1 and D 2 , respectively, and the open receptor form O [35].In addition, variable T represents neurotransmitter concentration.Mass action kinetics gives the following system of ODEs for the AMPA neuroreceptor model: State transition rates were assigned as follows: and O were set to 0 M [33], and initial values for C 0 and T were computed from a nonlinear least squares fit of the model to the whole cell recording data in Destexhe et al. [35].
2) GABA: GABA is the most abundant inhibitory neurotransmitter in the human brain [36].Like glutamate, GABA concentrations are altered by neurological disease and treatment.In STN DBS, for example, increased glutamate to the GPe increases GABA secretion to the GPi and SNr, resulting in greater GABA neuroreceptor binding in these regions [24].
There are two main categories of GABA neuroreceptors.Metabotropic GABA B receptors are slow-responding due to the secondary messenger biochemical network cascade necessary for ion channel activation.On the contrary, ligand-gated GABA A receptors are fast-responding due to their expedient ion channel opening rates.GABA A receptors are selective to chlorine with a reversal potential of approximately -70 mV.In addition, this receptor has two bound forms that can both trigger channel activation [35].
Figure 1b displays the kinetic binding model for the GABA A receptor that was utilized in this paper [26].In this model, there is the unbound receptor form C 0 , singly and doubly bound receptor forms C 1 and C 2 , slow and fast desensitized states D s and D f , and singly open and doubly open receptor forms O 1 and O 2 .This model incorporates the minimal forms needed to accurately reproduce GABA A R kinetics [37].Mass action kinetics gives the following ODE system for the GABA A neuroreceptor model: Transition rates for the GABA A R ODE system were assigned as follows: and O 2 were set 0 M, and C 0 and T were assigned the values 1 x 10 −6 M and 4096 x 10 −6 M, respectively [26].

B. Stiff ordinary differential equations
The stiffness ratio is defined as where λ i is the i th eigenvalue of the local Jacobian matrix [38], given by A general non-linear ODE system is stiff when L 1.For each neuroreceptor model, we estimated the eigenvalues numerically; a local Jacobian matrix is computed at each simulation time step using finite differences, and then its eigenvalues are computed using Matlab's eig function [39].For the AMPAR model L = 1.6 x 10 11 , and for the GABA A R model L = 3.5 x 10 11 .Thus, both of these systems are classified as stiff.

C. Implicit Runge-Kutta methods
Runge-Kutta methods are a family of numerical integrators that solve ODE systems with trial steps within the time step.These methods can be expressed with the following formulas: where ȳn is the current solution at time t n , h is the current time step, [a ij ] is the Runge-Kutta matrix, F is the ODE system, [c j ] represents intertime trial step nodes, [b j ] is the trial step solution weights, s is the number of stages, and ȳn+1 is the Biomath 4 (2015), 1512311, http://dx.doi.org/10.11145/j.biomath.2015.12.311 numerical solution at time t n+1 [28].A Runge-Kutta method can be fully defined with a Butcher table , i L-stable IRK methods are highly effective at solving stiff ODE systems [30]; these methods have no step size constraint to maintain numerical stability and quickly converge [41].Methods with second and third order accuracy were considered as these orders best match the numerical accuracy of fractional step algorithms typically employed with partial differential equation based multiscale models [16], [39].
The following L-stable IRK methods were selected for examination: SDIRK(2/1) [42], ES-DIRK23A [17], and RadauIIa(3/2) [30], [43].Each has demonstrated accuracy and computational efficiency when solving extremely stiff ODE systems.In addition, each provide an efficient local error estimator that enables error-based adaptive timestepping.For simplicity, these solvers will be referred to as SDIRK, ESDIRK, and Radau for the remainder of this paper.Butcher tables for these methods are displayed in Fig. 2. The SDIRK method is second order with an embedded first order formula for local error estimation.Each trial step, Zi , of the SDIRK solver can be solved for sequentially.Specifically, since a 12 = 0 (see Fig. 2a), the first stage of this method can be written as Z1 = h a 11 F (t n + c 1 h, ȳn + Z1 ) , and Z1 can be solved for first and used directly in the solution of Z2 = h a 21 F (t n + c 1 h, ȳn + Z1 ) + a 22 F (t n + c 2 h, ȳn + Z2 ) .
The Radau method has two stages like the SDIRK method (see Fig. 2b), but has third order accuracy with a second order error formula.This method's Runge-Kutta matrix is full, therefore the trial stages are solved as a coupled implicit system: Trial steps in the ESDIRK method are solved sequentially like the SDIRK method, after the initial explicit first stage (see Fig. 2c).This method is third order with an embedded second order formula for local error estimation, similar to the Radau solver.

D. Implementation
The three IRK methods were programmed in Matlab using principles specified in [30] and [44]; we refer these resources for a detailed explanation of Runge-Kutta method implementation and in this section provide just a brief overview of key aspects utilized in our implementations.
For each IRK method, Newton's method is used in solving system (3a).Typically, each inner Newton iteration involves computing the local Jacobian matrix and performing an LU factorization.To greatly decrease run-time, at each time step the Jacobian computation and LU factorization are performed just once on the first Newton iteration and retained for all remaining iterations.Execution time is further decreased by retaining the Jacobian in the subsequent time step if the IRK method converges with just one Newton iteration, or where k is the number of Biomath 4 (2015), 1512311, http://dx.doi.org/10.11145/j.biomath.2015.12.311 inner iterations for convergence and • is an error-normalized 2-norm [30], [45].
Newton iteration starting values are then given by: For each time step, local error is calculated and used for (i) step acceptance and (ii) subsequent step size prediction.The error at time step t n+1 can be computed by err = ŷn+1 − ȳn+1 , where The error calculations in the SDIRK and ES-DIRK methods are suitable for stiff systems [39], [41].For the Radau method, however, ŷn+1 − ȳn+1 will become unbounded and is therefore not appropriate for stiff systems [46].Instead, we use the formula err = (I − h b0 J) −1 (ŷ n+1 − ȳn+1 ) which is equivalent to where I is the identity matrix, J is the Jacobian, and b0 = √ 6 6 [46].
We can write ŷn+1 − ȳn+1 as follows [47]: To identify the coefficients e 1 and e 2 , we substitute Z1 and Z2 (3a) into (6): Collecting terms gives: From ( 5) and ( 7), we end up with the following system of equations: Using the Radau Butcher table (Fig. 2b) gives (e 1 , e 2 ) = b0 −9 2 , 1 2 .The error estimation is used to predict step size via the strategy proposed by Gustafsson [45].Further, step size following a rejected step due to excessive local error, namely err > 1, is 1 3 h.For large-scale simulations, e.g.multiscale applications, hundreds of millions of ODE systems may be solved at each time step.For these computationally intensive simulations, scripting languages such as Matlab are not ideal, and machinecompiled programs are generally necessary to achieve simulation results within reasonable computing time [48].Due to its superior accuracy in solving both the GABA A R and AMPAR models (see Sec. III), we selected the Radau method and configured a C++ implementation of it.Execution results of this version provide a measure of optimally expected computational performance.
We validated the implementation of each IRK method by comparing their GABA A R simulation results to those presented in Qazi et al. [37], and their AMPAR simulation results to whole cell recording data in Destexhe et al. [35].

E. Simulations
Numerical simulations were performed to assess the robustness of the IRK methods when solving the AMPAR and GABA A R models.Simulations were one second in duration, with rates and initial conditions as specified in Section II-A.Absolute and relative error tolerances were both set to 10 −8 , and initial step size, h, was set to 10 −4 .For each IRK method, the optimal number of maximum Newton iterations, k max , was identified by solving the AMPAR and GABA A R models with k max = 5, 6, ..., 20.For each value of k max , the mean execution time of five simulations was computed, and the value of k max that produced the lowest mean execution time was selected.Figure 3a displays the k max values selected for each model and method.
For each method, it was observed that a threshold value of k max exists, such that higher values do not result in faster simulations.Therefore, we selected the minimum k max value associated with the fastest execution speed.For example, for the Radau method solving the GABA A R model, simulation times begin to plateau for k max ≥ 10, and simulation times with k max ≥ 15 were the same (see Fig. 3b).Therefore, for this model and IRK method, k max = 15 was selected.
Figure 3b also shows that faster run times correlate with fewer solution time steps and LU factorizations, until a floor is reached; in the case of the Radau method solving the GABA A R model, this floor is 29 time steps and 30 LU factorizations.To a point, higher values of k max increase the probability of Newton method convergence, resulting in fewer time steps and fewer computationally expensive LU factorizations [30].For the Radau method solving the GABA A R model, values of k max ≥ 15 yield the fewest number of simulation times steps in addition to no steps where the Newton iteration fails to converge.Thus, when k max = 15, time steps and associated LU factorizations are minimized, yielding the fastest execution speeds.To evaluate the advantages that IRK methods have when solving fast-responding neuroreceptor models, we first compare the total number of simulation time steps and simulation step sizes of each IRK method to the following commonly used ERK methods: forward euler (FE), midpoint method (Mid), and 4 th order Runge-Kutta (RK4).Next, to compare each IRK method to one another and to the adaptive 5th order Dormand-Prince method (DP5) [49], metrics including local and global error, total simulation time steps, step sizes, execution times, and numbers of Jacobian computations and LU factorizations were evaluated.Absolute and relative error tolerances of the DP5 method were set to 10 −8 , matching the tolerances of the three implicit methods.
To more comprehensively assess performance differences amoung the IRK methods, workprecision diagrams using solution run times and scd values, where scd = -log 10 ( relative error at t = 1.0 sec ∞ ), were then generated [50].For the work-precision diagrams, relative error tolerances Biomath 4 (2015), 1512311, http://dx.doi.org/10.11145/j.biomath.2015.12.311 ) , m = 0, 1, ..., 25, absolute error tolerance was set to 10 −4 • rtol, and initial step size was 10 −4 .In addition, solution run times presented in these diagrams are the mean of five runs.For all accuracy calculations, solutions with a 5 th order adaptive time-stepping L-stable implicit Runge-Kutta method with a maximum step size of 10 −6 and both absolute and relative tolerances set to 10 −14 were used as true solutions.
Finally, the execution time of the Radau C++ implementation when solving both neuroreceptor models was assessed.All simulations were run on a Linux machine with an Intel i7 processor with a clock speed of 2.40 GHz.

A. GABA A R Model
Figure 4 presents the solution of the GABA A R model with the SDIRK method; ESDIRK and Radau solutions look identical.The sharp transition in the total open state concentration, O 1 (t) + O 2 (t), at the onset of neurotransmitter stimulus at t = 0 displays the necessity for smaller time steps in this region of the solution (Fig. 4a).Upon examining all receptor forms during the first 1.5 ms of the simulation, it is observed that both the unbound closed form, C 0 (t), and total bound closed form, C 1 (t) + C 2 (t), possess concentration transitions even greater than the open receptor form (Fig. 4b).These results show the stiffness possessed by the GABA A R system.
Table I displays simulation time step metrics for the three IRK methods and the FE, Mid, and RK4 ERK methods.The maximum step size of each explicit method was calculated with the GABA A R model stiffness index and the method's stability region [28], giving the largest step that can be taken while maintaining numerical stability.Then, the number of time steps required for each ERK method was computed by dividing the simulation duration by the maximum step size.The FE and Mid methods both require 2.1 x 10 4 time steps, and the RK4 method requires 1.5 x 10 4 , which is lower than the FE and Mid methods due to its larger stability region [30].On the contrary, each implicit method requires less than 30 simulation time steps.As displayed in Figure 4a, the majority of these time steps for the SDIRK method occur at the beginning of the simulation, within the region of rapid solution transition.
Similarly, the ESDIRK and Radau solvers demand noticeably more time steps at the onset of neurotransmitter stimulation (Fig. 5).Rejected steps, totalling three for the ESDIRK method (Fig.  ) and two for the Radau method (Fig. 5b) all occur at time t = 0; once the solution in this region has been accurately resolved, no further rejected steps occur.In addition, for all three IRK methods, all Newton iterations converged, which was facilitated by identifying optimal k max values (see Sec. II-E).Further, the smallest step sizes of the IRK methods, namely 1.4 x 10 −5 for the SDIRK and ESDIRK methods and 1.2 x 10 −5 for the Radau method, have the same order of magnitude as the largest stable step sizes of the ERK methods.
Next the accuracy and computational efficiency of the IRK methods were compared to one another and with the DP5 method (Table II).While the DP5 method possesses the lowest maximal local true solution deviation (3.2 x 10 −10 ), the 2-norm of its global error is one to two orders of magnitude higher than all three IRK methods.These results are explained by the fact that the solution of the DP5 solver oscillates around the true solution (Fig. 6).In addition, the DP5 method requires approximately 50,000 simulation time steps and takes 49.0 seconds to run.In comparison, the  [51], however, these approaches will result in even greater run times.
The Radau method has the greatest execution time of the three IRK methods, at 0.69 seconds.While the number of simulation time steps amoung the IRK methods are comparable, two factors contribute to the longer run time of the Radau method.First, this solver generally requires a greater number of iterations for Newton's method to converge (Fig. 3a).Second, the Radau method requires 30 Jacobian computations, versus just four for the SDIRK and ESDIRK methods.Despite its run time disadvantages amongst the IRK methods, the accuracy of the Radau method stands out as superior.It has the lowest global error 2-norm (8.7 x 10 −10 ), and its maximal deviation from the true solution (−3.7 x 10 −10 ) is comparable to that of the 5 th order DP5 method, the only IRK method examined where this is the case.Further, the Radau method has greater accuracy at every time step than both the SDIRK and ESDIRK methods.
These findings are reinforced by the workprecision diagram for the three IRK methods when solving the GABA A R model (Fig. 7).This diagram highlights the higher precisions attained by the third order methods, and in addition, also confirms the slower execution speeds achieved by the Radau method.However, when comparing graph points of similar relative tolerances, such as the symbols marked in yellow that represent rtol = 10 −6 , the Radau method is consistently more accurate.

B. AMPAR Model
Figure 8 presents solution results of the AMPAR model solved with the Radau method.Like the GABA A R model, the rapid transition in the open state concentration upon neurotransmitter stimulation demands a greater number of time steps (Fig. 8a).Specifically, the first 10% of the simulation (0.1 sec) encompasses approximately 96% of the simulation time steps.Once beyond this initial region, step size eventually increases by seven orders of magnitude (Fig. 8b).Similar to the GABA A R model, both unbound closed and bound closed forms contribute to the system's stiffness.
A noticeable difference, compared to the GABA A R simulation results, is the number of time  steps needed by the implicit methods to solve the AMPAR model.The Radau method, for example, requires 199 time steps (Fig. 8b), a 586% increase from the 29 steps needed to solve the GABA A R model.Similar increases are observed with the SDIRK and ESDIRK solvers, most notably the 531 steps required by the SDIRK method (Table III).In addition, the smallest step sizes of the IRK methods are two orders of magnitude lower with the AMPAR model (Fig. 8b), due to the stiffness index of the AMPAR system [27].Despite the elevated simulation time step counts, each IRK method still outperforms the explicit methods (Table III); maximum stable step sizes and simulation time steps for the explicit methods were again computed with their stiffness indices and stability regions [28].
While greater k max values eliminated nonconvergent Newton iterations in the GABA A R model, this is not the case with the AMPAR model.Each IRK method has two instances where Newton's method did not converge.In addition, the SDIRK method has four rejected steps, and the ESDIRK and Radau methods each have two, all occurring at time t = 0.    performed by these solvers.Figure 9a displays the Jacobian computations and LU factorizations of the SDIRK and ESDIRK methods.Each method has a near identical number of LU factorizations, however, the ESDIRK method requires 51 Jacobian computations, which is more than double the 24 performed by the SDIRK method.In addition, the Radau method requires 162 Jacobian computations.Therefore, despite having a lower number of simulation time steps, the computational advantages of the ESDIRK and Radau methods are diminished due to this elevated number of Jacobian computations.
Once again, the accuracy and computational performances of the IRK methods were compared to the DP5 method (Table IV).As observed with the GABA A R model, the DP5 method has inferior execution performance, requiring 1.1×10 5 simula-  The work-precision diagram for the AMPAR model (Fig. 10) again confirms the higher precision achieved by the third order ESDIRK and Radau solvers.More noticeable in this graph are the differences in the "slopes" of the curves, where "flatter" curves, i.e.ESDIRK and Radau, have more precision per unit CPU time [30].For the AMPAR model, the Radau method is slower than the ESDIRK method at all work-precision tolerances examined, yet at relative tolerances greater than 10 −6 , the Radau method becomes faster than the SDIRK method.Further, the Radau method is generally the most accurate of all three IRK methods.

C. C++ Radau Implementation
The Radau method consistently demonstrates the greatest accuracy of the methods examined, however, its main disadvantage is execution speed.For this reason, we selected the Radau method and configured a C++ implementation of it.Table VI displays execution times for the previous Radau Matlab implementation, as well as the new C++ version.
As expected, the C++ version is significantly faster.Specifically, the GABA A R model has a 99.6% decrease in execution time, and the AM-PAR model has a 99.7% decrease in execution time.Because the implementation algorithms between the two versions are the same, the C++ version maintains the accuracy of the Matlab prototype.

IV. CONCLUSIONS
Computational neurology is a valuable contributor in the diagnosis, treatment, and comprehension of neurological disease.To provide maximal utility to the scientific community, computational simulations should incorporate highlydetailed, neurotransmitter-based neuron models.Therefore, large-scale simulations involving populations of neurons will inevitably produce computational challenges.In this paper, we have shown that appropriate numerical solvers with efficient implementation strategies can alleviate computational difficulties.
Commonly used explicit methods are capable in solving a limited number of fast-responding ligand-gated neuroreceptor models.However, we have shown that poor stability properties make them non-ideal for large-scale applications.Rather, by addressing the stiffness possessed by these models, we show that implicit methods are highly advantageous.In particular, we demonstrate that L-stable implicit Runge-Kutta methods offer superior accuracy and run-time efficiency compared Biomath 4 (2015), 1512311, http://dx.doi.org/10.11145/j.biomath.2015.12.311 to their explicit siblings when solving biologicallybased AMPA and GABA A neuroreceptor models.To accelerate solutions, we utilize a range of strategies including embedded error estimators and simplified Newton iterations.In addition, we show that optimal execution times are achieved when costly Jacobian computations and LU factorizations are minimized.
The third order Radau IRK method demonstrates exceptional local and global accuracy compared to all other explicit and implicit methods examined.In addition, its numerical stability properties yield a relatively low number of simulation time steps and efficient step sizes when solving the AMPA and GABA A neuroreceptor models.Further, a C++ implementation of the Radau solver displays the computational faculty to enable largescale multi-cellular simulations.In future work, we plan to continue our investigation of numerical solvers for neurotransmitter-based neuron models by comparing the IRK methods to multi-step methods and exponential integrators.

Fig. 5 :
Fig. 5: Simulation step sizes for the GABA A R model.

Fig. 6 :
Fig. 6: GABA A R model open state concentration solution error.

TABLE I :
Simulation time steps results for the ERK and IRK methods when solving the GABA A R model.

TABLE II :
Accuracy and simulation run-time metrics of the DP5 and IRK methods when solving the GABA A R model.Boldface font denotes best results of each column.

TABLE IV :
Accuracy and simulation run-time metrics of the DP5 and IRK methods when solving the AMPAR model.Boldface font denotes best results of each column.

TABLE V :
Number of rejected and nonconvergent steps for each IRK method when solving the AMPAR model.error 2-norm (3.3 × 10 −8 ), which is again larger than those of the three IRK methods.The Radau method once again has the lowest global error 2-norm (1.6×10 −8 ) of all methods inspected.

TABLE VI :
Run times (seconds) for the Matlab and C++ Radau method when solving the GABA A R and AMPAR models.