Numerical simulation of solitons in the nerve axon using finite differences

Examensarbete 30 hp

September 2014

Numerical simulation of solitons

in the nerve axon using finite

differences

Teknisk- naturvetenskaplig fakultet UTH-enheten Besöksadress: Ångströmlaboratoriet Lägerhyddsvägen 1 Hus 4, Plan 0 Postadress: Box 536 751 21 Uppsala Telefon: 018 – 471 30 03 Telefax: 018 – 471 30 00 Hemsida: http://www.teknat.uu.se/student

Abstract

Numerical simulation of solitons in the nerve axon

using finite differences

Jonatan Werpers

A High-order accurate finite difference scheme is derived for a non-linear soliton model of nerve signal propagation in axons. Boundary conditions yielding well-posed problems are suggested and included in the scheme using a penalty technique. Stability is shown using the summation-by-parts framework for a frozen parameter version of the non-linear problem.

1 Introduction

Understanding the propagation of nerve signals is of great concern in neu-roscience. Specifically, the behavior of the signal in the axon of the nerve cell has been of great interest. In 1952 Hodgkin and Huxley introduced a mathematical model for this phenomena [4]. The model describes the nerve signal as a voltage pulse modeled by a set of coupled partial and ordinary differential equations. It builds on properties of specific proteins in the nerve cell. Since it was introduced, the model has provided insight into how neu-rons behave. However, it has also been made clear that the Hodgkin-Huxley model does not accurately describe certain phenomena observed in exper-iments. Some examples are the effects of anesthesia, temperature induced nerve pulses and thickness changes with pulse propagation. To try and ad-dress these discrepancies researchers are looking for new ways to describe the propagation of nerve signals. On approach has been to treat the signal as a thermodynamic phenomenon.[1] This has the potential of a more general description of what happens when a pulse propagates through a nerve cell, including the voltage potential described be the Hodgkin-Huxley model as well as several other phenomena show in experiments. Numerical simulation help in understanding and developing these models.

In the present study the soliton model in [2] will be examined. Its un-derlying principle is the wave equation for sound and it models the nerve pulse as a density wave in the lipid membrane of the nerve cell. To account for thermodynamic effects the model is non-linear and contains high order derivatives. Numerical simulation of such problems are known to require highly accurate and stable numerical methods to be useful. For accurate solutions high order accurate methods are required to avoid a prohibitively small grid spacing. The high orders of derivatives in the present model amplifies this in addition to complicating the boundary treatment. The finite difference method provided by the summation-by-parts framework ful-fill both of these requirements while also providing a flexible way to handle boundary closures. [6]

The SBP-SAT methods combines operators for discretizing the space dimension that satisfy certain conditions [5], with the simultaneous approx-imation term for implementing boundary conditions.[3] It does this in a way that help provide stable schemes. Typically the SBP-SAT method has been applied to problems involving first and second order differentiation [10, 7, 8]. Recently SBP operators was derived for third and fourth order derivatives [9]. The present study applies these new SBP operators to a non-linear problem as well as introduces a new kind of compatibility between SBP operators.

nu-merical schemes are developed and shown to be stable. Time integration is described in section 6 and convergence properties of the schemes are shown in 7. Section 8 show and discuss the results of a number of numerical simu-lations. The work is summarized in section 9.

2 The equation

The Soliton model describes the nerve signal as a partial differential equation in time and one dimension of space. Its basis is the wave equation for area density changes, ⇢A_, @2 @⌧2 ⇢ A₌ @ @z  c2 @ @z ⇢ A _{, (2.1)}

where ⌧ is time, z is position, c is the speed of sound and ⇢A_{is the density}

offset from the equilibrium. The sound velocity in the lipid membrane is considered a function of density, and is represented by a truncated power series,

c2 = c2₀+ p ⇢A+ q ⇢A 2, (2.2)

where p and q are parameters determined by measurements. To account for dispersion observed in experiments an extra term h@4

@z4 ⇢A is added. We

finally arrive at following mathematical model for the density of the lipid membrane in the axon,

@2 @⌧2 ⇢ A₌ @ @z ⇣ c2₀+ p ⇢A+ q ⇢A 2⌘ @ @z ⇢ A _h @4 @z4 ⇢ A_{, (2.3)} where • ⇢A_{= ⇢}A _⇢A 0

• ⇢A _{- Lateral density of the membrane}

• ⇢A

0 - Equilibrium density in the fluid phase

• c0 - Velocity of small amplitude sound

• p, q - Experimentally determined parameters determining the density dependence of sound velocity

2.1 Dimensionless version

By introducing the variable changes

u = 1 ⇢A 0 ⇢A, x = c0 hz, t = c2₀ p h⌧, B1= ⇢0 c2₀p, B2 = ⇢2₀ c2₀q, (2.4) we obtain the following dimension less version of (2.3),

utt = (B(u)ux)_x uxxxx, (2.5)

where

B(u) = 1 + B1u + B2u2, (2.6)

and subscripts denote differentiation. The parameters B1 and B2 are

deter-mined experimentally in [2] to the values B1 = 16.6and B2 = 79.5.

2.2 Analytical soliton solution

It is possible to find analytical soliton solutions to (2.5) in the infinite line as shown in [2] By introducing the change of variables

⇠ = x t (2.7)

and assuming solutions of the form u(x, t) = u(⇠) solutions are restricted to being of soliton form and one can arrive at an analytical soliton solution. One such soliton solution is given by

