Accurate numerical simulation of electrodiffusion and water movement in brain tissue
Ada J. Ellingsrud, Nicolas Boullé, Patrick E. Farrell, Marie E. Rognes
FFebruary 5, 2021 13:25 WSPC/INSTRUCTION FILE ws-m3as
Mathematical Models and Methods in Applied Sciences © World Scientific Publishing Company
Accurate numerical simulation of electrodiffusion and water movementin brain tissue
Ada Johanne Ellingsrud
Simula Research Laboratory, Oslo, [email protected]
Nicolas Boull´e
Mathematical Institute, University of Oxford, Oxford, [email protected]
Patrick E. Farrell
Mathematical Institute, University of Oxford, Oxford, [email protected]
Marie E. Rognes
Simula Research Laboratory, Oslo, [email protected]
Mathematical modelling of ionic electrodiffusion and water movement is emerging as apowerful avenue of investigation to provide new physiological insight into brain home-ostasis. However, in order to provide solid answers and resolve controversies, the ac-curacy of the predictions is essential. Ionic electrodiffusion models typically comprisenon-trivial systems of non-linear and highly coupled partial and ordinary differentialequations that govern phenomena on disparate time scales. Here, we study numericalchallenges related to approximating these systems. We consider a homogenized modelfor electrodiffusion and osmosis in brain tissue and present and evaluate different as-sociated finite element-based splitting schemes in terms of their numerical properties,including accuracy, convergence, and computational efficiency for both idealized scenar-ios and for the physiologically relevant setting of cortical spreading depression (CSD).We find that the schemes display optimal convergence rates in space for problems withsmooth manufactured solutions. However, the physiological CSD setting is challenging:we find that the accurate computation of CSD wave characteristics (wave speed andwave width) requires a very fine spatial and fine temporal resolution.
Keywords : electrodiffusion, osmosis, brain electrophysiology and mechanics, finite ele-ment method, splitting scheme, numerical convergenceAMS Subject Classification: 22E46, 53C35, 57S201 a r X i v : . [ m a t h . NA ] F e b ebruary 5, 2021 13:25 WSPC/INSTRUCTION FILE ws-m3as A. J. Ellingsrud, N. Boull´e, P. E. Farrell, and M. E. Rognes
1. Introduction
The movement of ions and molecules in and between cellular compartments isfundamental for brain function, and importantly for neuronal excitability andactivity , , . Vital processes such as action potential firing, transmitter release,and synaptic transmission are all driven by ionic gradients across the neuronalmembrane. Regulation of the extracellular volume by cellular swelling is closelyrelated to ionic dynamics, including potassium buffering , . Concurrently, sev-eral pathologies are associated with disruption to ionic homeostasis in the brain,e.g. Huntington’s disease , multiple sclerosis , migraine , epilepsy , Alzheimer’sdisease , and cortical spreading depression . Recent research efforts indicate thation concentrations in the extracellular space are not static, but vary across statessuch as locomotion and the sleep cycle .In spite of these aspects, mathematical and numerical models for describ-ing dynamics in brain tissue traditionally assume that the ion concentrationsare constant in time and space. Although such models have provided valu-able insight into the mechanisms underlying excitable cells, they fail to repre-sent essential dynamics related to altered ion concentrations in brain tissue. Re-cently, several mathematical models also including electrodiffusive effects have beenpresented , , , , , , , , , , . However, to date little attention has been paidto the numerical solution of these models.In this paper, we consider a mathematical framework proposed by Mori , con-sisting of a system of partial differential equations (PDEs) governing ionic elec-trodiffusion and water flow in biological tissue, coupled to a system of ordinarydifferential equations (ODEs) describing the temporal evolution of ionic membranemechanisms. The system predicts the dynamics of volume fractions, ion concen-trations, electrical potentials and mechanical pressure in an arbitrary number ofcellular compartments and in the extracellular space (ECS). The cellular compart-ments can communicate with the ECS via transmembrane ion and water fluxes. Thismathematical model extends on the celebrated bidomain model , , , and both rep-resent the tissue in a homogenized manner. Homogenized models are coarse-grained,and hence well suited for simulating phenomena on the tissue scale (mm). Impor-tantly, the two models differ in that the classical bidomain model only predictselectrical potentials, whereas the model for ionic electrodiffusion and water flowtakes into account how the movement of ions affect the excitable tissue, both interms of electrochemical and mechanical effects.Previously, the electrodiffusive model (in its zero flow limit) has been usedto study dynamics in brain tissue, and in particular cortical spreading depression(CSD) , . CSD is a slowly propagating wave of depolarization of brain cells, char-acterized by elevated levels of extracellular potassium, calcium, and glutamate,cellular swelling and pronounced ECS shrinkage , . Importantly, CSD is a funda-mental pattern of brain signalling that challenges ionic homeostasis mechanisms inthe brain. As such, a better understanding of the sequence of events in CSD has theebruary 5, 2021 13:25 WSPC/INSTRUCTION FILE ws-m3as Numerics for brain electrochemomechanics potential to provide new insight into underlying processes both in cerebral physi-ology and pathology . The aforementioned computational studies have focused onproviding new insight into the role of glial cells in CSD , and the role of glutamatedynamics in CSD . However, in order to provide true physiological insight, theaccuracy of the numerical and computational predictions is key; the difference be-tween conflicting experimental observations may very well be within the numericalerror of underresolved models.Previously considered numerical schemes for the electrodiffusive model are basedon finite difference , or finite volume discretizations in space and a backwardEuler scheme (with explicit treatment of the active membrane flux) in time. Thespatio-temporal discretization sizes of these schemes are reported to be on the orderof ∆ x ≈ . − . t ≈
10 ms. By applying a (comparable) finite elementscheme in space and a similar discretization scheme in time, we find that the CSDwave properties change substantially during spatial and temporal refinement. Inparticular, we observe that the wave speed and the width of the wave increase withdecreasing time resolution, and decrease with decreasing spatial resolution (Fig-ure 1). As the Mori framework comprises a system of non-linear and highly coupledpartial and ordinary differential equations, governing phenomena on disparate timescales (both fast electrotonic effects and much slower effects mediated by diffusion),theoretical analysis of the full equations is challenging. In this paper, we take a nat-ural first step towards accurate and efficient approximation schemes for the systemby comparing different schemes numerically, with an emphasis on convergence andperformance for simulating CSD waves.In the first instance, we consider a low-order finite element scheme in space forthe electrodiffusive model in the zero flow limit, and numerically study the effect of:(i) different operator splitting schemes, (ii) choice of time stepping schemes for thePDEs, (iii) choice of time stepping schemes for the ODEs, and (iv) a higher orderspatial discretization scheme. All schemes display optimal convergence rates duringrefinement in space and time for problems with smooth manufactured solutions.However, we find that the accurate computation of CSD wave characteristics (suchas wave speed and wave width) requires a very fine spatial and fine temporal reso-lution for all schemes tested. In particular, the different splitting schemes and PDEtime stepping give comparable results in terms of accuracy. We observe that higherorder ODE time stepping schemes (ESDIRK4, RK4) yield slightly faster conver-gence than the lower order backward Euler scheme. Finally, we find that applyinga higher order finite element scheme with a coarser spatial resolution gives compa-rable results to the lower order finite element scheme in terms of accuracy, but at ahigher cost, justifying the use of a low order scheme with fine spatial resolution.We then turn to consider the full mathematical model, where the electrochem-ical and mechanical dynamics described in the zero flow limit are coupled withmicroscopic fluid dynamics. To the best of our knowledge, the only previous studyinvolving this model was presented by O’Connell , studying the effects of mechani-cal pressures and compartmental fluid flow on CSD wave characteristics in a settingebruary 5, 2021 13:25 WSPC/INSTRUCTION FILE ws-m3as A. J. Ellingsrud, N. Boull´e, P. E. Farrell, and M. E. Rognes A Temporal refinement Sp a t i a l r e fin e m e n t B N ∆ t v CSD v CSD – 0.178 0.154 0.128 0.085 0.049 0.024
Fig. 1. Wave properties during refinement in space (N) and time (∆t, ms) in a 1D domain oflength 10 mm at t = 50 s. The PDEs are discretized in time by the previously presented first orderscheme from Mori, O’Connell, and Tuttle et al. , , and the ODEs are solved using backwardEuler. A : Neuron potential φ n ( x,
50) (mV) versus x ∈ Ω (mm). B : CSD mean wave speed ¯ v CSD (mm/min) and difference ∆¯ v CSD between consecutive refinements. with two compartments (neurons and ECS). Here, we present a low order finiteelement scheme and simulation results for the full model in three compartments(neurons, glial and ECS) and study the convergence of the scheme. The schemeagain displays optimal convergence rates during refinement in space and time for aproblem with smooth manufactured solutions, and meets similar challenges as thezero flow limit model for the CSD case.The paper is organized as follows. The Mori framework is summarized in Sec-tion 2. In Section 3, we present new numerical schemes for the zero flow limitversion of the model, along with numerical convergence and performance studiesfor the suggested schemes in Section 4 and Section 5. The next sections considerthe full model. We present a numerical scheme in Section 6 along with numericalconvergence studies in Section 7, while in Section 8 we present results from a fullmodel simulation of CSD including microscopic fluid dynamics. Finally, Section 9ebruary 5, 2021 13:25 WSPC/INSTRUCTION FILE ws-m3as
Numerics for brain electrochemomechanics contains a discussion and concluding remarks. The code used to obtain the sim-ulation results presented within this work is based on FEniCS , and is publiclyavailable .
2. Mathematical model
The tissue of interest is represented as a domain Ω ⊂ R d , with d ∈ { , , } . More-over, the tissue is composed of R compartments indexed by r = 1 , . . . , R . We assumethat r = R always denotes the ECS, and that other compartments communicatewith the ECS only. Let time t ∈ (0 , T ]. For a self-contained exposition, we summa-rize the Mori framework below. Governing equations
We will consider the following system of coupled, time-dependent, nonlinear partialdifferential equations. Find the volume fractions α r : Ω × (0 , T ] → [0 ,
1) such thatfor each t ∈ (0 , T ]: ∂α r ∂t + ∇ · ( α r u r ) = − γ rR w rR , for r = 1 , . . . , R − , (2.1a) ∂α R ∂t + ∇ · ( α R u R ) = R − (cid:88) r =1 γ rR w rR , (2.1b)where u r : Ω × (0 , T ] → R d is the fluid velocity field of compartment r (m/s).The transmembrane water flux w rR between compartment r and the ECS is drivenby osmotic and oncotic pressure, and will be discussed further in Section 2.2. Thecoefficient γ rR (1/m) represents the area of cell membrane between compartment r and the ECS per unit volume of tissue. By definition, the volume fractions sum to1, and hence: α R = 1 − R − (cid:88) r =1 α r . (2.2)where α r ≥
0. Further, for each ion species k ∈ K , find the ion concentrations [ k ] r : Ω × (0 , T ] → R (mol/m ) and the electrical potentials φ r : Ω × (0 , T ] → R (V)such that for each t ∈ (0 , T ]: ∂ ( α r [ k ] r ) ∂t + ∇ · J kr = − γ rR J krR , r = 1 , . . . , R − , (2.3a) ∂ ( α R [ k ] R ) ∂t + ∇ · J kR = R − (cid:88) r =1 γ rR J krR , (2.3b)where J kr : Ω × (0 , T ] → R d is the ion flux density (mol/(m s)) for each ion species k . Modelling the transmembrane ion flux density (mol/(m s)) J krR for each ionspecies k ∈ K will be discussed further in Section 2.2. Note that (2.3) follow fromfirst principles and express conservation of ion concentrations for the bulk of eachebruary 5, 2021 13:25 WSPC/INSTRUCTION FILE ws-m3as A. J. Ellingsrud, N. Boull´e, P. E. Farrell, and M. E. Rognes region. Moreover, we assume that the ion flux densities (and a fortiori the ionconcentrations and electrical potentials) satisfy: γ rR C rR φ rR = z F a r + α r F (cid:88) k z k [ k ] r , r = 1 , . . . , R − , (2.4a) − R − (cid:88) r γ rR C rR φ rR = z F a R + α R F (cid:88) k z k [ k ] R , (2.4b)where z k (dimensionless) is the valence of ion species k , a r (mol/m ) is the amountof immobile ions , F (C/mol) denotes Faraday’s constant , and C rR (F/m ) is the membrane capacitance of the membrane between compartment r and the ECS.In this paper, we assume that the ion flux densities J kr can be expressed interms of the ion concentrations, the electrical potentials, the volume fractions andthe fluid velocity field as: J kr = − D kr ∇ [ k ] r − D kr z k ψ [ k ] r ∇ φ r + α r u r [ k ] r , r = 1 , . . . , R. (2.5)Here, D kr (m /s) denotes the effective diffusion coefficient of ion species k in theregion r and may be a given constant or e.g. a spatially varying scalar field dependentof α r . The constant ψ = RT F − combines Faraday’s constant F , the absolutetemperature T (K), and the gas constant R (J/(K mol)). We assume that the systemis isothermal, i.e. T is constant. The ion flux density, i.e. the flow rate of ions per unitarea, is thus modelled as the sum of three terms: (i) the diffusive movement of ionsdue to ionic gradients − D kr ∇ [ k ] r , (ii) the ion concentrations that are transportedvia electrical potential gradients, i.e. the ion migration − D kr z k ψ − [ k ] r ∇ φ r where D kr ψ − is the electrochemical mobility and (iii) the convective movement α r u r [ k ] r of ions.We now turn to the dynamics of the fluid velocity field u r and the mechanicalpressure p r : Ω × (0 , T ] → R (Pa). Let the compartmental fluid velocity u r (m/s)be expressed as: u r = − κ r (cid:32) ∇ ˜ p r + F (cid:88) k z k [ k ] r ∇ φ r (cid:33) , ˜ p r = p r − RT a r α r , r = 1 , . . . , R, (2.6)where κ r (m /(N s)) is the water permeability in compartment r . The fluid velocityfield is thus modelled as the sum of three terms: (i) the mechanical pressure gradient − κ r ∇ p r , (ii) the oncotic pressure gradient κ r ∇ ( RT a r α r ), and (iii) the electrostaticforces − κ r F (cid:80) k z k [ k ] r ∇ φ r . The mechanical pressure p r in compartment r and inthe ECS p R may be related as: p r − p R = τ r ( α r ) , r = 1 , . . . , R − , (2.7)where the mechanical tension per unit area of the membrane τ r = τ r ( α r ) is to bemodelled. A simple relation could be τ r = S r ( α r − α r ), were S r (Pa/m ) denotesthe stiffness of the membrane between compartment r and the ECS, and α r is theebruary 5, 2021 13:25 WSPC/INSTRUCTION FILE ws-m3as Numerics for brain electrochemomechanics volume fraction at which the membrane has no mechanical tension. Furthermore,we assume that the volume–fraction weighted velocity is divergence free, that is: ∇ · (cid:32) R (cid:88) r =1 α r u r (cid:33) = 0 . (2.8)Upon inserting (2.6)–(2.7) into (2.8), we obtain the following equation for the ex-tracellular mechanical pressure p R : Ω × (0 , T ] → R (Pa), ∇· (cid:32) R (cid:88) r =1 κ r α r (cid:32) −∇ p R + ∇ RT a r α r − F (cid:88) k z k [ k ] r ∇ φ r (cid:33) − R − (cid:88) r =1 κ r α r ∇ τ r (cid:33) = 0 . (2.9)The combination of (2.1a) and (2.2)–(2.4), together with the insertion of (2.5)–(2.6) and (2.9), defines a system of | R || K | + 2 | R | + 1 equations for the | R || K | +2 | R | + 1 unknown scalar fields. Appropriate initial conditions, boundary conditions,and importantly membrane mechanisms close the system. Membrane mechanisms
The transmembrane water flux w rR is driven by a combination of mechanical andosmotic pressure, and can be expressed as: w rR = η r (cid:32) p r − p R + RT (cid:32) a R α R + (cid:88) k [ k ] R − a r α r − (cid:88) k [ k ] r (cid:33)(cid:33) , (2.10)where η r (m /(mol s)) denotes the water permeability . The compartmental andextracellular potentials across the membrane are coupled by defining φ rR as thejump in the electrical potential across the membrane between compartment r andthe ECS (c.f. (2.4)), φ rR = φ r − φ R . (2.11)The transmembrane ion flux J krR between compartment r and the ECS of ion species k is subject to modelling, and will typically take the form: J krR = a krR ( φ rR , [ k ] r ) + p krR ( φ rR , [ k ] r , s , . . . , s M ) . (2.12)Here, a krR represents active membrane mechanisms (e.g. ionic pumps) and p krR de-notes passive membrane mechanisms (e.g. leak or voltage gated ion channels, co-transporters). The passive membrane mechanisms typically depend on (unitless)gating variables s m = s m ( φ rR ) for m ∈ , . . . , M governed by an ODE system ofthe form: ∂s m ∂t = α m (1 − s m ) − β m s m , (2.13)where α m = α m ( φ rR ) (1/s) and β m = β m ( φ rR ) (1/s) are rate coefficients .ebruary 5, 2021 13:25 WSPC/INSTRUCTION FILE ws-m3as A. J. Ellingsrud, N. Boull´e, P. E. Farrell, and M. E. Rognes
Boundary conditions
The boundary conditions will strongly depend on the problem of interest. If nototherwise stated, we assume that no ion flux or fluid leaves on the boundary ∂ Ω,that is, J kr ( x, t ) = 0 and u r ( x, t ) = 0 , on ∂ Ω , (2.14)for r = 1 , . . . , R . Effective diffusion coefficients
We model the effective diffusion coefficients D kr for each ion species k ∈ K by: D kr = α r χ r D k and D kR = α R D k , (2.15)for r = 1 , . . . , R −
1, where D k (m /s) denotes the diffusion coefficient in water and χ r (dimensionless) reflects the cellular gap junction connectivity . Zero flow limit
We first consider a simplified version of the mathematical model presented in Sec-tion 2, where the compartmental fluid velocity u r for r = 1 , . . . , R is assumed to bezero. Thus, the advective terms in (2.1a) and (2.5) vanish and (2.9) is decoupledfrom the rest of the system. The remaining equations, (2.1a) and (2.2)–(2.4) withthe insertion of of (2.5), describe the dynamics of the volume fractions α r , the ionconcentrations [ k ] r and the electrical potentials φ r , for r ∈ { , . . . , R } and for each k ∈ K .
3. Numerical schemes for the zero flow limit
Below, we present two numerical schemes for the zero flow limit model based onthe finite element method.
Spatial discretization
To obtain a variational formulation, we multiply (2.1a) and (2.3)–(2.4) with suit-able test functions, integrate over the domain Ω, integrate terms with higher or-der derivatives by parts, and insert the boundary condition (2.14). Below, we let (cid:104) u, v (cid:105) = (cid:82) Ω uv d x . Let S r ⊂ L (Ω), V kr ⊂ H (Ω), V kR ⊂ H (Ω), T r ⊂ H (Ω), and T R ⊂ H (Ω) for r = 1 , . . . , R − k = 1 , . . . , | K | . The re-sulting system reads: find α r ∈ S r , [ k ] r ∈ V kr , φ r ∈ T r ( r = 1 , . . . , R − k ] R ∈ V kR ,ebruary 5, 2021 13:25 WSPC/INSTRUCTION FILE ws-m3as Numerics for brain electrochemomechanics and φ R ∈ T R such that: (cid:104) ∂α r ∂t , s r (cid:105) + γ rR (cid:104) w rR , s r (cid:105) = 0 , (3.1a) (cid:104) ∂α r [ k ] r ∂t , v kr (cid:105) − (cid:104) J kr , ∇ v kr (cid:105) + γ rR (cid:104) J krR , v kr (cid:105) = 0 , (3.1b) (cid:104) ∂α R [ k ] R ∂t , v kR (cid:105) − (cid:104) J kR , ∇ v kR (cid:105) − (cid:88) r γ rR (cid:104) J krR , v kR (cid:105) = 0 , (3.1c) (cid:104) γ rR C rR φ rR , t r (cid:105) − (cid:104) z F a r , t r (cid:105) − (cid:104) F α r (cid:88) k z k [ k ] r , t r (cid:105) = 0 , (3.1d) − (cid:88) r (cid:104) γ rR C rR φ rR , t R (cid:105) − (cid:104) z F a R , t R (cid:105) − (cid:104) F α R (cid:88) k z k [ k ] R , t R (cid:105) = 0 , (3.1e)for all s r ∈ S r , v kr ∈ V kr , v kR ∈ V kR , t r ∈ T r , t R ∈ T R . The compartmental ionflux J kr is given by (2.5), the transmembrane water flux w rR is given by (2.10),while the extracellular volume fraction α R is given by (2.2). The transmembraneion fluxes J krR will depend on the membrane mechanisms of interest. Note that inthe above formulation, the potentials φ r for r = 1 , . . . , R are only determined up toa constant.Next, we discretize the domain Ω by a simplicial mesh T h , where h denotesthe mesh size. We introduce separate finite element spaces for approximating theunknown fields in the weak formulation (3.1). Specifically, we approximate the vol-ume fractions α r using piecewise constants ( P ), and the ion concentrations [ k ] r ,[ k ] R and potentials φ r , φ R using continuous piecewise linear polynomials P c definedrelative to the mesh T h . Temporal PDE discretization: BDF2 and CN
We apply a second order Backward Differentiation Formula (BDF2) in time, re-sulting in the following system: given α nr , [ k ] nr , and [ k ] nR at time level n and α n − r ,[ k ] n − r , and [ k ] n − R at time level n −
1, find at time level n + 1 the volume fractions α r ∈ S r , the ion concentrations [ k ] r ∈ V kr , [ k ] R ∈ V kR and the potentials φ r ∈ T r ,ebruary 5, 2021 13:25 WSPC/INSTRUCTION FILE ws-m3as A. J. Ellingsrud, N. Boull´e, P. E. Farrell, and M. E. Rognes φ R ∈ T R such that: 1∆ t (cid:104) α r − α nr + 13 α n − r , s r (cid:105) + 23 γ rR (cid:104) w rR , s r (cid:105) = 0 , (3.2a)1∆ t (cid:104) α r [ k ] r − α nr [ k ] nr + 13 α n − r [ k ] n − r , v kr (cid:105) − (cid:104) J k,nr , ∇ v kr (cid:105) + 23 γ rR (cid:104) J k,nrR , v kr (cid:105) = 0 , (3.2b)1∆ t (cid:104) α R [ k ] R − α nR [ k ] nR + 13 α n − R [ k ] n − R , v kR (cid:105) − (cid:104) J k,nR , ∇ v kR (cid:105)− (cid:88) r γ rR (cid:104) J k,nrR , v kR (cid:105) = 0 , (3.2c) (cid:104) γ rR C rR φ rR , t r (cid:105) − (cid:104) z F a r , t r (cid:105) − (cid:104) F α r (cid:88) k z k [ k ] r , t r (cid:105) = 0 , (3.2d) − (cid:88) r (cid:104) γ rR C rR φ rR , t r (cid:105) − (cid:104) z F a R , t R (cid:105) − (cid:104) F α R (cid:88) k z k [ k ] R , t R (cid:105) = 0 , (3.2e)for all s r ∈ S r , v kr ∈ V kr , v kR ∈ V kR , t r ∈ T r , t R ∈ T R . The solutions at time level 1and 0 are given by a backward Euler step and the initial conditions, respectively.Note that the passive part of the membrane flux J k,nrR is treated implicitly, whereasthe active part is treated explicitly, both in the BDF2 and the BE time steppingschemes (see (2.12)).We will also compare with the following Crank-Nicholson (CN) scheme in time:given α nr , [ k ] nr , and [ k ] nR at time level n , find at time level n + 1 the volume fractions α r ∈ S r , the ion concentrations [ k ] r ∈ V kr , [ k ] R ∈ V kR , and the potentials φ r ∈ T r , φ R ∈ T R , such that 1∆ t (cid:104) α r − α nr , s r (cid:105) + γ rR (cid:104) w n +1 / rR , s r (cid:105) = 0 , (3.3a)1∆ t (cid:104) α r [ k ] r − α nr [ k ] nr , v kr (cid:105) − (cid:104) J k,n +1 / r , ∇ v kr (cid:105) + γ rR (cid:104) J k,n +1 / rR , v kr (cid:105) = 0 , (3.3b)1∆ t (cid:104) α R [ k ] R − α nR [ k ] nR , v kR (cid:105) − (cid:104) J k,n +1 / R , ∇ v kR (cid:105) − (cid:88) r γ rR (cid:104) J k,n +1 / rR , v kR (cid:105) = 0 , (3.3c)together with (3.2d) and (3.2e), for all s r ∈ S r , v kr ∈ V kr , v kR ∈ V kR , t r ∈ T r , t R ∈ T R Here, f n +1 / denotes the solution at time level n + 1 / f n + f n +1 ) /
2, for f ∈ { w rR , J krR , J kr } . Strang splitting scheme
We use a second order Strang splitting scheme where we solve the coupled systemsof ODEs and PDEs step-wise in the following manner:(1) Insert the previous solution of φ rR , [ k ] r , [ k ] R into the system ofODEs (2.13), and solve ODEs for a half timestep ∆ t/ s , s , . . . , s M from step (1) for one time step ∆ t ;ebruary 5, 2021 13:25 WSPC/INSTRUCTION FILE ws-m3as Numerics for brain electrochemomechanics (3) Insert solution for φ rR , [ k ] r , [ k ] R from step (2) into the system ofODEs (2.13), and solve the ODEs for a half timestep ∆ t/ for a detailed description the Godunov method.. ODE solvers
We consider three different schemes for the ODE time-stepping: a fourth orderRunge-Kutta (RK4) method, a fourth order explicit singly diagonal implicit Runge-Kutta (ESDIRK4) method or a second order backward Euler (BE) method. TheODEs are solved with the RK4 method unless otherwise stated. See e.g. Langtangenand Linge for details.
4. Numerical convergence study: smooth analytical solution
To evaluate the numerical accuracy of the various schemes presented above, we beginby constructing a smooth analytical solution using the method of manufacturedsolutions . Problem description
We consider a two-compartment version of the model in the zero flow limit anda neuronal ( n , r = 1) and an extracellular ( e , r = 2) compartment. We use 1 , n, e interchangeably for subscripts of our variables and model parameters. Ineach compartment, we model the movement of potassium (K + ), sodium (Na + ) andchloride (Cl − ). In the numerical experiments of this test, we consider a one dimen-sional domain Ω = [0 ,
1] uniformly meshed with N ∈ { , , , , } elements.We initially set ∆ t = 10 − s, and then halve the timestep with each spatial re-finement. The errors are evaluated at t = 2 × − s. Further, we assume that thetransmembrane ion flux J kne depends on the gating variables m , h and g governedby the following ODEs: ∂m∂t = φ ne , ∂h∂t = φ ne , ∂g∂t = φ ne , (4.1)where φ ne = φ n − φ e is the membrane potential. The analytical solution to the PDEsystem is given by: α n = 0 . − . πx ) exp( − t ) , [Na + ] n = 0 . . πx ) exp( − t ) , [Na + ] e = 1 . . πx ) exp( − t ) , [K + ] n = 0 . . πx ) exp( − t ) , [K + ] e = 1 . . πx ) exp( − t ) , [Cl − ] n = 1 . . πx ) exp( − t ) , [Cl − ] e = 2 . . πx ) exp( − t ) ,φ n = sin(2 πx ) exp( − t ) , φ e = sin(2 πx )(1 + exp( − t )) , (4.2)ebruary 5, 2021 13:25 WSPC/INSTRUCTION FILE ws-m3as A. J. Ellingsrud, N. Boull´e, P. E. Farrell, and M. E. Rognes and the solution to the system of ODEs (4.1) is m = cos( t ) cos( πx ) , h = cos( t ) cos( πx ) , g = cos( t ) cos( πx ) . (4.3)Parameter values are given in Table 1. Initial and boundary conditions are governedby the exact solutions (4.2).Parameter Symbol Value Unit Ref.Temperature T
310 K –Faraday’s constant F R . γ ne . × Membr. capacitance neuron C ne . × − F/m
Membr. water permeability neuron η ne . × − m /(mol s) Gap junction connectivity neuron χ n Diffusion coefficient Na + D Na . × − m /s Diffusion coefficient K + D K . × − m /s Diffusion coefficient Cl + D Cl . × − m /s Valence Na + z Na + z K − z Cl − Table 1. Physical model parameters. We use SI base units, that is, Kelvin (K), Coulomb (C), mole(mol), meter (m), second (s), and Joule (J). The values are collected from Hille et al. , Kager etal. , Mori , and O’Connell et al. . – indicates that a standard value is used. Convergence and convergence rates under refinement
Based on the approximation spaces and the time discretization, we expect the opti-mal theoretical rate of convergence to be 1 in the H -norm and 2 in the L -norm forthe concentrations [ k ] n , [ k ] e , the potentials φ n , φ e and the gating variables m, h, g ,and 1 in the L -norm for the volume fraction α n . Our numerical observations arein agreement with these optimal rates, both for the BDF2 scheme (Table 2A) andfor the CN scheme (Table 2B). We observe second order convergence in the L -norm for the approximation of the extracellular and intracellular concentrationsand potentials, and first order convergence in the H -norm. The neuron potentialapproximated by the CN scheme displays superconvergence between N = 64 and N = 128 (Table 2B). For the volume fractions, we observe a convergence rate of 1in the L -norm.ebruary 5, 2021 13:25 WSPC/INSTRUCTION FILE ws-m3as Numerics for brain electrochemomechanics A N (cid:107) [Na] e − [Na] eh (cid:107) L (cid:107) φ n − φ nh (cid:107) L (cid:107) α n − α nh (cid:107) L (cid:107) m − m h (cid:107) L N (cid:107) [Na] e − [Na] eh (cid:107) H (cid:107) φ n − φ nh (cid:107) H B N (cid:107) [Na] e − [Na] eh (cid:107) L (cid:107) φ n − φ nh (cid:107) L (cid:107) α n − α nh (cid:107) L (cid:107) m − m h (cid:107) L N (cid:107) [Na] e − [Na] eh (cid:107) H (cid:107) φ n − φ nh (cid:107) H Table 2. Selected L (upper panel) and H -errors (lower panel) and convergence rates (in paren-thesis) for the BDF2 ( A ) and CN ( B ) schemes at time t = 0 .
002 s. The test was run on theunit interval, and we initially let ∆ t = 1 × − s, and then halve the timestep with each meshrefinement. The spatial discretization consists of N intervals.
5. Numerical convergence study: physiological CSD wave
Next, we consider the simulation of cortical spreading depression with a sharp wavefront. This is a more challenging problem, with characteristics quite different fromthe previous smooth MMS case. We numerically study the effect of splitting scheme,time-discretization of the PDEs, and discretization of the ODEs.ebruary 5, 2021 13:25 WSPC/INSTRUCTION FILE ws-m3as A. J. Ellingsrud, N. Boull´e, P. E. Farrell, and M. E. Rognes
Problem description
Parameter Symbol Value Unit Ref.Na + leak conductance neuron g Na n, leak . K + leak conductance neuron g K n, leak . Cl − leak conductance neuron g Cl n, leak . Maximum pump rate neuron ˆ I n . Threshold for pump [Na + ] r m Na . Threshold for pump [K + ] R m K . Table 3. Physical parameters for the neuron membrane mechanisms. We use SI base units, i.e. meter(m), mole (mol), Siemens (S) and ampere (A). The values are collected from Kager et al. ,O’Connell et al. , and Yao et al. . We define a more physiological version of the mathematical model considered inSection 4.1 (two compartments, zero flow limit, and a neuronal ( n , r = 1) and anextracellular compartment ( e , r = 2)). In each compartment, we again model themovement of potassium (K + ), sodium (Na + ) and chloride (Cl − ). We consider a 1Ddomain of length 0 .
01 m (10 mm), and now apply physiologically relevant neuronalmembrane mechanisms as described below, notably including a system of ODEsdescribing the gating variables of voltage-gated sodium and potassium channels.The domain and the transmembrane ion flux densities are taken from the originalMori study .Concretely, the transmembrane ion flux densities J kne for k = { Na + , K + , Cl − } in (2.3) are modelled as: J Na ne = 1 F z Na (cid:0) I Na n, leak + I NaP + 3 I n, ATP + I Naex (cid:1) , (5.1a) J K ne = 1 F z K (cid:0) I K n, leak + I KDR + I KA − I n, ATP + I Kex (cid:1) , (5.1b) J Cl ne = 1 F z Cl (cid:0) I Cl n, leak + I Clex (cid:1) , (5.1c)where I Na n, leak , I K n, leak , I Cl n, leak , I NaP , I KDR , I KA , and I n, ATP denote the sodium leakcurrent, the potassium leak current, the chloride leak current, the persistent sodiumcurrent, the potassium delayed rectifier current, the transient potassium current,and the Na/K/ATPase current, respectively. Further, I Naex , I Kex , and I Clex are excita-tory fluxes used to trigger a cortical spreading depression wave and whose expres-sions are given by (5.5) in Section 5.1.1 below. Note that the currents (A/m ) areconverted to ion fluxes (mol/(m s)) by dividing by Faraday’s constant F times thevalence z k .The leak currents (A/m ) of ion species k over the membrane between compart-ebruary 5, 2021 13:25 WSPC/INSTRUCTION FILE ws-m3as Numerics for brain electrochemomechanics ment r (here with r = n ) and R are modelled as: I kr, leak = g kr, leak ( φ rR − E kr ) , E kr = RTF z k ln (cid:18) [ k ] R [ k ] r (cid:19) , (5.2)where E rk (V) denotes the Nernst potential. The values for the neuronal leak con-ductances g kn, leak are listed in Table 3.The current-voltage relation for the voltage-gated currents I NaP , I KDR and I KDR (A/m ) are described by the Goldman-Hodgkin-Katz (GHK) current equation: I k GHK = g k m p h q F µ ([ k ] e − [ k ] n e − µ )1 − e − µ . Here, µ = F φ ne / ( RT ) is dimensionless, while g k (m/s) denotes the product ofmembrane permeability and conductance. The gating variables m and h describethe proportion of open voltage–gated ion channels, and are governed by the followingODEs: ∂m∂t = α m ( φ ne )(1 − m ) − β m ( φ ne ) m, (5.3a) ∂h∂t = α h ( φ ne )(1 − h ) − β h ( φ ne ) h, (5.3b)where the activation rate functions α m : R → R and β m : R → R , and the inacti-vation rate functions α h : R → R and β h : R → R , are specified for each current inTable 4. The initial conditions are given in Table 12A (Supplementary Tables).Current Permeability Gates Voltage dependent rate constants(A/m ) (m/s) I NaP . × − m h α m = 1 / (6 + 6 exp( − . φ ne − . β m = 1 / − α m α h = 5 . × − exp( − . φ ne − . β h = 1 . × − / (1 + exp( − . φ ne − I KDR . × − m α m = 0 . − φ ne − . − . φ ne − . − β m = 0 .
25 exp( − . φ ne − . I KA . × − m h α m = 0 . − φ ne − . − . φ ne − . − β m = 0 . φ ne +29 . . φ ne +29 . − α h = 0 .
016 exp( − . φ ne − . β h = 0 . / (exp( − . φ ne − .
98) + 1)
Table 4. Permeability, gates and voltage dependent expressions for the activation rates ( α ) andthe inactivation rates ( β ) for the persistent sodium current ( I NaP ), the potassium delayed rectifiercurrent ( I KDR ) and the transient potassium current ( I KA ) in the neuron membrane. The valuesare collected from Kager et al. and Yao et al. . ebruary 5, 2021 13:25 WSPC/INSTRUCTION FILE ws-m3as A. J. Ellingsrud, N. Boull´e, P. E. Farrell, and M. E. Rognes
The Na/K/ATPase pump exchanges 2 potassium ions for 3 sodium ions, andthe pump current I r, ATP (A/m ) over the membrane between compartment r (herewith r = n ) and R is modelled as: I r, ATP = ˆ I r (1 + m K [K + ] R ) (1 + m Na [Na + ] r ) , (5.4)where ˆ I r is the maximum pump rate, and m Na and m K denote the sodium andpotassium pump threshold, respectively. The values for the neuron pump parame-ters are listed in Table 3.5.1.1. Initiation of the CSD wave
Following the original Mori study , we initiate a CSD wave by adding excitatoryfluxes to the transmembrane fluxes defined by (5.1) in the following manner: I k ex = G ex ( E kr − φ nR ) ,G ex ( x, t ) = (cid:40) G max cos ( πx/ L ex ) sin( πt/T ex ) if x ≤ L ex and t ≤ T ex ,0 otherwise, (5.5)for k = { Na + , K + , Cl − } . We set L ex = 2 . × − m, T ex = 2 s, and G max = 5 . . CSD wave characteristics
The excitatory flux stimulation leads to a wave of neuronal depolarization, ionicconcentration changes and neuronal swelling spreading through the tissue domain(Figure 2). We observe a dramatic depolarization of the neuron potential: from −
71 mV to − . − . . The neuronal depolarization wave isfollowed by a complete break-down of the ionic homeostatis: a substantial increase inthe concentrations of extracellular potassium, and decreases in extracellular sodiumand chloride. In response to the ionic shifts, the neurons swell with an increase involume fraction of up to 10%, while the extracellular space shrinks correspondingly.We observe that although we have not enforced positivity of α r explicitly, α r ≥ • The point of x peak , of peak neuron potential φ n at t = 50 s. • The mean CSD wave speed ¯ v CSD is computed by the distance between the peakneuronal potential at two points (as long as φ n > −
20 mV) divided by the timeelapsed to cover that distance.ebruary 5, 2021 13:25 WSPC/INSTRUCTION FILE ws-m3as
Numerics for brain electrochemomechanics + ),potassium (K + ) and chloride (Cl − ) concentrations ( A , B , C ), neuron and extracellular potentials( D , E , F ), and neuron and extracellular changes in volume fractions ( G , H , I ) – at, from left toright, t = 10 , ,
50 s. Numerical scheme and resolutions: BDF2, Strang, RK4 with N = 32000and ∆ t = 0 .
195 s. • The (temporal) duration d CSD of the CSD wave in terms of elevated extracellularpotassium levels at x = 1 mm, where T = { t | [K + ] R (1 , t ) > k thres } d CSD = max T − min T with k thres = 10 mM for t ∈ (0 , T ].ebruary 5, 2021 13:25 WSPC/INSTRUCTION FILE ws-m3as A. J. Ellingsrud, N. Boull´e, P. E. Farrell, and M. E. Rognes • The (spatial) width w CSD of the extracellular potassium wave at t = 50 s, where X = { x | [K + ] R ( x, > k thres } w CSD = max X − min X with k thres = 10 mM for x ∈ Ω. Convergence of the CSD wave characteristics
We begin by considering a reference scheme – based on Strang splitting, a BDF2method for the PDE time-discretization, and ESDIRK4 for the ODE time-stepping– to compute the mean speed, the (spatial) width, and the (temporal) duration ofthe CSD wave (cf. Section 5.2) during temporal and spatial refinement, before dis-cussing how variations in terms of splitting scheme, time-stepping and higher orderspatial discretization affect accuracy and convergence. Specifically, we apply the ref-erence scheme and calculate the quantities of interest for different mesh resolutionsand time steps: ∆ x N = 10 /N mm for N = 1000 , , , , , t i = 12 . /i s for i = 1 , , , , , ,
64. The results are presented in Figure 3.Wave speeds are converted from the native m/s to mm/min for interpretability.In general, the computed mean wave speed and width increases with decreasing∆ t , and decreases with decreasing ∆ x : the smaller the time step, the faster andwider the wave, while the smaller the mesh size, the slower and narrower the wave(Figure 3A). Regarding the mean wave speed (Figure 3B), we observe that the com-puted values vary substantially, ranging from 5 .
063 to 10 .
090 mm/min. We observethat the spatial errors, estimated by proxy by the difference between consecutivespatial refinements, dominate the temporal errors/differences. In particular, the dif-ferences ∆¯ v CSD in wave speed between the coarsest mesh sizes N = 1000 and 2000are in the range 1 . − .
45 mm/min. Conversely, the differences in wave speedbetween the coarsest time steps ∆ t = 12 . .
25 are in the range 0 . − . .
012 and 0 .
008 mm/min, respectively. Finally, we observe that the wavespeed seems to converge in space and time with an estimated wave speed of 5 . v CSD reduce as ∆ x and ∆ t is reduced. There is howeverno clear rate of convergence.Similarly, we observe large variations in the computed (spatial) CSD wave width,ranging from 1 .
121 to 4 .
420 mm (Figure 3C). The difference (∆ w CSD ) varies in therange 0 . − .
07 mm and 0 . − .
280 mm between respectively the coarsestmesh sizes N = 1000 and 2000 and the coarsest time steps ∆ t = 12 . .
25 s.For the finest time and mesh discretizations, we observe differences of 0 .
005 mmand 0 .
002 mm, respectively. Moreover, the differences ∆ w CSD reduce as the meshsize and time step are reduced. As with the mean wave speed, there is no clearrate of convergence. In contrast, the (temporal) duration d CSD of the CSD wavedoes not change substantially during refinement in space and time: we observe aduration of elevated extracellular potassium levels of 17 −
18 s for all the spatialebruary 5, 2021 13:25 WSPC/INSTRUCTION FILE ws-m3as
Numerics for brain electrochemomechanics and temporal resolutions considered (results not shown). Finally, we observe thatthis implicit higher-order (reference) scheme behaves qualitatively similar to theimplicit lower-order scheme – based on Godunov splitting, a BE method for thePDE time-discretization, and BE for the ODE time-stepping as shown in Figure 1. Choice of discretizations
We turn to evaluate the effect of discretization choices in terms of splitting scheme,time-stepping, and higher order spatial discretization.
Splitting scheme
To evaluate the second order Strang splitting scheme, we compared the computedCSD wave characteristics, peak neuron potentials, and wave speeds with those com-puted using a first order Godunov splitting scheme (cf. Section 3.3). The computedwave speeds are comparable at given resolutions, and we observe a similar differ-ence decay (∆¯ v ) for the Strang and Godunov splitting schemes (see Table 10 inSupplementary Tables). Time stepping
To assess how the choice of PDE time stepping affects the convergence of the nu-merical schemes, we repeat the convergence study presented above in Section 5.3replacing the BDF2 scheme by a Crank-Nicolson (CN) scheme and ESDIRK4 byRK4 for the ODE time stepping (Table 5). We note that choosing an explicit ODEtime stepping scheme (RK4) here is based on the observation that CN for the PDEtime stepping in combination with implicit ODE time stepping schemes results ina diverging (non-linear) ODE solver.Importantly, we note that the PDE Newton solver fails to converge for ∆ t ≥ .
563 ms for this scheme. Again, we observe that the computed speed and width ofthe wave decreases with decreasing ∆ x , and is essentially constant with decreasing∆ t (likely due to the already fine timestep required for convergence of the PDENewton solver). The spatial errors, estimated by proxy by the difference betweenconsecutive spatial refinements, are comparable to those reported for BDF2, andagain we obtain an estimated wave speed of 5 . t ≥ .
125 ms (see Supplementary Tables, Ta-ble 11A). Also for the BE scheme, the computed wave speed values are similar,ebruary 5, 2021 13:25 WSPC/INSTRUCTION FILE ws-m3as A. J. Ellingsrud, N. Boull´e, P. E. Farrell, and M. E. Rognes A Temporal refinement Sp a t i a l r e fin e m e n t B N ∆ t v CSD v CSD – 0.144 0.073 0.038 0.10 0.013 0.008 C N ∆ t w CSD w CSD – 0.062 0.032 0.016 0.005 0.006 0.002
Fig. 3. CSD wave characteristics and quantities of interest under refinement in space (rows) andtime (columns). A : Neuron potential φ n ( x,
50) (mV) versus x ∈ Ω (mm), where green, yellow, andred represent wave speeds differing by respectively ± ± ±
15% from ourestimated wave speed of 5 . B : CSD mean wave speed ¯ v CSD (mm/min) and difference∆¯ v CSD between consecutive refinements. C : Wave width w CSD (mm) at t = 50 s and difference∆ w CSD between consecutive refinements. Numerical scheme: Strang splitting, BDF2, ESDIRK4. ebruary 5, 2021 13:25 WSPC/INSTRUCTION FILE ws-m3as
Numerics for brain electrochemomechanics A N ∆ t v CSD ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ v CSD – – – – – 0.000 0.001 B N ∆ t w CSD ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ w CSD – – – – – 0.000 0.000
Table 5. CSD wave characteristics and quantities of interest under refinement in space (rows)and time (columns). A : CSD mean wave speed ¯ v CSD (mm/min) and difference ∆¯ v CSD betweenconsecutive refinements. B : Wave width w CSD (mm) at t = 50 s and difference ∆ w CSD betweenconsecutive refinements. Numerical scheme: Strang splitting, CN, RK4. ∗ indicates that the solverfailed to converge. but with slightly bigger differences in terms of differences between resolutions (seeSupplementary Tables, Table 11B). Higher order spatial discretization
To assess how the polynomial degree of the finite elements affects the convergenceof the numerical schemes, we repeat the convergence study presented in Section 5.3replacing the lowest order finite element pairings by higher order pairings. We con-sider the model described in Section 2 discretized with discontinuous piecewiselinear polynomials ( P ) for the volume fractions α r and continuous piecewise lin-ear polynomials of degree 2 ( P c ) for the ion concentrations and potentials. Thefinite element software used within this work (FEniCS ) has automated function-ality for solving coupled PDE-ODE systems via the PointIntegralSolver class. Theebruary 5, 2021 13:25 WSPC/INSTRUCTION FILE ws-m3as A. J. Ellingsrud, N. Boull´e, P. E. Farrell, and M. E. Rognes
PointIntegralSolver solves the ODEs at the vertices of the elements and only sup-ports the use of elements with degrees of freedom located at the vertices. Higherorder elements are thus not supported. To assess the accuracy and convergence ofthe higher order scheme, we here consider an alternative implementation where theODEs are solved at each degree of freedom of the spatial discretization using abackward Euler scheme and a Newton solver. A P - P c P - P c N ¯ v CSD ∆¯ v CSD ¯ v CSD ∆¯ v CSD B P - P c P - P c N w
CSD ∆ w CSD w CSD ∆ w CSD
Table 6. Comparison of CSD wave characteristics and quantities of interest under refinement inspace between the P - P c spatial discretization described in Section 3.1 and a P - P c discretization(both solved using the implementation allowing for higher order elements). A : CSD mean wavespeed ¯ v CSD (mm/min) and difference ∆¯ v CSD between consecutive refinements. B : Wave width w CSD (mm) at t = 50 s and difference ∆ w CSD between consecutive refinements. Numerical scheme:Strang splitting, BDF2, BE using ∆ t = 3 .
125 ms.
In Table 6, we compare the CSD wave characteristics (wave speed ¯ v CSD andwave width w CSD ) computed after solving the system with the two different spa-tial discretizations: P - P c and P - P c . The numerical scheme consists of a secondorder Strang splitting, together with a BDF2 time-stepping scheme for the PDEand a backward Euler (BE) scheme for the system of ODEs, with a timestep of∆ t = 3 .
125 ms. We observe that the wave speed and wave width computed us-ing a second order P - P c spatial discretization with a mesh of N elements (whereebruary 5, 2021 13:25 WSPC/INSTRUCTION FILE ws-m3as Numerics for brain electrochemomechanics φ N ( m V ) x (mm) -4 -2 Degree of polynomial M ag n i t ud e o f c o e ffi c i e n t Fig. 4. Left: Neuron potential φ n (blue) together with its polynomial approximation of degree 100(red) computed using the Chebfun software . Right: Magnitude of the Chebyshev coefficients ofthe polynomial approximation to φ n . N ∈ { , , , , } ) are similar to the quantities obtained by a P - P c discretization with 2 N elements (i.e. the same number of degrees of freedom).We also study the approximation of the neuron potential by arbitrary higher-order polynomials to compare low-degree finite element spatial discretization againsthigher order discretizations or global spectral methods . First, the system of equa-tions is solved using a P - P c finite element scheme with a large spatial discretizationof N = 32000. We then use the Chebfun software to approximate the neuron po-tential with Chebyshev polynomials , .In Figure 4 (left), we display the original function φ n (in blue), together with itspolynomial approximation of degree 100 (in red). We see that the approximationhas large oscillating errors near the wave front located at x = 4 mm, which is likelydue to Wilbraham–Gibbs phenomenon . On the right panel of Figure 4, we reportthe coefficients of the polynomial approximation to φ n of degree 100 constructed byChebfun. The magnitude of the Chebyshev coefficients seems to decay algebraicallyat a rate of − .
2, suggesting that a very large polynomial degree (of order 10 − )would be required to approximate the solution to (2.1a) and (2.3)–(2.4) accurately. Numerical performance for the zero flow limit
So far, we have investigated the convergence of the different schemes via studyingkey functionals of the solution. To evaluate the performance and scalability of theimplementation of the schemes for the model in the zero flow limit, we consideran additional set of experiments measuring the memory usage and CPU timings a .We evaluate both the standard implementation and the implementation allowingfor higher order finite elements. For the standard implementation, we consider sim-ulations with BDF2 and ESDIRK4 time stepping for the PDEs and the ODEs, a Timings were performed on a Lenovo ThinkPad 2.70GHz x 4 Intel Core i7-7500U CPU usingFEniCS 2019.1.0 without parallelization. ebruary 5, 2021 13:25 WSPC/INSTRUCTION FILE ws-m3as A. J. Ellingsrud, N. Boull´e, P. E. Farrell, and M. E. Rognes respectively, and compare the performance of Godunov and Strang splitting. Forthe implementation allowing for higher order finite elements we consider simula-tions with Strang splitting, BDF2 time PDE stepping, BE time ODE stepping, andcompare performance of using P - P c and P - P c elements.For both the Godunov and Strang splitting schemes and the standard imple-mentation, we observe that the memory usage increases linearly with the number ofdegrees of freedom: doubling the number of degrees of freedom leads to an increasein memory of a factor 2 (see Table 7A and B). We observe that the CPU time forthe simulations grows linearly with the number of degrees of freedom: doubling thenumber of degrees of freedom leads to an increase in total simulation time of a factor2. In the Strang splitting scheme, the ODEs are solved twice for each PDE step andwe thus expect the total ODE stepping time to be greater than for the Godunovsplitting scheme (where the ODEs are only solved once per PDE step). Indeed, thetotal ODE stepping time is about twice as large for Strang splitting as for Godunovsplitting (Table 7). Conversely, the time required for finite element assembly andLU solves is comparable for Godunov splitting and Strang splitting. In total, thesimulation time is higher for Strang splitting than Godunov (11 − P - P c discretization (with the imple-mentation allowing for higher order schemes) is reported in Table 7D. For a givennumber of mesh elements N ∈ { , , } , we compare the timings withthe P - P c discretization with a mesh of 2 N elements (Table 7C) so that the under-lying systems have the same number of degrees of freedom as well as approximatelythe same accuracy, following the linear convergence of the P - P c discretization asobserved in Section 5.4. We then see that the total time needed to run the full simu-lation is approximately 21 −
40% higher for the P - P c discretization. This differenceis explained by the higher density of the finite element matrices resulting from the P - P c discretization compared to the choice of P - P c .In conclusion, Strang and Godunov splitting yield comparable accuracy andmemory usage. We observe that Strang splitting yields higher total CPU time thanGodunov. When varying (ODE and PDE) time stepping schemes we observe minorvariations in terms of accuracy. The higher order element scheme ( P - P c ) has both ahigher total CPU time and memory usage than the lower order spatial scheme ( P - P ), whereas the accuracy is comparable. We find that the accurate computation ofCSD wave characteristics (wave speed and wave width) requires a very fine spatialand fine temporal resolution for all schemes tested.ebruary 5, 2021 13:25 WSPC/INSTRUCTION FILE ws-m3as Numerics for brain electrochemomechanics A N Dofs M (MiB) T A (s) T LU (s) T PDE (s) T
ODE (s) T tot (s)8000 72008 117 329 350 679 129 81316000 144008 199 659 707 1366 258 163232000 288008 396 1438 1520 2958 559 3526 B N Dofs M (MiB) T A (s) T LU (s) T PDE (s) T
ODE (s) T tot (s)8000 72008 118 333 352 685 259 95416000 144008 199 677 719 1396 536 193232000 288008 396 1367 1451 2818 1086 3916 C N Dofs M (MiB) T A (s) T LU (s) T PDE (s) T
ODE (s) T tot (s)8000 72008 130 – – 696 253 95316000 144008 223 – – 1411 506 192232000 288008 446 – – 3076 1115 4198 D N Dofs M (MiB) T A (s) T LU (s) T PDE (s) T
ODE (s) T tot (s)4000 72008 165 – – 1017 309 13318000 144008 285 – – 1693 511 220916000 288008 567 – – 3902 1182 5091
Table 7. CPU timings and memory usage for approximating solutions in the zero flow limit.Dofs: number of degrees of freedom in the linear (PDE) system, M: Maximal memory usage ofsimulation relative to baseline. T A : CPU time for finite element assembly, T LU : CPU time forLU solver, T ODE : CPU time for ODE stepping, T
PDE : CPU time for PDE stepping (in A and B this equals the sum of T A and T LU ), T tot : Total CPU time for simulation. T A and T LU arenot reported for C and D . All simulations have ∆ t = 3 .
125 ms and final time T = 5 s (i.e. 1600timesteps). Results from simulations with the standard implementation with BDF2, ESDIRK4, P - P c , and either A : Godunov splitting; or B : Strang splitting. Results from simulations with theimplementation allowing for higher order elements with BDF2, BE, Strang splitting, and either C : P - P c elements; or D : P - P c elements. ebruary 5, 2021 13:25 WSPC/INSTRUCTION FILE ws-m3as A. J. Ellingsrud, N. Boull´e, P. E. Farrell, and M. E. Rognes
6. Numerical solution of model including fluid dynamics
The previous schemes and experiments considered only the zero flow limit. Here, wepresent a numerical scheme for the full mathematical model. The variational formu-lation (6.1) is obtained by multiplying (2.1)–(2.4) and (2.8) by suitable test func-tions, integrating over the domain Ω, integrating terms with higher order derivativesby parts, and inserting the boundary conditions (2.14). Further, we use a backwardEuler scheme for the PDE time discretization. Let S r ⊂ H (Ω), V kr ⊂ H (Ω), V kR ⊂ H (Ω), T r ⊂ H (Ω), T R ⊂ H (Ω), and Q ⊂ H (Ω) be spaces of functionsfor r = 1 , . . . , R − k = 1 , . . . , | K | . Given α nr , [ k ] nr , and [ k ] nR at time level n ,at each time level n + 1 find the volume fractions α r ∈ S r , the ion concentrations[ k ] r ∈ V kr , [ k ] R ∈ V kR , the potentials φ r ∈ T r , φ R ∈ T R , and the extracellularmechanical pressure p R ∈ Q such that:1∆ t (cid:104) α r − α nr , s r (cid:105) − (cid:104) α r u r , ∇ s r (cid:105) + γ rR (cid:104) w rR , s r (cid:105) = 0 , (6.1a)1∆ t (cid:104) α r [ k ] r − α nr [ k ] nr , v kr (cid:105) − (cid:104) J kr , ∇ v kr (cid:105) + γ rR (cid:104) J krR , v kr (cid:105) = 0 , (6.1b)1∆ t (cid:104) α R [ k ] R − α nR [ k ] nR , v kR (cid:105) − (cid:104) J kR , ∇ v kR (cid:105) − R − (cid:88) r γ rR (cid:104) J krR , v kR (cid:105) = 0 , (6.1c) (cid:104) γ rR C rR φ rR , t r (cid:105) − (cid:104) z F a r , t r (cid:105) − (cid:104) F α r (cid:88) k z k [ k ] r , t r (cid:105) = 0 , (6.1d) − R − (cid:88) r (cid:104) γ rR C rR φ rR , t R (cid:105) − (cid:104) z F a R , t R (cid:105) − (cid:104) F α R (cid:88) k z k [ k ] R , t R (cid:105) = 0 , (6.1e) −(cid:104) R (cid:88) r α r u r , ∇ q (cid:105) = 0 , (6.1f)for all s r ∈ S r v kr ∈ V kr , v kR ∈ V kR , t r ∈ T r , t R ∈ T R , q ∈ Q . The compartmentalion flux J kr is given by (2.5), the compartmental fluid velocity u r is given by (2.6),the transmembrane water flux w rR is given by (2.10), while the transmembrane ionfluxes J krR are subject to modelling. As before, the potentials φ r for r = 1 , . . . , R andthe extracellular mechanical pressure p R are only determined up to a constant. Weemploy continuous piecewise linear elements for all variables. Note that α r ∈ L (Ω)is sufficient for the weak formulation of the zero flow limit model to be well-defined,whereas in the full model, the gradient of α r appears via the expression for thecompartmental fluid velocities u r in (2.6), thus suggesting α r ∈ H (Ω). For the ODEtime stepping, we apply an ESDIRK4 scheme and first order Godunov splitting.
7. Numerical convergence study with fluid dynamics: smoothanalytical solution
To evaluate the numerical accuracy of the scheme presented above in Section 6, weconstruct an analytical solution using the method of manufactured solutions . Weebruary 5, 2021 13:25 WSPC/INSTRUCTION FILE ws-m3as Numerics for brain electrochemomechanics consider the full model with three compartments, namely a neuronal ( n , r = 1),a glial ( g , r = 2) and an extracellular compartment ( e , r = 3). We use 1 , , n, g, e interchangeably for subscripts of our variables and model parameters.In each compartment, we model the movement of potassium (K + ), sodium (Na + )and chloride (Cl − ). The transmembrane ion fluxes J krR are taken to be passive leakfluxes, i.e.: J krR = 1 F z k I kr, leak , (7.1)where I kr, leak is defined by (5.2) for compartment r and ion species k . The neuronalleak conductances g kn, leak are given in Table 3, and the glial leak conductances g kg, leak for Na + and Cl − are given in Table 8, whereas the glial K + conductance is givenby (8.1). We consider a one dimensional domain Ω = [0 ,
1] uniformly meshed with N ∈ { , , , , } elements. We initially set ∆ t = 10 − s, and then halve thetimestep with each spatial refinement. The errors are evaluated at t = 2 × − s. Theanalytical solutions are given by (4.2)–(4.3) for the neuronal and extracellular tissuevariables, and by the following for the glial tissue variables and the extracellularmechanical pressure: α g = 0 . − . πx ) exp( − t ) , [Na] g = 0 . . πx ) exp( − t ) , [K] g = 0 . . πx ) exp( − t ) , [Cl] g = 1 . . πx ) exp( − t ) ,φ g = sin(2 πx ) exp( − t ) , p R = sin(2 πx ) exp( − t ) . (7.2)Parameters values are given in Table 1 and Table 8. Initial and boundary conditionsare governed by the exact solutions (4.2) and (7.2).Based on properties of the approximation spaces and the time discretization, weexpect the optimal rate of convergence to be 1 in the H -norm and 2 in the L -norm for the volume fractions α n , α g , α e , the ion concentrations [ k ] n , [ k ] g , [ k ] e , thepotentials φ n , φ g , φ e and the mechanical pressure p e . Our numerical observationsare in agreement with the theoretically optimal rates (Table 9).
8. Numerical convergence study: physiological CSD model withmicroscopic fluid mechanics8.1.
Problem description
Finally, we consider the full model with a neuronal ( n , r = 1), glial ( g , r = 2) andextracellular compartment ( e , r = 3) and three ion species, namely potassium (K + ),sodium (Na + ) and chloride (Cl − ). Cortical spreading depression is triggered byapplying excitatory fluxes to the neurons as described in (5.5) in the one-dimensionaldomain of length 10 mm. The physiological parameters values are given in Table 1and Table 8, whereas the initial conditions are given in Table 12B (SupplementaryTables).The neuronal membrane mechanisms are as described in Section 5.1. For theglial membrane mechanisms we follow O’Connell and consider leak currents mod-elled as in (5.2) (with r = g ) for sodium (Na + ) and chloride (Cl − ), and a potassiumebruary 5, 2021 13:25 WSPC/INSTRUCTION FILE ws-m3as A. J. Ellingsrud, N. Boull´e, P. E. Farrell, and M. E. Rognes
Parameter Symbol Value Unit Ref.Na + leak conductance glial g Na g, leak .
072 S/m Cl − leak conductance glial g Cl g, leak . KIR resting conductance glial g . Maximum NaKCl rate glial g NaKCl . × − A/m
Maximum pump rate glial ˆ I g . Membr. area-to-volume glial γ ge . × Membr. capacitance glial C ge . × − F/m
Membr. water permeability glial η ge . × − m /(mol s) Membrane stiffness neuron S ne . × Pa/m –Membrane stiffness glial S ge . × Pa/m –Gap junction connectivity glial χ g .
05 –Neuronal water permeability κ n /(N s) *Glial water permeability κ g . × − m /(N s) –ECS water permeability κ e . × − m /(N s) – Table 8. Physical parameters for the glial membrane and mechanical pressure. We use SI baseunits, that is, meter (m), mole (mol), Siemens (S) and ampere (A). The values are collected fromO’Connell et al. , Østby et al. , Steinberg et al. , and Yao et al. . The symbol – indicatesthat the value is chosen by the authors, as we could not find any relevant values in the literature.* The neuronal water permeability is set to zero as we assume no gap junctions connecting theneurons. inward rectifier current ( I KIR , A/m ). Following Steinberg et al. , the KIR con-ductance g KIR (S/m ) is modelled as: g KIR = g (cid:115) [K + ] R . . )1 + exp( φ ge − E g K +18 . . ) 1 + exp( − . − . . )1 + exp( − . φ ge . ) , (8.1)where g (S/m ) is the resting membrane conductance, and corresponds to theconductance when φ ge = E g K and [K + ] e = [K + ] e . The Na/K/ATPase pump (ATP)occurs in both the neuron and glial membrane, and is modelled as in (5.4). Finally,the current through the NaKCl cotransporter I NaKCl (A/m ) is modelled as: I NaKCl = g NaKCl ln (cid:32) [Na + ] g [K + ] g [Cl − ] g [Na + ] e [K + ] e [Cl − ] e (cid:33) . (8.2)In summary, the total currents over the glial membrane are modelled by (8.3)(with the currents (A/m ) are converted to ion fluxes (mol/(m s)) by dividing byebruary 5, 2021 13:25 WSPC/INSTRUCTION FILE ws-m3as Numerics for brain electrochemomechanics N (cid:107) [Na] e − [Na] eh (cid:107) L (cid:107) φ n − φ nh (cid:107) L (cid:107) α n − α nh (cid:107) L (cid:107) p e − p eh (cid:107) L n (cid:107) [Na] e − [Na] eh (cid:107) H (cid:107) φ n − φ nh (cid:107) H Table 9. Selected L -errors (upper panel) and H -errors (lower panel) and convergence rates (inparenthesis) for the full scheme at time t = 0 .
002 s. The test was run on the unit interval, andwe initially let ∆ t = 0 .
001 s, and then half the timestep in each series. The spatial discretizationconsists of N intervals. Faraday’s constant F times the valence z k ): J Na g = 1 F z Na (cid:0) I Na g, leak + 3 I g, ATP + I NaKCl + I Naex (cid:1) , (8.3a) J K g = 1 F z K (cid:0) I KIR − I g, ATP + I NaKCl + I Kex (cid:1) , (8.3b) J Cl g = 1 F z Cl (cid:0) I Cl g, leak + 2 I NaKCl + I Clex (cid:1) . (8.3c)Here, I Naex , I Kex , and I Clex are excitatory fluxes used to trigger a cortical spreadingdepression wave, see (5.5) in Section 5.1.1.
CSD wave characteristics
As in the zero flow limit scenario, excitatory flux stimulation leads to a CSD wavetraveling through the tissue: we observe neuronal depolarization, neuronal and ECSionic concentration changes, and neuronal swelling (Figure 5). Moreover, we observethat the glial potential depolarizes: from −
81 to −
31 mV, accompanied by a smalldrop in the extracellular potential from 0 to − .
5% and6 . A. J. Ellingsrud, N. Boull´e, P. E. Farrell, and M. E. Rognes compartments (neurons, ECS) in the zero flow limit (c.f. Figure 2), which is inaccordance with the (numerical) findings reported by O’Connell et al .The CSD wave is accompanied by a decrease in the mechanical pressures, fromthe baseline of 0 kPa down to − −
334 and −
389 kPa in the neuronal, glialand extracellular compartments, respectively (Figure 5 F). The pressure gradients(mechanical and osmotic) drive microscopic fluid flow within the glial and the ex-tracellular compartments. We observe flow rates of up to 1 . − . µ m/s in theglial and extracellular compartments, respectively; i.e. fluid flows in opposite direc-tions (Figure 6 B,C). During glial swelling, water moves across the glial membranefrom the ECS into the glial cells. In response, water within the glial cell networkwill be pushed away from this area. Indeed, we observe a positive flow rate to theright of the swelling, and a negative flow rate left of the swelling in the glial cellsand vice versa in the ECS. We observe no flow in the neuronal compartment, whichis expected as the neuronal water permeability is set to zero (Figure 6A).Next, we (quantitatively) evaluate the numerical convergence of the full modelscheme presented in Section 6. To this end, we consider the quantities of interestdefined in Section 5.2 in addition to the (spatial) width w CSD , p of the extracellularmechanical pressure wave at t = 50 s, given by X = { x | p R ( x, > p thres } w CSD , p = max X − min X with p thres = −
10 kPa.
Convergence of CSD wave characteristics during refinement
Drawing on the findings in Section 6, here, we consider the implicit lower-orderscheme based on Godunov splitting, a BE method for the PDE time-discretization,and ESDIRK4 for the ODE time-stepping – to compute the mean speed (cf. Sec-tion 5.2) and the (spatial) width of the extracellular mechanical pressure wave(cf. Section 8.2) for different mesh resolutions and time steps: ∆ x N = 10 /N mm for N = 1000 , , , , , t i = 12 . /i s for i =1 , , , , , ,
64. The results are presented in Figure 7. Wave speeds are convertedfrom the native m/s to mm/min for interpretability.As in the zero flow limit, the computed mean wave speed and the extracellularmechanical pressure wave width increases with decreasing ∆ t , and decreases withdecreasing ∆ x : the smaller the time step, the faster and wider the wave, whilethe smaller the mesh size, the slower and narrower the wave (Figure 7A–C). Thebehavior of the mean wave speed is qualitatively similar to that in the schemes forthe zero flow limit (cf. 1). The computed extracellular mechanical pressure wavewidth vary substantially, ranging from 2 .
409 to 4 .
67 mm.The differences ∆ w CSD , p in the extracellular mechanical pressure wave widthbetween the coarsest mesh sizes N = 1000 and 2000 are in the range 2 . − . w CSD , p between the coarsest time steps ∆ t = 12 . .
25 are in theebruary 5, 2021 13:25 WSPC/INSTRUCTION FILE ws-m3as
Numerics for brain electrochemomechanics t = 50 s of neuronal, glial and extracel-lular ion concentrations ( A , B , C ), electrical potentials ( D ), change in volume fractions ( E ) andmechanical pressure ( F ).Fig. 6. Full model simulation of a CSD wave: fluid velocities. Snapshots at t = 50 s of compart-mental fluid velocities in the neurons ( A ), the glial cells ( B ), and the ECS ( C ). range 2 . − .
910 mm. For the finest time and mesh resolutions, we observe that∆ w CSD , p = 0 .
005 and 0 .
344 mm, respectively; thus the spatial error continues todominate. Finally, we observe that ∆ w CSD , p decreases as we refine the discretiza-tions in time and space. There is however no clear convergence rate (Figure 7A, B).ebruary 5, 2021 13:25 WSPC/INSTRUCTION FILE ws-m3as A. J. Ellingsrud, N. Boull´e, P. E. Farrell, and M. E. Rognes A Temporal refinement Sp a t i a l r e fin e m e n t B N ∆ t w CSD , p w CSD , p C N ∆ t v CSD v CSD
Fig. 7. Full model simulation: CSD wave properties during refinement in space (N) and time(∆t, ms) in a 1D domain of length 10 mm at t = 50 s. A : Extracellular mechanical pressure p R ( x,
50) (kPa) versus x ∈ Ω (mm). B : Extracellular mechanical pressure wave width (mm) anddifference ∆ w CSD , p between consecutive refinements. C : CSD mean wave speed ¯ v CSD (mm/min)and difference ∆¯ v CSD between consecutive refinements. ebruary 5, 2021 13:25 WSPC/INSTRUCTION FILE ws-m3as
Numerics for brain electrochemomechanics
9. Discussion and concluding remarks
We have presented and analyzed finite element-based splitting schemes for a math-ematical framework modelling ionic electrodiffusion and water movement in biolog-ical tissue in general, and brain tissue in particular. We have evaluated the schemesin terms of their numerical properties, including accuracy, convergence, and compu-tational efficiency, for idealized scenarios as well as for challenging, physiologicallyrealistic problem settings.The schemes display optimal convergence rates in space for problems withsmooth manufactured solutions. However, the physiological CSD setting is chal-lenging: we find that the accurate computation of CSD wave characteristics (wavespeed and wave width) requires a very fine spatial and fine temporal resolutionfor all schemes tested. Indeed, different splitting and time stepping schemes andlower and higher order finite element schemes give comparable results in terms ofaccuracy. Overall, the error associated with the spatial discretization dominates.Explicit PDE and/or ODE time stepping schemes easily fail to converge even foronly moderately coarse time steps, but yield accurate results for very fine timesteps.In light of the long time scale associated with CSD (seconds to minutes), the smalltime steps imposed by the explicit schemes (less than a millisecond) represent asevere restriction.The mathematical framework studied here was presented by Mori in 2015 , andhas been used to simulate cortical spreading depression in a three-compartmentsetting (including neurons, glial cells and extracellular space) , and in multiplespatial dimensions . However, little has been reported on numerical properties ofdiscretizations of this model. The aforementioned studies , have used time stepsof the order 10 ms and mesh sizes of the order 0.156–0.02 mm (corresponding to N = 64 −
500 cf. Figure 1). Our findings indicate that high resolution is requiredto accurately compute CSD wave properties and that low-to-moderate resolutionscan substantially overestimate (or, but more rarely, underestimate) the CSD wavespeed. We expect our finite element findings to extend also to comparable finitedifference or finite volume discretizations.In terms of limitations, we have here compared different numerical schemes interms of accuracy, with less emphasis on computational complexity or cost. Weconsider these numerical investigatons as a starting point and guide for future the-oretical studies. Another research direction would be the extension of this study tothe two-dimensional CSD model studied by O’Connell and Tuttle , where thedevelopment of multigrid solvers seems crucial to reach the high spatial resolutionneeded to obtain accurate solutions.This paper focuses on numerical challenges related to approximating systemsfor ionic electrodiffusion and microscopic water movement. We remark that the fullmodel simulation, including extracellular mechanical pressure, yields pressure differ-ences far greater than what one might expect in this setting ( ∼ A. J. Ellingsrud, N. Boull´e, P. E. Farrell, and M. E. Rognes velocity model best represents the physiology, in particular the fluid velocity com-ponent driven by electrostatic forces. On the other hand, it is well-established thatlarge osmotic pressure gradients indeed are present in the brain environment.In conclusion, our findings show that numerical simulation of ionic electrodiffu-sion and water movement in brain tissue is feasible, but requires care numericallyand substantial computational resources. Numerical schemes or solution approachesthat retain accuracy at a lower computational expense would enable the study of awide array of phenomena in brain physiology, including in the context of patholog-ical conditions.
Appendix A. Supplementary Tables N ∆ t v CSD v CSD – 0.180 0.150 0.114 0.069 0.038 0.017
Table 10. CSD mean wave speed ¯ v CSD (mm/min) and difference in CSD mean wave speed ∆¯ v CSD between consecutive refinements in space (rows) and time (columns). Numerical scheme: Godunovsplitting, BDF2, ESDIRK4. ebruary 5, 2021 13:25 WSPC/INSTRUCTION FILE ws-m3as
Numerics for brain electrochemomechanics A N ∆ t v CSD ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ v CSD – – – – – 0.038 0.017 B N ∆ t v CSD v CSD – 0.195 0.158 0.132 0.081 0.046 0.022
Table 11. CSD mean wave speed ¯ v CSD (mm/min) and difference in CSD mean wave speed ∆¯ v CSD between consecutive refinements in space (rows) and time (columns). Numerical scheme: Strang,BDF2, and RK4 ( A ) or BE ( B ). ∗ indicates that the solver failed to converge. ebruary 5, 2021 13:25 WSPC/INSTRUCTION FILE ws-m3as A. J. Ellingsrud, N. Boull´e, P. E. Farrell, and M. E. Rognes A Parameter Symbol Value Unitneuron volume fraction α . + concentration neuron [Na] . K + concentration neuron [K]
132 mol/m Cl − concentration neuron [Cl] . Na + concentration ECS [Na]
137 mol/m K + concentration ECS [K] Cl − concentration ECS [Cl]
114 mol/m potential neuron φ − .
070 Vpotential ECS φ . B Parameter Symbol Value Unitneuron volume fraction α n . α g . + concentration neuron [Na] n . K + concentration neuron [K] n
132 mol/m Cl − concentration neuron [Cl] n . Na + concentration glial [Na] g
13 mol/m K + concentration glial [K] g
128 mol/m Cl − concentration glial [Cl] g . Na + concentration ECS [Na] e
137 mol/m K + concentration ECS [K] e Cl − concentration ECS [Cl] e
114 mol/m potential neuron φ n − .
070 Vpotential glial φ g − .
082 Vpotential ECS φ e . p e . Table 12. Initial values for state variables in the zero flow limit ( A ) and in the full model ( B ). Weuse SI base units; that is, meter (m), and mole (mol). ebruary 5, 2021 13:25 WSPC/INSTRUCTION FILE ws-m3as Numerics for brain electrochemomechanics Acknowledgments
The authors thank Didrik Bakke Dukefoss, Rune Enger, Geir Halnes, Erlend A.Nagelhus, and Klas Pettersen for valuable and constructive discussions on corticalspreading depression and brain electrophysiology.AJE and MER have received support the European Research Council (ERC)under the European Union’s Horizon 2020 research and innovation programme un-der grant agreement 714892. PEF acknowledges support from the Engineering andPhysical Sciences Research Council (grants EP/R029423/1 and EP/V001493/1).NB acknowledges support from the Engineering and Physical Sciences ResearchCouncil Centre for Doctoral Training in Industrially Focused Mathematical Mod-elling (grant EP/L015803/1) in collaboration with Simula Research Laboratory.
References
1. P. Aitken and G. Somjen, The sources of extracellular potassium accumulation in theCA1 region of hippocampal slices,
Brain Res. (1986) 163–167.2. F. Brocard, N. A. Shevtsova, M. Bouhadfane, S. Tazerart, U. Heinemann, I. A. Ry-bak and L. Vinay, Activity-dependent changes in extracellular Ca2+ and K+ revealpacemakers in the spinal locomotor-related network,
Neuron (2013) 1047–1054.3. A. C. Charles and S. M. Baca, Cortical spreading depression and migraine, Nat. Rev.Neurol. (2013) 637.4. J. R. Cressman, G. Ullah, J. Ziburkus, S. J. Schiff and E. Barreto, The influence ofsodium and potassium dynamics on excitability, seizures, and the stability of persistentstates: I. Single neuron dynamics, J. Comput. Neurosci. (2009) 159–170.5. F. Ding, J. O’Donnell, Q. Xu, N. Kang, N. Goldman and M. Nedergaard, Changesin the composition of brain interstitial ions control the sleep-wake cycle, Science (2016) 550–555.6. T. A. Driscoll, N. Hale and L. N. Trefethen,
Chebfun Guide (Pafnuty Publications,2014).7. R. Eisenberg, V. Barcilon and R. Mathias, Electrical properties of spherical syncytia,
Biophys. J. (1979) 151–180.8. R. S. Eisenberg and E. A. Johnson, Three-dimensional electrical field problems inphysiology, Prog. Biophys. Mol. Biol. (1970) 1–65.9. A. J. Ellingsrud and N. Boull´e, Supplementary code: Accurate numerical simulationof electrodiffu-sion and osmotic water movement in brain tissue, https://bitbucket.org/adaje/supplementary-code-accurate-numerical-simulation-of/src/master/ , 2020.10. A. J. Ellingsrud, A. Solbr˚a, G. T. Einevoll, G. Halnes and M. E. Rognes, Finiteelement simulation of ionic electrodiffusion in cellular geometries, Front. Neuroinform. (2020) 11.11. R. Enger, D. B. Dukefoss, W. Tang, K. H. Pettersen, D. M. Bjørnstad, P. J. Helm,V. Jensen, R. Sprengel, K. Vervaeke, O. P. Ottersen et al., Deletion of aquaporin-4curtails extracellular glutamate elevation in cortical spreading depression in awakemice, Cereb. Cortex (2017) 24–33.12. P. E. Farrell, J. E. Hake, S. W. Funke and M. E. Rognes, Automated adjoints ofcoupled PDE-ODE systems, SIAM J. Sci. Comput. (2019) C219–C244.13. M. S. Floater and K. Hormann, Barycentric rational interpolation with no poles andhigh rates of approximation, Numer. Math. (2007) 315–331. ebruary 5, 2021 13:25 WSPC/INSTRUCTION FILE ws-m3as A. J. Ellingsrud, N. Boull´e, P. E. Farrell, and M. E. Rognes
14. D. Gottlieb and S. A. Orszag,
Numerical analysis of spectral methods: theory andapplications (SIAM, 1977).15. G. Halnes, T. M¨aki-Marttunen, D. Keller, K. H. Pettersen, O. A. Andreassen andG. T. Einevoll, Effect of ionic diffusion on extracellular potentials in neural tissue,
PLoS Comput. Biol. (2016) e1005193.16. L. Hertz, Possible role of neuroglia: a potassium-mediated neuronal–neuroglial–neuronal impulse transmission system, Nature (1965) 1091–1094.17. B. Hille et al.,
Ion channels of excitable membranes , volume 507 (Sinauer Sunderland,MA, 2001).18. N. H¨ubel and M. A. Dahlem, Dynamics from seconds to hours in Hodgkin–Huxleymodel with time-dependent ion concentrations and buffer reservoirs,
PLoS Comput.Biol. .19. H. Kager, W. J. Wadman and G. G. Somjen, Simulated seizures and spreading de-pression in a neuron model incorporating interstitial space and ion concentrations, J.Neurophysiol. (2000) 495–512.20. R. K¨ohling and J. Wolfart, Potassium channels in epilepsy, Cold Spring Harb. Perspect.Med. (2016) a022871.21. S. Kuffler, J. Nicholls and R. Orkand, Physiological properties of glial cells in thecentral nervous system of amphibia., J. Neurophysiol. (1966) 768–787.22. H. P. Langtangen and S. Linge, Finite difference computing with PDEs: a modernsoftware approach (Springer Nature, 2017).23. A. Logg, K.-A. Mardal and G. Wells,
Automated solution of differential equations bythe finite element method: The FEniCS book , volume 84 (Springer Science & BusinessMedia, 2012).24. Y. Mori, A multidomain model for ionic electrodiffusion and osmosis with an appli-cation to cortical spreading depression,
Physica D (2015) 94–108.25. C. Nicholson, G. Ten Bruggencate, H. Stockle and R. Steinberg, Calcium and potas-sium changes in extracellular microenvironment of cat cerebellar cortex,
J. Neuro-physiol. (1978) 1026–1039.26. W. Noh, G. Choi, S. Pak, S. Yang and S. Yang, Transient potassium channels: Ther-apeutic targets for brain disorders, Front. Cell. Neurosci. (2019) 265.27. R. O’Connell, A computational study of cortical spreading depression, Ph.D. thesis,University of Minnesota, 2016.28. R. O’Connell and Y. Mori, Effects of glia in a triphasic continuum model of corticalspreading depression, Bull. Math. Biol. (2016) 1943–1967.29. I. Østby, L. Øyehaug, G. T. Einevoll, E. A. Nagelhus, E. Plahte, T. Zeuthen, C. M.Lloyd, O. P. Ottersen and S. W. Omholt, Astrocytic mechanisms explaining neural-activity-induced shrinkage of extraneuronal space, PLoS Comput. Biol. (2009)e1000272.30. L. Øyehaug, I. Østby, C. M. Lloyd, S. W. Omholt and G. T. Einevoll, Dependence ofspontaneous neuronal firing and depolarisation block on astroglial membrane trans-port mechanisms, J. Comput. Neurosci. (2012) 147–165.31. D. Pietrobon and M. A. Moskowitz, Chaos and commotion in the wake of corticalspreading depression and spreading depolarizations, Nat. Rev. Neurosci. (2014)379–393.32. J. Pods, J. Sch¨onke and P. Bastian, Electrodiffusion models of neurons and extracel-lular space using the Poisson-Nernst-Planck equations — numerical simulation of theintra-and extracellular potential for an axon model, Biophys. J. (2013) 242–254.33. P. J. Roache,
Verification and validation in computational science and engineering (Hermosa, 1998). ebruary 5, 2021 13:25 WSPC/INSTRUCTION FILE ws-m3as
Numerics for brain electrochemomechanics
34. M. J. Sætra, G. T. Einevoll and G. Halnes, An electrodiffusive, ion conserving Pinsky-Rinzel model with homeostatic mechanisms,
PLoS Comput. Biol. (2020) e1007661.35. G. Somjen, H. Kager and W. Wadman, Computer simulations of neuron-glia interac-tions mediated by ion flux, J. Comput. Neurosci. (2008) 349–365.36. R. Srivastava, M. Aslam, S. R. Kalluri, L. Schirmer, D. Buck, B. Tackenberg, V. Roth-hammer, A. Chan, R. Gold, A. Berthele et al., Potassium channel KIR4. 1 as animmune target in multiple sclerosis, N. Eng. J. Med. (2012) 115–123.37. C. Staehr, R. Rajanathan and V. V. Matchkov, Involvement of the Na+, K+-ATPaseisoforms in control of cerebral perfusion,
Exp. Physiol. (2019) 1023–1028.38. B. Steinberg, Y. Wang, H. Huang and R. Miura, Spatial buffering mechanism: Math-ematical model and computer simulations,
Math. Biosci. Eng (2005) 675–702.39. D. Sterratt, B. Graham, A. Gillies and D. Willshaw, Principles of computationalmodelling in neuroscience (Cambridge University Press, 2011).40. J. Sundnes, G. T. Lines, X. Cai, B. F. Nielsen, K.-A. Mardal and A. Tveito, Solvingsystems of ODEs, in
Computing the Electrical Activity in the Heart (Springer, 2006),pp. 149–173.41. X. Tong, Y. Ao, G. C. Faas, S. E. Nwaobi, J. Xu, M. D. Haustein, M. A. Anderson,I. Mody, M. L. Olsen, M. V. Sofroniew et al., Astrocyte Kir4. 1 ion channel deficitscontribute to neuronal dysfunction in Huntington’s disease model mice,
Nat. Neurosci. (2014) 694.42. L. N. Trefethen, Approximation Theory and Approximation Practice , volume 164 of
Other Titles in Applied Mathematics (SIAM, 2019).43. L. Tung, A bi-domain model for describing ischemic myocardial D-C potentials., Ph.D.thesis, Massachusetts Institute of Technology, 1978.44. A. Tuttle, Modeling Regional Variation of Cortical Spreading Depression: A Compu-tational Study, Ph.D. thesis, University of Minnesota, 2019.45. A. Tuttle, J. R. Diaz and Y. Mori, A computational study on the role of glutamate andNMDA receptors on cortical spreading depression using a multidomain electrodiffusionmodel,
PLoS Comput. Biol. .46. G. Ullah, J. R. Cressman Jr, E. Barreto and S. J. Schiff, The influence of sodium andpotassium dynamics on excitability, seizures, and the stability of persistent states: II.Network and glial dynamics, J. Comput. Neurosci. (2009) 171–183.47. D. Utzschneider, J. Kocsis and M. Devor, Mutual excitation among dorsal root gan-glion neurons in the rat, Neurosci. Lett. (1992) 53–56.48. W. Yao, H. Huang and R. M. Miura, A continuum neuronal model for the instigationand propagation of cortical spreading depression,
Bull. Math. Biol.73