Towards mesoscopic modeling of firing neurons: a feasibility study

UPTEC F 14011

Examensarbete 30 hp April 2014

Towards mesoscopic modeling of firing neurons: a feasibility study

Emil Berwald

Teknisk- naturvetenskaplig fakultet UTH-enheten

Besöksadress:

Ångströmlaboratoriet Lägerhyddsvägen 1 Hus 4, Plan 0 Postadress:

Box 536 751 21 Uppsala Telefon:

018 – 471 30 03 Telefax:

018 – 471 30 00 Hemsida:

http://www.teknat.uu.se/student

Abstract

Towards mesoscopic modeling of firing neurons: a feasibility study

Emil Berwald

Ion channel models are related to non-equilibrium statistical physics, fluid mechanics and electromagnetism. Some classes of ordinary differential equations that model ion channels can be seen as a limit of finite state-space continuous-time Markov chains.

The purpose of this thesis is to qualitatively investigate the numerical results of systems of equations that incorporate ion channels modeled by such Markov chains and an electrical circuit model of a single neuron with isopotential extracellular space.

This may be useful for making more detailed micro-physical simulations of neurons.

A subset of the Rallpack benchmarks is conducted in order to evaluate the accuracy of the electrical circuit model of the transmembrane voltage propagation. In order to test the tau-leap method employed to simulate the Markov-chain based ion channel models a cylindrical geometry is implemented. Convergence properties are presented in terms of mean interspike intervals of the transmembrane voltages for different time- and spatial discretisations. Accuracy of the tau-leap method is presented in relation to the deterministic versions of the ion channel models.

The results show that the method used to simulate the transmembrane voltages is accurate and that while the tau-leap method is convergent in the mean interspike interval sense, it is not conclusive how accurate it is compared to the corresponding ordinary differential equations or how efficient it is.

Examinator: Tomas Nyberg

Ämnesgranskare: Per Lötstedt

Handledare: Stefan Engblom

Popul¨ arvetenskaplig sammanfattning

Nervcellernas huvudfunktion i kroppen ¨ ar att skapa informationsfl¨ oden. P˚ a makroskopisk niv˚ a kopplas olika delar av hj¨ arnan och det perifera nervsys- temet samman av nerver. Enskilda nervceller kan kommunicera med varan- dra via synapser. Synapser reglerar nervsignalers fortskridning genom kemiska reaktioner med neurotransmittorer, som till exempel f¨ orts dit inuti vesiklar. I nervceller r¨ or sig nervsignaler i form av lokala elektriska sp¨ anningsskillnader

¨

over cellmembranet. Sp¨ anningssignalers r¨ orelse styrs till st¨ orsta delen av mem- branprotein som kallas jonkanaler. De kan ha v¨ aldigt komplex struktur som till˚ ater dem att sl¨ appa igenom eller st¨ anga ute vissa sorters joner fr˚ an att transporteras ¨ over membranet, och d¨ armed i f¨ orl¨ angningen p˚ averka membranets sp¨ anningsskillnader.

I sp¨ anningsstyrda jonkanaler regleras genomsl¨ appligheten via den lokala mem- bransp¨ anningen. Detta kan modelleras p˚ a flera s¨ att. Ett av de mest detaljerade s¨ atten ¨ ar via stokastiska Markovkedjor. Ett annat alternativ ¨ ar via determinis- tiska ordin¨ ara differentialekvationer (ODEer). Denna uppsats anv¨ ander sig av en teknik som konstruerar en Markovkedjemodell utifr˚ an en vanlig typ av ODE i detta sammanhang. Via denna konstruktion j¨ amf¨ ors beteendet f¨ or jonkanaler i deterministisk och stokastisk tappning i denna uppsats.

Nervcellen representeras h¨ ar matematiskt som en elektrisk krets vars dy- namik beskrivs med ett system av ODEer. Jonkanalernas ¨ oppnande och st¨ angande styrs i det stokastiska fallet av den lokala membransp¨ anningen och Markovked- jemodellerna, som simuleras med hj¨ alp av tau-leap metoden. Tau-leap metoden uppskattar hur jonkanalerna beter sig ¨ over sm˚ a tidssteg, likt forward Euler inom finita differenser. I det deterministiska fallet styrs jonkanalerna av ODEer som integreras ¨ over tidssteget. F¨ or att kontrollera att den deterministiska vari- anten av kretsmodellen ¨ ar korrekt genomf¨ ors en delm¨ angd av Rallpack-testerna (Bhalla et al., 1992). Sedan analyseras konvergens, noggrannhet och effektivitet hos den stokastiska varianten genom att tiden fr˚ an topp till topp i sp¨ anningssig- nalen hos en exempelmodell j¨ amf¨ ors i olika tids- och rumsdiskretiseringar.

Resultaten ¨ ar att den deterministiska modellen ¨ overenst¨ ammer med refer- ensl¨ osningarna i Rallpack-testerna och d¨ arf¨ or ¨ ar noggrann. I exempelmodellen noteras att noggrannare diskretiseringar leder till konvergenta resultat. Det ¨ ar dock n˚ agot oklart hur noggrann den stokastiska modellen ¨ ar, eftersom den inte helt ¨ overenst¨ ammer med den deterministiska varianten. Detta beh¨ over dock inte tyda p˚ a n˚ agot fel hos den stokastiska modellen. Effektiviteten hos tau-leap metoden ¨ ar sv˚ ar att uppskatta med tanke p˚ a att implementeringen av den ibland

¨

overg˚ ar till en mer noggrann metod. F¨ orslagsvis kan en annan exempelmodell

konstrueras som testar fler aspekter av tau-leap metoden.

Contents

1 Purpose 1

2 Biological background 1

2.1 Macroscopic level . . . . 1

2.2 Cellular level . . . . 1

2.3 Microscopic level . . . . 3

3 Mathematical background 4 3.1 Notation . . . . 4

3.2 Exponential and Poisson random variables . . . . 4

3.3 Homogeneous finite state space continuous-time Markov chains . 5 3.3.1 Lumpability . . . . 6

4 Physics background 6 4.1 Transport theory . . . . 6

4.2 Electromagnetism . . . . 7

4.2.1 Continuity equation . . . . 7

4.2.2 Maxwell’s equations . . . . 7

4.2.3 Response fields D and H . . . . 8

4.2.4 Movement of charges . . . . 8

4.2.5 Approximations of the Maxwell equations . . . . 9

5 Neuron model 11 5.1 Geometry and physical parameters . . . . 11

5.2 Cable equation . . . . 11

5.2.1 Total axial current flux . . . . 12

5.2.2 Total membrane current flux . . . . 14

5.2.3 Capacitive effects . . . . 15

5.2.4 Matrix model . . . . 15

