Analysis and implementation of an efficient solver for large-scale simulations of neuronal systems

(1)

Analysis and implementation of an

efficient solver for large-scale

simulations of neuronal systems

M A T T H E W T H E

(2)

(3)

Analysis and implementation of an

efficient solver for large-scale

simulations of neuronal systems

M A T T H E W T H E

Master’s Thesis in Numerical Analysis (30 ECTS credits) Master Programme in Computer simulation

for Science and Engineering (120 credits) Royal Institute of Technology year 2013

Supervisor at CSC was Michael Hanke Examiner was Michael Hanke TRITA-MAT-E 2013:33 ISRN-KTH/MAT/E--13/33--SE

Royal Institute of Technology School of Engineering Sciences

KTH SCI

SE-100 44 Stockholm, Sweden

(4)

(5)

Abstract

(6)

(7)

Referat

Analys och implementering av en effektiv

lösare for storskaliga simuleringar av

neuron-system

(8)

(9)

Introduction

1.1 Background and formulation of the problem

The human brain is built up of approximately 1011_{(100 billion) neurons,}

communi-cating using about 1014_{(100 trillion) synapses. Modelling and simulation of (parts}

of) this huge system could provide fundamental insights in brain functioning and could eventually lead to treatments for brain diseases.

As more and more computing power becomes available in these modern times, larger and larger systems can be analysed. To maximize this potential, efficient ways to solve the equations arising from neuron modelling are necessary, while at the same time keeping a good eye on the stability and accuracy of the methods.

Several software packages have been made available to assist in neuron mod-elling and simulation, e.g. NEURON[15], NEST[14], MOOSE[13] and GENESIS[6]. Several numerical integration methods are used in these packages, but these might not always be the optimal choice for simulation.

In the Master’s thesis of Maya Brandi[3], critical limitations in step size were met when using Hines’ method[10], a popular numerical integrator in neuron simu-lation. Its major drawback is that it is a constant step size solver, which can result in extremely long computation times if a tiny step size is necessary to maintain stability.

NEURON currently can make use of the BDF method as a variable step size solver. Even though this method has reasonable stability characteristics, it is not specially tailored to handle neuron simulations. The question remains if numerical solvers are available that make better use of the structure of the equations arising in neuron simulation.

1.2 Aim of the project

(12)

Using simple, but representative neuron models, these methods will be imple-mented in Matlab for analysis and a first comparison to other available solvers. Optimization steps will be taken regarding computational complexity, accuracy and numerical stability.

The final goal is to implement an efficient solver in one of the neuron simulation software packages and compare the results with existing solvers by using models from the respective database.

1.3 Scope and limitations

(13)

Chapter 2

Theory

2.1 Neuronal system simulation

A brain consist of a huge number of cells called neurons, connected to each other in enormous networks. The neurons communicate by sending electrical signals to each other through axons (sending signals away from the cell body) and dendrites (receiving signals from other neurons). The signals come in the form of so-called "ac-tion potentials" or "spikes", short electric pulses lasting in the range of milliseconds. Under certain conditions a neuron which receives this signal will then propagate this signal to other neurons. Sometimes a neuron spikes without stimulation from other cells and other times signals from multiple neurons at once are needed to make a single neuron spike.

In the study of the dynamics of neural networks two crucial issues arise: • which model should be used to describe the spiking dynamics of each neuron? • how should the neurons be connected into a network?

Concerning the first issue, many different kinds of spiking behaviour have been observed over the years and different models have been proposed to capture these behaviours. Two of the most notable models are the Hodgkin-Huxley model[11] and the integrate-and-fire model (and its variations). While integrate-and-fire models are much simpler and cheaper to simulate than Hodgkin-Huxley models, they lack a certain complexity to capture all different spiking behaviours. Therefore, for the most realistic simulation, Hodgkin-Huxley models are still the first choice [12]. The Hodgkin-Huxley model also has the advantage that it has a more direct relation to biophysical quantities (e.g. membrane conductance), which makes it easier to create and interpret models as well as interpreting the simulation results.

(14)

2.1.1 Hodgkin and Huxley model

Hodgkin and Huxley presented a model in the 1950s, representing neurons by its electrical properties and formulating an equivalent electrical circuit (see Fig. 2.1) [11]. Using the well established theory of electrical circuits, an ODE for the mem-brane potential V can be formulated:

Cm

dV

dt = gN a(EN a− V) + gK(EK− V) + gleak(Eleak− V) + Iinj (2.1)

Since currents sum linearly, this model can easily be extended to include other types of ions.

gK gNa gleak

Cm

EK ENa Eleak

Iinj

Figure 2.1. The electrical circuit equivalent to a neuron. The membrane is

repre-sented by a capacitance Cm, the ion channels are represented by conductances gN a

and gK, the gradients of ion concentration which drive the flow through the ion

chan-nels are represented by voltage sources EN a and EK. The other ion channels are

collected into a leak channel with a fixed conductance gleakand gradient Eleak. Not

depicted here is an injected current which can stimulate the cell to generate spikes.

Hodgkin and Huxley observed that the ion conductances gN a and gK depend

on the membrane potential. More specifically, they could be represented by a com-bination of "gates" that could be opened and closed. Some gates opened at high membrane potentials and closed at low membrane potentials (activation gates) and other gates did the exact opposite (inactivation gates). By fitting it to experimental data Hodgkin and Huxley came up with the following equations:

gN a = ¯gN am3h, gK = ¯gKn4

where m and n are activation gates and h is an inactivation gate and ¯gN aand ¯gK are

the maximum conductances for respectively the sodium and potassium channels. Gate variables are handled by simple first order kinetics (p = n, m, h):

dp

dt = αp(V )(1 − p) − βp(V )p (2.2)

where p is the probability that the gate is open, αp is the opening rate and βp is

the closing rate. αp and βp can be determined by fitting to experimental data.

2.1.2 Cable equation

(15)

2.2. NUMERICAL METHODS FOR NEURON SIMULATION

axons in finite time, an extension of this model is needed to allow for spatial prop-agation of currents. As a simplification we can consider all the components of the neurons as 1 dimensional cables with a certain electrical resistance. Using Ohms law a PDE can be derived for the membrane potential, better known as the cable equation: 1 Ra ∂2V ∂x2 = Im = − V Rm + Cm ∂V ∂t

with Im the membrane current, Ra the resistance of the intracellular fluid, Rm the

membrane resistance and Cm the membrane capacitance. _RV_m must be replaced by

the membrane currents (the right hand side of (2.1)).

A spatial discretisation can be performed on ∂2_V

∂x2 to transform this PDE into an

ODE: Cm dVi dt = Vi−1−2Vi_{+ V}i+1 Ra +gN a(EN a

−V)+gK(EK−V)+gleak(Eleak−V)+Iinj

(2.3) or allowing for varying axial resistance:

Cm dVi dt = Vi−1− Vi Ri−1a +Vi+1− Vi Ri a

+gN a(EN a−V)+gK(EK−V)+gleak(Eleak−V)+Iinj (2.4) By combining the Hodgkin-Huxley model with the cable equations in this way, one can describe both the temporal and spatial kinetics of one or several neurons.

2.2 Numerical methods for neuron simulation

The following mathematical characteristics were observed for the ODE system of (2.4) combined with gate equations (2.2):

