13.8. Stochastic Calcium Burst model with GHK currents

This chapter is based on the publication “Anwar H, Hepburn I, Nedelescu H, Chen W, De Schutter E (2013) Stochastic Calcium Mechanisms Cause Dendritic Calcium Spike Variability. The Journal of Neuroscience, 33(40): 15848-15867”. The model file is StochasticCaburst.py which is available at http://senselab.med.yale.edu/modeldb/ShowModel.asp?model=150635 and should be run within the complete modelling package because it depends on other files such as /extra/constants.py.

The simulation script described in this chapter is also available at STEPS_Example repository.

This chapter builds on previous chapters by simulating a model that includes both a reaction-diffusion component as well as electrical excitability. As described in 1, the two are closely coupled and this model contains ion channels where activation is both voltage-dependent and calcium-dependent. In addition, the calcium ions form an important part of the current across the membrane, further coupling the reaction-diffusion component with the electrical excitability. This chapter introduces an important new object: the GHK Current object, which is described in some detail in section P-type Calcium channel.

As in previous chapters we will go through the script, looking in some depth at new concepts, but only brief explanations will be offered of things that have been described in previous chapters.

13.8.1. Modelling solution

At the top of the script, as usual, we import some modules, including STEPS modules and some self-written helper modules. This import includes the extra/constants.py file, which includes all important parameters for the module such as physical constants, membrane properties, kinetic properties of the channels, initial conditions and so on:

import math
import time
from random import *
import sys
import os

import steps.model as smodel
import steps.geom as sgeom
import steps.rng as srng
import steps.utilities.meshio as meshio
import steps.solver as ssolver

import meshes.gettets as gettets
from extra.constants import *

13.8.1.1. extra/constants.py file

Although not part of this model script, for easy reference to many important parameters used in the model this is the contents of the imported extra.constants module. What these parameters are and how they are used will become apparent as we go through the model. All unit are S.I. units, with the exception of concentration where units are molar:

# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

import math

# # # # # # # # # # # # # # # # SIMULATION CONTROLS # # # # # # # # # # # # #

EF_DT = 2.0e-5          # The EField dt
NTIMEPOINTS =  25001

TIMECONVERTER =  2.0e-5

NITER = 1

############################ PARAMETERS ################################

init_pot = -60e-3

TEMPERATURE = 34.0

Q10 = 3

# Faraday constant: unit of FARADAY is C/mol
# Source: http://physics.nist.gov/cgi-bin/cuu/Value?f 24/2/2012
FARADAY = 96485.3365

# Molar Gas Constant: unit of R is J/mol K
# Source: http://physics.nist.gov/cgi-bin/cuu/Value?r 24/2/2012
R = 8.3144621

# Avogadro constant: unit of AVOGADRO is /mol
# Source: http://physics.nist.gov/cgi-bin/cuu/Value?na 24/2/2012
AVOGADRO = 6.02214129e23

# Elementary charge: unit of E_CHARGE is C
# Source: http://physics.nist.gov/cgi-bin/cuu/Value?e 24/2/2012
E_CHARGE = 1.602176565e-19


#FOR MSLO, THERE IS A NEW VALUE FOR Qt wrt to 25 degC
Qt = math.pow(Q10, ((TEMPERATURE-23)/10))
Qt_mslo = math.pow(Q10, ((TEMPERATURE-25)/10))

########## BULK RESISTIVITY ##########

Ra = 235.7*1.0e-2

########## MEMBRANE CAPACITANCE ##########

memb_capac = 1.5e-2


########## CaP channels density & permiability per channel ##########

# CaP_P is permiability per channel (m3/s)
# CaP_ro is channel/surface area (/m2)
# P in Ca Dynamics model is 0.95e-4 cm/s --> 0.95e-6 m/s

CaP_P = 2.5e-20
CaP_ro = 3.8e13

##########CaP channel parameters ####################

#Units (mV)
vhalfm = -29.458
cvm = 8.429

def minf_cap(V):
    #Units (mV)
    vhalfm = -29.458
    cvm = 8.429
    vshift = 0.0

    return (1.0/(1.0 + math.exp(-(V-vhalfm-vshift)/cvm)))

def tau_cap(V):
    vshift = 0.0
    if (V-vshift) >= -40:
        return (0.2702 + 1.1622 * math.exp(-(V+26.798-vshift)*(V+26.798-vshift)/164.19))
    else:
        return (0.6923 * math.exp((V-vshift)/1089.372))

def alpha_cap(V):
    return (minf_cap(V)/tau_cap(V))

def beta_cap(V):
    return ((1.0-minf_cap(V))/tau_cap(V))


## Intitial conditions

CaP_m0_p = 0.92402
CaP_m1_p = 0.073988
CaP_m2_p = 0.0019748
CaP_m3_p = 1.7569e-05


########## CaT channels density & permiability per channel ##########

# CaT_P is permiability per channel (m3/s)
# CaT_ro is channel/surface area (/m2)
# P in Ca Dynamics model is 6.2e-6 cm/s -->6.2e-8 m/s

CaT_P = 1.65e-20
CaT_ro = 3.7576e12


def minf_cat(V):
    #Units (mV)
    vhalfm = -52.0
    cvm = -5.0
    vshift = 0.0

    return (1.0/(1.0 + math.exp((V-vhalfm-vshift)/cvm)))

def taum_cat(V):
    vshift = 0.0
    if V > -90.0:
        return (1.0 + 1.0 / (math.exp((V+40.0-vshift)/9.0) + math.exp(-(V+102.0-vshift)/18.0)))
    else:
        return 1.0

def hinf_cat(V):
    vhalfh = -72.0
    cvh = 7.0
    vshift = 0.0
    return (1.0/(1.0 + math.exp((V-vhalfh-vshift)/cvh)))

def tauh_cat(V):
    vshift = 0.0
    return (15.0 + 1.0 / (math.exp((V+32.0-vshift)/7.0)))

def alpham_cat(V):
    return (minf_cat(V)/taum_cat(V))

def betam_cat(V):
    return ((1-minf_cat(V))/taum_cat(V))

def alphah_cat(V):
    return (hinf_cat(V)/tauh_cat(V))

def betah_cat(V):
    return ((1-hinf_cat(V))/tauh_cat(V))

## Initial conditions

CaT_m0h0_p = 0.58661
CaT_m1h0_p = 0.23687
CaT_m2h0_p = 0.023912
CaT_m0h1_p = 0.10564
CaT_m1h1_p = 0.042658
CaT_m2h1_p = 0.0043063

########## BK channels density & conductance per channel ##########

# Total conductance = BK_G (conductance/channel) * BK_ro (channel/surface area)
# BK in Ca Dynamics model is 4.25e-2 S/cm2 --> 4.25e2 S/m2


BK_G = 2.1e-10
BK_ro = 2.0238e12
BK_rev = -77e-3

######### BK channel parameters ######################

#Units (1)
Qo = 0.73
Qc = -0.67

#Units (/s)
pf0 = 2.39
pf1 = 5.4918
pf2 = 24.6205
pf3 = 142.4546
pf4 = 211.0220

pb0 = 3936
pb1 = 687.3251
pb2 = 234.5875
pb3 = 103.2204
pb4 = 11.6581

#Units(/M)
k1 = 1.0e6

#Units(/s)
onoffrate = 1.0e3

L0 = 1806

#Units (M)
Kc = 8.63e-6
Ko = 0.6563e-6


c_01 = 4.*k1*onoffrate*Qt_mslo
c_12 = 3.*k1*onoffrate*Qt_mslo
c_23 = 2.*k1*onoffrate*Qt_mslo
c_34 = 1.*k1*onoffrate*Qt_mslo
o_01 = 4.*k1*onoffrate*Qt_mslo
o_12 = 3.*k1*onoffrate*Qt_mslo
o_23 = 2.*k1*onoffrate*Qt_mslo
o_34 = 1.*k1*onoffrate*Qt_mslo

c_10 = 1.*Kc*k1*onoffrate*Qt_mslo
c_21 = 2.*Kc*k1*onoffrate*Qt_mslo
c_32 = 3.*Kc*k1*onoffrate*Qt_mslo
c_43 = 4.*Kc*k1*onoffrate*Qt_mslo
o_10 = 1.*Ko*k1*onoffrate*Qt_mslo
o_21 = 2.*Ko*k1*onoffrate*Qt_mslo
o_32 = 3.*Ko*k1*onoffrate*Qt_mslo
o_43 = 4.*Ko*k1*onoffrate*Qt_mslo