6 Ion channel model 16 6.1 Gating-particle models . . . . 17

6.2 Many Markov chains together . . . . 20

7 Simulation algorithms 21 7.1 Simulating the Markov chains . . . . 21

7.2 Simulating the potential . . . . 23

8 Test models 24 8.1 Rallpack . . . . 25

8.2 Tau-leap . . . . 26

9 Results 30 9.1 Rallpack . . . . 30

9.2 Tau-leap . . . . 34

10 Conclusions and Discussion 36

10.1 High accuracy of the ODE-schemes . . . . 36

10.2 Convergence of mean interspike interval . . . . 37

10.3 Inconclusive accuracy and efficiency of tau-leap model . . . . 37

10.4 Discussion and future work . . . . 38

Appendices 39 A Algorithms 39 A.1 Direct stochastic simulation algorithm (SSA) . . . . 39

A.2 Tau-leap stochastic simulation algorithm . . . . 40

A.3 Voltage simulation . . . . 42

B Mathematics 43 B.1 Logic and equivalence relations . . . . 43

B.2 Probability . . . . 43

B.3 Hadamard product . . . . 44

B.4 Stokes’ theorem for grad, div and curl . . . . 44

List of Figures 1 Close up of pyramidal neuron . . . . 2

2 Nervous system . . . . 2

3 Cell membrane . . . . 3

4 Geometry of neuron and directed adjacency matrix . . . . 11

5 Compartments . . . . 12

6 Connected cylinders . . . . 13

7 Ion channel circuit model . . . . 14

8 Cable equation circuit with isopotential outside (ground) . . . . 16

9 Sodium channel . . . . 19

10 SSA condition K-channel . . . . 28

11 SSA condition Na-channel . . . . 29

12 Rallpack 1 . . . . 31

13 Rallpack 2 . . . . 32

14 Rallpack 3 . . . . 33

15 Example run . . . . 34

16 Tau-leap mean interspike interval . . . . 35

17 Gating-particle mean interspike interval . . . . 36

List of Tables 1 Approximate constants and dimensions for neurons . . . . 10

2 Channel data . . . . 25

3 Rallpack output . . . . 30

4 Mean simulation time . . . . 36

1 Purpose

The purpose of this thesis is to qualitatively investigate the numerical results of systems of equations that model a neuron with stochastic ion channels and an isopotential extracellular space. The electromagnetic effects will be included through an electrical circuit model based on a system of ordinary differen- tial equations (ODEs). The stochastic effects will be included through Markov chains. The questions to answer are:

1. Is the implementation of the electrical circuit model accurate?

2. Is the used Markov chain method accurate, convergent and efficient?

The accuracy of the electrical circuit model will be tested through a subset of the Rallpack tests [1] [2] with ion channels modeled by coupled ODEs. To evaluate the accuracy of the Markov chain method the results will be compared to the results of the corresponding ODE-systems. The convergence and efficiency of the Markov chain method will be evaluated by analyzing the results of a parameter test.

2 Biological background

In this section some biological aspects of neurons will be covered. The mecha- nisms used to transmit signals via neurons are fascinating and complex. There are several levels of communication to consider, and the intricate connections be- tween them are not fully explored. Poorly functioning neurons underlie various neurodegenerative diseases such as Parkinson’s, Alzheimer’s and Huntington’s.

Assistive technologies such as hearing aids and biomechatronic prosthetics de- pend on the plasticity and ability of neurons to interpret signals. In Figure 1 one type of neuron is depicted.

2.1 Macroscopic level

Neurons are responsible for encoding information in the central nervous system, see Figure 2. Lower level functions are many times gathered at specific parts of the brain, processing information received from neurons throughout the body, often by means of signalling pathways called nerves which are bundles of ax- ons that reach out from the neurons. Higher level functions depend on large networks of neurons with various feedback loops.

2.2 Cellular level

Inside neurons information is transmitted by electric potentials as they propa-

gate throughout the cell. A process that underlies propagation is the regulation

of ion concentrations in both the intracellular cytoplasm and the extracellular

space caused by integral membrane proteins called ion channel gates and ion

pumps. Connections to other neurons are made possible at synapses where

neurotransmitters are extruded thanks to vesicles into the extracellular space

between the neurons. The neurotransmitters, in turn, are parts of a chemical

process which may start a potential wave in the postsynaptic neuron known as

an action potential.

Figure 1: ”SMI32 - immunoreactive pyramidal neuron in medial prefrontal cortex of macaque.” Brainmaps.org

Figure 2: The nervous system. theEmirr / Wikimedia Commons (CC-BY-3.0)

2.3 Microscopic level

There are many variants of ion channel proteins. Structural-functional studies and x-ray crystallography can show in very fine detail what the components of the ion channels are [3]. If the conformal state of the channel changes into an open state, they allow the passage of certain ions through pathways called pores which links the intra- and extracellular space, see Figure 3.

Figure 3: Picture of the cell membrane.

Wikimedia Commons

If the pores are closed, the ions are blocked at the membrane and cannot easily pass through. Some channels open due to ligands binding to receptor subunits of the channel protein, which can happen both on the inside of the membrane (second messengers) or outside of the membrane (ligand-gated), while others can open in an electrical voltage dependent way. Differently gated ion channels exists in some neurons such as heat sensitive nociceptors or photon sensitive neurons in the eye and so on.

Measuring some form of electrically induced response in a neuron using tech- niques such as the voltage clamp, current clamp or other experimental setup biologists can formulate mathematical models of the neuron. By conducting experiments on voltage dependent ion channels (for example using toxins to turn off certain ion channels), it is common to construct a mathematical model in the form of a system of ordinary differential equations. Based on the time rates in these equation systems a more detailed stochastic model in the form of a finite state space, continuous-time Markov chain can be constructed. In one study based on dynamical clamp experiments that link in silico models together with in vitro neurons the deterministic models were unable to phase lock their activity to real neurons, in contrast to the more detailed stochastic models [4, p. 331].

There are different levels of neuron simulation. In a network setting it will

be harder to create a computer model with appropriate detail, accuracy, validity

and time and data complexity. There is also the problem of not knowing the

exact physical configuration of the biological system. In this thesis a computa-

tionally tractable model of a single neuron will be suggested. It is driven by a

detailed biophysical description on the mesoscopic level.

3 Mathematical background

This section includes a summary of the mathematical notation used in the thesis as well as the underlying stochastic theory needed to simulate the ion channels with Markov chains.

3.1 Notation

The Iverson bracket explained in Appendix B.1 is used whenever a logical ex- pression is enclosed in brackets (e.g. rλ

i

¡ rs). Matrix equations will either be in component form (e.g. δ

ij

 ri  js), or in matrix form, for which they will be denoted by either boldface (e.g. A) or with indices enclosed by brackets (e.g. rA

ij

s). Integral theorems are summarised in Appendix B.4, note that the oriented line, area and volume elements will be denoted ˆ tdL