u(⇠) = 2a+a (a++ a ) + (a+ a ) cosh(⇠ p 1 2₎ (2.8) where a_±= B1 B2 1_± s 2 2 0 1 2₀ ! (2.9) and where satisfies 0 < | | < 1 where 0 = 0.649851. In this report

soliton solutions of this kind will be called

(x, t) (2.10)

3 Wellposedness

Well posedness of a problem requires the existence of a unique solution and that there is a bound on the growth of the solution. Since it is in general very hard to show existence and uniqueness of a solution for a non-linear problems we will skip this step and assume that there is a unique solution. This is reasonable since we will be using the minimum required number of boundary conditions. In this section two problems with two kinds of boundary condi-tions will be posed and shown to have a bounded energy growth. The highest order derivative determines the number of boundary conditions necessary to get a well posed problem. In this case the highest order derivative is 4, and so we need two boundary conditions on each boundary. The first of the two problems is the Dirichlet-Neumann problem. It can be consider the most physically natural condition on the boundary as it corresponds directly to physical properties of the nerve. The Dirichlet-Neumann boundary condition amounts to setting both u and ux on the boundary. Let g0, h0, g1, h1 2 R

be constants or functions of time, and f1, f2 be functions of x, then the full

continuous problem becomes 8 > > > > > > < > > > > > > : utt = (B(u)ux)x uxxxx, x2 [0, 1], t > 0 u = g0, ux = h0 x = 0 u = g1, ux = h1 x = 1 u = f1 x2 [0, 1], t = 0 ut= f2 x2 [0, 1], t = 0 . (3.1)

The second of the two problems is the Characteristic problem where charac-teristic boundary conditions are set. This type of boundary condition can be considered analytically natural as it arises from the analysis of the boundary interaction. In this case the full problem becomes

8 > > > > > > > > > > > < > > > > > > > > > > > : utt= (B(u)ux)_x uxxxx, x2 [0, 1], t > 0 utx uxx= g0 x = 0 ut B(u)ux+ uxxx = h0 x = 0 utx+ uxx= g1 x = 1 ut+ B(u)ux uxxx = h1 x = 1 u = f1 x2 [0, 1], t = 0 ut= f2 x2 [0, 1], t = 0 , (3.2)

where =p1 + B2_{(u). The non-linearities of these models not only makes}

To establish well posedness the PDE, equation (2.5), has to be associ-ated with an energy measure of the solution. This is done using the energy method. First, introduce the inner product

(f, g) = Z 1

0

f (x)g(x)dx (3.3)

and the norm

kfk2 = (f, f ) (3.4)

Further, notice that

d dtkfk

2_{= (f, f}

t) + (ft, f ) (3.5)

We multiply (2.5) with ut and integrate over the domain [0, 1] to get

(ut, utt) = (ut, (B(u)ux)x) (ut, uxxxx) (3.6)

By using integration by parts we may rewrite the two terms on the right hand side. For the first one we have

(ut, (B(u)ux)x) = utB(u)ux|10 (utx, B(u)ux) , (3.7)

and for the second,

(ut, uxxxx) = utuxxx|1₀ utx, uxx|1₀+ (utxx, uxx) . (3.8)

Utilizing these rewrites while adding (3.6) to itself and applying (3.5) yields @ @t ⇥ kutk2+kuxxk2 ⇤ =BT 2 (utx, B(u)ux) , (3.9)

in which the boundary terms are collected into BT, where BT = 2 utB(u)ux|1₀ 2 utuxxx|1₀+ 2 utxuxx|1₀ = 2 [ut(B(u)ux uxxx) + utxuxx]10. (3.10) An alternative version of (3.9) is @ @t ⇥ kutk2+kuxxk2⇤=BT Z 1 0 B(u)@ @tu 2 xdx. (3.11)

As shown in Appendix A B(u) will always be larger than zero so assuming a coefficient frozen to a function independent of t B(u) = b(x) the last term can be viewed as a weighted norm, written

Finally, terms can be rearranged to @

@tEb =BT, (3.13)

where Eb is called the energy associated with (2.5) and is defined as

Eb =kutk2+kuxk2b +kuxxk2. (3.14)

As shown in appendix A b(x) will always be larger than 0, so Eb will always

be a valid energy.

Following now are two theorems regarding the boundedness of the two problems described above

Theorem 1. Problem (3.1) conserves energy.

Proof. Taking the homogenouse boundary conditions from (3.1), differenti-ating with respect to t and inserting them into (3.13) gives

@

@tEb = 0. (3.15)

Theorem 2. The frozen coefficient version of problem (3.2) is strongly well posed

Proof. By introducing the vector

w =⇥ut ux utx uxx uxxx⇤ (3.16)

and freezing the coefficient B(u) we may write the energy estimate in (3.13) as @ @tEb= w >_Aw 1 0 (3.17) where A = 2 6 6 6 6 4 0 B(u) 0 0 1 B(u) 0 0 0 0 0 0 0 1 0 0 0 1 0 0 1 0 0 0 0 3 7 7 7 7 5 (3.18)

The eigenvalues and eigenvectors of A are

By introducing ˜

w = T>w ⇤ = T>AT (3.20)

where T is the transformation matrix from the basis of eigenvectors we may rewrite the energy estimate as

@ @tEb = ˜w >_{⇤ ˜w}1 0 = ⇥ 1w˜21+ 2w˜22+ 4w˜42+ 4w˜52 ⇤1 0 (3.21)

we now see that by setting eigenvectors corresponding to negative eigenvalues on the left boundary and eigenvectors corresponding to positive eigenvalues on the right boundary, we achieve the strong energy estimate,

@ @tEb = w˜ 2 1 1 w˜2 12 w˜24+0 w˜ 2 5 0+ g02+ h20+ g21+ h21, (3.22)

and hence problem 3.2 with B(u) frozen is strongly well posed.

Remark 1. The method used to obtain the energy estimate for the charac-teristic problem is in it self the way these boundary conditions were found. In this way these boundary conditions are analytically natural and may be viewed as a first order approximation of non-reflective boundary conditions. Remark 2. The problems are stated on the interval x 2 [0, 1] but the analysis holds for any finite interval of R.

4 The finite difference method

This section describes the Summation-by-parts framework that is used to obtain the numerical solutions to the problems in the previous section. It will also introduce the notions of a compatible set of SBP operators and borrowing, which will be used in showing stability. The summation-by-parts framework consists of carefully crafting central finite difference operators with associated boundary closures to mimic the integration-by-parts formula.

4.1 Defnitions

Matrices ei, Si,S2i and S3i are column vectors approximating values of the

function, the derivative, second derivative and third derivative of the solution at the boundary point i when multiplied with the solution vector v. The SBP-operators can be written as follows

where H is a diagonal matrix, Q and R are anti-symmetric matrices and M and N are symmetric and positive definite matrices.

Definition 1. A D1 and a D4 operator are said to be compatible if they

satisfy

e>_i D1= Si>, i = 1, m (4.5)

Remark 3. The requirement that e>

i D1 = Si> is not trivial, the D1 operator

must be specifically derived to have this property.

Lemma 3. Assuming compatible D1 and D4 operators, the combined

op-erator D1BD1, where B is a positive definite diagonal matrix, is an SBP

operator of the form D1BD1 = H 1

⇣

D₁>HBD1+ R2B

⌘

, R2B = b1e1S1>+ bmemSm> (4.6)

where b1 and bm are the first and last elements on the diagonal of B.

Proof. Using the definition of D1 and the properties Q = Q>, R1 = R>1

and the commutation of diagonal matrices we get

H 1D1BD1 = D1>HBD1+ H 12R1BD1 (4.7) further, 2R1BD1 = ( e1e>1 + eme>m)BD1 = b1e1S1>+ bmemSm>. (4.8) Then D1BD1= H 1 ⇣ D₁>BD1+ 2R1BD1+ R2B ⌘ (4.9) where D>

1BD1 is a positive definite matrix.

To simplify the proof of stability for the Dirichlet-Neumann problem we now the concept of borrowing from positive definite matrices. Borrowing is a way of splitting a positive definite matrix in such a way that we end up with another positive definite as well as some second matrix. This will be used combined with the energy method to prove stability of the semi-discretized form of the Dirichlet-Neumann problem.

Definition 2. (↵, S)is called a valid borrowing from an operator D = H 1(N + R) , N = N>

if the matrix

N ↵S

The amount available for borrowing, ↵ depends on the operator and what matrix,S, is being borrowed. Table 1 shows a set of numerically derived borrowings from the compatible D4 operator.

Matrix 4th order 6th order

S21S2>1 + S2mS2>m 0.4210h 0.06925h

S31S3>1 + S3mS3>m 0.7080h3 0.05128 h3

Table 1: Borrowing constants for the compatible D4operator. h refers to the grid size.

Borrowing multiple times from the same operator is summarized in the following lemma.

Lemma 4. Given two valid borrowings (↵, S), ( , T ) from an operator D, the borrowing (1, a↵S + b T ) is a valid borrowing from D if a, b 0 and a + b < 1

5 Semi-Discretization and stability

In this section semi-discretizations of the problems in section 3 will be pro-vided using the finite difference method described in the previous section. The semi-discretizations will be shown to have a bounded energy growth.

5.1 Dirichlet-Neumann problem

The problem given in (3.1) can be discretized using the SBP operators. It

then becomes ₈ > > > > > > < > > > > > > : vtt = D1B(v)D1v D4v, t > 0 e>₁v = g0, S>1v = h0, t > 0 e>_mv = g1, Sm>v = h1, t > 0 v = f1, t = 0 vt= f2, t = 0 (5.1)

where B(v) is a matrix with the values of the wave-speed function B(u) in-serted on the diagonal. Using the SAT technique to implement the boundary condition we get the numerical scheme

vtt= D1B(v)D1v D4v + H 1⌧0(e>1v g0) + H 1 0(S1>v h0)

+ H 1⌧1(e>mv g1) + H 1 1(Sm>v h1)

(5.2) Where ⌧0, ⌧1, 0 and 1 are penalty parameters on vector form that need to

version of this scheme. By multiplying (5.2) with v>

t H and adding the

transpose of the result to itself we arrive at @ @tEH = 2v > t e1e>1BD1v + 2v>t eme>mBD1v + 2v_t>e1S3>1v 2v>t emS3>mv 2v_t>S1S2>1v + 2v>t SmS2>mv + 2v_t>⌧0e>1v + 2v>t 0S>1v + 2v_t>⌧1e>mv + 2v>t 1S>mv (5.3)

where we have define the energy

EH =kv2k2H +kD1vk2HB+kvk2N (5.4)

Under the assumption that the SBP operators D1 and D4 are compatible

they satisfy e>

i D1 = Si> for i = 1, m and we may write

@ @tEH = 2b1v > t e1S1>v + 2bmv>t emSm>v + 2v>_t e1S3>1v 2v>t emS3>mv 2v>_t S1S2>1v + 2v>t SmS2>mv + 2v>_t ⌧0e>1v + 2v>t 0S1>v + 2v>_t ⌧1e>mv + 2v>t 1Sm>v (5.5)

where b1 and bm are the first and last diagonal element of B(v) This

ex-pression may be simplified by choosing the penalty parameters the be linear combinations of ei, Si,S2i and S3i

⌧i= ↵(1)i ei+ ↵(2)i Si+ ↵(3)i S2i+ ↵(4)i S3i (5.6) i= _i(1)ei+ _i(2)Si+ _i(3)S2i+ _i(4)S3i (5.7)

and by introducing the vectors

w₁>=⇥e>₁v S₁>v S2>₁v S3>₁v⇤ (5.8) w_m> =⇥e>_mv S_m>v S2>_mv S3>_mv⇤ (5.9) amd we may write (5.5) as

Going even further we can write @ @t h EH w>1A1w1 w>mAmwm i = 0, (5.12)

This energy mimics the energy from the continuous case with the exception of the terms due to A1 and Am. The soundness of the scheme is dependent

on the properties of this energy and these properties are affected by the matrices A1 and Am, which in turn are affected by the choice of penalty

parameters. By choosing these properly the scheme will be stable. Theorem 5 provides conditions on these parameters for the scheme to be stable.t Theorem 5. Assume the D1 and D4 operators are compatible. Let

⌧1= ↵1e1+ S31, 1= b1e1+ 1S1 S21,

⌧m= ↵mem S3m,

m= bmem+ mSm+ S2m,

(5.13)

where ↵1, ↵m, 1, m, ⇣1 and ⇣m are parameters. Under the condition that

these paramters satisfy

↵1, ↵m < 1 , 1, < 1 + ⇣ 2 1 ↵1 1 , m, < 1 + ⇣ 2 m ↵m 1 , ⇣1, ⇣m 2 R, (5.14)

where and are chosen so that (1, S21S2>1 + S2mS21> S31S3>1 + S3mS3>1

is a valid borrowing for D4, the numerical scheme given in (5.2) is stable.

For this proof homogeneous boundary data may be assumed with out any loss of generality.

Proof. To show stability we need to show (5.12) is a valid energy estimate, namely that the energy defined is positive definite. By introducing

↵(1)₁ = ↵1, 1(2) = 1

↵(2)₁ = b1+ ₁(1)= ⇣1

↵(3)₁ = ₁(4)= 0 ↵(4)₁ = ₁(3)= 1

and ↵(1)_m = ↵m, m(2)= m ↵(2)_m = bm+ m(1) = ⇣m ↵(3)_m = _m(4) = 0 (3) m = ↵(4)m = 1 (5.16)

which gives us the matrices A1 = 2 6 6 4 ↵1 ⇣1 0 1 ⇣1 1 1 0 0 1 0 0 1 0 0 0 3 7 7 5 Am= 2 6 6 4 ↵m ⇣m 0 1 ⇣m m 1 0 0 1 0 0 1 0 0 0 3 7 7 5 (5.17)

we enforce symmetricalness of the two matrices and thereby give them real eigenvalues. According to our and assumption and Definition 2 we now introduce the borrowing

(1, ⇣S21S2>1 + S2mS2>1

⌘ ⇣

S31S3>1 + S3mS3>1

⌘

(5.18) from D4. This allows us to write

N = ˜N ⇣S21S2>1 + S2mS2>1

⌘ ⇣

S31S3>1 + S3mS3>1

⌘

(5.19) , where ˜N is a positive definite matrix. Using ˜N instead of N in the energy yields

˜

EH =kvtk2H +kD1vk2HB+kvk2_N˜, (5.20)

which is an alternative energy. Using the above to rewrite (5.12) we get @ @t h ˜ EH w1>A˜1w1 wm>A˜mwm i = 0 (5.21) where ˜ A1= 2 6 6 4 ↵1 ⇣1 0 1 ⇣1 1 1 0 0 1 0 1 0 0 3 7 7 5 A˜m= 2 6 6 4 ↵m ⇣m 0 1 ⇣m m 1 0 0 1 0 1 0 0 3 7 7 5 (5.22)

We may now derive conditions on ↵1, 1,⇣1, ↵m, m,⇣m, and for ˜A1 and

˜

Am to be negative definite using Sylvester’s criterion. We get

Putting everything together, by choosing

⌧1= ↵1e1+ ⇣1S1+ S31 (5.24)

1= (b1+ ⇣1)e1+ 1S1 S21 (5.25)

⌧m= ↵mem+ ⇣mSm S3m (5.26)

m= ( bm+ ⇣m)em+ mSm+ S2m (5.27)

we get the energy estimate @ @t h ˜ EH w1>A˜1w1 wm>A˜mwm i = 0 (5.28)

where the negative definiteness of ˜A1 and ˜Am ensures that we have a valid

energy and so the scheme is stable.

Remark 4. The values of the parameters in the scheme can be chosen to minimize the spectral radius of the final scheme operator. In particular, after some numerical tests it seems choosing ⇣1 = ⇣m = 0both simplifies the

scheme and provides a smaller spectral radius than other choices.

5.2 Characteristic problem

The semi-discretization of (3.2) using SBP operators is 8 > > > > > > > > > > > > < > > > > > > > > > > > > : vtt= D1B(v)D1v D4v, t > 0 S₁>vt S2>1v = g0 1e>1vt b1S1>v + S3>1v = h0 S_m>vt+ S2>mv = g1 me>mvt+ bmSm>v S3>mv = h1 v = f1 t = 0 vt= f2 t = 0 (5.29)

where i and bi are and B(u) evaluated for vi respectively. Using the SAT

technique we get the following scheme vtt = D1B(v)D1v D4v+H 1⌧0 ⇣ S₁>vt S2>1v g0 ⌘ +H 1 0 ⇣ 1e>1vt b1S1>v + S3>1v h0 ⌘ +H 1⌧1 ⇣ S_m>vt+ S2>mv g1 ⌘ +H 1 1 ⇣ me>mvt+ bmSm>v S3>mv h1 ⌘ (5.30)

SBP operators we get @ @tEH = 2v > t R2Bv 2vt>R4v+2vt>⌧0 ⇣ S₁>vt S2>1v g0 ⌘ +2v_t> 0 ⇣ 1e>1vt b1S1>v + S3>1v h0 ⌘ +2v_t>⌧1 ⇣ S_m>vt+ S2>mv g1 ⌘ +2v_t> 1 ⇣ me>mvt+ bmSm>v S3>mv h1 ⌘ , (5.31) where EH =kvtk2H +kD1bk2HB+kvk2N (5.32)

is the energy exactly mimicking the energy from the continuous case.

Theorem 6. Let

⌧0= S1 0= e1

⌧1= Sm 1= em

. (5.33)

Then the semi-discretization in (5.30) is stable.

Proof. Plugging the chosen penalty parameters into (5.31) we see that most of the terms cancel and we end up with

@ @tEH = 2v > t S1 ⇣ S₁>vt g0 ⌘ 2v_t>e1 ⇣ 1e>1vt h0 ⌘ 2v_t>Sm ⇣ S_m>vt g1 ⌘ 2v_t>em ⇣ me>mvt h1 ⌘ , , @ @tEH = 2 ⇣ S₁>vt g0 ⌘2 + 2g₀2 2 1 ✓ e>₁vt h0 1 ◆2 + 2h 2 0 1 2⇣S_m>vt g1 ⌘2 + 2g2₁ 2 m ✓ e>_mvt h1 m ◆2 + 2h 2 1 m (5.34)

which shows strong stability.

6 Time integration

The focus of this work is the discretization of the equation in space and enforcing conditions on the boundaries of the domain. Therefore, to avoid complications brought by using implicit methods, two simpler methods were derived. The first, writing the system on first order form and using the Runge-Kutta 4th order method. The second applying 2nd order central finite differences. Due to extremely stiff nature of the equation coming from the fourth derivative both methods require very small time steps.

As we now focus our attention on the time aspects of the problem, for this analysis, we write the scheme for the Dirichlet-Neumann problem given in (3.1) as

vtt= Dpv + Sp (6.1)

and the scheme for the characteristic boundary condition problem given in (3.2) as vtt= Dcv + Ecvt+ Sc (6.2) where Dp = D1B(v)D1 D4+ H 1 h (↵1e1+ S31) e>1 + (b1e1+ 1S1 S21) S1> i + H 1h(↵mem S3m) e>m+ ( bmem+ mSm+ S2m) Sm> i (6.3) Sp = H 1[(↵1e1+ S31) g0+ (b1e1+ 1S1 S21) h0] + H 1[(↵mem S3m) g1+ ( bmem+ mSm+ S2m) h1] (6.4) and Dc = D1B(v)D1 D4+ H 1 ⇣ +S1S2>1 + b1e1S1> e1S3>1 ⌘ + H 1⇣ SmS2>m bmemSm> emS3>m ⌘ (6.5) Ec = H 1 ⇣ S1S1>+ 1e1e>1 + SmSm>+ meme>m ⌘ (6.6) Sc= H 1(S1g0+ e1h0+ Smg1+ emh1) (6.7)

for the respective problems. As shown in appendix B the eigenvalues scaling with respect to h, the grid space, is related to the way Dp, Dc and Ec scale

with respect to h. As we see from the expressions of these matrices we have

Dp ⇠ 1 h4 Dc ⇠ 1 h4 Ec⇠ 1 h3 (6.11)

for small h. For the first order in time systems we then have Mp⇠

↵p

h2 Mc⇠

↵c

h3, (6.12)

where ↵p and ↵c are constants. Using Matlab, the matrices of the different

problems we studied for different h to determine these constants. The results are shown in table 2.

Order

4 6

Dirichelt-Neumann 11.2060 31.2879

Characteristic 33.8129 26.1479

Table 2: Scaling constants of the matrices for schemes of orders 4 and 6 corresponding to the different boundary conditions.

6.1 CFL for Runge-Kutta

The eigenvalues of the operators are in this case purely real or imaginary. This simplifies the stability constraints for the Runge-Kutta method. The stability region in the complex plane for the 4th order Runge-Kutta method is shown in figure 1. By numerically determining where its boundary intersects the real and imaginary axis we get

2.7853 < k _Re> 0 _{|k Im| < 2.8284,} (6.13)

where k is the time step and Reand Im are real and imaginary eigenvalues.

Even though our schemes are non-linear, in which case you can’t talk about eigenvalues, the largest factor in determining the eigenvalues are the 4th derivative SBP operator D4. Thus it is still useful to examine the schemes

Order

4 6

Dirichlet-Neumann k < 2.8284

11.2060h2 k < 31.28792.8284h2

Characteristic k < _33.81292.7853h3 k < _26.14792.7853h3

Table 3: CFL conditions for the two first order in the schemes of orders 4 and 6.

Notice that the stability condition for the characteristic problem scales with h3_{, this unexpected result menans that the stiffness problem is even}

worse for the characteristic problem. This is due the boundary treatment of this scheme introduces off-diagonal elements to the first order in time operator. The Runge-Kutta method is therefore rendered completely useless for the characteristic problem.

4 3 2 1 0 1 2 3 2 1 0 1 2 3 Re Im

4th order Runge-Kutta stability region

Figure 1: The stability region for the 4th order Runge-Kutta method.

6.2 Centered finite difference scheme

Due to the bad scaling of the eigenvalues for the characteristic problem the central finite differences method was also applied. This method avoids the scaling behavior of the vt term by handling it in a implicit manner. Using

get utt⇡ 1 k2 (un+1 2un+ un 1) ut⇡ 1 2k(un+1 un 1) Given a scheme vtt= Dv + Evt+ S (6.14)

using these centered finite differences this can be written as

vn+1 = A 1(Bvn+ Cvn 1+ S) , (6.15) where A = 1 k2I + 1 2kE B = 2 k2I + D C = 1 k2I 1 2kE,

where I is the identity matrix. Since E only affect points close to the bound-aries A is a block diagonal matrix which inverse is also a block diagonal matrix making this feasible.

A CFL condtion may be arrived at by studying the method applied to the scalar equation

utt = aut+ bu. (6.16)

Applying the central finite difference approximations and introducing vn = n _{we get}

2_{+ 2}bh2+ 2

ah2 ₂

ah2+ 2

ah2 ₂ = 0, (6.17)

which has the solutions

= bh 2_{+ 2} ah2 ₂ ± s✓ bh2_{+ 2} ah2 ₂ ◆2 + ah2+ 2 ah2 ₂. (6.18)

From this we see that the behavior of is ⇠ bk_ak22. Stability requires that

< 1. For our scheme we then get the CFL

⇢(D)k2< K (6.19)

where ⇢(D) ⇠ 1

h4 is the spectral radius of D. Numerical experiments showed

7 Convergence

To test the convergence of the schemes the analytical solitons were used as reference solutions. Initial and boundary data were set to be that of the analytical solution. And the domain was chosen as x 2 [0, 10]. The system was then numerically evolved to time t = 1.0s where the global error was calculated according to e(h) = s hX 8i (vi u(xi, t))2, (7.1)

where v is the numerical solution, u is the analytical solution and i = 1, ..., m are the indecies of the grid points. The error was calculated for a number of grids and the convergence number was calculated according to

q(h) = loge(hprev) e(h) / log

h

hprev. (7.2)

The analytical solution used was the soliton described in Section 2.2 with = 0.8 for fourth and sixth order of accuracy in space. The central finite difference scheme was used for the time integration. The results are shown in Table 4 and Table 5.

The global error from our 2nd, 4th and 6th order schemes are expected to decrease as order 2, 4, and 5. These numbers depend on the highest order derivative present in the problem, the boundary conditions set as well as how the SBP-operator is built. [12]

2nd order 4th order 6th order

N log10(l2) q log10(l2) q log10(l2) q

21 -2.39 -3.54 -3.99 41 -2.86 -1.54 -4.89 -4.50 -5.78 -5.93 81 -3.24 -1.26 -5.99 -3.63 -7.26 -4.91 161 -3.56 -1.08 -7.05 -3.53 -8.80 -5.12 321 -3.87 -1.04 -8.41 -4.52 -8.84 -0.14 641 -4.18 -1.03 -8.79 -1.26 -8.73 0.38

10 2 ₁₀ 1 ₁₀0 ₁₀1 ₁₀2 ₁₀3 10 9 10 8 10 7 10 6 10 5 10 4 10 3 10 2 runtime e

Efficiency: Analytical Soliton with Dirichelt-Neumann BC

Order 2 Order 4 Order 6

Figure 2: Efficiency of computation with respect to error at t = 1 and compu-tation time. Analytical soliton with Dirichlet-Neumann boundary conditions.

The error and convergence calculations for the Dirichlet-Neumann prob-lem shown in Table 4. shows better than expected convergence for the first couple of refinements. However, for order 4 and 6, the global error quickly hits an lower limit that prevents further grid refinements from providing a more exact solution. We have reached grid convergence. Beyond a certain refinement the error becomes constant at around 1e 8.83. This is to be considered a fairly high error.

2nd order 4th order 6th order

N log10(l2) q log10(l2) q log10(l2) q

21 -3.94 -3.94 -4.24 41 -4.59 -2.15 -5.15 -3.99 -5.86 -5.38 81 -5.03 -1.48 -6.22 -3.55 -7.39 -5.09 161 -5.32 -0.96 -7.10 -2.95 -8.59 -3.98 321 -5.60 -0.92 -7.97 -2.88 -8.72 -0.43 641 -5.88 -0.95 -8.60 -2.11 -8.72 0.00

10 2 ₁₀ 1 ₁₀0 ₁₀1 ₁₀2 ₁₀3 10 9 10 8 10 7 10 6 10 5 10 4 10 3 runtime e

Efficiency: Analytical Soliton with characteristic BC

Order 2 Order 4 Order 6

Figure 3: Efficiency of computation with respect to error at t = 1 and com-putation time. Analytical soliton with Characteristic boundary conditions.

Table 5 shows the error and convergence calculations for the Character-istic problem. A similar pattern shows up here as for the Dirichlet-Neumann problem, with a better than expected convergence up to about 100 grid points and then a stagnation of the error.

The discrepancy between the expected convergence and the convergence seen in this study and in particular in the case of the Characteristic problem might be explained by the fact that we have not explicitly verified that the assumptions, for example regarding the boundary conditions, in [12] are fulfilled. We do see close to expected convergence for the Dirichlet-Neumann problem which suggests that the operators are working as expected and that the problem may instead lie in the specific boundary conditions chosen.

In figures 2 and 3 the efficiency of the schemes at different grids is shown. Here we clearly see why in some cases a higher order method is preferred. If we require an error no less than 10 7 _{for the characteristic problem, we}

see that the 6th order accurate scheme is about 10 times faster than the 4th order scheme. The 2nd order scheme never seem to get close to being exact enough.

ex-planation for the high bound might be the use of a longer interval that usual, [0, 10], instead of [0, 1]. Another explanation might be that the stiffness of the problem amplifies rounding errors made in the floating point computa-tions. Possible means of lowering this bound on the error is tuning of the penalty parameters to minimize the spectral radius of the schemes operator and use higher precision arithmetic.

8 Experiments

8.1 Comparison of a soliton and a non-soliton pulse

Shown here are two simulations of two different pulses traveling to the right in the domain. In Figure 4 we see a soliton of the form 0.65. As expected

it travels without being distorted despite the nonlinearity of the equation. The numerical scheme successfully captures the balancing of focusing and dispersive effects on the pulse. In Figure 5 we see a pulse with a similar shape to a soliton but with it’s amplitude changed so that it is no longer a soliton solution of the PDE, 2 0.76. The focusing and dispersive effects on

the pulse no longer cancel and smaller pulses are spreading off of the main one. Simulations were done using the 4th order accurate SBP scheme and the 4th order Runge-Kutta method. The number of grid points were set to 400 and the time step was chosen to achieve stability according to Section 6.

0 20 40 60 80 0 20 40 60 80 100 x t Regular Soliton 0 5· 10 2 0.1 0.15 0.2 u

0 20 40 60 80 0 20 40 60 80 100 x t Amplified soliton 0 5· 10 2 0.1 0.15 0.2 u

Figure 5: Non-Soliton pulse obtained by increasing the amplitude of a soliton with = 0.76 by a factor of 2. The pulse breaks up into several pulses as it travels.

8.2 Colliding solitons

0 50 100 150 200 0 50 100 150 200 250 300 x t

Colliding asymmetric solitons

0 5_{· 10} 2 0.1 0.15 0.2 0.25 u

Figure 6: Two asymmetric solitons colliding.

0 50 100 150 200 0 50 100 150 200 250 x t

Colliding symmetric solitons

0 5_{· 10} 2

0.1 0.15 u

8.3 Soliton Interacting with boundary

Here are shown how a single soliton interacts with the two different bound-ary conditions. In Figure 8 a soliton is shown bouncing on a homogeneous Dirichlet-Neumann condition. We notice that the soliton breaks down upon reaching the boundary and after reflection dispersion splits the pulse into several smaller pulses. However, the whole pulse is reflected back into the domain. 0 20 40 60 80 100 0 20 40 60 80 100 120 x t

Soliton Boundary Interaction

4 2 0 2 4 6 8 ·10u 2

Figure 8: A soliton bouncing off a homogeneous Dirichlet-Neumann boundary condition.

0 20 40 60 80 100 0 20 40 60 80 100 120 x t

Soliton Boundary Interaction

0 2 4 6 8 ·10u 2

Figure 9: A soliton interacting with a homogeneous Characteristic boundary condition.

9 Conclusions and future work

The goal of designing numerical schemes for the Soliton model of nerve sig-nal propagation was achieved. The schemes successfully enforce two different kind of boundary conditions and were shown to be stable. The convergence of the schemes was verified to be what can be expected. This work shows that the SBP-SAT finite difference method is a viable and efficient way to ob-tain numerical solutions to the Soliton model. It also shows that the method flexibly implements the given boundary conditions in a stable manner. How-ever, the final error after grid refinement seems to be higher than expected. The reason for this was not found and could be the subject of future studies. Interesting topics for future studies of this problem is for example how nu-merically treat boundary interfaces used when modeling axon junctions or soma and synapse junctions.

Acknowledgments

A Wave speed function analysis

The wave speed function in the Thermodynamic model is approximated by a power series truncated after the quadratic term and is written as

B(u) = 1 + B1u + B2u2. (A.1)

For any second degree polynomial on the form

f (x) = x2+ 2ax + b (A.2)

we know that

• b > a2 _{) f(x) > 0}

• b = a2 _{) f(x) 0}

• b < a2 _{) f(x) has a minimum.}

