INF-5610, Matematiske modeller i medisin Eksamen Forelesere: Det - - PowerPoint PPT Presentation

INF-5610, Matematiske modeller i medisin

Forelesere: Glenn Terje Lines (glennli@ifi.uio.no) Joakim Sundnes (sundnes@ifi.uio.no) Temaer: Hjertecellenes egenskaper Strømflyt i hjertet og i kroppen, EKG Modeller for disse fenomener Numeriske metoder for disse modellene To obligatorisk oppgaver

– p. 1

Eksamen

Det blir oppgitt seks temaer, senest to uker før eksamen Det skal forberedes 20 min foredrag for hvert tema På eksamen trekkes et av disse temaene og holde foredraget Det vil også bli stilt spørsmål fra andre temaer

– p. 2

Forelesningsplan, del I

Anatomi, om celler og hjertet. Campell, KS Grunnleggende biophys proseeser. KS Kap. 2

• Ionekanaler. KS Kap. 3

Eksiterbarhet og progagering av signal. KS Kap. 4 Nevroner og cell-koplinger KS Kap. 7&8 Calcium dynamikk KS Kap. 5

– p. 3

Forelesningsplan, del II

Bidomene modell. KS Kap. 11. Presentasjon av oblig. Presentasjon av en anvendelse Modeller for mekanikk. Numeriske metoder: fra modell til software, splitting. Numeriske metoder: prinsipp bak ODE og PDE løsere.

– p. 4

Levels of modeling

Body Organ Tissue Cell Organelles Proteins

– p. 5

The Cell Membrane

Consist of a bilipid layer Embedded proteins for transport control Selectively permeable Maintains concentration gradients Has a transmembrane potential

– p. 6

The Cell Membrane

– p. 7

The Cell Membrane

– p. 8

Two types of transmembrane flow

Passive: Diffusion along the concentration gradient Through the membrane (H2O, O2, CO2) Through specialized channels (Na+, K+, Cl−) Carrier mediated transport Active: Energy driven flow against the gradients ATP driven pumps (Na+ − K+, Ca2+) Exchangers driven by concentration gradients (Na+ − Ca2+)

– p. 9

Transmembrane flow

– p. 10

Active Transport

– p. 11

Mathematical models of chemical reactions

– p. 12

The Law of Mass Action

Chemical A and B react to produce chemical C: A + B

k

− → C The rate constant k determines the rate of the reaction. It can be interpreted as the probability that a collision between the reactants produces the end results. If we model the probability of a collision with the product [A] [B] we get the law of mass action: d[C] dt = k[A][B]

– p. 13

A two way reaction

The reverse reaction may also take place: A + B

k+

− → ← −

k−

C The production rate is then: d[C] dt = k+[A][B] − k−[C] At equilibrium when d[C]/dt = 0 we have: k−[C] = k+[A][B] (1)

– p. 14

If A + B

k

− → C is the only reaction involving A and C then d[A]/dt = −d[C]/dt so that [A] + [C] = A0 (2) Substituting (2) into (1) yields: [C] = A0 [B] Keq + [B] where Keq = k−/k+. Notice that [B] = Keq = ⇒ [C] = A0/2 and [B] → ∞ = ⇒ [C] → A0

– p. 15

Enzyme Kinetics

Characteristics of enzymes: Made of proteins Acts as catalysts for biochemical reactions Speeds up reactions by a factor > 107 Highly specific Often part of a complex regulation system

– p. 16

Reaction model of enzymatic reaction

S + E

k1

− → ← −

k−1

C

k2

− → P + E with S: Substrate E: Enzyme C: Complex P: Product

– p. 17

Mathematical model of enzymatic reaction

Applying the law of mass action to each compound yields: d[S] dt = k−1[C] − k1[S][E] + JS d[E] dt = (k−1 + k2)[C] − k1[S][E] d[C] dt = k1[S][E] − (k2 + k−1)[C] d[P] dt = k2[C] − JP Here we also supply the substrate at rate JS and the product is removed at rate JP.

– p. 18

Equilibrium

Note that In equilibrium d[S]/dt = d[E]/dt = d[C]/dt = d[P]/dt = 0 it follows that that JS = JP. Production rate: J = JP = k2[C]

– p. 19

In equilibrium we have d[E] dt = 0 that is (k−1 + k2)[C] = k1[S][E] Since the amount of enzyme is constant we have [E] = E0 − [C] This yields [C] = E0[S] Km + [S] with Km = k−1+k2

k1

and E0 is the total enzyme concentration. Production rate: d[P]

dt = k2[C] = Vmax [S] Km+[S], where Vmax = k2E0.

– p. 20

Cooperativity

S + E

k1

− → ← −

k−1

C1

k2

− → E + P S + C1

k3

− → ← −

k−3

C2

k4

− → C1 + P with S: Substrate E: Enzyme C1: Complex with one S C1: Complex with two S P: Product

– p. 21

Mathematical model of cooperativ reaction

Applying the law of mass action to each compound yields: ds dt = −k1se + k−1c1 − k3sc1 + k−3c2 dc1 dt = k1se − (k−1 + k2)c1 − k3sc1 + (k4 + k−3)c2 dc2 dt = k3sc1 − (k4 + k−3)c2

– p. 22

Equilibrium

Set dc1

dt = dc2 dt = 0, and use e0 = e + c1 + c2,

c1 = K2e0s K1K2 + K2s + s2 c2 = e0s2 K1K2 + K2s + s2 where K1 = k−1+k2

k1

, K2 = k4+k−3

k3

Reaction speed: V = k2c1 + k4c2 = (k2K2 + k4s)e0s K1K2 + K2s + s2

– p. 23

Case 1: No cooperation

The binding sites operate independently, with the same rates k+ and k−. k1 and k3 are associated with events that can happen in two ways, thus: k1 = 2k3 = k+ k−3 = 2k−1 = k− Which gives this reaction speed: V = 2 k2e0s K + s, K = k− + k2 k+ Note that this is the same as the reaction speed for twice the amount of an enzyme with a single binding site.

– p. 24

Case 2: Strong cooperation

The first binding is unlikely, but the next is highly likely, i.e. k1 is small, and k3 is large. We go to the limit: k1 → 0, k3 → ∞, k1k3 = const so K2 → 0, K1 → ∞, K1K2 = const In this case the rection speed becomes: V = k4e0s2 K2

m + s2 = Vmax

s2 K2

m + s2

with K2

m = K1K2, and Vmax = k4e0

– p. 25

The Hill equation

In general with n binding sites, the reaction rate in the limit will be: V = Vmax sn Kn

m + sn

This model is often used when the intermediate steps are unknown, but cooperativity suspected. The parameters Vmax, Km and n are usually determined experimentally.

– p. 26

Carrier-Mediated Transport

Some substances can not pass the membrane on their own, but are helped by a carrier protein. Types of transport: Uniport: Transport of single substance Symport: Transport of several substances in same direction Antiport: Transport of several substances in opposite directions With symport and antiport the carrier molecule as several binding sites.

– p. 27

Uniport

Substrate S combines with a carrier protein C to form a complex P . The protein has two conformal states. Model: Si + Ci

k+

− → ← −

k−

Pi

k

− → ← −

k

Pe

k−

− → ← −

k+

Se + Ce Ci

k

− → ← −

k

Ce

– p. 28

Model for Carrier Mediated Transport, Uniport

Applying the law of mass action: d[Si] dt = k−[Pi] − k+[Si][Ci] − J d[Se] dt = k−[Pe] − k+[Se][Ce] + J d[Pi] dt = k[Pe] − k[Pi] + k+[Si][Ci] − k−[Pi] d[Pe] dt = k[Pi] − k[Pe] + k+[Se][Ce] − k−[Pe] d[Ci] dt = k[Ce] − k[Ci] + k−[Pi] − k+[Si][Ci] d[Ce] dt = k[Ci] − k[Ce] + k−[Pe] − k+[Se][Ce] Here J is the influx of the glucose molecules (S).