• the system is stiff, i.e. transitions occur on different time scales • there exist eigenvalues close to the imaginary axis

• the voltage derivatives are linear in the voltage variables and the gate deriva-tives are linear in the gate variables:

˙Y = f1(Z)Y + f2(Z) (2.5)

˙Z = g1(Y ) + g2(Y )Z (2.6)

where Y is the vector containing the voltages and Z is the vector containing the gates. This also holds for most other HH-type models.

(16)

Method Stability Order reduction

stiff system Linear systems perstep Implicit RK methods (Gauss,

Radau) A-stable 2s/2s − 1 → s Newton iteration withsystem of size ns

Diagonal Implicit RK or

ROW methods A-stable (for right choice ofparameters)

s+1 → unknown s systems of size n

BDF methods unbounded stability domain

no A-stability for s > 2 none (stays s) Newton iteration onsystems of size n Implicit Peer methods A-stable (for right choice of

parameters) none (stays s−1) s Newton iteration onsystems of size n Linearly Implicit Peer

meth-ods A-stable (for right choice ofparameters) none (stays s−1) s systems of size n

Table 2.1. Overview of implicit numerical integrators and their characteristics.

Peer methods are not as well known as Runge-Kutta and BDF methods and a short introduction will be given in section 2.2.2.

Some nuances have to be given to Table 2.1:

• Even though there is no order limit for A-stability for (linear) implicit peer methods ([20], p.348), the methods used in practice are sometimes not A-stable, because other conditions (e.g. consistency) have to be fulfilled first. • According to [2] linearly implicit peer methods fail for certain problems with

crude tolerances (10−2_,₁₀−3_).

• The s Newton iterations on systems of size n for implicit peer methods might not be much worse than the s linear systems of size n of the linear variant. However, according to [2] there is not much difference in number of steps and accuracy (at least for finer tolerances, see previous comment) and therefore the linearly implicit peer method is preferred.

The two most promising methods are the BDF method and the linearly implicit peer method. The BDF method is already implemented in the NEURON environ-ment and an impleenviron-mentation of the linearly implicit peer method was programmed in NEURON for this project as well to provide a comparison.

2.2.1 BDF methods

Introduced by Curtiss and Hirschfelder in 1952 [5], backwards differentiation formu-las (BDF) methods are currently the most popular methods to solve stiff differential equations. The short overview here is based on the treatment in the books of Hairer, Nørsett and Wanner [8], [9].

(17)

here yn+1 is the unknown, fm = f(tm, ym) and αj and βj are coefficients of a

particular linear multistep method.

For BDF, one sets βj = 0 for j 6= k, i.e. only the function value at the (unknown)

yn+1 is used. Since fn+1 = f(tn+1, yn+1) ≈ y0(tn+1), we can approximate this

function value by numerical differentiation.

For constant step size this approximation can be done using backward differ-ences: k X j=1 1 j∇ j_y n+1= hfn+1 (2.8)

where ∇j _{is the backward difference operator, recursively defined by}

∇0_y

n+1= yn+1, ∇j+1yn+1 = ∇jyn+1− ∇jyn. (2.9)

Variable step sizes lead to step size dependent coefficients αj.

In both cases the BDF method leads to an implicit equation in yn+1 (nonlinear

if f(t, y) is nonlinear), which is usually solved by a Newton-like procedure. Often the Jacobians needed for these Newton iterations are reused over multiple steps to reduce computation cost.

If LU decompositions are used with variable step size, the computation cost can be reduced even more by forcing the step size to be constant over multiple steps (quasi-constant step size). This allows the reuse of the LU decomposition.

From the numerical differentiation, one can see that if the yn−k+1 up to yn+1

are used, the method is consistent of order k (k + 1 nodes, therefore an error of

O(hk+1_{)). These are then called the BDF method of order k or simply BDF-k}

It can be seen that for k ≥ 7, the method is unstable and that for k > 3 the method cannot be A-stable (second Dahlquist barrier prevents all linear multistep methods to be A-stable for k > 2). The BDF methods for k ≤ 6 are however A(α)-stable with

k 1 2 3 4 5 6

α 90◦ 90◦ 86.03◦ 73.35◦ 51.84◦ 17.84◦ (2.10)

Because of the small angle α for BDF6, it is usually not used. The stability domains for the BDF methods are given in Fig. 2.2. The higher order methods have a relatively small part of the imaginary axis included in their stability domain and therefore have problems with eigenvalues close to the imaginary axis.

(18)

−20 −10 0 10 20 −20i −15i −10i −5i 0 5i 10i 15i 20i BDF1 BDF2 BDF3 BDF4 BDF5 BDF6

Figure 2.2. The stability domains of BDF1-6. Only BDF1 and BDF2 are A-stable,

the A(α) angle decreases for higher order BDF methods.

2.2.2 Linearly implicit peer methods

Peer methods are a relatively new family of numerical integrators (introduced in 2004, see [18]) that in a way combine the ideas of Runge-Kutta methods with linear multistep methods. The overview below is based on the treatment in [20] and [2].

Peer methods use s stages between tn and tn+1 where an approximate solution

is calculated, comparable to a Runge-Kutta method. Instead of throwing this in-formation away for the next step (as is the case in RK methods), they are used to

calculate the approximate solution for the next step between tn+1 and tn+2 as well.

Since there is no real difference in treatment of the different stages (they are treated as "peers"), all the stage solutions are of the same order, which prevents order reduction for stiff problems, as seen in RK-methods.

The equations for the stage variables of an s step implicit peer method can be written as: Ym,i= s X j=1 bijYm−1,j+ hm i X j=1 gijFm,j (2.11)

(19)

A convenient way to write (2.11) is by using the Kronecker/tensor product ⊗:

Ym = (Bm⊗ I)Ym−1+ hm(Gm⊗ I)F (tm, Ym) (2.12)

where Ym= (Ym,1, Ym,2, . . . , Ym,s)T and F (tm, Ym) similarly formed with the f(tm+

hmci, Ym,i). Note that Bm= (bij)si,j=1 and Gm = (gij)si,j=1 may depend on the step size ratio σm = hm/hm−1, and that Gm= (gij)si,j=1 is lower triangular.

For the simplest notation, we consider a scalar autonomous ODE, which can easily be generalized back to the tensor notation:

Ym= BmYm−1+ hmGmF(Ym) (2.13)

For the normal (also called "full") implicit peer method, this equation would be solved with s Newton iterations (1 for each stage). For variable step sizes, the maximum consistency order that can be reached with s stages is p = s − 1. These

consistency conditions fix the matrix Bmin terms of the matrix Gm. The coefficient

matrix Gm and the stage locations vector c can be chosen freely. s(s − 1)/2 degrees

of freedom in the matrix Gm are used to satisfy zero-stability, since zero-stability

and consistency of order p guarantee convergence of order p. The remaining degrees

of freedom in Gm and c can be used to obtain superconvergence p = s at constant

step size, large angles α for A(α)-stability and/or a small norm of the local error term.

In the singly-implicit peer method, the diagonal of Gm is taken to be a

con-stant γ. This reduces the degrees of freedom, but will allow the reuse of the LU decomposition in the different stages.

If we take the starting vector for the Newton iterations Y0

m,i to be estimates of