_g

, ˆ ndA

_g

and dV

_g

.

3.2 Exponential and Poisson random variables

If a random variable T : Ω Ñ R belongs to the exponential distribution with parameter lambda, T  Exppλq, its probability density function and cumulative distribution function is [5, pp. 109,126]

f

_T

ptq  λ expp
λ  tqrt ¡ 0s

F

T

ptq  PpT ¤ tq  r1
expp
λ  tqsrt ¥ 0s

In [6, p. 72] we have that for a countable sequence of independent random variables pT

j

q

_jPJ

, T

_j

 Exppq

j

q, 0   °

jPJ

q

_j

  8, there is an exponentially distributed random variable

T  inf

iPJ

T

_i

 Exp

 ¸

jPJ

q

_j



(3.1)

where the index of the random variable which has the infimum value has a distribution

PpT

k

 T q  ° q

k jPJ

q

j

(3.2)

independent of T , in other words

PpT

k

 T ^ T ¥ tq  q

_k

°

jPJ

q

j

exp p
¸

jPJ

q

_j

t q

A random variable K : Ω Ñ N

^¥0

belongs to the Poisson distribution with parameter lambda, K  P oissonpλq, if the probability mass function and cu- mulative distribution function is [7]:

p

T

pkq  PpK  kq  Pptω P Ω: Kpωq  kuq  expp
λq λ

^k

k!

F

T

pkq  PpK ¤ kq 

¸

tku i0

p

T

piq  1 tku!

»

₈

λ

y

^tku

exp p
λqdy

3.3 Homogeneous finite state space continuous-time Markov chains

For a Markov chain with state space S, the probability that it is in state i at time t is given by the probability distribution λ

i

ptq, °

iPS

λ

i

ptq  1. A Markov chain possesses the Markov property: the conditional probability of being in any allowed subset of the state space at any point in the future is the same whether given the entire history of visited states or given only the state cur- rently occupied [8, p. 167]. The Markov chain is homogeneous if the transition probabilities between any two states only depend on the time difference and not the current time [8, p. 168]. A closed class is a subset of the state space such that no state outside the subset is accessible from states inside it. If the only closed class of the Markov chain is its entire state space it is called irreducible [8, p. 32]. A state is recurrent if it almost surely is accessible to itself infinitely many times [6, p. 114]. The rate of the transition from state i to j is given by the infinitesimal generator Q

ij

which is related to the transition probability P

ij

p∆tq by:

Q

_ij

 lim

∆tÑ0

P

_ij

p∆tq
ri  js

∆t rP

ij

sp∆tq  exppQ∆tq  ¸

⁸

n0

rQ

ij

s

ⁿ

∆t

ⁿ

n!

From [6, pp. 108,109] it is found that the infinitesimal generator of a finite state space Markov chain satisfies:

@i:
8   Q

ii


¸

ji

Q

ij

¤ 0, @j  i: 0 ¤ Q

ij

¤ 8

The jump times J

_{n 1}

, holding times S

_{n 1}

and embedded discrete-time jump chain Y

_n

of a finite Markov chain [6, pp. 69,90] are:

S

n 1

 J

n 1

J

n

, J

n 1

 inftt ¥ J

n

: X

t

 X

J_n

u, J

0

 0 Y

n

 X

J_n

Given the current state i the holding time S

_{n 1}

is an exponentially distributed random variable which is the infimum of exponentially distributed random vari- ables T

_ij

 ExppQ

ij

q over j  i [ 6, p. 89] (compare with (3.1)):

T

i

 inf

ji

T

ij

 Exp

 ¸

ji

Q

ij



the embedded jump matrix represents the probability that the next transition will be from i to j (compare with (3.2)):

PpX

J_{n 1}

 j | X

J₀

, . . . , X

J_n

 iq  Π

ij



$ ' ' '

&

' ' ' %

Q

_ij

°

ji

Q

ij

if j  i ^ °

ji

Q

_ij

 0

0 if j  i ^ °

ji

Q

ij

 0 rj  is if °

ji

Q

ij

 0

and is independent of the holding time [6, pp. 87-89,109][8, p. 180].

If Q is irreducible and recurrent there exists a non-negative left kernel which is unique up to scalar multiples [6, p. 118]. The limiting probability distribution is the unique stationary distribution satisfying global balance and conservation of probability, that is [4, p. 32]:

0  λ

^J

Q ô 0

^J

 Q

^J

λ @i: 0 ¤ λ

i

¤ 1, ¸

i

λ

i

 1 (3.3)

3.3.1 Lumpability