We may rewrite (2.6) into B(u) = B2 ✓ u2+B1 B2 u + 1 B2 ◆ , (A.3)

and we now see that 4 B2 > B 2 1 B2 2 ^ B2 > 0) B(u) > 0. (A.4) Alternatively 4B2 > B12^ B2> 0) B(u) > 0.

Using the experimentally derived values for the constants, B1 = 16.6 and

B2 = 79.5, we see that this holds.

B Determinants of block matrices

Two identities regarding block determinants are useful when examining eigen-value problems for block matrices. They are listed below.

If A,B, C and D are matrices of the same order and C and D commute we have that det ✓ A B C D ◆ = det (AD BC) (B.1)

If A,B, C and D are matrices of the same order and A is invertible we have that det ✓ A B C D ◆

We may use these relations to infer some things about the eigenvalues of operators for second order in time ODE written on first order form. For an ODE of the form vtt = Dv + S we may write

 v v_{t t}= ✓ 0 I D 0 ◆  v vt +  0 S . (B.3)

The eigenvalue problem of this operator is det ✓ I I D I ◆ = 0. (B.4)

Since I commutes with every matrix we may use the first relation to obtain

det 2I D = 0 (B.5)

which means that the eigenvalues of the first order operator are ±pµ where µare the eigenvalues of D.