the Ym,i of order s − 1, one Newton iteration is enough to keep consistency of order

p= s − 1 [16]. This is known as the linearly implicit peer method.

Since the linearly implicit peer method is our method of choice, some extra notes should be made:

• For the linearly implicit peer method, it is not necessary to use the exact Jacobian in the (single) Newton iteration of (2.13) to maintain consistency of order p = s − 1 [16]. Comparable to Rosenbrock W-methods (ROW) the

Jacobian at the last stage of the previous step J(Ym−1,s) is used for all stages

in one step. For even less computation, the Jacobian can be reused during multiple steps, which resembles the normal W-methods.

• It would be possible to also use the Fm−1,j in (2.11), at the cost of extra

(20)

• In [2] it was noted that the linearly implicit peer method performs about the same in number of steps as the fully implicit peer method, but is faster because of its simplicity. However, for some problems the linearly implicit peer method failed for crude tolerances.

Applying the test equation y0 _{= λy to a (linearly) implicit peer method, the}

stability matrix is given by

M(z) = (I − zGm)−1Bm (2.14)

and the domain of absolute stability is given by

S = {z ∈ C : ρ(M(z)) < 1} (2.15)

where ρ(M(z)) is the spectral radius of M(z). The second Dahlquist barrier does not apply to peer methods and it is possible to construct A-stable peer methods of high order, as given in [20]. In practice, other characteristics might be more important, such as a low error constant, and giving up A-stability can be a good trade-off. Since the domain of absolute stability depends on the chosen parameters, even for constant step size, there is no unique domain shape. For our implemented version of the linearly implicit peer method, the domain will be presented in Sec. 4.1.2.

Similar to linear multistep methods, multiple initial values are necessary, usually generated by some 1-step method (BDF1, ROW).

2.2.3 System decoupling

As mentioned before, neuronal equation systems often have a quasi-linear form:

˙Y = f1(Z)Y + f2(Z) (2.16)

˙Z = g1(Y ) + g2(Y )Z (2.17)

where Y contains the voltage variables and Z contains the gate variables. The Jacobian then becomes:

J = " f1(Z) ∂f1_∂Z(Z)Y +∂f_∂Z2(Z) ∂g1(Y )Z ∂Y + ∂g2(Y ) ∂Y g2(Y ) # (2.18) Implicit numerical integrators use the Jacobian in the linear system that needs to be solved, in the form

(I − γhJ)x = b. (2.19)

Since we often work with a large number of voltage and gate variables, solving the linear system can quickly become the bottleneck of fast computation. For dense

systems solving the system is of order O(n3_{) and for sparse systems this can often}

(21)

Approximating the Jacobian by the block diagonal matrix

J∗ = " f1(Z) 0 0 g2(Y ) # (2.20) would reduce solving a linear system of size n to 2 smaller systems of roughly size

n/2. Winning a factor of 8 for dense solvers.

According to first experiments, the BDF method has a significant increase in the number of steps and Jacobian evaluations using this approximation, whereas the linearly implicit peer method performs equal to working with the exact Jacobian.

For the linearly implicit peer method, using this block diagonal Jacobian can be written as: I − hmγ " f1(Zm−1,s) 0 0 g2(Ym−1,s) #! " ∆Ym,i ∆Zm,i # = − " Ym,i Zm,i # + s X j=1 bij " Ym−1,j Zm−1,j # + hm i−1 X j=1 gij "

f1(Zm,i)Ym,i+ f2(Zm,i)

g1(Ym,i) + g2(Ym,i)Zm,i

#

(2.21) i.e. in split form:

(I − hmγf1(Zm−1,s))∆Ym,i= −Ym,i+ s X j=1 bijYm−1,j+ hm i−1 X j=1

gij(f1(Zm,i)Ym,i+ f2(Zm,i))

(I − hmγg2(Ym−1,s))∆Zm,i= −Zm,i+ s X j=1 bijZm−1,j+ hm i−1 X j=1

(22)

(23)

Chapter 3

Models

In order to compare the BDF method with linearly implicit peer methods, some models were selected and implemented in Matlab and NEURON. We tried to select representative models for different classes of neurons with its distinct behaviours.

3.1 Matlab

3.1.1 Hodgkin-Huxley model with 3 compartments

In the Master’s thesis of Maya Brandi [3] a 3 compartment model was proposed to model a single neuron. The 3 compartments represent the soma, dendrite and spine. Current is injected at the soma, resulting in an action potential travelling through the dendrite and arriving at the spine. At the soma, sodium and potassium ion channels are present, and at the spine, calcium ion channels exist.

Together, the current and gate equations produce a nonlinear ODE system of 9 equations. The parameters were regulated to produce a certain realistic calcium concentration in the spine, but besides this, the model does not try, nor claim, to model an actual neuron accurately.

The model presented below is taken ad verbatim from [3], the values of the parameters can also be found there. The model parameters are given in Table 3.1:

The HH-equations combined with the cable equation are:

(24)

C_mi , Ri_m, V_mi capacitance, resistance and membrane potential of the respective

compartments

m, n, h, s, r probability of the respective gate being open

KCa, K_Ca∗ number of respectively inactive and active potassium channels

Eion Nernst equilibrium potential for the respective ion

¯gion maximum conductance of the respective ion channel

Iinj injected current at the soma

αp, βp voltage dependent opening and closing rates of the respective gate

Table 3.1. Model parameter descriptions of the 3 compartment model in [3]

where the currents are given by

IN a = (EN a− Vm1)¯gN am3h (3.5) IK = (EK− Vm1)¯gKn4 (3.6) ICa= (ECa− Vm3)¯gCas2r (3.7) IKCa = (EK− Vm3)¯gKCa K_Ca∗ KCa+ K_Ca∗ (3.8)

Note that this model fits the quasi-linearity as mentioned in (2.17). The membrane

potential over time for the soma (V1

m) is shown in Fig. 3.1 0 0.02 0.04 0.06 0.08 0.1 −0.08 −0.06 −0.04 −0.02 0 0.02 0.04 V1 (V) time (s)

Figure 3.1. Membrane potential vs. time at the soma (Vm1), current injection is

(25)

3.1. MATLAB

3.1.2 Minimal Hodgkin-Huxley models for different neuron classes

In [17] the authors tried to capture the behaviour of the 4 most prominent classes of neurons (fast spiking, regular spiking, intrinsically bursting and low-threshold spik-ing) in a minimal number of HH-type equations. The models were initially imple-mented in NEURON by the authors (however, without variable step size support), but because of its simple nature, we implemented them in Matlab for comparison.

The model parameters were fitted to a large set of experimental data, both in vitro as in vivo, of different animals and different parts of the brain. Only single neurons were modelled every time using just 1 compartment, but these provide the essential building blocks for large network simulations.

The model presented below is taken ad verbatim from [17] (with a different but equivalent description of the p gate), the values of the parameters can also be found there. The model parameters are given in Table 3.2:

Cm, V capacitance and membrane potential of the neuron

m, n, h, p, q, r, u, s(V ) probability of the respective gate being open

gleak, Eleak resting membrane conductance and potential

Eion Nernst equilibrium potential for the respective ion ¯gchannel maximum conductance of the respective channel

IN a, IKd sodium and potassium currents responsible for action potentials

IM slow potassium current responsible for spike-frequency adaptation

