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 obtained from NeuroMorpho.Org and initialize based on path distance from the soma.

!wget -N https://raw.githubusercontent.com/neuronsimulator/resources/8b1290d5c8ab748dd6251be5bd46a4e3794d742f/notebooks/rxd/c91662.CNG.swc

--2022-12-21 10:19:45--  https://raw.githubusercontent.com/neuronsimulator/resources/8b1290d5c8ab748dd6251be5bd46a4e3794d742f/notebooks/rxd/c91662.CNG.swc
Resolving raw.githubusercontent.com (raw.githubusercontent.com)... 185.199.108.133, 185.199.109.133, 185.199.110.133, ...
Connecting to raw.githubusercontent.com (raw.githubusercontent.com)|185.199.108.133|:443... connected.
HTTP request sent, awaiting response... 200 OK
Length: 55074 (54K) [text/plain]
Saving to: ‘c91662.CNG.swc’

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

Last-modified header missing -- time-stamps turned off.
2022-12-21 10:19:46 (6.13 MB/s) - ‘c91662.CNG.swc’ saved [55074/55074]

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.0)


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.0)


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")