For an ODE of the form vtt = Dv + Evt+ S we may write

 v v_{t t}= ✓ 0 I D E ◆  v vt +  0 S . (B.6)

The eigenvalue problem of this operator is det ✓ I I D I E ◆ = 0. (B.7)

Since I is its own inverse we may use the second relation to obtain det ( I) det ✓ I E D1I ◆ = 0, (B.8) det 2I E D = 0 (B.9)

From this we see that the eigenvalues of the first order operator scales as the square root of the eigenvalues of D and as the eigenvalues of E.

C Diagonal-norm SBP operators

Below we list the first- and fourth-derivative SBP operators for the second-, fourth- and sixth-order cases. The first-derivative SBP operator is given by,

D1= H 1 ✓ Q 1 2e1e T 1 + 1 2eme T m ◆ . And the fourth-derivative SBP operator is given by,

D4 = H 1 N e1d3;1+ emd3;m+ dT1;1d2;1 dT1;md2;m .

The boundary closures (the coefficients) are presented for the H, Q and N matrices. Notice that the coefficients in the norm should be multiplied by the grid-step h and the coefficients in N should be multiplied by 1

C.1 Second order case.

The boundary derivative operators are given by d1;1v = v1_h+v2, d1;mv = vm _hvm 1

d2;1v = v1 2v_h22+1v3, d2;mv = vm 2vm 1_h2 +vm 2