IL, IT respectively high and low threshold calcium current

Iinj injected current at the soma

αz(V ), βz(V ) voltage dependent opening and closing rates of the respective gate

Table 3.2. Model parameter descriptions of the minimal HH-type models in [17]

The HH-equations are given by

Cm

dV

dt = −gleak(V − Eleak) − IN a− IKd− IM − IT − IL− Iinj (3.9)

dz

dt = αz(V )(1 − z) − βz(V )z, z ∈ {m, n, h, p, q, r, u} (3.10)

where the currents are given by

IN a = (V − EN a)¯gN am3h (3.11)

IKd = (V − EK)¯gKdn4 (3.12)

IM = (V − EK)¯gMp (3.13)

IL= (V − ECa)¯gLq2r (3.14)

IT = (V − ECa)¯gTs(V )2u. (3.15)

Note that this model does not completely fit the quasi-linearity of (2.17), because

of the nonlinear s(V ) term in IT. However, most of the neuron classes have ¯gT = 0,

(26)

The different neuron classes can be expressed by this equations by including and excluding certain currents:

Neuron class IN a IKd IM IL IT

Regular spiking x x x

Fast spiking x x

Intrinsically bursting x x x x

Regular spiking x x x x

Table 3.3. Present currents for the different neuron classes, an "x" indicates that

the particular current is present.

In Fig. 3.2, 3.3 and 3.4 the spiking behaviour of the different classes are shown.

0 200 400 600 800 1000 −80 −60 −40 −20 0 20 40 60 Regular Spiking V (mV) time (ms) 0 200 400 600 800 1000 −80 −60 −40 −20 0 20 40 60 Fast Spiking V (mV) time (ms)

Figure 3.2. Left: Voltage over time for the regular spiking neuron model, current

(27)

3.1. MATLAB 0 500 1000 1500 2000 2500 3000 −100 −50 0 50 Intrinsically Bursting V (mV) time (ms) 0 500 1000 1500 2000 2500 3000 −100 −50 0 50

Intrinsically Bursting (Repetitive bursting)

V (mV)

time (ms)

Figure 3.3. Voltage over time for the intrinsically bursting neuron model,

cur-rent injection takes place between 500 and 2500 ms. Left: the membrane potential shows an initial burst, followed by a regular spiking pattern. Right: larger L-type conductance causes multiple bursts.

0 200 400 600 800 1000 −100 −50 0 50 V (mV) time (ms)

Low Threshold Spiking (Hyperpolarised)

0 200 400 600 800 1000

−100 −50 0 50

Low Threshold Spiking (Depolarised)

V (mV)

time (ms)

Figure 3.4. Voltage over time for the low threshold spiking neuron model, current

(28)

3.2 NEURON

One of the cornerstones of the NEURON environment is an easy way to build up physical models and adapt these to the accuracy needs of the user. For example, it is easy to define different compartments of the neuron (e.g. soma, dendrites, spines) and connect them. If more accuracy is needed within one of these compartments, it is possible to split it up into multiple compartments, called segments, without having to redefine all the parameters (they are automatically recomputed by NEURON). This segmentation occurs very often in NEURON models and it is this discretisation mechanism that can quickly lead to large ODE systems.

3.2.1 Minimal Hodgkin-Huxley

In the web based NEURON tutorial [7] the basic Hodgkin-Huxley equations are modelled in NEURON. We only used the simplest single compartment neuron with sodium and potassium channels, resulting in a system of 4 ODEs (1 voltage and 3 gate variables). The model reproduces the basic spiking behaviour, but has not yet adapted its parameters to physically relevant values. The spiking behaviour is shown in Fig. 3.5. 0 50 100 150 200 250 300 −80 −60 −40 −20 0 20 40 60 time (ms) voltage (mV)

Figure 3.5. Membrane potential vs. time, current injection is present between 100

and 200 ms.

3.2.2 Pyramidal neuron

(29)

3.2. NEURON

apical dendrite. The soma and basal/apical dendrites each consist of 3 segments, resulting in a total of 10 segments. The spiking behaviour is shown in Fig. 3.6.

It is difficult to say how many variables are actually used, since some of the variables seem to be unused because of certain conditions. The total number of variables is 222, but only about 70 components in the solution vector are non-zero.

0 200 400 600 800 1000 −70 −60 −50 −40 −30 −20 −10 0 10 time (ms) voltage (mV)

Figure 3.6. Membrane potential vs. time at the soma, short synaptic pulses are

given at 200, 300, 400 and 500 ms. Only the first pulse leads to an action potential

3.2.3 Reduced models for pyramidal neurons

In [1] an attempt was made to fit experimental data of layer 5 pyramidal neurons to models with only a few compartments. The idea behind it is that these reduced models, compared to fully detailed compartmental models, allow faster fitting to ex-perimental data and also reduce computation time when simulating large networks. They produced 10 sets of parameters that all managed to fit the experimental data accurately. There was a significant difference between these 10 parameter sets, even though they all managed to fit the same data, leading to the hypothesis that it is mainly the ratio between certain parameters that cause a certain behaviour.

There are 7 compartments with in total 20 segments, and 86 membrane variables, leading to a total of 106 equations.

Three different behaviours of the pyramidal neurons were simulated: • Regular spiking by current injection to the soma (Fig. 3.7)

• Calcium spiking by a small current injection, combined with dendritic stimu-lation (Fig. 3.8, left)

(30)

0 200 400 600 800 −80 −60 −40 −20 0 20 40 60 time (ms) voltage (mV) Regular Spiking

Figure 3.7. Membrane potential vs. time at the soma, current injection takes place

between 100 and 600 ms. 0 50 100 150 200 250 300 −80 −60 −40 −20 0 20 40 60 time (ms) voltage (mV) Calcium Spiking I 0 50 100 150 200 250 300 −80 −60 −40 −20 0 20 40 60 time (ms) voltage (mV) Calcium Spiking II