f_0 = lambda mV: pf0*Qt_mslo*(math.exp((Qo* FARADAY* mV) / (R* (TEMPERATURE + 273.15))))
f_1 = lambda mV: pf1*Qt_mslo*(math.exp((Qo* FARADAY* mV) / (R* (TEMPERATURE + 273.15))))
f_2 = lambda mV: pf2*Qt_mslo*(math.exp((Qo* FARADAY* mV) / (R* (TEMPERATURE + 273.15))))
f_3 = lambda mV: pf3*Qt_mslo*(math.exp((Qo* FARADAY* mV) / (R* (TEMPERATURE + 273.15))))
f_4 = lambda mV: pf4*Qt_mslo*(math.exp((Qo* FARADAY* mV) / (R* (TEMPERATURE + 273.15))))

b_0 = lambda mV: pb0*Qt_mslo*(math.exp((Qc* FARADAY* mV) / (R* (TEMPERATURE + 273.15))))
b_1 = lambda mV: pb1*Qt_mslo*(math.exp((Qc* FARADAY* mV) / (R* (TEMPERATURE + 273.15))))
b_2 = lambda mV: pb2*Qt_mslo*(math.exp((Qc* FARADAY* mV) / (R* (TEMPERATURE + 273.15))))
b_3 = lambda mV: pb3*Qt_mslo*(math.exp((Qc* FARADAY* mV) / (R* (TEMPERATURE + 273.15))))
b_4 = lambda mV: pb4*Qt_mslo*(math.exp((Qc* FARADAY* mV) / (R* (TEMPERATURE + 273.15))))


# Initial conditions
BK_C0_p= 0.99997
BK_C1_p= 4.3619e-07
BK_C2_p= 4.1713e-09
BK_C3_p= 4.4449e-11
BK_C4_p= 6.3132e-14

BK_O0_p= 2.5202e-05
BK_O1_p= 1.1765e-06
BK_O2_p= 6.6148e-08
BK_O3_p= 2.4392e-09
BK_O4_p= 4.0981e-11

########## SK channel density & conductance per channel #############

# Total conductance = SK_G (conductance/channel) * SK_ro (channel/surface area)
# SK in Ca Dynamics model is 3.1e-4 S/cm2 --> 3.1 S/m2


SK_G = 1.0e-11
SK_ro = 31.0e10

SK_rev = -77e-3

######### SK channel parameters ###################

#Units (/s)
invc1 = 80
invc2 = 80
invc3 = 200

invo1 = 1000
invo2 = 100

diro1 = 160
diro2 = 1200

#Units ( /s M)

dirc2 = 200e6
dirc3 = 160e6
dirc4 = 80e6

invc1_t = invc1*Qt
invc2_t = invc2*Qt
invc3_t = invc3*Qt

invo1_t = invo1*Qt
invo2_t = invo2*Qt

diro1_t = diro1*Qt
diro2_t = diro2*Qt

dirc2_t = dirc2*Qt/3.0
dirc3_t = dirc3*Qt/3.0
dirc4_t = dirc4*Qt/3.0


# Intital conditions
SK_C1_p= 0.96256
SK_C2_p= 0.036096
SK_C3_p= 0.0010829
SK_C4_p= 6.4973e-06

SK_O1_p= 0.00017326
SK_O2_p= 7.7967e-05


######### leak current channel density & conductance per channel ########
# Total conductance = 1e-6 S/cm2 --> 1e-2 S/m2

L_G = 4.0e-14
L_ro = 25.0e10

L_rev = -61e-3


######### Pump parameters ###################

P_f_kcst = 3e9
P_b_kcst = 1.75e4
P_k_kcst = 7.255e4


############################CALCIUM BUFFERING MODEL################################

########## Ca concentrations #########

Ca_oconc = 2e-3
Ca_iconc = 45e-9

########## Mg concentrations #########

Mg_conc = 590e-6

########## Buffer concentrations #############

iCBsf_conc = 27.704e-6
iCBCaf_conc = 2.6372e-6
iCBsCa_conc= 1.5148e-6
iCBCaCa_conc= 0.14420e-6

CBsf_conc= 110.82e-6
CBCaf_conc= 10.549e-6
CBsCa_conc= 6.0595e-6
CBCaCa_conc= 0.57682e-6

PV_conc= 3.2066e-6
PVCa_conc= 16.252e-6
PVMg_conc= 60.541e-6

# Diffusion constant of Calcium
DCST = 0.223e-9
# Diffusion constant of Calbindin (CB)
DCB = 0.028e-9
# Diffusion constant of Parvalbumin (PV)
DPV = 0.043e-9

#iCBsf-fast
iCBsf1_f_kcst = 4.35e7
iCBsf1_b_kcst = 35.8

#iCBsCa
iCBsCa_f_kcst = 0.55e7
iCBsCa_b_kcst = 2.6

#iCBsf_slow
iCBsf2_f_kcst = 0.55e7
iCBsf2_b_kcst = 2.6

#iCBCaf
iCBCaf_f_kcst = 4.35e7
iCBCaf_b_kcst = 35.8

#CBsf-fast
CBsf1_f_kcst = 4.35e7
CBsf1_b_kcst = 35.8

#CBsCa
CBsCa_f_kcst = 0.55e7
CBsCa_b_kcst = 2.6

#CBsf_slow
CBsf2_f_kcst = 0.55e7
CBsf2_b_kcst = 2.6

#CBCaf
CBCaf_f_kcst = 4.35e7
CBCaf_b_kcst = 35.8

#PVca
PVca_f_kcst = 10.7e7
PVca_b_kcst = 0.95

#PVmg
PVmg_f_kcst = 0.8e6
PVmg_b_kcst = 25

# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

13.8.1.2. Command line execution

Back to the StochasticCaburst.py script, and next we take a slightly new approach to previous models, in which we utilise command line arguments (sys.argv). The ‘zeroth’ argument (sys.argv[0]) is the script pathname, then other arguments are command line arguments which are read as strings. In this sense it is intended that the StochasticCaburst.py script is run from the command line with a statement such as:

$ python StochasticCaburst.py Cylinder2_dia2um_L10um_outer0_3um_0.3shell_0.3size_19156tets_adaptive.inp ~/stochcasims/ 1

and, in the present form, this script can not be run interactively as previous examples can. In the script, we read command line arguments thusly:

meshfile_ab, root, iter_n = sys.argv[1], sys.argv[2], sys.argv[3]

Looking at the above example, sys.argv[1] would be the string ‘Cylinder2_dia2um_L10um_outer0_3um_0.3shell_0.3size_19156tets_adaptive.inp’ and will be stored as variable meshfile_ab, sys.argv[2] will be the string ‘~/stochcasims/’ and stored as variable root which, as we will see, defines where to store simulation output, and sys.argv[3] will be the string ‘1’ and stored as variable iter_n which is involved in random number initialisation (after conversion to an integer) and data storage. We will look at these variables later in the script as they are used.

The last thing to be done before moving onto the biochemical model description is to set a flag to detect whether we are dealing with the 160um mesh (available meshes that are intended for this simulation are cylinders all of diameter 2um and lengths of 10um, 20um, 40um, 80um and 160um). This flag is necessary because the 160um mesh has slightly different properties than the other meshes in that it contains no extracellular tetrahedral compartment, and so there will be some different behaviours throughout the script depending on whether we are using the 160um mesh or not:

if meshfile_ab == 'Cylinder2_dia2um_L160um_outer0_0.3shell_0.3size_279152tets_adaptive.inp': cyl160=True
else: cyl160=False

13.8.2. Model specification

Since this is a relatively large model we will split its description up into two sections, after first creating the parent model container object, one volume system and one surface system:

mdl = smodel.Model()

# Vol/surface systems
vsys = smodel.Volsys('vsys', mdl)
ssys = smodel.Surfsys('ssys', mdl)

13.8.2.1. Calcium dynamics

The following lines of code describe the calcium and calcium buffer reactions and diffusion. Since these are ‘ordinary’ dynamics with no voltage-dependence we will not look look at this part in detail. A more detailed explanation is offered in 1 and 2. Most parameters come from the extra/constants.py file:

# Calcium
Ca = smodel.Spec('Ca', mdl)
Ca.setValence(2)

# Pump
Pump = smodel.Spec('Pump', mdl)
# CaPump
CaPump = smodel.Spec('CaPump', mdl)

# iCBsf
iCBsf = smodel.Spec('iCBsf', mdl)
# iCBsCa
iCBsCa = smodel.Spec('iCBsCa', mdl)
# iCBCaf
iCBCaf = smodel.Spec('iCBCaf', mdl)
# iCBCaCa
iCBCaCa = smodel.Spec('iCBCaCa', mdl)

# CBsf
CBsf = smodel.Spec('CBsf', mdl)
# CBsCa
CBsCa = smodel.Spec('CBsCa', mdl)
# CBCaf
CBCaf = smodel.Spec('CBCaf', mdl)
# CBCaCa
CBCaCa = smodel.Spec('CBCaCa', mdl)

# PV
PV = smodel.Spec('PV', mdl)
# PVMg
PVMg = smodel.Spec('PVMg', mdl)
# PVCa
PVCa = smodel.Spec('PVCa', mdl)
# Mg
Mg = smodel.Spec('Mg', mdl)

diff_Ca = smodel.Diff('diff_Ca', vsys, Ca)
diff_Ca.setDcst(DCST)
diff_CBsf = smodel.Diff('diff_CBsf', vsys, CBsf)
diff_CBsf.setDcst(DCB)
diff_CBsCa = smodel.Diff('diff_CBsCa', vsys, CBsCa)
diff_CBsCa.setDcst(DCB)
diff_CBCaf = smodel.Diff('diff_CBCaf', vsys, CBCaf)
diff_CBCaf.setDcst(DCB)
diff_CBCaCa = smodel.Diff('diff_CBCaCa', vsys, CBCaCa)
diff_CBCaCa.setDcst(DCB)
diff_PV = smodel.Diff('diff_PV', vsys, PV)
diff_PV.setDcst(DPV)
diff_PVCa = smodel.Diff('diff_PVCa', vsys, PVCa)
diff_PVCa.setDcst(DPV)
diff_PVMg = smodel.Diff('diff_PVMg', vsys, PVMg)
diff_PVMg.setDcst(DPV)

#Pump
PumpD_f = smodel.SReac('PumpD_f', ssys, ilhs=[Ca], slhs=[Pump], srhs=[CaPump])
PumpD_f.setKcst(P_f_kcst)

PumpD_b = smodel.SReac('PumpD_b', ssys, slhs=[CaPump], irhs=[Ca], srhs=[Pump])
PumpD_b.setKcst(P_b_kcst)

PumpD_k = smodel.SReac('PumpD_k', ssys, slhs=[CaPump], srhs=[Pump])
PumpD_k.setKcst(P_k_kcst)

#iCBsf-fast
iCBsf1_f = smodel.Reac('iCBsf1_f', vsys, lhs=[Ca,iCBsf], rhs=[iCBsCa], kcst = iCBsf1_f_kcst)
iCBsf1_b = smodel.Reac('iCBsf1_b', vsys, lhs=[iCBsCa], rhs=[Ca,iCBsf], kcst = iCBsf1_b_kcst)

#iCBsCa
iCBsCa_f = smodel.Reac('iCBsCa_f', vsys, lhs=[Ca,iCBsCa], rhs=[iCBCaCa], kcst = iCBsCa_f_kcst)
iCBsCa_b = smodel.Reac('iCBsCa_b', vsys, lhs=[iCBCaCa], rhs=[Ca,iCBsCa], kcst = iCBsCa_b_kcst)

#iCBsf_slow
iCBsf2_f = smodel.Reac('iCBsf2_f', vsys, lhs=[Ca,iCBsf], rhs=[iCBCaf], kcst = iCBsf2_f_kcst)
iCBsf2_b = smodel.Reac('iCBsf2_b', vsys, lhs=[iCBCaf], rhs=[Ca,iCBsf], kcst = iCBsf2_b_kcst)

#iCBCaf
iCBCaf_f = smodel.Reac('iCBCaf_f', vsys, lhs=[Ca,iCBCaf], rhs=[iCBCaCa], kcst = iCBCaf_f_kcst)
iCBCaf_b = smodel.Reac('iCBCaf_b', vsys, lhs=[iCBCaCa], rhs=[Ca,iCBCaf], kcst = iCBCaf_b_kcst)

#CBsf-fast
CBsf1_f = smodel.Reac('CBsf1_f', vsys, lhs=[Ca,CBsf], rhs=[CBsCa], kcst = CBsf1_f_kcst)
CBsf1_b = smodel.Reac('CBsf1_b', vsys, lhs=[CBsCa], rhs=[Ca,CBsf], kcst = CBsf1_b_kcst)

#CBsCa
CBsCa_f = smodel.Reac('CBsCa_f', vsys, lhs=[Ca,CBsCa], rhs=[CBCaCa], kcst = CBsCa_f_kcst)
CBsCa_b = smodel.Reac('CBsCa_b', vsys, lhs=[CBCaCa], rhs=[Ca,CBsCa], kcst = CBsCa_b_kcst)

#CBsf_slow
CBsf2_f = smodel.Reac('CBsf2_f', vsys, lhs=[Ca,CBsf], rhs=[CBCaf], kcst = CBsf2_f_kcst)
CBsf2_b = smodel.Reac('CBsf2_b', vsys, lhs=[CBCaf], rhs=[Ca,CBsf], kcst = CBsf2_b_kcst)

#CBCaf
CBCaf_f = smodel.Reac('CBCaf_f', vsys, lhs=[Ca,CBCaf], rhs=[CBCaCa], kcst = CBCaf_f_kcst)
CBCaf_b = smodel.Reac('CBCaf_b', vsys, lhs=[CBCaCa], rhs=[Ca,CBCaf], kcst = CBCaf_b_kcst)

#PVca
PVca_f = smodel.Reac('PVca_f', vsys, lhs=[Ca,PV], rhs=[PVCa], kcst = PVca_f_kcst)
PVca_b = smodel.Reac('PVca_b', vsys, lhs=[PVCa], rhs=[Ca,PV], kcst = PVca_b_kcst)

#PVmg
PVmg_f = smodel.Reac('PVmg_f', vsys, lhs=[Mg,PV], rhs=[PVMg], kcst = PVmg_f_kcst)
PVmg_b = smodel.Reac('PVmg_b', vsys, lhs=[PVMg], rhs=[Mg,PV], kcst = PVmg_b_kcst)

13.8.2.2. P-type Calcium channel

The P-type calcium channel is a different type of ion channel to those we have seen before. In previous chapters we saw Hodgkin-Huxley sodium and potassium channels that conducted an Ohmic current. The sodium and potassium ions in that situation were not explicitly simulated, which was reasonable because those ions were not involved in other processes we were interested in, and we could assume their concentrations inside and outside the cell were not altered significantly during their conduction. However, with calcium we need a different approach. Here calcium is involved in intracellular processes such as potassium channel-activation (as we will see), buffering and diffusion, and so we must simulate the influx of calcium through these P-type channels. Furthermore, the Ohmic approximation is no longer sufficient for our purposes. The large differences between intracellular and extracellular concentration along with large changes in intracellular concentration mean that, in effect, channel conductance has some voltage and concentration dependence and is described much better by the GHK flux equation. The GHK flux equation itself is derived under certain simplifying assumptions that are good approximations for many ion channels, specifically those where channel occupancy and competition are negligible. Please see 3 for further discussion on the use of the GHK flux equation and the behaviour of the GHK current object in STEPS. It is worth noting that use of the GHK flux equation means that (instead of conductance) we must specify the channel’s permeability, which can be more difficult to parameterize.

The P-type calcium channel kinetics are described in detail in 1. To create the channel we first describe the channel states:

CaPchan = smodel.Chan('CaPchan', mdl)

CaP_m0 = smodel.ChanState('CaP_m0', mdl, CaPchan)
CaP_m1 = smodel.ChanState('CaP_m1', mdl, CaPchan)
CaP_m2 = smodel.ChanState('CaP_m2', mdl, CaPchan)
CaP_m3 = smodel.ChanState('CaP_m3', mdl, CaPchan)