d3;1v = v1+3v_h233v3+v4, d3;mv = +vm 3vm 1+3v_h3 m 2 vm 3

The internal schemes of D1,4 are given by

(D1v)j = vj _{2 h}1+vj+1, (D4v)j = vj 2 4vj 1+6v_h4j 4vj+1+vj+2

The norm is the traditional second order norm, H = diag(1

2, 1, . . . , 112).

The upper part of Q, is given by Q1,2 = 1₂.

The upper part of N, N1,1= 4₅ N1,2= 7₅ N1,3 = 2₅ N1,4 = 1₅ N2,2 = 16₅ N2,3 = 11₅ N2,4= 2₅ N3,3= 21₅ N3,4 = 17₅ N4,4 = 29₅

C.2 Fourth order case.

The boundary derivative operators are given by

d1;1v = 11v1+18v_{6 h}2 9v3+2v4, d1;mv = +11vm 18vm 1_{6 h}+9vm 2 2vm 3

d2;1v = 2v1 5v2_h+4v2 3 v4, d2;mv = 2vm 5vm 1+4v_h2m 2 1vm 3

d3;1v = v1+3v2_h33v3+v4, d3;mv = +vm 3vm 1+3v_h3 m 2 vm 3

The internal schemes of D1,4 are given by

(D1v)j = vj 2 8vj 1_{12 h}+8vj+1 vj+2