Figure 3.8. Voltage over time at the soma due to dendritic stimulation (starting

(31)

3.2. NEURON

3.2.4 Signal propagation in a neural network

In [4] benchmarks were created to evaluate different neural network simulation tools. One benchmark was based on a neural network for signal propagation from [22]. The model consisted of a large network of 4000 neurons interconnected by 32000 synapses. Originally it was modelled with integrate-and-fire neurons but [4] con-verted this to a network with HH-type neurons.

A part of the neurons is randomly stimulated during the first 100 ms, after this the network becomes self sustaining with neurons that keep firing irregularly without additional stimulation (see Fig. 3.9).

Because the 4000 neuron network was computationally very expensive for dense solvers, we considered a scaled down network of 200 neurons with 600 synapses here instead. 200 was the smallest number of neurons consistently reproducing the same type of behaviour as the 4000 neuron network.

0 100 200 300 400 500 600 700 800 900 1000 0 20 40 60 80 100 120 140 160 180 200 time (ms) neuron #

Figure 3.9. Spike times of a 200 neuron network with 600 synapses, a few neurons

(32)

(33)

Chapter 4

Implementation

Both in the BDF method as well as the linearly implicit peer method certain strate-gies have to be chosen for:

• step size control • Jacobian reuse

• quasi-constant step size enforcement • order control (BDF)

• coefficient choice (linearly implicit peer method)

Some choices effect the efficiency of the method more than others, and sometimes a "suboptimal" choice has to be made as a result of constraints of the implementation environment.

4.1 Matlab

4.1.1 BDF method

In [19] the ODE suite of Matlab is extensively documented and the treatment below is just a summary of the relevant information presented there.

The BDF method is implemented in the function ode15s, which however defaults to a related family of formulas, the so-called Numerical Differentiation Formulas (NDF). The NDF methods try to address the poor stability of higher order BDF methods. The allowed step size is increased, but the A(α)-stability angle decreases (Table 4.1).

Another interesting aspect is that instead of the solution vectors yn−k, . . . , yk, the backward differences yn, ∇yn, . . . , ∇kynare saved. This leads to a Nordsieck-like scheme (see Sec. 4.2.1).

The advantage of this scheme is that it is very easy to recompute the solution

(34)

order 1 2 3 4 5

BDF stability angle 90◦ ₉₀◦ ₈₆◦ ₇₃◦ ₅₁◦

NDF stability angle 90◦ ₉₀◦ ₈₀◦ ₆₆◦ ₅₁◦

step ratio percent 26% 26% 26% 12% 0%

Table 4.1. Comparison of BDF and NDF methods used in ode15s. The NDF

methods trade a bigger step size ratio percent ((NDF step - BDF step)/BDF step) for a smaller stability angle

with old step size h. This allows us to use fixed coefficients with variable time steps, as if it was a constant step size scheme. Otherwise, it would be necessary to

recompute the coefficients αj from (2.7) in each step, which depend on the previous

and current step sizes hn−k, . . . , hn.

The strategies for step size control, order control and Jacobian reuse are intri-cately interwoven, also shown in Fig. 4.1:

• the step size is kept constant for kmax+2 steps, unless there is a problem with

convergence of the Newton iteration.

• if the convergence rate of the Newton iteration falls below 0.9 or the solution has not converged after the maximum number of iterations, the Jacobian is checked. If it is not current, the Jacobian is recomputed and the step is tried again. If it is current, the step size is reduced to 0.3 · h, or, if the maximum iterations was reached, the order can be reduced and integration continues with a newly computed step size.

• the maximum number of Newton iterations is 4.

The chosen strategy is based on the assumption that Jacobian computation is relatively cheap (because of the nature of problems that are usually computed with Matlab). If this is not the case, the method might take a long time to compute.

The quasi-constant step size enforcement could be less stable than a fully vari-able step size implementation, but according to [19], the quasi-constant step size implementation has been quite satisfactory in practice.

(35)

4.1. MATLAB

calculate predictor and set up linear system

do Newton iteration

converged?

conv. rate < 0.9 & max it. not

reached?

Jacobian current?

recompute Jacobian

h reused

< kmax+ 2 times?

calculate new step size and order

max it. reached?

multiply step size with 0.3 no yes no no yes no yes yes no yes

Figure 4.1. Step size control, Jacobian reuse, order control and quasi-constant step

(36)

4.1.2 Linearly implicit peer method

The implementation of the linearly implicit peer method is quite more straightfor-ward than the BDF implementation. We have adapted the code provided by Steffen Beck, Rüdiger Weiner, Helmut Podhaisky and Bernhard Schmitt from the Univer-sity of Halle, which has been used in [2]. For our experiments we have used the coefficients of the superconvergent singly-implicit peer method (SSIP) ([2], Table 3).

The A(α) stability of this coefficient choice is given in Table 4.2, where a com-parison with the BDF methods of the same order was made. For the same s, quite an improvement in the angle α is made. However the superconvergence p = s of the peer method only holds for constant step sizes and reduces to p = s − 1 for variable step sizes, therefore actually a comparison between s − 1 for the peer method and

sfor BDF should be made. In this case, the angles are quite comparable.

s 2 3 4 5

BDF stability angle 90.0◦ _86.0◦ _73.4◦ _51.8◦

SSIP stability angle 86.1◦ _83.2◦ _75.7◦

Table 4.2. Comparison of the A(α)-stability of the BDF (ode15s) and

superconver-gent singly-implicit peer methods (SSIP) [2]. It should be noted that the comparison is slightly misleading, since the convergence order of the peer methods here is p = s−1 for variable step sizes (p = s holds only for constant step sizes), whereas the BDF methods converge with order p = s also with variable step sizes.

Also a study of the actual stability domains was conducted, presented in Fig. 4.2, with more detailed pictures for specific s in Fig. 4.3, 4.4 and 4.5. The stability domains show a definite advantage of the SSIP methods compared to the BDF methods for eigenvalues close to the imaginary axis with negative real part.

−20 −10 0 10 20 −20 −15 −10 −5 0 5 10 15 20 SSIP3−5 SSIP3 SSIP4 SSIP5 −20 −10 0 10 20 −20 −15 −10 −5 0 5 10 15 20 BDF2−5 vs SSIP3−5 BDF2 BDF3 BDF4 BDF5 SSIP3 SSIP4 SSIP5

Figure 4.2. Left: The stability domains of SSIP-5. Right: The stability domains of

(37)

4.1. MATLAB −1 −0.5 0 0.5 1 −10 −5 0 5 10 BDF2 BDF3 SSIP3

Figure 4.3. Comparison of stability domains of BDF2, BDF3 and SSIP3. Note that

the axis are not equal. If variable step sizes are used, SSIP3 should be compared to BDF2, in which case SSIP is worse on every level. For constant step sizes SSIP3 should be compared to BDF3. In this case there is better stability for eigenvalues close to the imaginary axis with negative real part (better A(α)-stability).

−1 −0.5 0 0.5 1 −10 −5 0 5 10 BDF3 BDF4 SSIP4

(38)

−1 −0.5 0 0.5 1 −10 −5 0 5 10 BDF4 BDF5 SSIP5

the axis are not equal. If variable step sizes are used, SSIP5 should be compared to BDF4. For constant step sizes SSIP5 should be compared to BDF5. The same comments hold as for s = 4 in Fig. 4.4.

On recommendation of the authors, a special initial vector for the Newton iter-ation was used, of the form

Y_m,i(0) = ΘYm−1+ hmG0Fm, (4.1)

where Θ is a full matrix and G0 is strictly lower triangular.

For the estimation of the local error term, also on recommendation of the

au-thors, the solution at the final stage Ym,swas estimated by polynomial extrapolation

of order s − 1 by using the first s − 1 stage solutions:

∆m= Ym,s− Ps−1(tm,s) = y(tm,s) + O(hsm) − y(tm,s) − O(hs−1m ) = O(hs−1m ), (4.2)

where Ps−1(t) is the extrapolating polynomial. Then the WRMS norm is taken for

the local error estimation:

err= v u u t 1 n n X i=1 ∆m,i

Atoli+ |Ym−1,s,i|Rtoli

. (4.3)

If err < 1, the step was accepted, otherwise rejected.

As is common in variable step size methods, the new step size (either for the next step, or to retry the current step) is then estimated by

hnew= h · min{facmax, max{facmin, fac · (1/err)1/(s−1)}}. (4.4)

Here facmax = 1.5, facmin = 0.2, fac = 0.8 were chosen as safety factors.

(39)

4.2. NEURON ENVIRONMENT

The strategy for Jacobian reuse was very basic: update the Jacobian if it was used for more than x (e.g. 20) steps, or if the step was rejected. This seems rather arbitrary, but worked all right in practice.

The strategy for quasi-constant step size enforcement was also very basic: keep-ing the step size constant, if the estimated next time step was between h and h · T , where T is a certain threshold depending on the order of the method.

As initial solver both ode15s as well as an own implementation of a Rosenbrock W-method with dense output were used. numjac was used to calculate the numerical

Jacobian. The initial step was taken to be hinitial= min{

√

reltol,10−4}. This is a rather naive approach, trusting the step size control to quickly reach the optimal step size. However, since (in our Matlab models) this only happens at the start and no integration restarts occur, this works fine in practice.

4.2 NEURON environment

The NEURON code is written in C and C++. For variable step sizes it currently uses the SUNDIALS suite [21], which can handle both ODEs (CVODES package) and DAEs (IDA package).

Some DAEs arise, because zero-capacitance compartments are used. In many cases it is possible to rewrite these equations to ODEs and solve them with CVODES, but in other cases this is not possible and IDA has to be used.

The main difference with the Matlab implementation will be the computation of the Jacobian. Since there is no direct equivalent for the Matlab function numjac in C, the Jacobian will differ.

In fact, the default linear solver in NEURON does not form the Jacobian ex-plicitly, but instead makes use of the sparse pattern to only do the necessary com-putations. In this stage, the Jacobian is estimated in block diagonal form, similar to (2.20). There might, however, be some differences to the Jacobian of (2.20), because the Jacobian of the gate equations is estimated by a diagonal Jacobian by default. Furthermore, this strategy cannot make use of Jacobian nor LU decompo-sition reuse, but does make excellent use of the sparsity pattern.

It is possible to change the solver to use the dense solver CVDense, also included in the SUNDIALS suite. Choosing this allows Jacobian and LU decomposition reuse, but does not effectively make use of the sparse pattern. The numerical Jaco-bian calculation is less accurate than numjac, and it does not have the possibility to make use of the sparse pattern for a decreased number of function evaluations.

It is also possible to use CVDiag, which approximates the Jacobian by its diag-onal. However, this is generally not an efficient solver.

In general relatively crude tolerances are used, which seem to be sufficient for

its purposes. The default is abstol = 10−2 _{and no relative tolerance; in articles}

(40)

4.2.1 BDF method

The BDF method implemented in the NEURON environment is taken from the CVODES package from the SUNDIALS suite [21]. It uses the Nordsieck implemen-tation described in [20], p. 152 & p. 331.

This scheme saves the derivatives of y(tm) in a so-called Nordsieck vector:

z{m} = y(tm), hy0(tm), h2 2 y 00_(t m), . . . , hk k!y (k)_(t m) !T + O(hk+1_). _(4.5)

This reduces changing the step size, to multiplying with a diagonal matrix:

z_new[m+1]= diag 1,hnew

h , h2_new h2 , . . . , hk_new hk ! z[m+1] (4.6)

instead of having to recompute the coefficients, which can become quite expensive if the step size is changed often.

The strategies for step size control, Jacobian reuse, order control and quasi-constant step size enforcement are slightly different from the Matlab implementa-tion, mainly because Jacobian evaluations are assumed to be very expensive and will therefore be avoided as much as possible. Furthermore the (numerical) Jacobian considerations mentioned before should be taken into account.

The strategies for step size control, order control and Jacobian reuse are

intri-cately interwoven, also shown in Fig. 4.6 (γ = h · βk, where h is the step size and

βk as in (2.7)):

• the step size is updated if:

– more than 20 steps have been taken since the last update

– the change in γ is more than 0.3γ.

– the Newton iteration failed to converge

– the step size was rejected

• the Jacobian is updated if:

– more than 50 steps have been taken since the last evaluation

– Newton iteration failed to converge and the Jacobian was not up to date

and the change in γ was less than 0.2γ.

– Newton iteration failed to converge and forced a step size reduction.

• the order q is changed if:

– at least q + 1 steps have been taken and the change in step size with the

new order will be bigger than 1.5.

• the maximum number of Newton iterations is 3.

(41)

4.2. NEURON ENVIRONMENT

calculate predictor and set up linear system

do Newton iteration

converged?

conv. rate < 0.5 & max it. not

reached?

multiply step size with 0.25

Jacobian current?

recompute Jacobian

error test passed?

h reused > 20 times?

calculate new step size and order

jac reused > 50 times no yes yes no no no yes no no yes yes yes

Figure 4.6. Step size control, Jacobian reuse, order control and quasi-constant step

(42)

4.2.2 Linearly implicit peer method

The Matlab implementation was rewritten into C with some minor changes: • As initial solver only the Rosenbrock W-method (LSTAB, fourth order) with

dense output is used, since interaction with the existing BDF code was difficult due to the Nordsieck form.

• The initial step size calculation underwent a minor change. The NEURON environment restarts the integration every time current injection is started or stopped, and in certain configurations also whenever a spike is encountered. This becomes a problem with the naive initial step size implemented in Matlab (hinitial = min{

√

reltol,1e − 4}), especially for the integration restarts when

spikes are encountered. The most natural way to resolve this was to take the last step size before integration was restarted, which is reasonable for restarts due to spiking, but probably sub-optimal for restarts due to current injection. Furthermore, since we are actually providing an initial step for the Rosenbrock W-method, the previous step size might not be suitable at all. For future implementation, a scheme as in [8], p. 169 could be considered instead.

(43)

Chapter 5

Analysis and Results

5.1 Spectrum analysis

In order to be able to explain some of the simulation results, it is useful to consider the eigenvalues of the systems we will encounter. This will only be done for the small problems in Matlab, since this would be quite expensive for large systems and because a qualitative treatment of HH-type equations should provide sufficient information.

Large positive eigenvalues and high frequency oscillations are usually the cause of a forced step size reduction in variable time step methods, or instabilities in constant step size methods.

(44)

5.1.1 Hodgkin-Huxley model with 3 compartments

For the simulations, we have used: Abstol = Reltol = 10−8 _{and computed the}

Jacobians at regular intervals by interpolating the solution vectors. The following observations were made:

• The largest and smallest two eigenvalues (by absolute value) are real

– The largest is practically constant at −3.5 · 106_.

– The second largest eigenvalue varies between −6500 and −2800

– The two smallest eigenvalues are constant at −5.

• There is a real eigenvalue oscillates between −800 and −600

• The remaining 4 eigenvalues intertwine into complex pairs (see Fig. 5.1 and 5.2).

– 1 complex pair follows a big bow-like path in the complex plane

– 2 complex pairs follow small ellipses in the complex plane, both around

−130 + 0i.

– 1 of the eigenvalues attains a large positive real value (up to 1300).

−1000 −500 0 500 1000 1500 −1000 −800 −600 −400 −200 0 200 400 600 800 1000 real part imaginary part Complex Plane

Figure 5.1. The middle 5 eigenvalues in the complex plane over one period. There

(45)

5.1. SPECTRUM ANALYSIS 0.02 0.025 0.03 0.035 0.04 −1000 −500 0 500 1000 1500

Real Part vs. Time

time

real part

big bow−like pattern ellipses around −130+0i

Figure 5.2. The real part of the middle 5 eigenvalues over one period. They

intertwine and follow two ellipses in the complex plane around 22 and 32 ms, and follow a big bow-like pattern in the complex plane between 23 and 37 (see Fig. 5.1).

The huge difference between −3.5 · 106 _{and −5 characterizes the stiff behaviour,}

which is known to exist in these HH-type models. The crossing of the imaginary axis by the big bow-like path is also characteristic of HH-type models. The small ellipses in the complex plane will most likely only cause minor oscillations and are not particularly important.

These observations hold during current injection. In the periods before and after current injection, only real negative eigenvalues occur, which regulate the transition to equilibrium potential.

5.1.2 Minimal Hodgkin-Huxley models for different neuron classes

For the simulations, we have used: Abstol = Reltol = 10−8 _{and computed the}

Jacobians at regular intervals by interpolating the solution vectors. The following observations were made:

Regular spiking neurons

• There is a real eigenvalue which oscillates between −10 and −35 • There is a real eigenvalue which oscillates between −0.003 and −0.2

(46)

– 2 complex pairs follow ellipses around −2 and −4.

−25 −20 −15 −10 −5 0 5 10 15 −10 −5 0 5 10 Complex Plane real part imaginary part

Figure 5.3. The 3 intertwining eigenvalues of the RS neuron model in the complex

plane at the start of an action potential. There are 3 patterns with a non-zero imaginary part. Two ellipses, one around −2 and one around −4 and a large bow-like pattern, which crosses the imaginary axis in a very particular manner.

158 158.5 159 159.5 160 160.5 161 161.5 162 −25 −20 −15 −10 −5 0 5 10 15

time (ms)

real part

start of bow−like pattern end of bow−like pattern 2 ellipses

Figure 5.4. The real part of the 3 intertwining eigenvalues of the RS neuron model

at the start of an action potential. They intertwine and follow two ellipses in the complex plane around 159 ms, and follow a big bow-like pattern in the complex plane between 159.5 and 161 (see Fig. 5.3).

(47)

5.1. SPECTRUM ANALYSIS

different units are used here). Once again, we see the crossing of the imaginary axis with a big bow-like path and the minor ellipses.