– p. 29

Size of flux in equilibrium

The flow in equilibrium can be setting the derivatives to zero and solve for J. This yields a system of six eq. and seven unknowns. The amount of protein is conserved so we have: [Ci] + [Ce] + [Pi] + [Pe] = C0 Solving for J in equilibrium then gives: J = 1 2kKC0 [Se] − [Si] ([Si] + K + Kd)([Se] + K + Kd) − K2

d

with K = k−/k+ and Kd = k/k+.

– p. 30

Size of flux in equilibrium

J = 1 2kKC0 [Se] − [Si] ([Si] + K + Kd)([Se] + K + Kd) − K2

d

Factors affecting the flux: The amount of Carrier molecules C0 The rate constants Substrate gradient

– p. 31

Model for symport

Two different substances S and T are transported in the same

• direction. The carrier C has m binding sites for S and n for T:

mSi + nTi + Ci

k+

− → ← −

k−

Pi

kp

− → ← −

k−p

Pe

k−

− → ← −

k+

mSe + nTe + Ce Ci

k

− → ← −

k

Ce

– p. 32

Need to model mathematically the process mS + nT + C

k+

− → ← −

k−

P Consider the simpler reaction A + B + C

k+

− → ← −

k−

ABC If we assume that the reaction takes place in two steps A + B

k1

− → ← −

k−1

AB AB + C

k+

− → ← −

k−

ABC

– p. 33

cont. A + B

k1

− → ← −

k−1

AB AB + C

k+

− → ← −

k−

ABC If the intermediate step is fast, we can assume it to be in equilibrium: d[AB] dt = k1[A][B] − k−1[AB] = 0 ⇒ [AB] = k1/k−1[A][B] For the total reaction: d[ABC] dt = k+[AB][C] − k−[ABC] = k+ k1 k−1 [A][B][C] − k−[ABC]

– p. 34

Flux for symport

With repeated use of similar arguments d[P] dt = k+[S]m[T]n[C] − k−[P] The symport model will be identical to the uniport model by substituting [S] with [S]m[T]n. Flux: J = 1 2KdKk+C0 [Se]m[Te]n − [Si]m[Ti]n ([Si]m[Ti]n + K + Kd)([Se]m[Te]n + K + Kd) − K2

d

– p. 35

Antiport

In antiport the two substances travels in opposite direction (exchangers). Model: mSi + nTe + Ci

k+

− → ← −

k−

Pi

kp

− → ← −

k−p

Pe

k−

− → ← −

k+

mSe + nTi + Ce Mathematically almost the same flux, but with subscript of T toggled: J = 1 2KdKk+C0 [Se]m[Ti]n − [Si]m[Te]n ([Si]m[Te]n + K + Kd)([Se]m[Ti]n + K + Kd) − K2

d

– p. 36

Active Transport

The Na-K pump. Sodium is pumped out, potassium is pumped in, both against their concentration gradients. Fueled by ATP .

– p. 37

The Sodium-Potassium Pump

Model: Na+

i + C k1

− → ← −

k−1

NaC NaC + ATP

kp

− → ← −

k−p

NaCP + ADP NaCP

k2

− → ← −

k−2

Na+

e + CP

CP + K+

e k3

− → ← −

k−3

KCP

k4

− → ← −

k−4

P + KC KC

k5

− → ← −

k−5

K+

i + C

– p. 38

Flux for the Sodium-Potassium Pump

Solving for J in equilibrium: J = C0 [Na+

i ][Ke]K1K2 − [Na+ e ][K+ i ]K−1K−2[P]

([K+

e ]K2 + [K+ i ]K−2)Kn + ([Na+ i ]K1 + [Na+ e ]K−1)Kk

With K1 = k1k2kp, K−1 = k−1k−2k−p,K2 = k3k4k5, K−2 = k−3k−4k−5,Kn = k−1k−p + k2k−1 + k2kp and Kk = k−3k−4[P] + k−3k5 + k4k5 Assuming no ATP synthesis (k−p = 0): J = C0 [Na+

i ][Ke]K1K2

([K+

e ]K2 + [K+ i ]K−2)Kn + [Na+ i ]K1Kk

– p. 39

Nernst Equilibrium Potential

– p. 40

Diffusion

The conservation law for a compound with concentration c: rate change of c = local production + accumulation due to transport. Model: d dt

• Ω

c dV =

• Ω

p dV −

• ∂Ω

J · n dA Here p represents the production and J is the flux of c. The divergence theorem:

• ∂Ω

J · n dA =

• Ω

∇ · J dV The law is valid for every volume, thus: ∂c ∂t = p − ∇ · J Models for p and J are needed to compute c.

– p. 41

Fick’s Law

J = −D∇c The diffusion coefficient D depends upon the solute and the temperature of the embedding fluid: D = kT f T is the temperature measured on Kelvin, f is a frictional constant and k is the Boltzmann’s constant. The conservation law with this assumption is a reaction-diffusion equation: ∂c ∂t = ∇ · (D∇c) + p

– p. 42

1D Diffusion through a pore in the membrane

∂c ∂t = D ∂2c ∂2x Fixed intra and extra cellular concentration: c(0, t) = [C]i c(L, t) = [C]e At steady state: ∂c ∂t = 0 = ⇒ D ∂2c ∂2x = 0 = ⇒ ∂c ∂x = a = ⇒ c(x) = ax + b Taking the boundary condition into consideration yields: c(x) = [C]i + ([C]e − [C]i) x L and a constant flux: J = −D ∂c

∂x = D L ([C]i − [C]e)

– p. 43

Flow through a semi-preamble membrane

Consider two solutions: A: Contains 100mM Cl− ions and 100mM Na+ ions B: Contains 10mM Cl− ions and 10mM Na+ ions Both are neutral. If they are only separated by a membrane permeable to Cl− but not Na+, this will happen: Cl− will diffuse from A to B due the concentration gradient [Cl−]A will drop and [Cl−]B will increase [Na+]A and [Na+]B will remain fixed (no flow) A and B will no longer be neutral [Na+]A > [Cl−]A ⇒ A > 0, [Cl−]B > [Na+]B ⇒ B < 0. Cl− will be attracted to A and repelled from B

– p. 44

The Nernst Equilibrium Potential

We now have two forces driving Cl− across the membrane: Flow from A to B due to the concentration gradient Flow from B to A due to the charge gradient At some point an equilibrium is reached were the net flow is zero. The transmembrane potential at that point is called the Nernst Equilibrium Potential. An expression for this potential will now be derived

– p. 45

Plank’s equation

Models the ion-flux caused by an electrical field: J = −m z |z|c∇φ with m - mobility of the ions in the liquid z/|z| - sign of the charge of the ion c - the concentration of the ion ∇φ - the electrical field

– p. 46

Fick’s law: J = −D∇c Relationship between m and D: m = D|z|F RT here R is the gas constant and F is Faraday’s constant. Combined effect of concentration gradient and an electric field: J = −D(∇c + zF RT c∇φ)

– p. 47

Nernst Equilibrium Potential

Consider equilibrium in 1D flow: dc dx + zF RT cdφ dx = 0 1 c dc dx + zF RT dφ dx = 0 Integrating from inside (x=0) to outside (x=L) yields: ln(c)|c(L)

c(0) = − zF

RT (φ(L) − φ(0)) We define the transmembrane potential to be v = φi − φe The value of the transmembrane potential at zero flux is then veq = RT zF ln(ce ci ) (3) Equation (3) is referred to as the Nernst Equilibrium Potential.