(D4v)j = vj 3+12vj 2 39vj 1+56v_{6 h}4j 39vj+1+12vj+2 vj+3

The upper part of the norm H, H1,1 = ₁₁3 H2,2 = 2125516311_1311004640 H3,3= _1966506960278735189 H4,4= 285925927_163875580 H5,5= 1284335339_1966506960 H6,6= 4194024163_3933013920

The upper part of N, N1,1= 227176919517319_{94899692875680} N1,2= 15262605263734_{2965615402365} N1,3= 20205404771243_{6778549491120} N1,4= _{23724923218920}3998303664097 N1,5= _{94899692875680}1088305091927 N1,6= _{23724923218920}1686077077693 N2,2= 280494781164181_{23724923218920} N2,3= 46417445546261_{5931230804730} N2,4= 1705307929429_{1694637372780} N2,5= _{5931230804730}553547394061 N2,6= _{23724923218920}5615721694973 N3,3= 4135802350237_551742400440 N3,4= 4140981465247_{1078405600860} N3,5= 75538453067437_{47449846437840} N3,6= _{11862461609460}4778134936391 N4,4= 20760974175677_{2965615402365} N4,5= 138330689701889_{23724923218920} N4,6= 23711317526909_{11862461609460} N5,5= 120223780251937_{13557098982240} N5,6= 151383731537477_{23724923218920} N6,6= 220304030094121_{23724923218920}