With an equivalence relation over the state space „ S  S the quotient set S {  is the partition of the state space with respect to the equivalence relation;

it groups them according to equivalence (see Appendix B.1). A time continuous Markov chain with state space S and infinitesimal generator Q

ij

is lumpable (and the lumped chain is a time continuous Markov chain as well) for any initial distribution vector, if the row sums

Q

_risrjs

 ¸

jPrjsPS{

Q

ij

, i P ris P S{  (3.4)

are equal for all i P ris for all pairs ris, rjs as proved in [ 9, pp. 66,69, 10, p. 124].

4 Physics background

In this section the aim is to show the physics that is behind the scenes of the voltage propagation of neurons. Another purpose of this section is to motivate the simplification of the full electromagnetic theory into the theory of electric circuits, which is the setting the neuron simulation in this thesis will use.

4.1 Transport theory

From [11, pp. 216-220,454] there is the following transport theorem for a vector field X with flow F

t

:

_dt^d

F

t

pxq  XpF

t

pxqq and function fpx, tq:

d dt

»

FtpUq

f dV

g



»

FtpUq

rB

t

f ∇ pfXqsdV

g

Allow the time rate of the extensive physical property Ψ, with density ψ, within volume F

t

pUq to change by the two mechanisms:

a net flux:

¾

BFtpUq

n

ψ

 ˆndA

g

, and a net production:

»

F_tpUq

π

ψ

px, tqdV

g

where ˆ n is the unit outward pointing normal vector field and n

_ψ

px, tq is the flux density of the property, and π

ψ

px, tq is the temporal rate of production of the property per unit volume. Assuming the shape of the volume is such that Stokes’

theorem for differential forms applies the transport equation for the vector field v following the fluid element in some sense then gives this continuity equation:

d dt

»

F_tpUq

ψdV

_g



»

F_tpUq

r
∇pn

ψ

q π

ψ

sdV

g



»

F_tpUq

rB

t

ψ ∇ pψvqsdV

g

(4.1)

If one chooses F

_t

pUq  U and v  0 instead one gets the continuity equation for a fixed region.

4.2 Electromagnetism

This subsection will include a summary of the electromagnetic theory of Maxwell as well as some theory for the movement of particles due to electrical forces. To conclude the subsection arguments in favor of using the simplified electric circuit theory will be presented to motivate the use of them in the remaining parts of the thesis.

4.2.1 Continuity equation

The electric charge q is an extensive thermodynamic quantity. It can either be positive or negative. Examples are protons with positive charge and electrons with negative charge. The charge density ρ is an intensive thermodynamic quantity describing the amount of charge per volume (or area or length).

Denote the current density or flux density of the electric charge by J rather than n

_ρ

. As a thought model, with classical particles the contribution of a single particle to the net flow through an oriented differential element is its charge multiplied by its ”trajectory flux” during dt: positively counted if it crossed the boundary towards the outside and negatively counted if it crossed the boundary towards the inside. The current density associated with an area element is the volume current density, similarly there are lower dimensional current densities.

By electrical charge conservation we know that processes involving electri- cally charged particles leave the net electric charge unchanged in time and are local phenomena [13, p. 10]. Therefore there is no net production of electric charge within a volume so the continuity equation (4.1) for a fixed region be-

comes: »

V

B

t

ρ ∇  JdV

g

 0 (4.2)

4.2.2 Maxwell’s equations

In [14] a method to go from relativistic quantum electrodynamics to macroscopic electromagnetism is proposed, to fill the gap between empirical knowledge and theory. From [14] we have the following definitions of the microscopic Maxwell’s equations:

ρ pxq  ¸

l

q

_l

δ px
x

l

q, Jpxq  ¸

l

q

_l

dx

l

dt δ px
x

l

q (4.3)

Gauss’s laws ∇  E  

¹0

ρ, ∇  B  0 (4.4)

Maxwell-Faraday equation ∇  E 
B

t

B (4.5)

Amp` ere’s law with Maxwell’s correction ∇  B  µ

0

pJ 

0

B

t

E q (4.6) where the interface conditions are derived from their integral forms [15].

Using ∇ p∇pXqq  0, ∇p∇pfqq  0, ( 4.4), (4.5) and demanding that the electric and magnetic fields E, B are invariant of the choice of gauge potentials φ, A, and gauge transforms Λ gives the defining relations [16, p. 240]:

B  ∇ A A Ñ A

^{Λ 1}

: A

¹

 A ∇ Λ

E 
∇ φ
B

t

A φ Ñ φ

^{Λ 1}

: φ

¹

 φ
B

t

Λ

The approach where the charge is divided into free or true parts (subscript f or t) and polarised or bound parts (subscript P or b) and the current into free or conducting parts (subscript f or c), polarisation parts (subscript P ) and magnetization parts (subscript M ) is unfortunately not uniquely defined and problematic in other ways as well [14, pp. 12-16]. Since some of the references used in this thesis uses such a decomposition it will be included anyway [14]:

D  

0

E P, H  µ

¹0

B
M ρ  ρ

t

ρ

P

 ρ

t

∇ P

J  J

c

J

_P

J

_M

 J

c

B

t

P ∇  M

Gauss’s law for the electric field ∇  D  ρ

t

(4.7) Amp` ere’s law with Maxwell’s correction ∇  H  J

c

B

t

D (4.8) As an example, one of the problematic approaches groups particles into freely moving entities of one or more members, where the group members are consid- ered bound to each other (such as the particles in atoms, ions or molecules), where the bound particles in a group are referenced relative to the group posi- tion [18, p. 1191, 13, pp. 56–57, 16, p. 251]. It utilizes either ensemble averages or weighted spatial averages xy that includes ”many” particles [ 17, p. 23][18][16, p. 249][13, p. 215]. In [18, p. 1191] the particles are partitioned into groups m with one or more members in each group, using point x

m

to represent group m:

ρ

_f

pxq x ¸

mPM

p ¸

jPm

q

_j

qδpx
x

m

qy (4.9)

J

_f

pxq x ¸

mPM

p ¸

jPm

q

_j

q dx

_m

dt δ px
x

m

qy (4.10) where groups of only one member are called conduction charges.

4.2.3 Response fields D and H

The approach where the charge is divided into free and bound charge and current can give expressions of the response fields that unfortunately do not include quantum effects and large scale effects [18] [17] [19].

Usually macroscopic constitutive equations are used such as the permittivity

: D  pE, xq and the permeability µ: B  µpH, xq. In linear response media D  

0



r

E  

0

p1 χ

e

qE

B  µ

0

µ

r

H  µ

0

p1 χ

m

qH

with χ

constant in linear homogenous isotropic nondispersive media [14][15].

4.2.4 Movement of charges

From [20, pp. 346-351] we get Lorentz’ force law per unit volume:

f  f

^j

e

j

 ρE ρv  B, ρv  J

Note that infinitesimally f 

^Br_Bt

dt  ρE  vdt  E  Jdt. The total rate of work done by the fields in a finite volume V is [16, pp. 151,153,197,189]:

dW dt 

»

V

E  J

f

dV

_g



»

V

∇pE  Hq
H  B

t

B
E  B

t

DdV

_g

For most substances, the current is proportional to the force per unit charge as J  σ˜f where σ is the conductivity [ 20, p. 285]. When the velocities of the charges are sufficiently small and the force is electromagnetic this becomes Ohm’s law J

f

 σE [ 20, p. 393], as can alternatively be seen in the diffusion model below.

When neglecting magnetic fields, [21, pp. 40,63] gives for dilute electrolyte media and chemical species S the Nernst-Planck equation with convection:

J

f

 ¸

kPS

z

k

F n

k

 ¸

kPS

z

k

F pz

k

|z

k

|

¹

u

k

c

k

E
D

k

∇ c

k

c

k

v q (4.11)

 ¸

kPS

|z

k

|eu

k

n

k

E
z

k

{|z

k

|

¹

u

k

k

B

T ∇ n

k

z

k

en

k

v N

_A

 6.02214129p27q  10

²³

[mol

¹

], Avogadro’s number

e  1.602176565p35q  10

¹⁹

[A s], elementary charge F  eN

A

k

_B

 1.3806488p13q  10

²³

[J K

¹

], Boltzmann’s constant R  k

B

N

A

E  electrical field [V m

¹

= J (A s)

¹

m

¹

]

n

_k

 sum of migration, diffusion and convection flux of species k [mol m

²

s

¹

] z

k

 valence of species k

u

k

 mobility of species k [m s

¹

(V m

¹

)

¹

] n

_k

 number density of species k [m

³

], c

_k

 N

A
¹

n

_k

T  absolute temperature [K]

D

_k

 u

k

RT |z

k

|

¹

F

¹

 u

k

k

_B

T |z

k

|

¹

e

¹

[m

²

s

¹

], diffusion (Einstein relation) where v is a velocity field representing the movement of the fluid.

For a linear, isotropic medium with charge carrier concentrations and mo- bilities independent of E (Ohmic conductors), the electrical conductivity can be defined as [21, p. 42]:

σ  ¸

kPS

|z

k

|F u

k

c

k

 ¸

kPS

|z

k

|eu

k

n

k

Then Ohm’s law is arrived at when v and ∇ n

k

are neglected.

4.2.5 Approximations of the Maxwell equations

By introducing typical dimensional scales and parameters, see Table 1, an esti- mate of the error of the quasistatic approximation can be obtained, as in [22].

Ohm’s law J

f

 σE will be used in the approximation.

It is possible to rewrite Maxwell’s equations using the gauge potentials and

the d’Alembertian partial derivative operator. If the partial time derivative

Table 1: Approximate constants [21, pp. 45,46] and dimensions for neurons µ

₀

 p4π10

⁷

 10

⁶

q V s A

¹

m

¹

L  1 cm  10

²

m

c  2.99792458  10

⁸

m s

¹

T  1 µs  10

⁶

s rµ

0



0

s

¹

 c

²

ñ 

0

 10

¹¹

A s V

¹

m

¹

σ  1 S m

¹

µ

r

 1, 

r

 10

²

component of the d’Alembertian is small, electroquasistatics (EQS) and mag- netoquasistatics (MQS) may be applicable. [22]

p
∇

²

µ

₀



₀

B

t²

q ÞÑ L

²

p
1 µ

0



₀

L

²

T

²

q

µ

₀



₀

L

²

T

²

 10

^{6
11
4 12}

 10

⁹

! 1, µL

²

T

²

 10

⁷

! 1 EQS : B

t

B  0 MQS : B

t

D  0, B

t

ρ

f

 B

t

∇  D  ∇ B

t

D  0 To use EQS, one wishes that the effect of B

t

B is small.

(4.5) : ∇  E 
B

t

B ÞÑ L

¹

E

F

 T

¹

µH

(4.8) : ∇  H  J

f

B

t

D  σE B

t

E ÞÑ L

¹

H  pσ T

¹

 qE

A

E

F

E

A

 L

²

T

¹

µ pσ T

¹

 q  10

^{4 6
6}

p10

⁰

10

^{6
11 2}

q  10

⁴

! 1 ñ EQS To use MQS, one wishes that the effect of B

t

D is small.

∇  H

D

 B

t

D  B

t

E ÞÑ L

¹

H

D

 T

¹

E

A

H

_D

H  LT

¹

E

_A

H ¤ LT

¹

 pE

A

E

_F

q

H  LT

¹

 pL

¹

pσ T

¹

 q

¹

H LT

¹

µH q H

 pT

¹

pσ T

¹

 q

¹

L

²

T

²

µ q ¤ pT

¹

σ

¹

L

²

T

²

µ q

 10

^{11 2}

p10

^6
0

10

^{4 12
6}

q  10

³

! 1 ñ MQS

So within the regions where the approximations are legitimate the approximate equations are:

∇  E  0, ∇ D  ∇ E  ρ

f

, ∇  H  0, ∇ H  J

f

Poincar´ es Lemma then gives that locally around ∇  E  0 the vector field can be expressed as E 
∇ φ [ 11, pp. 421-424].

A capacitor stores electromagnetic energy. Capacitors are defined under electrostatic conditions, that is when the time derivatives of the fields in the Maxwell equations are zero. The net charge of a conductor and the potential difference with respect to another conductor can be related through the capac- itance C as Q

a

 Crφ

a

φ

b

s [ 20, pp. 103,104]. In harmonic systems where the frequencies are low, such as when the size of the system is smaller than the free-space wavelength λ  c{f or when energy dissipation is mostly due to ohmic losses and not radiation, the field definitions of impedance can be traced to the usual definitions of capacitance and inductance [16, pp. 264-267].

Hopefully this motivates the use of electric circuit theory instead of the full

set of Maxwell’s equations.

5 Neuron model

In this section the chosen representation of the neuron geometry will be pre- sented. Furthermore the system of ordinary differential equations corresponding to the electrical circuit model of the neuron will be derived.

5.1 Geometry and physical parameters

As a simple model of a neuron’s geometry, approximately cover it with a tempo- rally fixed set of cylinders called compartments and index them by t1, . . . , Mu.

If the neuron’s volume is connected in the intersection or region between two neighbouring compartments, encode that data in a unidirectional adjacency ma- trix dA of the graph with the compartments as nodes, see Figure 4. In this simple approach the cylinders’ side strips are thought to represent the membrane, and the end plates represent either endings or cross-sections of the neuron.

Figure 4: Geometry of neuron and directed adjacency matrix. Removing the root node (1) gives an ordering on the cylinders. Associate with the neuron a specific membrane leakage conductivity g

L

(Siemens per square meter) and/or leak reversal potential E

L

(volt), specific membrane capacitance c

m

(Farad per square meter), cytoplasm conduc- tivity %

c

(ohm meter) and resting potential E

r

(volt). Then extrapolate this for each cylinder into membrane leak conductance G

^a_L

(Siemens) and membrane capacitance C

m^a

(Farad). TREES toolbox [23]

Denote the types of ion channel species to be modeled as T  tT

¹

, . . . , T

^c

u.

Assign for each cylinder the number of ion channels of each type, N

_T^a

: a P t1, . . . , Mu ^ T P T , and denote for each cylinder the number of open ion channels of each type by O

^a_T

: a P t1, . . . , Mu ^ T P T . It will be assumed that they have uniformly distributed and fixed positions on the membrane and that the number of channels do not change.

5.2 Cable equation

The total current flux over an oriented surface S is the current I  ³

S

n ˆ  JdA

g

. Let compartment Ω

a

have branches to compartments tΩ

b

u

bPB

. Compartment Ω

a

then has axial currents going through the surfaces A

b

 BΩ

a

X BΩ

b

, b P B and membrane currents going through the surfaces M  BΩ

a

z ”

bPB

BΩ

b

. The

continuity equation (4.2) gives:

¾

BΩa

J  ˆndA

g


¸

b

¾

Ab

J  ˆndA

g

¾

M

J  ˆndA

g



»

Ωa

B

t

ρdV

g

I

^BΩa


¸

b

I

_b^a

I

m^a


I

axial^a

I

m^a



»

Ωa

B

t

ρdV

g

(5.1)

in Figure 5 an example geometry is visualised.

Figure 5: A stylistic cross-section of a neuron with several connected compartments B  tb

1

, b

2

, b

3

, b

4

u. The arrows are ˆn at different positions.

5.2.1 Total axial current flux

Using the approximations we have E 
∇φ. Integrating the electric field over a curve γ : r0, 1s Ñ R

³

going from a to b we have from Stokes’s theorem:

»

γ

E  ˆtdL

g



»

γ

∇φ  ˆtdL

g


rφpbq
φpaqs  φpaq
φpbq  V

b^a

where the potential difference is denoted V

_b^a

. Assuming the material obeys Ohm’s law we have E  %J

f

, %  σ

¹

:

»

γ

%J

f

 ˆtdL

g

 V

b^a

For an object with piecewise constant cross-sectional area along ˆ t:

A : t ÞÑ A

^a

r0 ¤ t   l

^a

s A

^b

rl

^a

  t ¤ l

^a

l

^b

s

homogeneous resistivity % pxq  % and current Jpxq  I

^a

A

¹

ˆ t along the straight curve γ : r0, 1s Ñ R

³

: t ÞÑ p

¹₂

l

^a

t

¹₂

pl

^a

l

^b

qqˆt we get:

»

γ

% pxqJ  ˆtdL

g

 %I

^a

»

γ

A

¹

ˆ t  ˆtdL

g

 %I

^a

rpA

^a

q

¹

»

l^aˆt

1 2l^aˆt

dL

g

pA

^b

q

¹

»

_pl^{a lb}_2qˆt

l^aˆt

dL

g

s

 %I

^a

r 1 2

l

^a

A

^a

1 2

l

^b

A

^b

s  V

b^a

ñ I

^a

 r%r 1 2

l

^a

A

^a

1 2

l

^b

A

^b

ss

¹

V

_b^a

 G

^ab

V

_b^a

This is equivalent to the expression of the conductance between two connected cylinders in Dayan and Abbott [24, p. 219], the cross-sectional area of a cylinder being A

^k

 πpd

k

{2q

²

where d

k

is the diameter. See Figure 6 for an example.

Figure 6: Example of two connected compartments of cylindrical shape with cross- sectional areas A

k

, lengths l

k

and assumed current direction ˆ t.

Using the above expression as an approximation for the conductance the total axial current for compartment Ω

a

with respect to connected compartments tΩ

b

u

bPB

is:

I

_axial^a

 ¸

b

G

^a_b

V

_b^a

Adding and subtracting the extracellular potentials φ popkqq for the compartment and its branches allows the total axial current to be written as:

I

_axial^a

 ¸

b

G

^a_b

rV

b^a

s  ¸

b

G

^a_b

rV

o^apaq

V

o^bpbq

V

_o^o_pbq^paq

s,



 ¸

kPta,bu

φ popkqq
φpopkqq  0 

where o pkq denotes the extracellular region of compartment k. Assuming the outside is isopotential turns the current flux into:

I

_axial^a

 ¸

b

G

^a_b

rV

o^apaq

V

o^bpbq

s  ¸

b

G

^a_b

V

_o^a_paq

¸

b

G

^a_b

V

_o^b_pbq

, (5.2)

which is the form that will be used in this thesis.

5.2.2 Total membrane current flux

In order to have livable intracellular conditions, ions are transported over the membrane creating certain preferred concentrations of ions on the inside and outside. This creates a preferred equilibrium potential difference V

_m^a

 φpinside aq

φ poutside aq  φpaq
φpopaqq  V

_o^a_paq

over the membrane.

In the time stationary equilibrium of the continuity equation for the concen- tration c

k

of species k (the flux is related to (4.11), see [21] for details), neglecting the fluid velocity and production term, assuming a Boltzmann distribution of the species (such that they are non-interacting and in thermal equilibrium), the equilibrium potential E

^k

of that species is [21, p. 90][25, p. 13]

V

m

|

equilibrium

 E

^k

 RT z

k

F ln

 pc

k

q

outside a

pc

k

q

a

which is also approximately valid when the net flux is much less than the migra- tion and diffusion fluxes [21, p. 91]. More complicated models need to be used for ion channels that lets through many species of ions, such as the Goldman- Hodking Katz voltage equation [25, p. 344]. In this thesis the reversal potentials will be given as constants.

A linear current-voltage relation for ion channels is [25, p. 15]:

I

T

pV

m

, E

^T

, t q  G

T

rV

m

E

^T

s

where E

^T

is the reversal potential for the channel type and G

T

is a conductance representing the single channel conductances γ

^T

in parallel [25, p. 8] and is thus the sum of them. See Figure 7.

Figure 7: Ion channel model with single channel conductances γ

^T

, total conductance G

T

, reversal potential E

^T

and potentials φ depicted.

Assuming that closed channels have zero conductivity and that open chan-

nels have the same single channel conductance γ

^T

, this sum equals multiplying

γ

^T

by the number of open channels O

^a_T

:

I

_T^a

pV

m^a

, E

^T

, t q  G

^aT

rV

m^a

E

^T

s  O

^aT

γ

^T

rV

m^a

E

^T

s  O

^a_T

N

_T^a

N

_T^a

γ

^T

rV

m^a

E

^T

s In this thesis all other membrane currents in each compartment will be lumped together into a leak current I

_L^a

pV

m^a

, E

_L^a

, t q  G

^aL

rV

m^a

E

L^a

s with G

^aL

, E

_L^a

constants. If g

L

is given then G

^a_L

 A

^a

g

L

and E

_L^a

must be computed. As a simplification the end plates, which are not connected to other cylinders in the geometry, will not give any current contribution. So the total membrane current flux becomes:

I

_m^a

 I

L^a

¸

T

I

_T^a

 pG

^aL

¸

T

G

^a_T

qV

m^a

pG

^aL

E

_L^a

¸

T

G

^a_T

E

^T

q (5.3)

5.2.3 Capacitive effects

Compartment Ω

_a

with the membrane as an isolator gives an expression for the time-rate of change of the total charge, assuming the transmembrane potential difference V

_o^a_paq

is equal throughout the compartment and the capacitance is constant over time:

»

Ω_a

ρ

e

px, tqdV

g

 Q

^a

 C

m^a

V

_o^a_paq

»

Ω_a

B

t

ρ

_e

px, tqdV

g

 r d

dt C

_m^a

sV

o^apaq

C

_m^a

d

dt V

_o^a_paq

 C

m^a

d

dt V

_o^a_paq

(5.4) 5.2.4 Matrix model

Using the cable equation (5.1), the simplified axial current flux (5.2), the total membrane current flux (5.3), the capacitative effects (5.4) and a source term from outside the system I

_inj^a

gives:

C

_m^a

d

dt V

_o^a_paq

 I

inj^a

I

^BΩa

 I

inj^a

I

m^a

I

axial^a

 I

inj^a

pI

L^a

¸

T

I

_T^a

q
I

axial^a

 I

inj^a

pG

^aL

E

^a_L

¸

T

G

^a_T

E

^T

q
pG

^aL

¸

T

G

^a_T

qV

o^apaq

¸

b

G

^a_b

V

_o^a_paq

¸

b

G

^a_b

V

_o^b_pbq

which, with V

_m^a

 V

_o^a_paq

and the expressions for the G

^a_T

conductances can be written in matrix form as:



  C

_m¹

.. . C

_m^M



  d d dt



  V

_m¹

.. . V

_m^M



  



  I

_inj¹

.. . I

_inj^M



 



  G

¹_L

E

_L¹

.. . G

^M_L

E

_L^M



 



 

O

¹_T1

. . . O

¹_Tc

.. . .. . O

^M_T1

. . . O

^M_Tc



 



  γ

^T1

E

^T1

.. . γ

^Tc

E

^Tc



 



 



 

O

¹_T1

. . . O

¹_Tc

.. . .. . O

^M_T1

. . . O

^M_Tc



 



  γ

^T1

.. . γ

^Tc



 



  G

¹_L

.. . G

^M_L



  rG

^ab

s



  1 .. . 1



  

 d



  V

_m¹

.. . V

_m^M



  rG

^ab

s



  V

_m¹

.. . V

_m^M



 

ô C

m

d d

dt V

_m

 I

inj

G

_L

d E

L

O pγ d Eq
pOγ G

L

G1 q d V

m

GV

_m

ô diagpC

m

q d

dt V

_m

 I

inj

G

_L

d E

L

O pγ d Eq
diagpOγ G

L

G1 qV

m

GV

_m

ô C d

dt V

m

 Apt, Oq BpOqV

m

(5.5)

where d is the Hadamard product Appendix B.3. A sketch of the circuit around a compartment a with branches to tb

1

, b

2

, b

3

u can be seen in Figure 8.

Figure 8: Cable equation circuit with isopotential outside (ground)

The resting potential E

_r

is the value of V

_m^a

such that dV

_m^a

{dt  0 when I

_inj^a

 0 for all a. Using this we get an equation relating G

^a_L

, E

_L^a

and E

_r

:

E

L



 1 G

^a_L



d rOpγ d Eq
pOγ G

L

G1 q d p1E

r

q Gp1E

r

qs

G

L



 1

E

_L^a

E

r



d rOpγ d Eq
pOγ G1q d p1E

r

q Gp1E

r

qs

where the number of open channels O

^a_T

are taken as their steady state values for the potential E

_r

. An alternative is to give both parameters G

^a_L

, E

_L^a

and checking that dV

_m^a

{dt  0.

6 Ion channel model

This section will present how the dynamics of channel opening and closing will

be modeled in this thesis. Since it is common for neurophysiologists to create

models of gating-particle types, which can be understood as finite state-space

continuous-time homogenous Markov chains [4, pp. 308-309], this approach will

be used here. A very detailed simulation algorithm is that of [26, 27], but it

is deemed too inefficient for the scope of this thesis. A faster approach stems

from the tau-leap approximation, which will be described later in this section.

Another way of increasing the efficiency is approximating the Markov chains as stochastic differential equation models, for example as in [28]. Another option is converting it into a system of ordinary differential equations.

6.1 Gating-particle models

In this thesis, each ion channel species is assumed to only have one open, con- ducting state. Any closed states are non-conducting. In a patch of membrane a with area A

^a

, let there be a voltage gated ion channel species T present. De- note by

^O_N^aTa

T

the ratio of the number of open T -channels over the total number of T -channels in the patch, %

^a_T

the density of that species in the patch, and γ

^T

the single channel conductance of one open T -channel, assumed equal for all channels. Hodgkin and Huxley theorised that there were some kind of gating- particles sitting on the membrane, deciding whether the ion channel would open or close depending on the voltage. An interpretation of them are as charged parts of a voltage sensitive subunit of the ion channel [25, p. 55]. Following [25, p. 46], the model can be understood as follows. Introduce #T types of p

_Ti

identical gating particles o

^a_T

i

, each representing the probability that its associ- ated particle is in the correct position to set up an open channel. Let V

_m^a

be the potential difference over the membrane, that is the potential of the inside of the membrane patch minus the potential of the outside of the membrane patch. The gating particles make transitions between the permissive/open and non-permissive/closed forms according to a first order reaction with voltage and time dependent rates α

^a_T

i

and β

_T^a

i

. Closed particles p1
o

^aT_i

q become open with rate α

^a_Ti

and open particles o

^a_Ti

become closed with rate β

_T^a_i

. In the form of a system of ordinary differential equations this becomes (see for example [4, p. 306]):

$ ' ' '

&

' ' ' %

I

_T^a



^O_N^aTa

T

%

^a_T

A

^a

γ

^T

rV

m^a

E

^T

s

O^a_T N_T^a



^#T

±

i1

ro

^aT_i

s

^pTi

do^a_Ti

dt

 α

^aT_i

 p1
o

^aT_i

q
β

^aT_i

 o

^aT_i

[(milli seconds)

¹

],

(6.1)

where α

^a_T

i

 α

Ti

pV

m^a

E

r

q: R [milli volt] Ñ R [(milli seconds)

¹

] β

^a_Ti

 α

T_i

pV

m^a

E

r

q: R [milli volt] Ñ R [(milli seconds)

¹

]

are the rate functions. By comparing the rates of [30, p. 346], [27, p. 70], and [31, p. 3014] it is apparent they are equal up to choices of E

_r

. The rates are probably fitted around the resting potential, thus needing it as an extra parameter besides V

m

. Often the rates are composed of exponentials, which need the arguments to be dimensionless. Therefore the factors scaling the potential inside the exponential are probably dimensional, requiring the input to be in mV . Similarly the output is usually in rmss

¹

. Gating particle models can (at the present date) be found online at [29].

Using combinatorial arguments [24, pp. 176-177] shows the equivalence of

the lumped Markov chain based on the rates given in the particle-model of the

single delayed-rectifier K channel and the corresponding system of ordinary differential equations. Generalising from [30, p. 384] and [25, pp. 485–490] the derivation of the lumped Markov chain can in some detail be done as follows, but first an example. The ODE-system for the sodium m

³

h

¹

channel is:

$ ' ' '

&

' ' ' %

I

N a



^O_N^{N a}_{N a}

%

N a

Aγ

^{N a}

rV

m

E

^{N a}

s

ON a

N_{N a}

 ro

N a₁

s

^{pN a1}

ro

N a₂

s

^{pN a2}

 m

³

h

¹

dm

dt

 α

m

 p1
mq
β

m

 m [(ms)

¹

]

dh

dt

 α

h

 p1
hq
β

h

 h [(ms)

¹

]

The number of particle types is #N a  2. Denote the first particle type by o

N a₁

 m, its power by p

N a₁

 p

m

 3 and its rates by α

m

and β

m

. Similarly let the second particle type be denoted o

N a₂

 h with associated power p

N a₂

 p

h

 1 and rates α

h

and β

h

. The state space of the channel is composed of two tuples of length p

m

and p

h

: ppm

1

, m

2

, m

3

q, phqq, where the elements are binary and represent whether a particular particle is open or closed. Create a Markov chain by connecting the states which only differ by one element being open or closed and associate with those transitions a rate that is α

_{N ai}

if it is an opening transition and β

_{N ai}

if it is a closing transition. Then make an equivalence relation that groups the tuples according to how many of the elements are open, so for example ppO, C, Cq, pOqq  ppC, O, Cq, pOqq  ppC, C, Oq, pOqq P pr1s, r1sq P S{ . Then create a lumped version of the Markov chain by using (3.4). In Figure 9 the procedure is shown graphically. As for the general case, first assume the transition rates between the open and closed states of the particles are instant. Also assume the state of the pore changes by one particle’s state at a time. For a channel with #T particle types and p

T_i

identical gating particles, let the state space S

0

be the Cartesian product of the family of p

T_i

- ary Cartesian products over tC, Ou indexed by T

¹

 tT

1

, . . . , T

#T

u; let the state space be S

0

 ±

T_iPT¹

±

k1...,pT

tC, Ou. The transition rates for s Ñ ˆs, ˆs  s are

Q

sˆs



$ '

&

' %

α

T_i

if D!pT

i

, k q: s

_pTi,kq

 ˆs

_pTi,kq

^ s

_pTi,kq

 C ^ ˆs

_pTi,kq

 O β

T_i

if D!pT

i

, k q: s

_pTi,kq

 ˆs

_pTi,kq

^ s

_pTi,kq

 O ^ ˆs

_pTi,kq

 C 0 otherwise

where the state element s

_pTi,kq

is the associated element of the projection map from S

₀

to the term indexed by particle type T

_i

in the family and k in the p

_Ti

- ary particle tuple. Each state element then represents whether the k

^th

particle of type T

_i

is open or closed.

Define an equivalence relation „ S

0

 S

0

by

s  ˆs ô @T

i

P T

¹

: |tk P t1, . . . , p

T_i

u: s

_pTi,kq

 Ou|  |tk P t1, . . . , p

T_i

u: ˆs

_pTi,kq

 Ou|

which partitions the state space into classes rss P S

0

{  with equal number of open state elements for all subunit types. Application of (3.4) gives the new transition rates between the classes, since the transitions only happen in steps of one opening or one closing in one type of particle, and those rates are equal.

In [32] a fast algorithm for lumping is described, but since the state spaces of

gating particle models are usually small this is probably unnecessary.

Figure 9: Sodium channel m

³

h

¹

. The opening rates α

T

are associated with the edges

ending with Ÿ and the closing rates β

T

with the edges ending with ž, for both particle

types, T P tm, hu. 0 means closed and 1 means open. The different particle types are

separated by commas in the tuples. The equivalence class lumping is seen below, the

edge labels are the scale factors associated with the sum of the rates from application

of (3.4). For example the rate from (000,1) to (001,1) is α

m

, from ([0],[1]) to ([1],[1])

is 3α

m

and the rates from ([1],[0]) to ([1],[1]) and (001,0) to (001,1) are both α

h

.

6.2 Many Markov chains together

Let there be N

_T^a

independent and identical Markov chains with infinitesimal generator Q

^{pT q}_ij

and state space S

^{pT q}

. Count the number of Markov chains that in the time frame rt, t τs are in state i and denote that by n

^{pa,T q}i

. Then the Markov chains are distributed over the states as the sequence

n

^{pa,T q}

 pn

^{pa,T q}i

q

iPS^{pT q}

: ¸

iPS^{pT q}

n

^{pa,T q}_i

 N

T^a

^ n

^{pa,T q}i

P N

^¥0

(6.2)

and have holding times from state i P S

^{pT q}

to j P S

^{pT q}

ztiu equal to exponential random variables T

_ijk^{pnpa,T qq}

 ExppQ

^{pT q}_ij

q, k P t1, . . . , n

^{pa,T q}_i

u. Taking the infi- mum of T

_ijk^{pnpa,T qq}

over pi, j  iq and k, the time of the first transition in all of the N

_T^a

Markov chains becomes

T  inf

pi,jiq

inf

k

T

_ijk^{pnpa,T qq}

 Exp

 ¸

i

¸

j

rj  is ¸

k

Q

^{pT q}_ij



 Exp

 ¸

i

¸

j

rj  isn

^{pa,T q}i

Q

^{pT q}_ij



(6.3)

as seen by applying (3.1) two times, and according to (3.2) the pi, j  iq that is chosen is independent of (6.3) and distributed as:

PpT  inf

k

T

_ijk^{pnpa,T qq}

|n

^{pa,T q}

, t q  rj  isn

^{pa,T q}i

Q

^{pT q}_ij

°

i

°

j

rj  isn

^{pa,T q}i

Q

^{pT q}_ij

(6.4)

assuming the denominator is nonzero. The probability density of (6.3) and the equation (6.4) are called p

1

and p

2

in [33, p. 1717].

Let there be an indexed set tN

a

u

aPt1,...,Mu

of distributions of different Markov chains n

^pa,Tbq

, T

^b

P T , all independent. Then according to ( 3.1) the time of the first state transition in all of the Markov chains is:

T  inf

pn^{pa,T bq}P—

aNa,i,jiq

inf

k

T

_i,j,k^{pnpa,T bqq}

 Exp



 ¸

n^{pa,T bq}P—

aNa

¸

i

¸

j

rj  is ¸

k

Q

^pT_ij^bq



 Exp



 ¸

n^{pa,T bq}P—

aNa

¸

i

¸

j

rj  isn

^pa,Ti ^bq

Q

^pT_ij^bq



(6.5)

where —

is the disjoint union

¹

. Furthermore the probability of the first transi- tion according to (3.2) is, assuming nonzero denominator, distributed as:

PpT  inf

k

T

_i,j,k^{pnpa,T bqq}

|tN

a

u, tq  rj  isn

^pa,T_i ^bq

Q

^pT_ij^bq

°

n^{pa,T bq}P—

aN_a

°

i

°

j

rj  isn

^pa,Ti ^bq

Q

^pT_ij^bq

(6.6)

1This means identical elements in different Naare still in the union as separate elements, indexed here by a.