Different from the 3 compartment model is that the bow-like shape only occurs in a tiny part of the period (1.5 ms in a 140 ms period), also the shape is very different.

Fast spiking neurons

Practically everything is similar to the case of the regular spiking neurons, except for the small real negative eigenvalue, which disappears (which relates to the slow potassium current).

Intrinsically bursting neurons

• There is a real eigenvalue which oscillates between −10 and −30

• 2 real eigenvalues are almost constant, one at −0.002 and another at −0.2. • The other 6 eigenvalues intertwine and form complex pairs around the start

of the action potential, but are negative and real during most of the period (see Fig. 5.5 and 5.6).

– about 4 or 5 complex pairs follow ellipses with centres between −1 and

−4.

The behaviour copies certain parts of the RS neurons, but much more com-plex pairs occur. They all have small imaginary parts however, resulting in slow oscillations lasting only for short intervals.

Low threshold spiking neurons

(48)

−25 −20 −15 −10 −5 0 5 10 15 −10 −5 0 5 10 Complex Plane real part imaginary part

Figure 5.5. The 6 intertwining eigenvalues of the IB neuron model in the complex

plane at the start of an action potential. There are many patterns with a non-zero imaginary part. Multiple ellipses, with centres between −1 and −4 and small imaginary part. Then there is a large bow-like pattern, which crosses the imaginary axis in a very particular manner.