and the voltage-dependent kinetics. Remember for each of these discrete channels this voltage will be read from the local voltage across the membrane triangle where the channel resides:

CaPm0m1 = smodel.VDepSReac('CaPm0m1', ssys, slhs = [CaP_m0], srhs = [CaP_m1], k= lambda V: 1.0e3 *3.* alpha_cap(V*1.0e3)* Qt)
CaPm1m2 = smodel.VDepSReac('CaPm1m2', ssys, slhs = [CaP_m1], srhs = [CaP_m2], k= lambda V: 1.0e3 *2.* alpha_cap(V*1.0e3)* Qt)
CaPm2m3 = smodel.VDepSReac('CaPm2m3', ssys, slhs = [CaP_m2], srhs = [CaP_m3], k= lambda V: 1.0e3 *1.* alpha_cap(V*1.0e3)* Qt)

CaPm3m2 = smodel.VDepSReac('CaPm3m2', ssys, slhs = [CaP_m3], srhs = [CaP_m2], k= lambda V: 1.0e3 *3.* beta_cap(V*1.0e3)* Qt)
CaPm2m1 = smodel.VDepSReac('CaPm2m1', ssys, slhs = [CaP_m2], srhs = [CaP_m1], k= lambda V: 1.0e3 *2.* beta_cap(V*1.0e3)* Qt)
CaPm1m0 = smodel.VDepSReac('CaPm1m0', ssys, slhs = [CaP_m1], srhs = [CaP_m0], k= lambda V: 1.0e3 *1.* beta_cap(V*1.0e3)* Qt)

We come to creating our GHK current object (steps.API_1.model.GHKcurr). This object will calculate single-channel current for a given channel state by the GHK flux equation:

(1)\[I_{s}=P_{s}z_{s}^{2}\frac{V_{m}F^{2}}{RT}\frac{[S]_{i}-[S]_{o}exp(-z_{s}V_{m}F/RT)}{1-exp(-z_{s}V_{m}F/RT)}\]

where \(I_{s}\) is the single-channel current (amps) of ion S, \(P_{s}\) is the single-channel permeability of ion S (\(m^{3}.s^{-1}\)), \(z_{s}\) is the valence of ion S, \(V_{m}\) is the membrane voltage (volts), F is the Faraday constant, R is the gas constant, T is temperature (Kelvin), \([S]_{i}\) is the intracellular concentration of ion S (\(mol.m^{-3}\)) and \([S]_{o}\) is the extracellular concentration of ion S (\(mol.m^{-3}\)).

When a GHK current is applied in STEPS it (optionally) results in movement of ions between the ‘outer’ and ‘inner’ compartments, the direction of which will depend on the sign of the current and the valence of the ions.

Many of the values required for calculating a GHK current are simulation variables, such as concentrations and voltage, simulation constants such as temperature, or fixed constants such as the Faraday constant and the gas constant. Such values are either known or can be found by STEPS during runtime and so are not part of object construction, with the exception of single-channel permeability which we will come to later. First let’s look at the required arguments to the object constructor, which are, in order: a string identifier, parent surface system, the channel state (reference to a steps.API_1.model.ChanState object) and the permeable ion (reference to a steps.API_1.model.Spec object). There are also optional keyword arguments (‘virtual_oconc’ and ‘computeflux’) and we’ll see that which of these optional arguments are used depends on whether the mesh has an extracellular ‘outer’ compartment available (e.g. the 10um, 20um, 40um and 80um meshes) or not (e.g. the 160um mesh):

if cyl160:
    OC_CaP = smodel.GHKcurr('OC_CaP', ssys, CaP_m3, Ca, virtual_oconc = Ca_oconc, computeflux = True)
else:
    OC_CaP = smodel.GHKcurr('OC_CaP', ssys, CaP_m3, Ca, computeflux = True)

First let’s look at the ‘virtual_oconc’ argument. This option allows us to not explicitly model the extracellular (‘outer’) concentration of the ion, useful because often the extracellular compartment is not modelled. This option, rather, allows a fixed ‘outer’ concentration for the ion to be specified and that number will be used in the GHK flux calculations. The value of the parameter Ca_oconc in the extra.constants module is 2mM, so when the 160um mesh is used (when the cyl160 flag is True) where there is no extracellular compartment, the extracellular concentration of Ca2+ in all GHK flux calculations will be 2mM.

The second optional argument is ‘computeflux’. This flag (which defaults to True) tells STEPS whether to model this GHK current process as ion transport or not. If ‘computeflux’ is True, then the calculated GHK current will result in transport of ions between the ‘outer’ and ‘inner’ compartments. For example, if over some 0.01ms time step, somewhere on the membrane a mean current of approximately 1.6pA is calculated through a membrane channel to which a GHK current is applied, then for an ion of valence 2+ this means that 50 ions moved from one compartment to the other. The direction of movement depends on the signs of the current and the ion valence. The movement only occurs between surface tetrahedrons surrounding the membrane triangles in which the channels reside and so, for ions where this kind of process occurs, for accuracy it is necessary to model diffusion of these ions at least within the inner compartment and often within both compartments. This can be an expensive computation, particularly where concentrations are in the millimolar range, which shows the value of the ‘computeflux’ flag- if the GHK flux is applied to an ion which does not have any other particularly important effects in the model other than its effect on membrane excitability (a possible example is potassium) then it may be a good labour-saver to clamp ‘inner’ and ‘outer’ concentrations of the ion and turn off the transport of ions as an approximation. However, in this model if we set ‘computeflux’ to False then the result would be no intracellular calcium, which is obviously not desirable, and so the ‘computeflux’ flag is set to True, as it usually will be for most ions in most models.

We might notice by equation (1) that there is some missing information that we did not supply to the steps.API_1.model.GHKcurr constructor, specifically the valence and the single-channel permeability. We will come to the latter soon, but first if we go back to the Calcium dynamics description, where we created the calcium species in the system we see this:

# Calcium
Ca = smodel.Spec('Ca', mdl)
Ca.setValence(2)

For calcium (and only for calcium) we used function steps.API_1.model.Spec.setValence() to specify a valence of 2. ‘Valence’ can be an ambiguous term, but here it means the net elementary electrical charge per ion, which in this example for Ca2+ is +2. Negative valences can of course be specified by using a negative number. It is essential that this function is called to set a valence for any ion that will be used for a GHK current in the simulation- if no valence is specified the result will be an error.

The last parameter we need to set is single-channel permeability. Because conductance is not constant for a GHK current (apart from under certain unusual conditions) one value for a conductance parameter does not suffice. However, since single-channel permeability is often rather a difficult parameter to define, STEPS does provide functionality for estimating the permeability. So we have two options for setting single-channel permeability: steps.API_1.model.GHKcurr.setP() and steps.API_1.model.GHKcurr.setPInfo(). The first is straightforward and simply means providing single-channel permeability in S.I. units of cubic metres / second. In this model the parameter can be found in the extra/constants.py file and takes the value 2.5e-20 cubic metres / second 4:

OC_CaP.setP(CaP_P)

The second option, the steps.API_1.model.GHKcurr.setPInfo() function, requires some explanation. In effect, the conductance of a channel that is modelled by the GHK flux equation varies with voltage (Figure 10.1) with a dependence on the ‘outer’ and ‘inner’ concentrations of the ion (in fact conductance is only constant with voltage when these concentrations are equal), as well as weakly on temperature.

Figure 10.1: A single-channel GHK flux in the physiological range for a typical monovalent cation compared to an Ohmic approximation. The GHK flux is calculated with single-channel permeability of 9e-20 cubic metres / second, fixed extracellular concentration of 4mM, fixed intracellular concentration of 155mM and temperature of 20 Celsius. The single-channel Ohmic conductance is 20pS with reversal potential -77mV.

STEPS is able to estimate single-channel permeability from single-channel conductance, but for STEPS to do so the user must supply information about the conditions under which the conductance was measured, and in theory this should be enough to find the single-channel permeability since it is assumed constant (although there are occasions when permeability too can have some weak voltage dependence 3, which is, however, currently not possible to model with STEPS). Specifically, the steps.API_1.model.GHKcurr.setPInfo() function requires arguments of: estimated single-channel conductance 5 (units: Siemens), one voltage within the range at which conductance was measured (Volts), temperature (Kelvin), ‘outer’ concentration of the ion (molar), and ‘inner’ concentration of the ion (molar). Since the valence of the ion is known it is not necessary to supply that information to the steps.API_1.model.GHKcurr.setPInfo() function. So, for example, for some GHKcurrent object called K_GHK, if we measured single-channel conductance as 20pS in a small voltage range around -22mV at 20 degrees Celsius (293.15 Kelvin) with an estimated extracellular ion concentration of 4mM and intracellular concentration of 155mM, then we would call the function like so:

K_GHK.setPInfo(g = 20e-12, V = -22e-3, T = 293.15, oconc = 4e-3, iconc = 155e-3)

and the single-channel permeability would be set to approximately 9e-20 cubic metres / second. The behaviour of such a channel is shown in Figure 10.1.

We are now familiar, through aspects discussed so far in this chapter and other chapters, with most of the concepts applied for this model, so a very detailed description is not necessary for most remaining parts of the model. We move on to our other three ion channels in the model.

13.8.2.3. T-type Calcium channel

Like the P-type Calcium channel, transitions between channel states of the T-type Calcium channel are voltage-dependent and we model the calcium current as a GHK current:

CaTchan = smodel.Chan('CaTchan', mdl)

CaT_m0h0 = smodel.ChanState('CaT_m0h0', mdl, CaTchan)
CaT_m0h1 = smodel.ChanState('CaT_m0h1', mdl, CaTchan)
CaT_m1h0 = smodel.ChanState('CaT_m1h0', mdl, CaTchan)
CaT_m1h1 = smodel.ChanState('CaT_m1h1', mdl, CaTchan)
CaT_m2h0 = smodel.ChanState('CaT_m2h0', mdl, CaTchan)
CaT_m2h1 = smodel.ChanState('CaT_m2h1', mdl, CaTchan)


CaTm0h0_m1h0 = smodel.VDepSReac('CaTm0h0_m1h0', ssys, slhs = [CaT_m0h0], srhs = [CaT_m1h0], k= lambda V: 1.0e3 *2.* alpham_cat(V*1.0e3))
CaTm1h0_m2h0 = smodel.VDepSReac('CaTm1h0_m2h0', ssys, slhs = [CaT_m1h0], srhs = [CaT_m2h0], k= lambda V: 1.0e3 *1.* alpham_cat(V*1.0e3))

CaTm2h0_m1h0 = smodel.VDepSReac('CaTm2h0_m1h0', ssys, slhs = [CaT_m2h0], srhs = [CaT_m1h0], k= lambda V: 1.0e3 *2.* betam_cat(V*1.0e3))
CaTm1h0_m0h0 = smodel.VDepSReac('CaTm1h0_m0h0', ssys, slhs = [CaT_m1h0], srhs = [CaT_m0h0], k= lambda V: 1.0e3 *1.* betam_cat(V*1.0e3))

CaTm0h1_m1h1 = smodel.VDepSReac('CaTm0h1_m1h1', ssys, slhs = [CaT_m0h1], srhs = [CaT_m1h1], k= lambda V: 1.0e3 *2.* alpham_cat(V*1.0e3))
CaTm1h1_m2h1 = smodel.VDepSReac('CaTm1h1_m2h1', ssys, slhs = [CaT_m1h1], srhs = [CaT_m2h1], k= lambda V: 1.0e3 *1.* alpham_cat(V*1.0e3))

CaTm2h1_m1h1 = smodel.VDepSReac('CaTm2h1_m1h1', ssys, slhs = [CaT_m2h1], srhs = [CaT_m1h1], k= lambda V: 1.0e3 *2.* betam_cat(V*1.0e3))
CaTm1h1_m0h1 = smodel.VDepSReac('CaTm1h1_m0h1', ssys, slhs = [CaT_m1h1], srhs = [CaT_m0h1], k= lambda V: 1.0e3 *1.* betam_cat(V*1.0e3))


CaTm0h0_m0h1 = smodel.VDepSReac('CaTm0h0_m0h1', ssys, slhs = [CaT_m0h0], srhs = [CaT_m0h1], k= lambda V: 1.0e3 *1.* alphah_cat(V*1.0e3))
CaTm1h0_m1h1 = smodel.VDepSReac('CaTm1h0_m1h1', ssys, slhs = [CaT_m1h0], srhs = [CaT_m1h1], k= lambda V: 1.0e3 *1.* alphah_cat(V*1.0e3))
CaTm2h0_m2h1 = smodel.VDepSReac('CaTm2h0_m2h1', ssys, slhs = [CaT_m2h0], srhs = [CaT_m2h1], k= lambda V: 1.0e3 *1.* alphah_cat(V*1.0e3))

CaTm2h1_m2h0 = smodel.VDepSReac('CaTm2h1_m2h0', ssys, slhs = [CaT_m2h1], srhs = [CaT_m2h0], k= lambda V: 1.0e3 *1.* betah_cat(V*1.0e3))
CaTm1h1_m1h0 = smodel.VDepSReac('CaTm1h1_m1h0', ssys, slhs = [CaT_m1h1], srhs = [CaT_m1h0], k= lambda V: 1.0e3 *1.* betah_cat(V*1.0e3))
CaTm0h1_m0h0 = smodel.VDepSReac('CaTm0h1_m0h0', ssys, slhs = [CaT_m0h1], srhs = [CaT_m0h0], k= lambda V: 1.0e3 *1.* betah_cat(V*1.0e3))

if cyl160:
    OC_CaT = smodel.GHKcurr('OC_CaT', ssys, CaT_m2h1, Ca, virtual_oconc = Ca_oconc, computeflux = True)
else:
    OC_CaT = smodel.GHKcurr('OC_CaT', ssys, CaT_m2h1, Ca, computeflux = True)

OC_CaT.setP(CaT_P)

13.8.2.4. BK-type Calcium-activated Potassium channel

The BK channel in the model undergoes both voltage-dependent and non-voltage dependent processes and so its Channel State transitions are described by both steps.API_1.model.SReac and steps.API_1.model.VDepSReac objects. This is an example of the same functionality for Channel State objects as for Species objects from which they are derived (as described in 3) : Channel State objects can be used interchangeably anywhere a Species object can be used, and so they can interact with other Channel States and Species through Surface Reactions, they may diffuse on the surface, or even diffuse in volumes and undergo volume reactions. Here we will notice that Channel States (e.g. BK_C0) appear alongside Species (Ca) in steps.API_1.model.SReac constructors:

BKchan = smodel.Chan('BKchan', mdl)

BK_C0 = smodel.ChanState('BK_C0', mdl, BKchan)
BK_C1 = smodel.ChanState('BK_C1', mdl, BKchan)
BK_C2 = smodel.ChanState('BK_C2', mdl, BKchan)
BK_C3 = smodel.ChanState('BK_C3', mdl, BKchan)
BK_C4 = smodel.ChanState('BK_C4', mdl, BKchan)
BK_O0 = smodel.ChanState('BK_O0', mdl, BKchan)
BK_O1 = smodel.ChanState('BK_O1', mdl, BKchan)
BK_O2 = smodel.ChanState('BK_O2', mdl, BKchan)
BK_O3 = smodel.ChanState('BK_O3', mdl, BKchan)
BK_O4 = smodel.ChanState('BK_O4', mdl, BKchan)


BKCAC0 = smodel.SReac('BKCAC0', ssys, slhs = [BK_C0], ilhs = [Ca], srhs = [BK_C1], kcst = c_01)
BKCAC1 = smodel.SReac('BKCAC1', ssys, slhs = [BK_C1], ilhs = [Ca], srhs = [BK_C2], kcst = c_12)
BKCAC2 = smodel.SReac('BKCAC2', ssys, slhs = [BK_C2], ilhs = [Ca], srhs = [BK_C3], kcst = c_23)
BKCAC3 = smodel.SReac('BKCAC3', ssys, slhs = [BK_C3], ilhs = [Ca], srhs = [BK_C4], kcst = c_34)

BKC0 = smodel.SReac('BKC0', ssys, slhs = [BK_C1], srhs = [BK_C0], irhs = [Ca], kcst = c_10)
BKC1 = smodel.SReac('BKC1', ssys, slhs = [BK_C2], srhs = [BK_C1], irhs = [Ca], kcst = c_21)
BKC2 = smodel.SReac('BKC2', ssys, slhs = [BK_C3], srhs = [BK_C2], irhs = [Ca], kcst = c_32)
BKC3 = smodel.SReac('BKC3', ssys, slhs = [BK_C4], srhs = [BK_C3], irhs = [Ca], kcst = c_43)