– p. 48

Membrane currents

– p. 49

Accumulation around the membrane

The membrane has properties similar to a capacitor Consists of two conducting medias (intra- and extra cellular space) These are separated by an insulating material (the membrane) The potential over a capacitor is proportional to the separated charge (q): v = q/c The factor c is called the capacitance of the capacitor.

– p. 50

The cell membrane modeled as a leaky capacitor

As any real capacitor the membrane conducts some current. The flux of ions (Iion) will cause a change in q and thus v. Consider the change over a time interval ∆t. It follows that

∆v ∆t = 1 c ∆q ∆t and in the limit we get:

dv dt = 1 c dq dt The term dq

dt is called the capacitive current and is denoted Ic.

– p. 51

Electrical circuit model of the cell membrane

Intracellular space Extracellular space

Iion Ic

The membrane behaves like resistor and capacitor in parallel: It = Iion + Ic If the loop is closed then It = 0. In that case all the ions passing the membrane accumulate and change the membrane potential accordingly.

– p. 52

Ionic currents

For passive ionic channels the flow through it must obey the equilibrium potential, i.e. be zero when v = veq. An number of models exists, two common are: Linear: I(v) = g(v − veq) Here g is the conductance of the channel. Goldman-Hodgkin-Katz: I(v) = gvci − cee

−zvF RT

1 − e

−zvF RT

Derived from Nernst-Planck equation with assuming a constant (non-zero) field strength.

– p. 53

Nernst-Plack equation: J = −D(∇c + zF RT c∇φ) Consider 1D flow through a channel and assume ∇φ is constant in space and that c and φ are in steady-state. dφ dx = ∆φ ∆x = φ(L) − φ(0) L − 0 = φe − φi L = −v/L The equation is reduced to an ODE: J/D = − dc dx − zF RT c(−v/L) = − dc dx + kc where k = zFv

RTL

– p. 54

The equation J/D = − dc dx + kc is solved by e−kxc = ci + J D 1 k(e−kx − 1) We determine J by using c(L) = ce: J = Dkci − c(L)e−kL 1 − e−kL = D zFv RTL ci − cee

−zvF RT

1 − e

−zvF RT

J has dimension moles per area per time, an expression for current is given by I = zFJ = D L z2F 2 RT vci − cee

−zvF RT

1 − e

−zvF RT

This is the Goldman-Hodgkin-Katz current equation.

– p. 55

Channel gating

– p. 56

Channel gating

The conductance of a channel varies with time and with transmembrane potential. Model for current per membrane area: I(V, t) = g(V, t)φ(V ) (4) Current through a single open channel is φ(V ) and the amount of

• pen channels per membrane area is g(V, t).

– p. 57

Two State K+-channel

Assumes that the channel can exist in two states, closed(C) and

• pen(O):

C

α(v)

− → ← −

β(v)

O Applying law of mass action: d[0] dt = α(v)[C] − β(v)[O] Dividing by the total amount of channels ([C]+[O]) yields dg dt = α(v)(1 − g) − β(v)g where g is the portion of open channel ( [O]/([C]+[O])).

– p. 58

Multiple sub-units

For some channels it is more appropriate to model the gate as series of several sub-gates. Example with two gates:

S00

α

→ ←

β

S10

α ↓↑ β α ↓↑ β

S01

α

→ ←

β

S11

Using law of mass action we get a system of four equation. Will try to reduce this number to one!

– p. 59

Combine the states S01 and S10 into S1 = S01 + S10 :

S01 dt

= αS00 + βS11 − (α + β)S01 +

S10 dt

= αS00 + βS11 − (α + β)S10 =

S1 dt

= 2αS00 + 2βS11 − (α + β)S1 Define S0 = S00 and S2 = S11, we can then write:

S0

2α

→ ←

β S1 α

→ ←

2β S2

– p. 60

Only two independent variables since S0 + S1 + S2 = ST,

• constant. Define xi = Si/ST. Claim:

x2 = n2, with dn dt = α(1 − n) − βn

– p. 61

Sodium

Behavior of the Sodium conductance can not be described by a chain of identical gates. Two subunits of type m and one of type h.

S00

2α

→ ←

β

S10

α

→ ←

2β

S20

γ ↓↑ δ γ ↓↑ δ γ ↓↑ δ

S01

2α

→ ←

β

S11

α

→ ←

2β

S21

Here Sij denotes i open m gates and j open h gates. Arguments similar to the one used above leads to these equations for m and h: dm dt = α(1 − m) − βm, dh dt = γ(1 − h) − δh

– p. 62

Excitability

– p. 63

Excitable Cells

Unlike other cells, excitable cells can be triggered to set off an action potential. During the action potential the transmembrane potential departs from its resting potential, reaches a peak potential and returns to the resting potential after some time. Nerve cells and cardiac cells uses the action potential as a signal to neighboring cells. The trigger must be of a certain size, if below the threshold the cell will not “fire”. As long as the trigger is above the threshold, the shape of the

– p. 64

The Hodgkin-Huxley Model

Developed to study the action potential of the squid nerve cells. Assumed three different current INa, IK and IL Assumed also linear current-voltage relationship: −Cm dv dt = Iion = gNa(v − vNa) + gK(v − vK) + gL(v − vL)

– p. 65

Can collect the current terms due to linearity: Cm dv dt = −geff(v − veq) where geff = gNa + gK + gL and veq = gNavNa + gKvK + gLvL geff veq is a weighted average of the individual equilibrium potentials. The weighing factors are time and voltage dependent.

– p. 66

A steady applied current Iapp moves the membrane potential to different equilibrium. Cm dv dt = −geff(v − veq) + Iapp = 0 Implies v = veq + 1 Cm geff Iapp The applied current will be compensated by an ionic current going the opposite way, thus the net current will be zero. For a sufficiently large Iapp, v will pass the threshold potential and an action potential is triggered. The conductivities will vary greatly.

– p. 67

Voltage Clamp measurements

The transmembrane potential is forced by an applied current to a fixed value. Since Iion = −Iapp for a fixed v, we can measure Iion as a function of time for a given level of v. Since v is fixed the observed variations must be due to temporal variation in the conductivities.

– p. 68

Total membrane current for different steps

– p. 69

From measurements to models

Initially, Hodgkin and Huxley assumed Iion = INa + IK. Two kind

• f experiments conducted:

1: Normal concentrations 2: [Na]e replaced by cohline ⇒ affects INa but not IK. Assumed further: Initially IK = 0 I1

Na/I2 Na = C, constant

I1

K = I2 K

Once I1

ion and I2 ion is recorded we can determine C from the first

and the second assumptions.

– p. 70

Expressions for the currents in terms of measurable quantities can now be obtained: I1

Na =

C C − 1(I1

ion − I2 ion)

IK = 1 1 − C (I1

ion − KI2 ion)

Assuming linear current-voltage relationships we get expressions for the conductivities: gNa = INa V − VNa , gK = IK V − VK For each pair of voltage clamp experiment (with a given voltage step), we now have a time course for gNa and gK.

– p. 71

Potassium and Sodium conductance

– p. 72

Model for the Potassium conductance

Assumed

dgK dt = f(v, t).

Ended up with introducing a second variable: gK = gKn4, with dn dt = α(v)(1 − n) − β(v)n and g is the maximum conductance. Forth power was chosen to get the correct shape of the solution.

– p. 73

The solution of τn dn dt = n∞ − n with constant coefficients is n(t) = n∞ + (n(0) − n∞)e−t/τn If we assume that n∞(0) = 0 a step from from 0 to v yields: n(t) = n∞(v) + (n∞(0) − n∞(v))e−t/τn(v) = n∞(v)(1 − e−t/τn(v)) A step in the other direction gives: n(t) = n∞(0) + (n∞(v) − n∞(0))e−t/τn(v) = n∞(v)e−t/τn(v)