160.5 161 161.5 162 162.5 163 −25 −20 −15 −10 −5 0 5 10 15

time

real part

Figure 5.6. The real part of the 6 intertwining eigenvalues of the IB neuron model at

(49)

5.2. LINEARLY IMPLICIT PEER METHODS

5.2 Linearly implicit peer methods

In order to make the linearly implicit peer method competitive with existing stiff solvers, strategies to reuse Jacobian evaluations and LU decompositions were needed. Practically all BDF method implementations contain these kind of strategies, but these are sadly not easily transferable to the linearly implicit peer methods, since there is no Newton iteration where we can check the convergence rate.

We will try out different strategies and optimize the parameters for the models implemented in Matlab, in the hope that these will transfer reasonably well to other problems in this family. More weight was given to the 3 compartment model of Maya Brandi, since the models of Pospischil did not contain multiple compartments, which takes away the characteristics of the cable equation.

All simulations were run on a Dell Optiplex 780 equipped with a Intel Core2 Quad Q9550 processor running at 2.83GHz and 3.7 GiB RAM, running Ubuntu 12.04.

5.2.1 Fully implicit vs. linearly implicit

The first justification we have to make is to replace the Newton iteration using convergence checks, with a single iteration using only a local error check at the end. In order to do this, we compared the fully implicit implementation from [2]

peer_028_exact to our linearly implicit implementation peerli.

Chosen parameters and strategies:

• superconvergent singly implicit peer method coefficients (SSIP)

• tolerances 1e−k_{, k}_{= 2, . . . , 6 (both absolute and relative tolerance were used:}

abstol= reltol = T OL)

• full Jacobian

• no Jacobian reuse (and therefore no LU decomposition reuse) • ode15s as initial solver

The simulation statistics for the 3 compartment model can be found in Tables 5.1 and 5.2. The function evaluation count fcn includes function evaluations needed for Jacobian evaluations, unless mentioned otherwise.

(50)

## peer_028_exact s3

# 3 Compartment model, Full Jacobian, StSol ode15s, 0x Jacobian reuse #---# TOL err time nacc rej fcn idid Jac LU #---1.00E-02 3.15E-02 0.1354 59 14 1194 1 59 73 1.00E-03 2.79E-03 0.2226 123 33 1817 1 123 156 1.00E-04 2.01E-04 0.4737 322 15 3872 1 322 337 1.00E-05 5.38E-06 1.3349 988 11 10730 1 988 999 1.00E-06 1.54E-07 3.7545 3113 1 28796 1 3113 3114 ## peerli s3

# 3 Compartment model, Full Jacobian, StSol ode15s, 0x Jacobian reuse #---# TOL err time nacc rej fcn idid Jac LU #---1.00E-02 9.56E-02 0.0837 63 12 629 1 64 75 1.00E-03 1.58E-03 0.1542 120 28 1190 1 121 148 1.00E-04 1.38E-04 0.3687 321 14 2962 1 322 335 1.00E-05 4.79E-06 1.1002 987 12 8953 1 988 999 1.00E-06 1.53E-07 3.4362 3112 1 28058 1 3113 3113 =========================================================================== ## peer_028_exact s4

# 3 Compartment model, Full Jacobian, StSol ode15s, 0x Jacobian reuse #---# TOL err time nacc rej fcn idid Jac LU #---1.00E-02 2.30E-02 0.0892 60 14 682 1 61 74 1.00E-03 8.25E-04 0.1631 109 33 1266 1 113 142 1.00E-04 5.95E-05 0.2956 212 43 2323 1 213 255 1.00E-05 2.64E-06 0.5480 424 27 4382 1 425 451 1.00E-06 1.43E-07 1.1126 896 14 9063 1 897 910

Table 5.1. Simulation statistics for the 3 compartment model of fully implicit peer

(51)

## peer_028_exact s5

Table 5.2. Simulation statistics for the 3 compartment model of fully implicit peer