BKCAO0 = smodel.SReac('BKCAO0', ssys, slhs = [BK_O0], ilhs = [Ca], srhs = [BK_O1], kcst = o_01)
BKCAO1 = smodel.SReac('BKCAO1', ssys, slhs = [BK_O1], ilhs = [Ca], srhs = [BK_O2], kcst = o_12)
BKCAO2 = smodel.SReac('BKCAO2', ssys, slhs = [BK_O2], ilhs = [Ca], srhs = [BK_O3], kcst = o_23)
BKCAO3 = smodel.SReac('BKCAO3', ssys, slhs = [BK_O3], ilhs = [Ca], srhs = [BK_O4], kcst = o_34)

BKO0 = smodel.SReac('BKO0', ssys, slhs = [BK_O1], srhs = [BK_O0], irhs = [Ca], kcst = o_10)
BKO1 = smodel.SReac('BKO1', ssys, slhs = [BK_O2], srhs = [BK_O1], irhs = [Ca], kcst = o_21)
BKO2 = smodel.SReac('BKO2', ssys, slhs = [BK_O3], srhs = [BK_O2], irhs = [Ca], kcst = o_32)
BKO3 = smodel.SReac('BKO3', ssys, slhs = [BK_O4], srhs = [BK_O3], irhs = [Ca], kcst = o_43)

BKC0O0 = smodel.VDepSReac('BKC0O0', ssys, slhs = [BK_C0], srhs = [BK_O0], k=lambda V: f_0(V))
BKC1O1 = smodel.VDepSReac('BKC1O1', ssys, slhs = [BK_C1], srhs = [BK_O1], k=lambda V: f_1(V))
BKC2O2 = smodel.VDepSReac('BKC2O2', ssys, slhs = [BK_C2], srhs = [BK_O2], k=lambda V: f_2(V))
BKC3O3 = smodel.VDepSReac('BKC3O3', ssys, slhs = [BK_C3], srhs = [BK_O3], k=lambda V: f_3(V))
BKC4O4 = smodel.VDepSReac('BKC4O4', ssys, slhs = [BK_C4], srhs = [BK_O4], k=lambda V: f_4(V))

BKO0C0 = smodel.VDepSReac('BKO0C0', ssys, slhs = [BK_O0], srhs = [BK_C0], k=lambda V: b_0(V))
BKO1C1 = smodel.VDepSReac('BKO1C1', ssys, slhs = [BK_O1], srhs = [BK_C1], k=lambda V: b_1(V))
BKO2C2 = smodel.VDepSReac('BKO2C2', ssys, slhs = [BK_O2], srhs = [BK_C2], k=lambda V: b_2(V))
BKO3C3 = smodel.VDepSReac('BKO3C3', ssys, slhs = [BK_O3], srhs = [BK_C3], k=lambda V: b_3(V))
BKO4C4 = smodel.VDepSReac('BKO4C4', ssys, slhs = [BK_O4], srhs = [BK_C4], k=lambda V: b_4(V))

steps.API_1.model.OhmicCurr objects are applied to 5 different channel states, demonstrating the support for multiple conducting/permeable states for a channel:

OC_BK0 = smodel.OhmicCurr('OC_BK0', ssys, chanstate = BK_O0, erev = BK_rev, g = BK_G )
OC_BK1 = smodel.OhmicCurr('OC_BK1', ssys, chanstate = BK_O1, erev = BK_rev, g = BK_G )
OC_BK2 = smodel.OhmicCurr('OC_BK2', ssys, chanstate = BK_O2, erev = BK_rev, g = BK_G )
OC_BK3 = smodel.OhmicCurr('OC_BK3', ssys, chanstate = BK_O3, erev = BK_rev, g = BK_G )
OC_BK4 = smodel.OhmicCurr('OC_BK4', ssys, chanstate = BK_O4, erev = BK_rev, g = BK_G )

13.8.2.5. SK-type Calcium-activated Potassium channel

The SK channel does not have any voltage dependence, and contains two conducting states:

SKchan = smodel.Chan('SKchan', mdl)

SK_C1 = smodel.ChanState('SK_C1', mdl, SKchan)
SK_C2 = smodel.ChanState('SK_C2', mdl, SKchan)
SK_C3 = smodel.ChanState('SK_C3', mdl, SKchan)
SK_C4 = smodel.ChanState('SK_C4', mdl, SKchan)
SK_O1 = smodel.ChanState('SK_O1', mdl, SKchan)
SK_O2 = smodel.ChanState('SK_O2', mdl, SKchan)


SKCAC1 = smodel.SReac('SKCAC1', ssys, slhs = [SK_C1], ilhs = [Ca], srhs = [SK_C2], kcst = dirc2_t)
SKCAC2 = smodel.SReac('SKCAC2', ssys, slhs = [SK_C2], ilhs = [Ca], srhs = [SK_C3], kcst = dirc3_t)
SKCAC3 = smodel.SReac('SKCAC3', ssys, slhs = [SK_C3], ilhs = [Ca], srhs = [SK_C4], kcst = dirc4_t)

SKC1 = smodel.SReac('SKC1', ssys, slhs = [SK_C2], srhs = [SK_C1], irhs = [Ca], kcst = invc1_t)
SKC2 = smodel.SReac('SKC2', ssys, slhs = [SK_C3], srhs = [SK_C2], irhs = [Ca], kcst = invc2_t)
SKC3 = smodel.SReac('SKC3', ssys, slhs = [SK_C4], srhs = [SK_C3], irhs = [Ca], kcst = invc3_t)

SKC3O1 = smodel.SReac('SKC3O1', ssys, slhs = [SK_C3], srhs = [SK_O1], kcst = diro1_t)
SKC4O2 = smodel.SReac('SKC4O2', ssys, slhs = [SK_C4], srhs = [SK_O2], kcst = diro2_t)

SKO1C3 = smodel.SReac('SKO1C3', ssys, slhs = [SK_O1], srhs = [SK_C3], kcst = invo1_t)
SKO2C4 = smodel.SReac('SKO2C4', ssys, slhs = [SK_O2], srhs = [SK_C4], kcst = invo2_t)

OC1_SK = smodel.OhmicCurr('OC1_SK', ssys, chanstate = SK_O1, erev = SK_rev, g = SK_G )
OC2_SK = smodel.OhmicCurr('OC2_SK', ssys, chanstate = SK_O2, erev = SK_rev, g = SK_G )

13.8.2.6. Leak channel

The leak conductance is described as a leak channel, though another option for setting the leak would have been to (later) use function steps.API_1.solver.Tetexact.setMembRes() (also supported in TetODE solver: steps.API_1.solver.TetODE.setMembRes()):

L = smodel.Chan('L', mdl)
Leak = smodel.ChanState('Leak', mdl, L)

OC_L = smodel.OhmicCurr('OC_L', ssys, chanstate = Leak, erev = L_rev, g = L_G)

13.8.3. Geometry specification

This model is set up for relatively simple geometry- cylinders of diameter 2um and varying lengths from 10um to 160um. As discussed previously, the 160um cylinder does not include an extracellular compartment within the mesh whereas the 10um, 20um, 40um and 80um cylinders do, so the initialisation is slightly different for the 160um mesh compared to the others.

First we separate the tetrahedrons into the ‘inner’ tetrahedrons, which will form the cytosolic compartment, and the ‘outer’ tetrahedrons, which will form the extracellular compartment. We do that by finding the tetrahedrons within a 1um radius along the cylinder axis (which is along the z axis) to form the inner compartment, and exclude those tetrahedrons from the complete set to find the outer compartment:

mesh = meshio.loadMesh('./meshes/'+meshfile_ab)[0]

outer_tets = range(mesh.ntets)

# Will return all tetrahedrons within a 1um radius along the z-axis
inner_tets = gettets.getcyl(mesh, 1e-6, -200e-6, 200e-6)[0]

for i in inner_tets: outer_tets.remove(i)
assert(outer_tets.__len__() + inner_tets.__len__() == mesh.ntets)