– p. 74

Gating variable raised to different powers

2 4 6 8 10 0.2 0.4 0.6 0.8 1 n1 2 4 6 8 10 0.2 0.4 0.6 0.8 1 n2 2 4 6 8 10 0.2 0.4 0.6 0.8 1 n3 2 4 6 8 10 0.2 0.4 0.6 0.8 1 n4

– p. 75

Sodium conductance model

H&H realized that two different sub units were at work. Ended up with dgNa dt = gNam3h Values for mτ, m∞, hτ and h∞ obtained by fitting the solution to plots of gNa.

– p. 76

The Hodgkin-Huxley model

Introduces a third current, not time dependent: Cm dv dt = −gKn4(v − vK) − gNam3h(v − vNa) − gL(v − vL) with dg dt = αg(v)(1 − g) − βg(v)g, g = m, h, n Model based on voltage clamp measurement. How will it behave under normal conditions? The model will predict the action potential.

– p. 77

Analysis of the Hodgkin-Huxley model

– p. 78

Qualitative analysis

Would like to reduce the number of state variables to simplify analysis. One way is to treat the slowest variables as constants. Of the three gating variables m has the fastest dynamics. (Controls the activation of the Na-current). Reduced model: Cm dv dt = −gKn4

0(v − vK) − gNam3h0(v − vNa) − gL(v − vL)

– p. 79

Equilibria in the reduced HH-model

The nullclines dv

dt = 0 and dm dt = 0 form curves in the (v, m)-plane.

Their intersections are the equilibria. Initially three steady states vr, vs and ve. vr and ve are stable and vs unstable. As n0 and h0 changes, the dv

dt = 0 line will shift. ve will decrease,

coincide with vs and disappear. vr will become the only stable equilibrium.

– p. 80

Phase plot for the fast sub-system

– p. 81

Alternative reduction: m is very fast, almost in equilibrium: m = m∞(v) h + n almost constant: h = 0.8 − n We then have Cm dv dt = −gKn4(v −vK)−gNam3

∞(v) h

• (0.8 − n)(v −vNa)−gL(v −vL)

Equilibria found by looking at the crossing of the nullclines dv

dt = 0

and dn

dt = 0 in the (v, n)-plane.

– p. 82

Phase plot for the fast-slow reduced system

– p. 83

Properties of the phase plot

dv dt = 0 cubic form, with two stable and one unstable branch dn dt = 0 sigmoid form

One crossing with default parameters Trajectories horizontal due faster dynamic of v Starting points to the left of the unstable branch converges to equilibrium without crossing the unstable branch Starting points to the right of the unstable branch crosses this branch, reaches the rightmost branch, follows this branch and the trajectory continues to rise until dn

dt = 0 is

• crossed. The trajectory finally reaches the leftmost branch

and follows it to the equilibrium points.

– p. 84

Simulations with different initial conditions

– p. 85

Modified model

The point (0,0) is no longer a stable equilibrium.

– p. 86

Other models of the action potential

– p. 87

The FitzHugh-Nagumo model

Purpose: Keep the qualitative behavior of the Hodgkin-Huxley system, but in a simplified form. Derivation based on a an electrical circuit model. On dimensionless form: ǫ dv

dt

= f(v) − w − w0

dw dt

= v − γw − v0 The variable w is called the recovery variable.

– p. 88

Typically f is chosen to be “cubic”, i.e. with three zeros, f(0) = f(α) = f(1) and 0 < α < 1. Some choices: f(v) = Av(v − α)(1 − v) f(v) =

• −v,

v < α 1 − v, v > α f(v) =      −v, v < α/2 v − α, α/2 < v < (1 + α)/2 1 − v, v > (1 + α)/2

– p. 89

Cardiac cells

Excitable like neurons, display great variability SA node cells: Pace maker cells, controls the heart rate, self depolarizing AV node cells: Transmit signal from atria to ventricles with a delay Purkinje cells: Very high conductivity, propagate signal from AV out to the ventricles. Myocardial cells: Muscle cells (can contract) These cells have different action potentials. The HH-model was based on neurons. Other models necessary for cardiac cells.

– p. 90

The Beeler-Reuter model

A model for ventricular cells, includes three currents, six gates and one ionic concentration. −Cm dV dt = INa(V, m, h, j) + IK(V, x) + IS(V, f, g, [Ca]i) Here m, h, j, x, f, g are gating variables and [Ca]i is the intra cellular Calcium concentration The action potential is much longer then for HH. Early repolarization (notch).

– p. 91

Action potential produced by the Beeler-Reuter

– p. 92

Currents of the Beeler-Reuter model

Sodium current: Third gating variable included to model the slow recovery (long refractory period). The model also include an ungated “leakage” current: INa = (4m3hj + 0.003)(V − 50) Potassium: One singled gated (with x) and one ungated component: IK = f(v) + xg(v)

– p. 93

Calcium: To gates, d activates, f inactivates: IS = 0.09fg(V − VCa) In addition the [Ca]i is updated: dc dt = 0.07(1 − c) − IS where c = 107[Ca]i

– p. 94

The cable equation A.K.A. the monodomain model

– p. 95

Neurons

– p. 96

Electric flow in neurons

The neuron consists of three parts: Dendrite-tree, the “input stage” of the neuron, converges on the soma. Soma, the cell body, contain the “normal” cellular functions Axon, the output of the neuron, may be branched. Synapses at the ends are connected to neighboring dendrites. The axon has an excitable membrane, gives rise to active

• conduction. Will first look at conduction in the dentrites, passive

conduction.

– p. 97

The cable equation

The cell typically has a potential gradient along its length. Radial components will be ignored. Notation: Vi and Ve are intra- and extra cellular potential Ii and Ie are intra- and extra cellular (axial) current ri and re are intra- and extra cellular resistance per unit length

– p. 98

Discrete cable

– p. 99

Ohmic resistance assumed: Vi(x + ∆x) − Vi(x) = −Ii(x)ri∆x Ve(x + ∆x) − Ve(x) = −Ie(x)re∆x In the limit: Ii = − 1 ri ∂Vi ∂x and Ie = − 1 re ∂Ve ∂x

– p. 100

Conservation of current yields: Ii(x) − Ii(x + ∆x) = −(Ie(x) − Ie(x + ∆x)) = It∆x (5) where It is transmembrane current, per unit length. In the limit (5) yields: It = −∂Ii ∂x = ∂Ie ∂x We would like to express It in terms of V . 1 re ∂2Ve ∂x2 = − 1 ri ∂2Vi ∂x2 = − 1 ri (∂2V ∂x2 + ∂2Ve ∂x2 ) ( 1 re + 1 ri )∂2Ve ∂x2 = − 1 ri ∂2V ∂x2

– p. 101

cont. ( 1 re + 1 ri )∂2Ve ∂x2 = − 1 ri ∂2V ∂x2 ∂2Ve ∂x2 = −

1 ri 1 re + 1 ri

∂2V ∂x2 = − re re + ri ∂2V ∂x2 so It = ∂Ie ∂x = − 1 re ∂2Ve ∂x2 = 1 re + ri ∂2V ∂x2

– p. 102

From the membrane model previously derived we have It = p(Cm ∂V ∂t + Iion) where p is the circumference of the cable. The total expression will be in Ampere/meter. The total 1D cable model is then: p(Cm ∂V ∂t + Iion(V )) = ( 1 re + ri ∂2V ∂x2 )

– p. 103

Dimensionless form

