Initialization Strategies

A version of this notebook may be run online via Google Colab at https://tinyurl.com/rxd-basic-initialization (make a copy or open in playground mode).

The time series of a chemical concentration necessarily depends on its initial conditions; i.e. the concentration at time 0. An analogous statement is true for gating variables, etc. How do we specify this?

Option 1: NEURON and NMODL defaults

If the species corresponds to one with initial conditions specified by NMODL (or in the case of sodium, potassium, or calcium with meaningful NEURON defaults), then omitting the initial argument will tell NEURON to use those rules. e.g.

from neuron import h, rxd
from neuron.units import mV

soma = h.Section(name='soma')
cyt = rxd.Region(soma.wholetree(), name='cyt', nrn_region='i')

ca = rxd.Species(cyt, name='ca', charge=2, atolscale=1e-6)
na = rxd.Species(cyt, name='na', charge=1)
k = rxd.Species(cyt, name='k', charge=1)
unknown = rxd.Species(cyt, name='unknown', charge=-1)

h.finitialize(-65 * mV)

print(f'ca: {ca.nodes[0].concentration} mM')
print(f'na: {na.nodes[0].concentration} mM')
print(f'k: {k.nodes[0].concentration} mM')
print(f'unknown: {unknown.nodes[0].concentration} mM')

ca: 5e-05 mM
na: 10.0 mM
k: 54.4 mM
unknown: 1.0 mM

As shown here, unknown ions/proteins are by default assigned a concentration by NEURON of 1 mM. The atolscale value for calcium has no effect on the initialized value, but is included here as an example of best practice for working with low concentrations.

Importantly, the NEURON/NMODL rules only apply if there is a corresponding classical NEURON state variable. That is, nrn_region must be set and the Species must have a name assigned.

Running what is otherwise the same code without the nrn_region assigned causes everything to default to 0 µM:

from neuron import h, rxd
from neuron.units import mV

soma = h.Section(name='soma')
cyt = rxd.Region(soma.wholetree(), name='cyt')

ca = rxd.Species(cyt, name='ca', charge=2)
na = rxd.Species(cyt, name='na', charge=1)
k = rxd.Species(cyt, name='k', charge=1)
unknown = rxd.Species(cyt, name='unknown', charge=-1)

h.finitialize(-65 * mV)

print(f'ca: {ca.nodes[0].concentration} mM')
print(f'na: {na.nodes[0].concentration} mM')
print(f'k: {k.nodes[0].concentration} mM')
print(f'unknown: {unknown.nodes[0].concentration} mM')

ca: 0.0 mM
na: 0.0 mM
k: 0.0 mM
unknown: 0.0 mM

## get rid of previous model
soma = ca = na = k = unknown = None

For extracellular species, there is no equivalent traditional NEURON state variable (as those only exist within and along the cell), however NEURON’s constant initialization parameters for the nrn_region=‘o’ space are used if available; e.g.

from neuron import h, rxd
from neuron.units import mV

ecs = rxd.Extracellular(-100, -100, -100,
                        100, 100, 100,
                        dx=20, volume_fraction=0.2, tortuosity=1.6)

## defining calcium on both intra- and extracellular regions
ca = rxd.Species(ecs, name='ca', charge=2)

## global initialization for NEURON extracellular calcium
h.cao0_ca_ion = 0.42

h.finitialize(-65 * mV)

print(f'ca: {ca.nodes[0].concentration} mM')

ca: 0.42 mM

We could do something similar using cai0_ca_ion to set the global initial intracellular calcium concentration.

## get rid of previous model
soma = ca = ecs = None

Option 2: Uniform initial concentration

Setting initial= to a Species or State assigns that value every time the system reinitializes. e.g.

from neuron import h, rxd
from neuron.units import mV

soma = h.Section(name='soma')

cyt = rxd.Region([soma], name='cyt')
m = rxd.State(cyt, initial=0.47)

h.finitialize(-65 * mV)
print(f'm = {m.nodes[0].value}')

m = 0.47

## get rid of previous model
m = cyt = soma = None

Option 3: Initializing to a function of position

The initial= keyword argument also accepts a callable (e.g. a function) that receives a node object. Nodes have certain properties that are useful for assinging based on position, including .segment (intracellular nodes only) and .x3d, .y3d, and .z3d. Segment-to-segment (or the segment containing a node) distances can be measured directly using h.distance.

Using h.distance:

Here we use the morphology c91662.CNG.swc from NeuroMorpho.Org and initialize based on path distance from the soma.

!wget http://neuromorpho.org/dableFiles/amaral/CNG%20version/c91662.CNG.swc

--2022-03-28 10:16:52--  http://neuromorpho.org/dableFiles/amaral/CNG%20version/c91662.CNG.swc
Resolving neuromorpho.org (neuromorpho.org)... 129.174.130.34
Connecting to neuromorpho.org (neuromorpho.org)|129.174.130.34|:80... connected.
HTTP request sent, awaiting response... 200 OK
Length: 54984 (54K) [text/plain]
Saving to: ‘c91662.CNG.swc’

c91662.CNG.swc      100%[===================>]  53.70K  --.-KB/s    in 0.04s

2022-03-28 10:16:52 (1.27 MB/s) - ‘c91662.CNG.swc’ saved [54984/54984]

from neuron import h, gui, rxd
from neuron.units import mV

h.load_file('stdrun.hoc')
h.load_file('import3d.hoc')

## load the morphology and instantiate at the top level (i.e. not in a class)
cell = h.Import3d_SWC_read()
cell.input('c91662.CNG.swc')
h.Import3d_GUI(cell, 0)
i3d = h.Import3d_GUI(cell, 0)
i3d.instantiate(None) # pass in a class to instantiate inside the class instead

## increase the number of segments
for sec in h.allsec():
    sec.nseg = 1 + 2 * int(sec.L / 20)

soma = h.soma[0]

def my_initial(node):
    # return a certain function of the distance
    return 2 * h.tanh(h.distance(soma(0.5), node) / 1000.)

cyt = rxd.Region(h.allsec(), name='cyt', nrn_region='i')
ca = rxd.Species(cyt, name='ca', charge=2, initial=my_initial)

h.finitialize(-65 * mV)

1.0

import plotly
ps = h.PlotShape(False)
ps.variable(ca[cyt])
ps.scale(0, 2)
ps.plot(plotly).show(renderer='notebook_connected')

Using position:

We continue the above example adding a new species, that is initialized based on the x-coordinate. This could happen, for example, on a platform with a nutrient or temperature gradient:

def my_initial2(node):
    # return a certain function of the x-coordinate
    return 1 + h.tanh(node.x3d / 100.)

alpha = rxd.Parameter(cyt, name='alpha', initial=my_initial2)

h.finitialize(-65 * mV)

1.0

import plotly
ps = h.PlotShape(False)
ps.variable(alpha[cyt])
ps.scale(0, 2)
ps.plot(plotly).show(renderer='notebook_connected')

Option 4: to steady state

Sometimes one might want to initialize a simulation to steady-state where e.g. diffusion, ion channel currents, and chemical reactions all balance each other out. There may be no such possible initial condition due to the interacting parts.

In principle, such initial conditions could be assigned using a variant of the option 3 approach above. In practice, however, it may be simpler to omit the initial= keyword argument, and use an h.FInitializeHandler to loop over locations, setting the values for all states at a given location at the same time. A full example is beyond the scope of this tutorial.