# Create an intracellular compartment, and extracellular compartment if not the 160um mesh
cyto = sgeom.TmComp('cyto', mesh, inner_tets)
cyto.addVolsys('vsys')
if not cyl160: outer = sgeom.TmComp('outer', mesh, outer_tets)

Then at this point we find the tetrahedron at the centre of the mesh, where we will record voltage from during simulation 6:

# Record voltage from the central tetrahedron
cent_tet = mesh.findTetByPoint([0.0,0.0,0.0])

Now we move on to find the triangles that form the cell membrane between the intracellular and extracellular compartments. The way to do that is different depending on whether we are using the 160um cylinder or not. The circular faces at each end of the cylinder are excluded from the membrane:

if cyl160:
    # Find points a small distance inside the circular boundaries: then exclude triangles outside these bounds
    LENGTH = mesh.getBoundMax()[2] - mesh.getBoundMin()[2]
    boundminz = mesh.getBoundMin()[2] + LENGTH/mesh.ntets
    boundmaxz = mesh.getBoundMax()[2] - LENGTH/mesh.ntets

    memb_tris = list(mesh.getSurfTris())
    minztris = []
    maxztris = []
    for tri in memb_tris:
        zminboundtri = True
        zmaxboundtri = True
        tritemp = mesh.getTri(tri)
        trizs = [0.0, 0.0, 0.0]
        trizs[0] = mesh.getVertex(tritemp[0])[2]
        trizs[1] = mesh.getVertex(tritemp[1])[2]
        trizs[2] = mesh.getVertex(tritemp[2])[2]
        for j in range(3):
            if (trizs[j]>boundminz): zminboundtri = False
        if (zminboundtri):
            minztris.append(tri)
            continue
        for j in range(3):
            if (trizs[j]< boundmaxz): zmaxboundtri = False
        if (zmaxboundtri):
            maxztris.append(tri)

    for t in minztris: memb_tris.remove(t)
    for t in maxztris: memb_tris.remove(t)

else:
    # Find common triangles between inner and outer compartments
    out_tris = set()
    for i in outer_tets:
            tritemp = mesh.getTetTriNeighb(i)
            for j in range(4): out_tris.add(tritemp[j])

    in_tris = set()
    for i in inner_tets:
            tritemp = mesh.getTetTriNeighb(i)
            for j in range(4): in_tris.add(tritemp[j])

    memb_tris = out_tris.intersection(in_tris)
    memb_tris = list(memb_tris)

We also find the submembrane tetrahedrons, that is all tetrahedrons connected to a membrane triangle from the intracellular side:

########## Find the submembrane tets

memb_tet_neighb = []
for i in memb_tris:
    tettemp = mesh.getTriTetNeighb(i)
    for j in tettemp:
        memb_tet_neighb.append(j)

submemb_tets = []
for i in memb_tet_neighb:
    if i in inner_tets:
        submemb_tets.append(i)

# Find the volume of this region
vol = 0.0
for i in submemb_tets:
    vol = vol + mesh.getTetVol(i)

And we are ready to create our membrane. That is a steps.API_1.geom.Memb object for which we are able to model electrical excitability by adding ion channels and solving potential across the membrane and within the intracellular conduction volume. For details of the method see 3. First we need to create a patch, named a little confusingly memb:

########## Create a membrane
if cyl160:
    memb = sgeom.TmPatch('memb', mesh, memb_tris, cyto)
else:
    memb = sgeom.TmPatch('memb', mesh, memb_tris, cyto, outer)

memb.addSurfsys('ssys')

And then we create the membrane. Here we take advantage of previously found and stored connectivity optimisation (by function steps.API_1.solver.Tetexact.saveMembOpt()), in files such as /meshes/Cylinder2_dia2um_L10um_outer0_3um_0.3shell_0.3size_19156tets_adaptive.inp_optimalidx. Connectivity optimisation is discussed in 3 and Simulating Membrane Potential:

membrane = sgeom.Memb('membrane', mesh, [memb], opt_file_name = './meshes/'+meshfile_ab+"_optimalidx")

13.8.4. Simulation with Tetexact

13.8.4.1. Initialization

Before simulation we create the random number generator, but this time initialise with a fixed number. The reason for doing this is to ensure that initial conditions (placement of ion channels etc) is the same for each simulation iteration, so that the stochastic effects observed are purely from stochastic kinetics, and not due to different arrangements of channels:

r = srng.create_mt19937(512)
r.initialize(7)

We create the spatial stochastic solver object, turning the voltage calculation on by setting the ‘calcMembPot’ flag to True. Recall that we need one (and only one) steps.API_1.geom.Memb object to exist in the geometry for this to work:

sim = ssolver.Tetexact(mdl, mesh, r, True)
sim.reset()

Next we see a new function, steps.API_1.solver.Tetexact.setTemp(), which sets the simulation temperature. Currently, this will only influence any GHK flux rates, and will have no influence on any other kinetics. The value for TEMPERATURE is 34 degrees Celsius and we need to set temperature in Kelvin, following the usual S.I. rule in STEPS:

sim.setTemp(TEMPERATURE+273.15)

Next we inject ions, buffers and channels. Most values appear in the extra/constants.py file:

if not cyl160:
    sim.setCompConc('outer', 'Ca', Ca_oconc)
    sim.setCompClamped('outer', 'Ca', True)

sim.setCompConc('cyto', 'Ca', Ca_iconc)

sim.setCompConc('cyto', 'Mg', Mg_conc)

surfarea = sim.getPatchArea('memb')

pumpnbs = 6.022141e12*surfarea
sim.setPatchCount('memb', 'Pump', round(pumpnbs))
sim.setPatchCount('memb', 'CaPump', 0)

sim.setCompConc('cyto', 'iCBsf', iCBsf_conc)
sim.setCompConc('cyto', 'iCBCaf', iCBCaf_conc)
sim.setCompConc('cyto', 'iCBsCa', iCBsCa_conc)
sim.setCompConc('cyto', 'iCBCaCa', iCBCaCa_conc)

sim.setCompConc('cyto', 'CBsf', CBsf_conc)
sim.setCompConc('cyto', 'CBCaf', CBCaf_conc)
sim.setCompConc('cyto', 'CBsCa', CBsCa_conc)
sim.setCompConc('cyto', 'CBCaCa', CBCaCa_conc)

sim.setCompConc('cyto', 'PV', PV_conc)
sim.setCompConc('cyto', 'PVCa', PVCa_conc)
sim.setCompConc('cyto', 'PVMg', PVMg_conc)

sim.setPatchCount('memb', 'CaP_m0' , round(CaP_ro*surfarea*CaP_m0_p))
sim.setPatchCount('memb', 'CaP_m1' , round(CaP_ro*surfarea*CaP_m1_p))
sim.setPatchCount('memb', 'CaP_m2' , round(CaP_ro*surfarea*CaP_m2_p))
sim.setPatchCount('memb', 'CaP_m3' , round(CaP_ro*surfarea*CaP_m3_p))

sim.setPatchCount('memb', 'CaT_m0h0' , round(CaT_ro*surfarea*CaT_m0h0_p))
sim.setPatchCount('memb', 'CaT_m1h0' , round(CaT_ro*surfarea*CaT_m1h0_p))
sim.setPatchCount('memb', 'CaT_m2h0' , round(CaT_ro*surfarea*CaT_m2h0_p))
sim.setPatchCount('memb', 'CaT_m0h1' , round(CaT_ro*surfarea*CaT_m0h1_p))
sim.setPatchCount('memb', 'CaT_m1h1' , round(CaT_ro*surfarea*CaT_m1h1_p))
sim.setPatchCount('memb', 'CaT_m2h1' , round(CaT_ro*surfarea*CaT_m2h1_p))

sim.setPatchCount('memb', 'BK_C0' , round(BK_ro*surfarea*BK_C0_p))
sim.setPatchCount('memb', 'BK_C1' , round(BK_ro*surfarea*BK_C1_p))
sim.setPatchCount('memb', 'BK_C2' , round(BK_ro*surfarea*BK_C2_p))
sim.setPatchCount('memb', 'BK_C3' , round(BK_ro*surfarea*BK_C3_p))
sim.setPatchCount('memb', 'BK_C4' , round(BK_ro*surfarea*BK_C4_p))