We can scale the variables to reduce the number of parameters. Defines a membrane resistance: 1 Rm = ∆Iion ∆V (V0) where V0 is the resting potential. Multiplication with Rm CmRm ∂V ∂t + RmIion = Rm p(ri + re) ∂2V ∂x2 Here we have assumed ri and re constant. Defining f = −RmIion, τm = CmRm (time constant) and λ2

m = Rm/(p(ri + re)) (space constant squared) we can write

τm ∂V ∂t − f = λ2

m

∂2V ∂x2 (6)

– p. 104

Introduces the dimensionless variables: T = t/τm and X = x/λm We can then write: ∂V ∂T = f + ∂2V ∂X2 (7) A solution ˆ V (T, X) of (7) will imply that V (t, x) = ˆ V (t/τm, x/λm) will satisfy (6).

– p. 105

The reaction term

The form of f depends on the cell type we want to study. For the axon Iion(m, n, h, V ) of the HH-model is a good candidate. For the dendrite, which is non-excitable, a linear resistance model is good. Shift V so V = 0 is the resting potential: ∂V ∂T = ∂2V ∂X2 − V Need boundary and initial values. Initially at rest: V (X, 0) = 0

– p. 106

Boundary conditions

Types of boundary conditions: Dirichlet: V (xb, T) = Vb, voltage clamp. Neumann: ∂V

∂X = −riλmI, current injection.

Justification: ∂Vi ∂x = −riIi ⇒ ∂V ∂x − ∂Ve ∂x = −riIi

re=0

= ⇒ ∂V ∂x = −riIi

– p. 107

Gap junctions

– p. 108

Gap junctions

– p. 109

Gap junctions

Gap junctions are located between cell and ions may pass through them. The junctions have a high resistance to flow compared to the intra cellular environment.

– p. 110

Consider a 1D line of cell and assume that Fick’s law holds in the interior: J = −D dc dx Between cells we must have continuity of flow: −Ddc(x−

b )

dx = −Ddc(x+

b )

dx Here x−

b and x+ b indicates that the function is evaluated in the

limit from left and right, respectively. Furthermore we assume this flow to be proportional to the fall: J = F[c(x−

b ) − c(x+ b )]

where F is a permeability constant.

– p. 111

Would like relate F and D into an average, large scale, effective diffusion coefficient. Consider N cells of length L: J = −De ∆c ∆x = −De c1 − c0 NL Steady flux with fixed gradient. At steady state J is constant, so from J = −D(dc/dx) we have that c is linear in the interior. The solution c will be piecewise linear with jump at the cell boundaries. Continuity of flow over the gaps gives dc/dx = −λ for all

• interfaces. For the same reason the steps must all be equal,

c(x−

b ) − c(x+ b ) = ∆.

– p. 112

The size on ∆ and λ must fit the drop (c0 − c1): N intervals of length L: NLλ N jumps of size ∆: N∆ So in total we must have: NLλ + N∆ = c0 − c1 From the definition of F we have Dλ = F∆. In steady state the flux is the same on every scale: J = Dλ =

De NL(c0 − c1) = De NL(NLλ + N∆)

=

De NL(NLλ + N Dλ F ) = Deλ(1 + D FL)

So D = De(1 + D FL) ⇒ 1 De = 1 D + 1 FL

– p. 113

Calcium dynamics

– p. 114

Calcium dynamics

Calcium is an important ion in the biochemistry of cells Is used a signal carrier, i.e. causes contraction of muscle cells Is toxic at high levels The concentration is regulated through buffers and intra cellular compartments

– p. 115

Ryanodine Receptors

Sits at the surface of intra cellular calcium stores Endoplasmic Reticulum (ER) Sarcoplasmic Reticulum (SR) Sensitive to calcium. Both activation and inactivation. Upon stimulation calcium is released from the stores. To different pathways Triggering from action potential through extra cellular calcium inflow. Calcium oscillations observed in some neurons at fixed membrane potentials.

– p. 116

Compartments and fluxes in the model

– p. 117

Model equations

d[c] dt = JL1 − JP1 + JL2 − JP2 d[cs] dt = −JL2 + JP2 JL1 = k1(ce − c), Ca2+ entry JP1 = k2c, Ca2+ extrusion JL2 = k3(cs − c), Ca2+ release JP2 = k4c, Ca2+ uptake

– p. 118

The calcium sensitivity

Release modelled with Hill type dynamics: JL2 = k3(cs − c) = (κ1 + κ2cn Kn

d + cn)(cs − c)

– p. 119

Experiments and simulations

– p. 120

Good agreement between experiments and simulations Inactivation through calcium not included, but does not seem to be an important aspect

– p. 121

A more refined model

Inclusion of both activation and inactivation sites at the RyR Four states of the RyR

S00 No Ca ions attached, closed S10 Ca attached to activating site, open S01 Ca attached to inactivating site, closed S11 Ca attached both sites, closed

Define the fractions xi: x1 = S00/ST x2 = S10/ST x3 = S11/ST x4 = S01/ST = 1 − x1 − x2 − x3

– p. 122

The state transitions

The state transitions

S00

ck1

→ ←

k−1

S10

k−2 ↑↓ ck2 k−2 ↑↓ ck2

S01

ck1

→ ←

k−1

S11

Better models for the pumps JP = Vmax c2 K2 + c2

– p. 123

Model equations

dx1 dt = k−1x2 + k−2x4 − (k1 + k2)x1c dx2 dt = −k−1x2 + k−2x3 + (k1x1 − k2x2)c dx3 dt = (k2x2 + k1x4)c − (k−2 + k−1)x3 dc dt = vc(JL2 − JP2) + JL1 − JP1 dcs dt = −JL2 + JP2 JL1 = g2(ce − c) JL2 = (kfx2 + g1)(cs − c) JP1 = q1c2 q2

2 + c2

JP2 = p1c2 p2

2 + c2

– p. 124

Modeling cardiac propagation

– p. 125

Outline

Cardiac propagation. The bidomain concept. Derivation of the bidomain model from assumptions on the cell membrane and basic electromagnetic relations. A model for the surrounding body. Reduction to a monodomain model.

– p. 126

Cardiac propagation

Cardiac cells has two properties and corresponding function Excitable → Propagates the AP Contractive → Pumps blood Furthermore, the arrival of an AP triggers contraction. Cell to cell

• coupling. Two types:

Tight junctions: Transfer mechanical energy Gap junctions: Inter cellular channels where ions can flow

– p. 127

The conduction system

– p. 128

Cardiac propagation

Electrical signal initiated in the sinoatrial node (SA node). The action potential propagates the atria, which are insulated from the ventricles by a septum of non-excitable cells. The signal is conducted to the ventricles through the atrioventricular node (AV-node), located at the base of the atria. Conduction through the AV-node is fairly slow, but from there the action potential enters the so-called bundle of HIS, made up of fast-conducting Purkinje fibers.

– p. 129 – p. 130

The Purkinje fibers spread in a tree-like branching, ending

• n the endocardiac surface of the ventricles.

Muscle cells are stimulated at the end of the Purkinje fibers, and the action potential spreads through the muscle tissue. The electrical propagation in the heart is both

• ne-dimensional (along the purkinje fibers) and

three-dimensional (in the muscle tissue).

– p. 131

Modelling propagation in heart tissue

Because of the large number of cells, it is impossible to model the tissue by modeling each individual cell. The cells are connected, and the heart may be seen as consisting of two continuous spaces, the intracellular and the extracellular domain. The geometries of the two domains are too complex to be represented accurately.

– p. 132

The bidomain concept

Instead of accurately modeling the geometry of the two domains, they are assumed to be overlapping, both filling the complete volume of the heart muscle. Hence, every point in the myocardium lies in both the intracellular and the extracellular domain.

– p. 133