C.3 Sixth order case.

The boundary derivative operators are given by d1;1v = 12700800 7493827v1+18502332189925924v2+3979148944962962v3 3824325314987654v4+7609944744962962v5 3493607589925924v6 h , d1;mv = + 12700800 7493827vm 185023321 89925924vm 1 39791489 44962962vm 2+ 38243253 14987654vm 3 76099447 44962962vm 4+ 34936075 89925924vm 5 h , d2;1v = 35v1 104v2+114v_{12 h}23 56v4+11v5, d2;mv = 35vm 104vm 1+114v_{12 h}m 22 56vm 3+11vm 4 d3;1v = 5v1+18v2 _{2 h}24v33+14v4 3v5, d3;mv = +5vm 18vm 1+24v_{2 h}m 23 14vm 3+3vm 4

The internal schemes of D1,4 are given by

(D1v)j = vj 3+9vj 2 45vj_{60 h}1+45vj+1 9vj+2+vj+3

(D4v)j = 7vj 4 96vj 3+676vj 2 1952vj 1+2730v_{240 h}4j 1952vj+1+676vj+2 96vj+3+7vj+4

The upper part of Q, Q1,2= 26431903_43545600 Q1,3= _15240960039791489 Q1,4= 12747751_16934400 Q1,5= _15240960076099447 Q1,6= _121927681397443 Q1,7= 0 Q1,8= 0 Q2,3= 13847476213_19559232000 Q2,4= 35844843977_11735539200 Q2,5= 63413503537_23471078400 Q2,6 = 4764412871_3911846400 Q2,7= 1668252557_5867769600 Q2,8= _29338848000842644697 Q3,4= 73834802771_23471078400 Q3,5= 1802732209_325987200 Q3,6= 65514173_16299360 Q3,7= 79341409141_58677696000 Q3,8= 1282384321_7823692800 Q4,5= 5274106087_1173553920 Q4,6= 33743985841_5867769600 Q4,7= 6482602549_2607897600 Q4,8= 1506017269_3911846400 Q5,6= 7165829063_2607897600 Q5,7= 23903110999_11735539200 Q5,8= _117355392005346675911 Q6,7= _117355392001060918223 Q6,8= _3911846400628353989 Q7,8= 25889988599_39118464000