methods (peer_028_exact) and linearly implicit peer methods (peerli) using the SSIP parameters and s = 5. The same comments hold as for s = 3, 4 in Table 5.1.

For all s, the linearly implicit variant performs similar in terms of number of steps and global error, and performs significantly better in terms of number of function evaluations. It seems that the extra Newton iterations are not necessary to achieve the same stability and accuracy.

One note that has to be made, is that this result depends greatly on the

approx-imation of the solution vector that starts the Newton iteration Y(0)

m,i = ΘYm−1 +

hmG0Fm. Using a naive prediction, like Ym,i(0) = Ym,i−1, does lead to a significant increase in number of steps.

We can observe the order of convergence of p = s − 1 from these results as well:

for a 10× decrease of tolerance there is roughly an increase of 101/(s−1)_×_{in number}

of steps (3.2×, 2.2×, 1.8× for s = 3, 4, 5 respectively).

For the minimal HH-type models of Pospischil, similar behaviour was seen for

s = 3 and s = 4. For s = 5, there was however an increase of 40% in number of

steps for tol = 10−2 _{for RS and IB neurons, but the number of function evaluations}

stayed roughly the same. For finer tolerances, no significant increase in number of steps was seen however.

5.2.2 Jacobian reuse

(52)

peer methods, which is the Jacobian reuse over the different stages within one step. Jacobian reuse is also necessary to allow for LU decomposition reuse.

We have tested multiple values of maximum Jacobian reuse, with the other method parameters kept the same as in the previous section. Part of the simulation statistics for the 3 compartment model can be found in Table 5.3 and 5.7.

## peerli s3

Table 5.3. Simulation statistics for the 3 compartment model of linearly implicit

(53)

For the 3 compartment model with s = 3, Jacobian reuse works very well, managing a decrease in function evaluations of 40% for the crudest tolerance up to 65% for the finest tolerance. For s = 4 and s = 5 the decrease in number of function evaluations is less, mainly because there are less steps in total, but also because the rejection rate is higher for these methods. The decrease in number of function evaluations is shown in Fig. 5.7 (3 compartment model), 5.8 (Pospischil IB neurons) and 5.9 (Pospischil LTS neurons). The Pospischil RS neurons show similar behaviour to the 3 compartment model and the FS neurons to the IB neurons.

0 5 10 20 40 80 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 s = 3

max Jacobian reuse

ratio fcn wrt no Jac reuse

tol = 1e−2 tol = 1e−3 tol = 1e−4 tol = 1e−5 tol = 1e−6 0 5 10 20 40 80 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 s = 4

max Jacobian reuse

0 5 10 20 40 80 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 s = 5

max Jacobian reuse

Figure 5.7. Ratio of number of function evaluations for different maximum Jacobian

reuses w.r.t. no Jacobian reuse for the 3 compartment model (the RS neuron model from Pospischil shows similar results). Note that the x-axis scale is not linear. The lower the order and the finer the tolerance, the more profit from Jacobian reuse, because more steps are taken.

In all of the models, we can see that from around a maximum of 10 Jacobian reuses, the relative function evaluation decrease is rather low because the function evaluations needed for the stage evaluations are dominating the ones needed for Jacobian evaluation. This could however shift for larger systems, since the number of function evaluations needed for Jacobian evaluation is proportional to the system size.

(54)

0 20 50 100 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 s = 3

max Jacobian reuse

tol = 1e−2 tol = 1e−3 tol = 1e−4 tol = 1e−5 tol = 1e−6 0 20 50 100 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 s = 4

max Jacobian reuse

0 20 50 100 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 s = 5

max Jacobian reuse

reuses w.r.t. no Jacobian reuse for the Pospischil IB neuron model (FS neurons show similar results). Note that the scale on the x-axis is not linear. Low order and fine tolerance still seem to be the leading factor to decrease the function evaluations. For tol = 10−2 there is however an unusually big decrease, probably the result of problems in the simulation with no Jacobian reuse. Also, for s = 5 the number of function evaluations rises again if the maximum Jacobian reuse is increased, which is caused by an increase in the number of steps.

0 20 50 100 0 0.2 0.4 0.6 0.8 1 1.2 1.4 s = 3

max Jacobian reuse

tol = 1e−2 tol = 1e−3 tol = 1e−4 tol = 1e−5 tol = 1e−6 0 20 50 100 0 0.2 0.4 0.6 0.8 1 1.2 1.4 s = 4

max Jacobian reuse

0 20 50 100 0 0.2 0.4 0.6 0.8 1 1.2 1.4 s = 5

max Jacobian reuse

(55)

5.2.3 Quasi-constant step size enforcement

If we keep the step size constant over multiple steps, while the Jacobian is kept the same as well, the matrix of the linear system does not change and the LU decomposition can be reused (assuming these are used).

Keeping the step size constant when only a small increase in step size would be realized (lower than a threshold T ), is relatively risk free, since the local error will only decrease. The downside is that more steps might be needed. Also the step size could get stuck at a certain value (see Fig. 5.10) if the Jacobian is not renewed often enough and the threshold is too high.

0 0.02 0.04 0.06 0.08 0.1 −0.08 −0.06 −0.04 −0.02 0 0.02 0.04 time(s) Membrane potential (V) 0 0.02 0.04 0.06 0.08 0.10 1 2 3 4 5 x 10−4 Step size V1 step size

Figure 5.10. Step size plot of the 3 compartment model using the linearly implicit

peer method with s = 4 and T = 1.55 and a maximum of 100 Jacobian reuses, showing that the threshold was chosen too high. For a high threshold and many Jacobian reuses, the step size could get stuck at a low value, due to a fixed step size and/or Jacobian that are sub-optimal.

For the 3 compartment model with a 10 maximum Jacobian reuse, we have tested different thresholds. The results are displayed in Table 5.4 for s = 4. For

s= 3 and s = 5 similar behaviour is seen, but with different values of the threshold.

Higher values of maximum number of Jacobian reuses (20× and 40× were tested) seem of little influence to the threshold value.

The largest allowable threshold before the number of steps significantly increases is T = 1.4, 1.45, 1.5 for s = 3, 4, 5 respectively. For the Pospischil models the optimal thresholds were roughly the same for RS and FS neurons. For IB and LTS neurons, a small but significant increase in steps could be seen, especially for higher orders

and cruder tolerances for these values. Except for tol = 10−2_{, the number of LU}

Analysis and implementation of an efficient solver for large-scale simulations of neuronal systems

Analysis and implementation of an

efficient solver for large-scale

simulations of neuronal systems

M A T T H E W T H E

Analysis and implementation of an

efficient solver for large-scale

simulations of neuronal systems

M A T T H E W T H E

Abstract

Referat

Analys och implementering av en effektiv

lösare for storskaliga simuleringar av

neuron-system

Contents

Chapter 1

Introduction

1.1

Background and formulation of the problem

1.2

Aim of the project

1.3

Scope and limitations

Chapter 2

Theory

2.1

Neuronal system simulation

2.2

Numerical methods for neuron simulation

Chapter 3

Models

3.1

Matlab

3.2

NEURON

Chapter 4

Implementation

4.1

Matlab

4.2

NEURON environment

Chapter 5

Analysis and Results

5.1

Spectrum analysis

5.2

Linearly implicit peer methods