Basic equations

Maxwells equations states that ∇ × E + ˙ B = 0 where E and B are the strength of the electrical and magnetic field, respectively. Since ˙ B denotes the time derivative of the magnetic field, if the fields are stationary we have. ∇ × E = 0

– p. 134

The quasi-static condition

Although the electrical and magnetic fields resulting from cardiac activity are not stationary, the variations are fairly slow. We may therefore assume that the fields are stationary, an assumption called the “quasi-static condition”. As above, we have ∇ × E = 0, which means that the field E may be written as E = −∇u for some potential function u.

– p. 135

Current

In a conducting body, the electrical current is given by J = ME, where M is the conductivity of the medium. With the definition of E given above, the current is given by J = M∇u .

– p. 136

Following the bidomain concept introduced above, we introduce two electrical potentials: Intracellular potential ui Extracellular potential ue Since every point in the heart muscle lies in both the extracellular and intracellular domain, both ui and ue are defined in every point.

– p. 137

The currents in the two domains are given by Ji = −Mi∇ui, Je = −Me∇ue, and if we assume no accumulation of charge, the total current entering a small volume must equal the total current leaving the

• volume. This gives
• ∂V

(Ji + Je) · nds = 0 Since the volume V is arbitrarily chosen, this may be written as ∇ · (Ji + Je) = 0

– p. 138

Inserting the expressions for the currents, we get ∇ · (−Mi∇ui) + ∇ · (−Me∇ue) = 0 (This equation states that all current leaving one domain must enter the other domain.)

– p. 139

The two domains are separated by the cell membrane. Hence, all current going from one domain to the other must cross the cell membrane. We have −∇ · (−Mi∇ui) = ∇ · (−Me∇ue) = Im

– p. 140

We have previously modeled the membrane current Im as the sum of a capacitive current and an ionic current. However, that current was measured per membrane area, while we are now interested in the current per volume. The current per volume is achieved by multiplying with a scale factor χ, which is the ratio of cell membrane surface to cell volume. Im = χ(Cm ∂V ∂t + Iion), where V is the transmembrane potential.

– p. 141

To summarize, we now have two relations for the unknown potentials ui and ue: ∇ · (Mi∇ui) = χCm ∂V ∂t + χIion and ∇ · (Mi∇ui) + ∇ · (Me∇ue) = 0 We see that we have three unknown potentials ui, ue and V , and

• nly two equations. But V is defined as the difference between

the intracellular end the extracellular potential, and this may be used to eliminate one of the unknowns.

– p. 142

We have V = ui − ue, or ui = ue + V . Inserting this into the two equations, we get ∇ · (Mi∇(V + ue)) = χCm ∂V ∂t + χIion ∇ · (Mi∇(V + ue)) + ∇ · (Me∇ue) =

– p. 143

These equations may be rewritten as ∇ · (Mi∇V ) + ∇ · (Mi∇ue) = χCm ∂V ∂t + χIion ∇ · (Mi∇V ) + ∇ · ((Mi + Me)∇ue) = This is the standard formulation of the bidomain model.

– p. 144

Potential in the surrounding body

The tissue surrounding the heart is mostly non-excitable, meaning that the cells do not actively change their electric properties. The body surrounding the heart may hence be modeled as a passive conductor.

– p. 145

Extracardiac potential uo

Introducing the extracardiac potential uo, and using the arguments presented above, we derive the equation ∇ · (Mo∇uo) = 0, where Mo is the (averaged) conductivity of the tissue.

– p. 146

Boundary conditions

To complete the mathematical model, we need boundary conditions for V and ue on the heart surface, and for uo on the surface of the heart and the surface of the body. It is natural to assume that the body is insulated from its surroundings, implying that no current leaves the body. This gives the condition n · (Mo∇uo) = 0

• n the body surface.

– p. 147

Several choices of boundary conditions have been used on the heart surface. A common assumption is that the intracellular domain is insulated from the tissue surrounding the heart, while the extracellular domain connects directly to the surrounding tissue.

– p. 148

We get n · Ji = 0 ⇒ n · (Mi∇V + Mi∇ue) = 0 and ue = uo n · Je = n · Jo ⇒ n · (Me∇ue) = n · (Mo∇uo)

– p. 149

The complete model

∂s ∂t

= F(v, s) x ∈ H χCm ∂v

∂t + χIion(v, s)

= ∇ · (Mi∇v) + ∇ · (Mi∇ue) x ∈ H ∇ · ((Mi + Me)∇ue) = −∇ · (Mi∇v) x ∈ H Mi∇v · nT = x ∈ ∂H ue = uo x ∈ ∂H Me∇ue · nT = Mo∇uo · nT x ∈ ∂H ∇ · (Mo∇uo) = x ∈ T (Mo∇uo) · nT = x ∈ ∂T

– p. 150

Reduction to a monodomain model

The bidomain model is a very complex system of equations. Many (most) simulations are based on a simpler monodomain equation. The derivation of the monodomain model is based on the assumption of equal anisotropy ratios: Me = λMi

– p. 151

With this simplification, the bidomain equations may be written as χCm ∂v ∂t + χIion(v, s) = ∇ · (Mi∇v) + ∇ · (Mi∇ue) ∇ · (Mi(1 + λ)∇ue) = −∇ · (Mi∇v)

– p. 152

The second equation gives ∇ · (Mi∇ue) = − 1 1 + λ∇ · (Mi∇v), and if we insert this into the first equation we get χCm ∂v ∂t + χIion(v, s) = ∇ · (Mi∇v) − 1 1 + λ∇ · (Mi∇v)

– p. 153

Finally, we get χCm ∂v ∂t + χIion(v, s) = λ 1 + λ∇ · (Mi∇v) A problem with the monodomain model is that it can not be coupled directly to a surrounding body. Reproduction of ECG signals hence require additional computations.

– p. 154

Conclusions

Because of the excitability of cardiac cells, a simple volume-conductor model is not sufficient for modeling the heart muscle. By conceptually dividing the tissue into extracellular and intracellular domains, we are able to construct a mathematical model which describes signal propagation in the excitable tissue. By making a (non-physiological) assumption, the complex model may be reduced to a simpler monodomain model.

– p. 155

From PDE to Software: Numerical techniques for the bidomain model

– p. 156

Outline

The model Possible strategies. Operator splitting. Solving systems of ODEs. Solving the PDEs.

– p. 157

The model

∂s ∂t

= F(v, s) x ∈ H χCm ∂v

∂t + χIion(v, s)

= ∇ · (Mi∇v) + ∇ · (Mi∇ue) x ∈ H ∇ · ((Mi + Me)∇ue) = −∇ · (Mi∇v) x ∈ H Mi∇v · nT = x ∈ ∂H ue = uo x ∈ ∂H Me∇ue · nT = Mo∇uo · nT x ∈ ∂H ∇ · (Mo∇uo) = x ∈ T (Mo∇uo) · nT = x ∈ ∂T

– p. 158

Challenges

Highly complex system of PDEs and ODEs. Difficult to construct efficient solution algorithms. Simulation software tends to be complex and difficult to maintain.

– p. 159

Coping with the complexity

Several possible strategies exist for solving the bidomain equations. Using explicit schemes for the ODEs and the PDE involving the time derivative. Fully coupled approach, i.e. using a fully implicit time discretization to integrate the complete system simultaneously.

– p. 160

Explicit schemes

sn+1−sn ∆t

= F(vn, sn) x ∈ H χCm vn+1−vn

∆t

+ χIion(vn, sn) = ∇ · (Mi∇vn) + ∇ · (Mi∇un

e )

x ∈ H ∇ · ((Mi + Me)∇ue

n+1)