The upper part of N, N1,1= 600485868980522851_{274825314114120000} N1,2= 1421010223681841_{348984525859200} N1,3= _{1586293299360000}38908412970187 N1,4= 10224077451922837_{2243471951952000} N1,5= 7577302712815639_{1744922629296000} N1,6= 138091642084013_{59351109840000} N1,7= 3775041725375197_{4486943903904000} N1,8= _{61072292025360000}9907210230881393 N2,2= 3985852497808703_{407903991264000} N2,3= 90048788923861_{15579666333000} N2,4= 4312795866499_997098645312 N2,5= 4414634708891947_{448694390390400} N2,6= 886174803100459_{99709864531200} N2,7= 4333₁₀₀₀ N2,8= 13800578064893047_{15704303663664000} N3,3= 2071682582321887_{113306664240000} N3,4= 769471337294003_{41545776888000} N3,5= 112191585452033_{166183107552000} N3,6= 7204491902193671_{623186653320000} N3,7= 24847093554379_{3115933266600} N3,8= 943854037768721_{545288321655000} N4,4= 3086874339649421_{81580798252800} N4,5= 396009005312111_{16618310755200} N4,6= 348854811893087_{249274661328000} N4,7= 895954627955053_{224347195195200} N4,8= 184881685054543_{166183107552000} N5,5= 3774861828677557_{112173597597600} N5,6= 5693689108983593_{249274661328000} N5,7= 803944126167107_{99709864531200} N5,8= 19547569411550791_{15704303663664000} N6,6= 73965842628398389_{2492746613280000} N6,7= 2184472662036043_{124637330664000} N6,8= 46667₁₀₀₀₀ N7,7= 37593640125444199_{2243471951952000} N7,8= 37₄ N8,8= 12766926490502478779_{1099301256456480000}

References

[1] Søren SL Andersen, Andrew D Jackson, and Thomas Heimburg. To-wards a thermodynamic theory of nerve pulse propagation. Progress in neurobiology, 88(2):104–113, 2009.

potential in nerves. Advances in Planar Lipid Bilayers and Liposomas, 16:271–279, 2012.

[3] M. H. Carpenter, D. Gottlieb, and S. Abarbanel. Time-stable bound-ary conditions for finite-difference schemes solving hyperbolic systems: Methodology and application to high-order compact schemes. J. Com-put. Phys., 111(2), 1994.

[4] Alan L Hodgkin and Andrew F Huxley. A quantitative description of membrane current and its application to conduction and excitation in nerve. The Journal of physiology, 117(4):500, 1952.

[5] H.-O. Kreiss and G. Scherer. Finite element and finite difference meth-ods for hyperbolic partial differential equations. Mathematical Aspects of Finite Elements in Partial Differential Equations., Academic Press, Inc., 1974.

[6] Heinz-Otto Kreiss and Joseph Oliger. Comparison of accurate methods for the integration of hyperbolic equations. Tellus XXIV, 3, 1972. [7] K. Mattsson, M. Svärd, M.H. Carpenter, and J. Nordström. High-order

accurate computations for unsteady aerodynamics. Computers & Fluids, 36:636–649, 2006.

[8] K. Mattsson, M. Svärd, and M. Shoeybi. Stable and accurate

schemes for the compressible navier-stokes equations. J. Comput. Phys., 227(4):2293–2316, 2008.

[9] Ken Mattsson. Diagonal-norm summation by parts operators for finite difference approximations of third and fourth derivatives. Journal of Computational Physics, (0):–, 2014.

[10] J. Nordström and M. H. Carpenter. Boundary and interface conditions for high-order finite-difference methods applied to the Euler and Navier-Stokes equations. J. Comput. Phys., 148:341–365, 1999.

[11] G. Strang. Accurate partial difference methods II. non-linear problems. Num. Math., 6:37–46, 1964.