sim.setPatchCount('memb', 'BK_O0' , round(BK_ro*surfarea*BK_O0_p))
sim.setPatchCount('memb', 'BK_O1' , round(BK_ro*surfarea*BK_O1_p))
sim.setPatchCount('memb', 'BK_O2' , round(BK_ro*surfarea*BK_O2_p))
sim.setPatchCount('memb', 'BK_O3' , round(BK_ro*surfarea*BK_O3_p))
sim.setPatchCount('memb', 'BK_O4' , round(BK_ro*surfarea*BK_O4_p))

sim.setPatchCount('memb', 'SK_C1' , round(SK_ro*surfarea*SK_C1_p))
sim.setPatchCount('memb', 'SK_C2' , round(SK_ro*surfarea*SK_C2_p))
sim.setPatchCount('memb', 'SK_C3' , round(SK_ro*surfarea*SK_C3_p))
sim.setPatchCount('memb', 'SK_C4' , round(SK_ro*surfarea*SK_C4_p))

sim.setPatchCount('memb', 'SK_O1' , round(SK_ro*surfarea*SK_O1_p))
sim.setPatchCount('memb', 'SK_O2' , round(SK_ro*surfarea*SK_O2_p))

sim.setPatchCount('memb', 'Leak', int(L_ro * surfarea))

And finally set some parameters for the ‘E-Field’ voltage calculation using solver functions, in order: steps.API_1.solver.Tetexact.setEfieldDT() to set the communication time-step between the voltage calculation and the reaction-diffusion simulation (from extra/constants.py file: EF_DT = 0.02e-3 seconds), steps.API_1.solver.Tetexact.setMembPotential() to set the initial membrane potential (init_pot = -60e-3 V), steps.API_1.solver.Tetexact.setMembVolRes() to set resistivity of the conduction volume enclosed by the membrane (Ra = 2.357 ohm.m) and steps.API_1.solver.Tetexact.setMembCapac() to set membrane capacitance (memb_capac = 1.5e-2 F/m²):

sim.setEfieldDT(EF_DT)
sim.setMembPotential('membrane', init_pot)
sim.setMembVolRes('membrane', Ra)
sim.setMembCapac('membrane', memb_capac)

Next we create data storage location and files, applying some try-except tests which are there because the locations have been created previously. Three files will be created for data storage, currents.dat, voltage.dat and calcium.dat, the location of which will depend on the root directory read from the command line and the present system time:

c=time.ctime()

dc = c.split()[1]+c.split()[2]+'_'+c.split()[3]+'_'+c.split()[4]
dc= dc.replace(':', '_')

try: os.mkdir(root+'data')
except: pass
try: os.mkdir(root+'data/' +  'StochasticCaburst')
except: pass
try: os.mkdir(root+'data/' +  'StochasticCaburst/'+meshfile_ab)
except: pass
os.mkdir(root+'data/' +  'StochasticCaburst/'+meshfile_ab+'/'+iter_n+'__'+dc )

datfile =  open(root+'data/' +  'StochasticCaburst/'+meshfile_ab+'/'+iter_n+'__'+dc + '/currents.dat', 'w')
datfile2 = open(root+'data/' +  'StochasticCaburst/'+meshfile_ab+'/'+iter_n+'__'+dc + '/voltage.dat', 'w')
datfile3 = open(root+'data/' +  'StochasticCaburst/'+meshfile_ab+'/'+iter_n+'__'+dc + '/calcium.dat', 'w')

Before simulation we reinitialize the random number generator with a random seed (based on the iter_n number given from the command line). This is because initialization was from a fixed number earlier to ensure consistent initial conditions, and now we need to reinitialize to ensure each simulation iteration is unique:

r.initialize(100*int(iter_n))

13.8.4.2. Running the simulation

At last we are ready to run the simulation, which is achieved simply by calls to function steps.API_1.solver.Tetexact.run() within a simulation loop. The rest of the following code is simply used for recording data. Functions steps.API_1.solver.Tetexact.getTriGHKI() and steps.API_1.solver.Tetexact.getTriOhmicI() allow us to record the currents from all membrane triangles, which are then summed to record the total current for each channel type. Voltage is recorded from just one central location in the conduction volume with function steps.API_1.solver.Tetexact.getTetV() and calcium is recorded from both the submembrane region using function steps.API_1.solver.Tetexact.getTetCount() and summing over submembrane tets, and for the cytosol compartment as a whole with function steps.API_1.solver.Tetexact.getCompConc()

for l in range(NTIMEPOINTS):

    sim.run(TIMECONVERTER*l)

    # total P-type calcium current
    tcur_CaP = 0.0
    # total T-type calcium current
    tcur_CaT = 0.0
    # total BK-type potasium current
    tcur_BK = 0.0
    # total SK-type potassium current
    tcur_SK = 0.0

    # The number of calcium ions in the submembrane region
    tca_count = 0.0

    # Record calcium in the submembrane region
    for m in submemb_tets:
        tca_count = tca_count + sim.getTetCount(m,'Ca')

    # Record the membrane currents
    for m in memb_tris:
        tcur_CaP = tcur_CaP + sim.getTriGHKI(m,'OC_CaP')
        tcur_CaT = tcur_CaT + sim.getTriGHKI(m,'OC_CaT')
        tcur_BK = tcur_BK + sim.getTriOhmicI(m,'OC_BK0') \
            + sim.getTriOhmicI(m,'OC_BK1') \
            + sim.getTriOhmicI(m,'OC_BK2') \
            + sim.getTriOhmicI(m,'OC_BK3') \
            + sim.getTriOhmicI(m,'OC_BK4')
        tcur_SK = tcur_SK + sim.getTriOhmicI(m,'OC1_SK') + sim.getTriOhmicI(m,'OC2_SK')

    # Write the data
    datfile.write('%.6g' %(1.0e3*TIMECONVERTER*l) + ' ')
    datfile.write('%.6g' %((tcur_CaP*1.0e-1)/surfarea) + ' ')
    datfile.write('%.6g' %((tcur_CaT*1.0e-1)/surfarea) + ' ')
    datfile.write('%.6g' %((tcur_BK*1.0e-1)/surfarea) + ' ')
    datfile.write('%.6g' %((tcur_SK*1.0e-1)/surfarea) + ' ')
    datfile.write('\n')

    datfile2.write('%.6g' %(1.0e3*TIMECONVERTER*l) + ' ')
    datfile2.write('%.6g' %(sim.getTetV(cent_tet)*1.0e3) + ' ')
    datfile2.write('\n')

    datfile3.write('%.6g' %(1.0e3*TIMECONVERTER*l) + ' ')
    datfile3.write('%.6g' %(sim.getCompConc('cyto', 'Ca')*1.0e6) + ' ')
    datfile3.write('%.6g' %tca_count + ' ')
    datfile3.write('\n')

datfile.close()
datfile2.close()
datfile3.close()

Figure 10.2 shows five example runs of this script using the 80um mesh with different random number seeds. This model is analysed in much more depth in 1.

Figure 10.2: Voltage traces from five iterations on the 80um mesh.

Footnotes

1(1,2,3,4)

Anwar H, Hepburn I, Nedelescu H, Chen W, De Schutter E (2013) Stochastic Calcium Mechanisms Cause Dendritic Calcium Spike Variability. The Journal of Neuroscience, 33(40): 15848-15867, doi: 10.1523/​JNEUROSCI.1722-13.2013.

2

Anwar H, Hong S, De Schutter E (2012) Controlling Ca2+-activated K+ channels with models of Ca2+ buffering in Purkinje cells. Cerebellum, 11(3):681-93, doi: 10.1007/s12311-010-0224-3.

3(1,2,3,4,5)

Hepburn I, Cannon R and De Schutter E (2013) Efficient calculation of the quasi-static electrical potential on a tetrahedral mesh and its implementation in STEPS. Frontiers in Computational Neuroscience: 7:129, doi: 10.3389/fncom.2013.00129

4

The same considerations for converting membrane permeability to single-channel permeability apply as for conductance discussed in Simulating Membrane Potential, requiring some estimate of the channel density.

5

Since it is assumed that conductance is measured by estimating the slope of an I-V curve over some small voltage range, the conductance will be treated as a slope conductance for the purposes of single-channel permeability estimation.

6

We may record voltage from anywhere on the membrane surface or within the ‘conduction volume’ (here and in most models the conduction volume is the cytosolic compartment).