= −∇ · (Mi∇vn) x ∈ H Problem with the explicit approach: severe restrictions on the time step, and the second PDE can not be solved explicitly.

– p. 161

“Solution”: reduction to a monodomain model:

sn+1−sn ∆t

= F(vn, sn) x ∈ H χCm vn+1−vn

∆t

+ χIion(vn, sn) = ∇ · (Mi∇vn) x ∈ H Fully explicit solution scheme, but the approach still suffers from severe time step restrictions.

– p. 162

Fully implicit solution

sn+1−sn ∆t

= F(vn+1, sn+1) x ∈ χCm vn+1−vn

∆t

+ χIion(vn+1, sn+1) = ∇ · (Mi∇vn+1) + ∇ · (Mi∇ue

n+1)

x ∈ ∇ · ((Mi + Me)∇ue

n+1)

= −∇ · (Mi∇vn+1) x ∈ A fully implicit discretization avoids the strict stability restrictions

• n the time step, but requires solution of large systems of

non-linear equations for each time step.

– p. 163

Operator splitting

An attractive approach for handling the complex equations is to use some form of operator splitting. The complexity of the problem may be significantly reduced, without sacrificing too much in accuracy. Several approaches exist, but we will focus on a type of methods often referred to as fractional step methods.

– p. 164

A simplified problem

To introduce the ideas of operator splitting, we consider an initial value problem on the form u′(t) = (L1 + L2)u, u(0) = u0, (8) where L1 and L2 are differential operators.

– p. 165

For a given ∆t, an approximate solution at t = ∆t may be found by first solving v′(t) = L1v, 0 < t ≤ ∆t v(0) = u0, and then w′(t) = L2w, 0 < t < ∆t w0 = v(∆t), and finally set u(∆t) = w(∆t).

– p. 166

While it may seem that we have integrated the equation a step of length 2∆t, we have only included parts of the differential operator in each step. It is natural that the procedure of solving the equation in two steps introduces an error. This “splitting error” may be computed by comparing the Taylor expansion of the solution of the original equation and the solution obtained with the two-step approach.

– p. 167

In the Taylor expansion we use that ut = (L1 + L2)u . If L1 and L2 do not depend explicitly on t we can differentiate both sides of the equation to obtain utt = (L1 + L2)ut = (L1 + L2)2u In general we have ∂n

t u = (L1 + L2)nu

– p. 168

The solution of the original equation at time ∆t may now be written u(∆t) = u(0) + ∆t(L1 + L2)u(0) + 1 2∆t2(L1 + L2)2u(0) + . . . = (I + ∆t(L1 + L2) + 1 2∆t2(L1 + L2)2 + . . .)u(0)

– p. 169

Similarly, the solution obtained with the splitting method may be written as ˜ u(∆t) = (I + ∆tL2 + 1 2∆t2L2

2 + . . .)

(I + ∆tL1 + 1 2∆t2L2

1 + . . .)u(0)

= (I + ∆t(L1 + L2) + 1 2∆t2(L2

1 + 2L1L2 + L2 2) + . . .)u(0)

– p. 170

Subtracting ˜ u from u, we get ˜ u(∆t) − u(∆t) = 1 2∆t2(L1L2 − L2L1)u(0) + O(∆t3) We see that the splitting introduces an error proportional to ∆t2 for each time step, and after n steps this is expected to accumulate to an O(∆t) error.

– p. 171

The described technique gives a first order approximation to u at any given time step tn = n∆t. The technique is called Godunov splitting.

– p. 172

A slight modification of the Godunov algorithm is to solve the equation with the three-step algorithm v′(t) = L1v, 0 < t ≤ ∆t/2 v(0) = u0, w′(t) = L2w, 0 < t ≤ ∆t w(0) = v(∆t/2), v′(t) = L1v, ∆t/2 < t ≤ ∆t v(0) = w(∆t), and finally set u(∆t) = v(∆t).

– p. 173

By comparing the Taylor series of the solutions, it is now possible to show that the O(∆t2) terms now cancel, giving an error of O(∆t3) for each step. After n ≈ 1/∆t steps, the accumulated error is O(∆t2). This procedure is called Strang splitting.

– p. 174

• Op. splitting for the bidomain model

We consider the PDE containing the time derivative: χCm ∂v

∂t + χIion(v, s)

= ∇ · (Mi∇v) + ∇ · (Mi∇ue) x ∈ H In simplified form, this equation may be written as ∂v ∂t = (L1 + L2)v L1 is a linear PDE-operator, while L2 corresponds to the ionic current term Iion(s, v), and is a non-linear ODE-operator.

– p. 175

Inserting these expressions into the Strang splitting scheme, we get a three-step algorithm: Step 1. Solve ∂v ∂t = − 1 Cm Iion(v, s), tn < tn ≤ tn + ∆t/2. Step 2. Use the obtained approximation to v(tn + ∆t/2) as initial condition, and solve χCm ∂v ∂t = ∇ · (Mi∇v) + ∇ · (Mi∇ue), tn < t ≤ tn + ∆t . Step 3. Finally, use the new value of v as initial condition, and solve ∂v ∂t = − 1 Cm Iion(v, s), tn + ∆t/2 < t ≤ t + ∆t .

– p. 176

We see that the equation is solved by alternately solving an ODE and a PDE. The ODE may be integrated together with the cell model ODEs, and the PDE (which is now linear) may be solved simultaneously with the other PDEs of the model. Free choice of solvers for the two sub-problems (sub-problem solving accuracy must be at least 2 to get

• verall second order accuracy).

– p. 177

To simplify the operator splitting formulation for the complete model, we introduce the notation vn

∗ for the value of v after Step 1, i.e. after the ODE solver

has been applied. vn+1

∗

for the value after the second step, i.e. after the PDE system is solved but before the second application of the ODE solver sn+1/2 = s(tn + ∆t/2)

– p. 178

The complete model

We assume vn and sn are known: 1. Compute sn+1/2 and vn

∗ by solving the system ∂v ∂t = − 1 Cm Iion(v, s),

v(tn) = vn

∂s ∂t = F(v, s),

s(tn) = sn , tn < t < tn+1/2. 2. Compute vn+1

∗

(and new values of ue and uo) by solving the coupled PDE system χCm ∂v

∂t

= ∇ · (Mi∇v) + ∇ · (Mi∇ue) x ∈ H ∇ · ((Mi + Me)∇ue) = −∇ · (Mi∇v) x ∈ H ∇ · (Mo∇uo) = x ∈ T for tn < t ≤ tn+1, with vn

∗ as initial condition.

3. Compute sn+1 and vn+1 by solving the system

∂v ∂t = − 1 Cm Iion(v, s),

v(tn+1/2) = vn+1

∗ ∂s ∂t = F(v, s),

s(tn+1/2) = sn+1/2 , tn+1/2 < t ≤ tn+1.

– p. 179

In the continuous case all variables, including the cell model state vector s, are defined throughout the heart muscle. The solution of the equations is based on discretizing the computational domain in a finite number of points (a grid/mesh), and use a numerical method to determine all variables in each point. Applying this type of discretization to the ODEs in step 1 and 3 of the solution algorithm gives one ODE system for each point in the grid.

– p. 180

The computational work for each time step consists of two separate tasks:

• 1. Solving a large number of non-linear ODE-systems (Step 1

and 3).

• 2. Discretizing and solving a coupled PDE system (Step 2).

Both tasks contribute significantly to the total computational time, so it is important to choose efficient solvers for both sub-problems.

– p. 181

Solving the ODEs

The task in Step 1 and Step 3 is basically to solve a large number of initial value problems on the form dy dt = f(t, y) y(t0) = y0. The form and complexity of the ODE system depends on the chosen cell model.

– p. 182

Simpler cell models: Hodgkin-Huxley (4 ODEs), Beeler-Reuter (8 ODEs), Fitzhugh-Nagumo (2 ODEs). More realistic models: Luo-Rudy (13 ODEs), Winslow (31 ODEs) ++ The more realistic cardiac cell models typically have much faster dynamics than the simpler models, adding further to the difficulty

• f solving these equations.

– p. 183

The Fitzhugh-Nagumo model

100 200 300 400 −0.5 0.5 1 Original FitzHugh−Nagumo model

– p. 184

The Winslow model

100 200 300 400 −100 −50 50

– p. 185

Solvers for ODE systems

A huge variety of methods have been developed for solving systems of ODEs. Simple, well known methods are the forward and backward Euler methods. A (large) family of widely used methods are the so-called Runge-Kutta methods.

– p. 186

Runge-Kutta methods

Definition: Let bi, aij, ci (i, j = 1, . . . , s) be real numbers. The method ki = f  t0 + ci∆t, yn + ∆t

s

• j=1

aijkj   yn+1 = yn + ∆t

s

• i=1

biki is called an s-stage Runge-Kutta method.

– p. 187

Explicit methods

If aij = 0 for i ≥ j we have an explicit (ERK) method, that can be expressed by k1 = f(t0, y0) k2 = f(t0 + c2∆t, y0 + ∆ta21k1) k3 = f(t0 + c3∆t, y0 + ∆t(a31k1 + a32k2)) . . . ks = f(t0 + cs∆t, y0 + ∆t(as1k1 + · · · + as,s−1ks−1)) y1 = y0 + ∆t(b1k1 + · · · + bsks)

– p. 188

Explicit RK methods are the classical, most well-known RK methods. Efficient methods with good accuracy. Stability problems when applied to the stiff cell model equations.

– p. 189

Integrating the ODE-system in the Winslow cell model from t = 0 to t = 400. A good explicit RK-solver requires about 60000 timesteps.

50 100 150 200 250 300 350 400 −100 −50 50 – p. 190

A semi-implicit method

A class of implicit methods is commonly referred to as semi-implicit methods. An example such a method is the following 4-step scheme: k1 = f(tn, yn) k2 = f (tn + ∆tc2, yn + ∆t (a21k1 + γk2)) k3 = f

• tn + ∆tc3, yn + ∆t
• ˆ

b1k1 + ˆ b2k2 + γk3

• k4

= f (tn + ∆tc4, yn + ∆t (b1k1 + b2k2 + b3k3 + γk4)) yn+1 = y0 +

4

• i=1

biki en+1 =

4

• i=1

(bi − ˆ bi)ki.

– p. 191

Main differences between explicit methods and the semi-implicit method. A system of non-linear equations must be solved to determine each kj. For an ODE system of n equations, each kj is computed by solving a system of n algebraic equations. The stability is dramatically improved.

– p. 192

Application to the cell model ODEs

Numerical experiments have confirmed that implicit methods are much more well-suited for solving the cell-model ODE systems than explicit methods. The efficiency difference between different implicit methods (semi- and full) is much smaller, and depends on the accuracy requirements. Because of the coupling to the PDEs, the accuracy demands for our application are generally low.

– p. 193

The semi-implicit method

As an example, consider the 4-step semi-implicit method introduced above. Application of this method to the cell model ODEs involves the following tasks:

• 1. Compute the first stage derivative (explicitly) from

k1 = f(tn, yn)

• 2. For j = 2, 3, 4 use kj−1 as start vector and compute kj by

solving a non-linear equation with Newtons method.

• 3. Compute yn+1 from

yn+1 =

4

• i=1

biki

– p. 194

The amount of work for each time step is considerably larger than for explicit methods, because we need to solve non-linear algebraic equations. For stiff problems like the cell model ODEs, the extra work per time step is easily outweighed by the enhanced stability properties, allowing the equations to be integrated with much larger time steps. The required number of time steps to integrate the Winslow ODE system with an implicit method is approximately 150. (Compared to about 60000 for the explicit method).

– p. 195 50 100 150 200 250 300 350 400 −100 −50 50 – p. 196

Discretization of the PDE system

χCm ∂v

∂t

= ∇ · (Mi∇v) + ∇ · (Mi∇ue) x ∈ H ∇ · ((Mi + Me)∇ue) = −∇ · (Mi∇v) x ∈ H ∇ · (Mo∇uo) = x ∈ T

– p. 197

Outline of PDE discretization

To obtain the overall second order accuracy of the Strang splitting, the PDE system must be solved with second order

• accuracy. This is achieved with the Crank-Nicolson scheme.

To simplify the spatial discretization of the problem, the second and third equation of the system are combined to

• ne equation.

A standard finite element method is used for the spatial discretization of the resulting equations.

– p. 198

Time discretization

A widely used scheme for time discretization of PDEs is the Crank-Nicolson scheme ∂v ∂t = G(v), an approximate solution is found by vn+1 − vn ∆t = 1/2G(vn+1) + 1/2G(vn) .

– p. 199

Crank-Nicolson for coupled problems

For a more complex problem on the form ∂u ∂t = G(v) + H(u), we may apply a combination of Crank-Nicolson and themidpoint rule un+1 − un ∆t = 1/2G(vn+1) + 1/2G(vn) + H(un+1/2) . This also gives a second order accurate scheme.

– p. 200

Application to the bidomain PDEs

∇ · Mi∇(1/2vn+1

∗

+ 1/2vn

∗ )

+∇ · Mi∇un+1/2

e

= vn+1

∗

− vn

∗

∆t x ∈ H, (9) ∇ · Mi∇(1/2vn+1

∗

+ 1/2vn

∗ )

+∇ · (Mi + Me)∇un+1/2

e

= 0 x ∈ H, (10) ∇ · MT ∇un+1/2

T

= 0 x ∈ T, (11)

– p. 201

Notes on space discretization

To simplify the space discretization of the problem, the continuity conditions on the heart surface are used to reduce the three equations to two PDEs. The result is one equation which is defined only in the heart muscle, and one which is defined throughout the heart and the surrounding body. A standard finite element technique is used to discretize the resulting equations.

– p. 202

Notes on implementation

The main goal of the operator splitting technique was to handle the complexity of the problem. Although the solution algorithm looks complex, the splitting

• f the problem makes it much easier to implement the

algorithm. The modular nature of the discretized problem fits particularly well with object oriented programming.

– p. 203

The simulator software: Numerical intensive part implemented in compiled languages (C++/C/Fortran). Python layer on top calling the different components

– p. 204

Design advantages

Ease of debugging and maintenance. Easy to switch solvers (and models) for smaller parts of the problem.

– p. 205

Summary

Developing an efficient and reliable simulator based on the bidomain model is a challenging task, because of the complexity of the equations. The most commonly used techniques are different forms of explicit schemes. These have poor numerical properties for the problem, but are easy to construct and implement. Various operator splitting techniques are good alternatives for constructing numerically efficient schemes for the complex problem.

– p. 206

Strang splitting gives second order accuracy, as long as the sub-problems are solved to at least this accuracy. The splitting algorithm reduces the non-linear PDE problem to linear PDEs and non-linear ODEs. PDEs: A Crank-Nicolson scheme is used for time discretization. Standard finite element techniques can be used for the spatial discretization.

– p. 207

ODEs: Stiff, non-linear, fairly complex systems. Because of stiffness, explicit methods perform poorly. Implicit RK methods perform well for the present application. Implicit methods are more complicated to implement than explicit methods, and require more work per time

• step. This is outweighed by the improved stability

(larger time steps).

– p. 208

The modular nature of the splitting algorithms suits particularly well for object oriented programming, and makes it easier to develop and maintain the simulation software.

– p. 209