This page was generated from docs/intro_spatial.ipynb. Interactive online version: Colab badge Deepnote badge. Binder badge Download notebook.

NeuroML and spatially detailed cells

The “NeuroML basics” article introduced the NeuroML framework and data format for expressing spiking neural networks. In the process, it explained how simplified neuron models, such as LIF and Izhikevich types, can be specified, parameterised and used to construct networks.

These simplified models are also often called abstract cells or point neurons, because they reduce the detailed anatomy and biophysical specificity found in nature, into just a few scalar values that produce a (roughly) equivalent response.

There is also another, more “direct” approach to simulate biophysics, by taking into account the spatial extent of the neuron and the immediate electro-chemical reactions that are supposed to take place along it. The “spatial” or “physically-based” formulation is thus composed of the shape (called morphology in neuro-anatomy) of the neuron, and the biophysical mechanisms that are present across the morphology. Under this formulation, most cell-internal biophysics properties (such as ion channels, concentration models, and physical membrane and cytosol properties) are modelled as being finely distributed over spans of neurite; whereas points of interaction with the neuron such as synapses and input probes are modelled as non-spatial (“point”) processes that are located at a specific spot on the neuron.

Aside: Once the biophysical effects are captured and quantified well enough through in-vitro and/or in-silico experiments, these effects can be added back to simplified neuron models: the aim is to preserve the effect while simplifying the dynamics equations and their computational requirements. These improved simplified neurons can be expressed in NeuroML with the help of LEMS, as shown in the relevant article.

This article introduces the additional attributes that spatially-detailed cell models have compared to point neuron models, outlines how these attributes can be stated through NeuroML, and shows the usual steps to construct spatially detailed neuron models, run them, and display and assess the simulation’s results.

Putting together a cell type

In NeuroML, spatially-modelled neurons are described through the generic-sounding <cell> tag, one for each neuron type. This tag must “contain” two other tags: a <morphology> tag that describes the shape of the neuron type, and a <biophysicalProperties> tag that describes the biophysics going on throughout the shape of the cell. Each of these tags will be explained in adequate depth in the following sections.

  • Tip: If each neuron has a different shape, each neuron will need to be in its own <cell> type.

When different cell types share a morphology or set of biophysics, or it’s just more convenient that way, the <morphology> and <biophysicalProperties> tags may have unique id attributes and be defined independently in the top-level <neuroml> tag. In this case, the <cell> tag may, instead of containing a <morphology> or <biophysicalProperties> tag inside it, have a morphology or biophysicalProperties attribute that names the stand-alone tag with associated id.

  • Note: If different <cell> types point to the same <biophysicalProperties> this way, the contents of that tag will be duly interpreted according to each <cell> type’s <morphology> (including <segmentGroup>s that may differ in each case).

Let’s see how the attributes of a spatially-modelled neuron are described with NeuroML; first the shape of the neuron, and then the biophysics that are present on each part of the morphology.

Peeking at the neuron model’s structure

Before diving into the details, to provide some context on what we’re about to discuss, let’s visualise a spatially detailed neuron model first. For that we’ll use explain_cell that will be further explained in the following.

A typical neuron model with high spatial detail is 𝙽𝙼𝙻𝙲𝙻𝟶𝟶𝟶𝟼𝟸𝟻: A Neocortical L6 basket cell from the NeuroML-DB. Feel free to search for and upload other models as well; a simpler stick-figure model is also suggested at the end of this section.



nmldb_cell_name = 'NMLCL000625'; zip_file_name = f'{nmldb_cell_name}.zip'
# Download the zip file with the model
import urllib.request
# Because NMLDB is having some issues with HTTPS, don't demand HTTPS verification
import ssl; ssl._create_default_https_context = ssl._create_unverified_context
urllib.request.urlretrieve(f'http://neuroml-db.org/GetModelZip?modelID={nmldb_cell_name}&version=NeuroML', zip_file_name)
# and unpack it
from zipfile import ZipFile
with ZipFile(zip_file_name, 'r') as zipp: zipp.extractall(nmldb_cell_name+'/')
# Find the cell's .nml file
import os
cell_filenames = [ nmldb_cell_name+'/'+name for name in os.listdir(nmldb_cell_name) if name.endswith('.cell.nml') ]
assert len(cell_filenames) == 1, "Didn't find a lonely cell nml file!"

An important thing to understand and remember is that when simulating across the body of a neuron, this is not done in perfect detail, for practical reasons. The shape of the neuron is discretised (sort of “chopped”) into chunks called compartments (or other terms, as seen on the terminology table.

To illustrate, we’ll paint each compartment differently and show both the shape of the modelled neuron, as well as how it is discretised for simulation. Right after, we’ll see how this neuron was built in NeuroML.



# Get some information about the cell, like shape and no. of compartments
import eden_simulator
cells_info = eden_simulator.experimental.explain_cell(cell_filenames[0])
# one cell in the model file, get it
cell_info = list(cells_info.values())[0]
# count the per-compartment information that explain_cell provides
nCellComps = len(cell_info['comp_midpoint'])

# Create some different colours to show the different compartments
import numpy as np
rng = np.random.default_rng(seed=20240901)
comp_colors = rng.random([nCellComps, 3])



# Show a neuron, painted by compartment
from eden_simulator.display.spatial import k3d as ek3d; import k3d
plot = Comps_plot = ek3d.Plot(grid_visible=False, menu_visibility=False, camera_auto_fit=True);
plot += ek3d.plot_neuron(cell_info, comp_colors)
# TODO zoom a bit closer
plot.show_html()

Try to read the NeuroML description of the cell (check for code-folding arrows at the start of nested lines, on your web browser or code editor). Scroll to near the end of the file and find that these endless lines of <segment>s and <segmentGroup>s have given way to a <biophysicalProperties> tag with only a few <channelDensity>s and other properties inside it. Don’t push yourself to understand all of the file just yet; all of its intricacies will be explained in depth by this article.

In case the above file is unwieldy to use (or out of interest), see also a simpler <cell> model with a stylised (or “stick figure”) shape, for granule cells of the cerebellar cortex; scroll down to Granule_98_3D.cell.nml and click on ( view NeuroML ) to read the NML description. As you’ll see by studying its <morphology> in the following (and also by clicking on Morphology on that page), this cell has a small spherical soma linked to a very long axon, bifurcated into a T shape. This axon is called a parallel fiber and it goes in the same direction relative to the Cb cortex, for all granule cells. Likewise, other types of neuron in the brain have different distinctive anatomical features and shapes.

Let’s see what NeuroML <cell>s are made of then.

The morphology of the cell

In neuro-anatomy, what’s called the morphology of a cell is all the spatial information about it. This is described in NeuroML with a <morphology> tag, which contains <segment>s that the cell is built of and <segmentGroups> which define areas of particular interest on the cell.

See also the official NeuroML guide on the following concepts and moving data to/from other technologies.

Segments

The morphology of a neuron, in NeuroML and virtually all simulators at the modelling level of networks with multi-compartment neurons, is formulated as a tree structure of branching neurite tubes. The tubes can have a varying diameter over their length; they can be broken down into truncated cones (specifically non-oblique frusta) that are joined together to form the neuron’s neurite tree.

Here is a sketch of a NeuroML-style frustum for reference, from Wikipedia:

A frustum.

A <segment>’s geometry can be defined as a pair of ends, where each end is a circular disc with a centre in space, and a diameter (or equivalently radius). The planes on which the discs lie are perpendicular to the line segment joining the two centres, for non-oblique frusta.

Within a <morphology>, the <segment>s are arranged in a directed tree structure, which means that all except for one are connected to a “parent” segment. The one that remains is called the “root” and is the “ancestor” of all other segments. By convention, the root is considered to be the soma (from which all processes emerge, after all), so that each “parent” segment is in turn located a little closer to the soma. Hence, for each segment we name its end that’s near the parent as proximal and the one that’s further from the parent as distal. (For the root segment, we can name the two ends either way.)

Of course, some questions arise when it’s time to model the soma as a tree of frusta. The usual practices are:

  • When the soma can be modelled as one single compartment and looks bulbous, it is allowed to be modelled as a sphere instead:

    • The segment’s two ends have the same diameter, and the distance between them is either none (as a convention that’s accepted to represent a sphere) or equal to the diameter (in which case, the resulting cylinder has the same (membrane) surface area as a sphere with the same diameter).

  • Alternatively, for an e.g. elongated soma, a modeller could slice it like a tuber along the longest axis, to better preserve the original shape.

  • When modelling a bipolar cell, modellers may opt to have neurites grow from the “proximal” end of a segment that represents the soma. More commonly for any type of cell, neurites may sprout from a segment near the middle of the soma.

Note

It is worth noting that the soma is not that nicely captured in this formulation, which may impact simulation accuracy within the soma. To correct for this, the biophysical factors over the soma are usually empirically adjusted to yield the correct “effective” values.

The neuron’s shape is then described inside the <morphology> as a list of <segment> tags like this:

<morphology>
    ...
    <segment id="0" name="Soma_0">
        <proximal x="0" y="0" z="0" diameter="20"/>
        <distal x="10" y="0" z="0" diameter="20"/>
    </segment>
    <segment id="10" name="Axon_1">
        <parent segment="0"/>
        <proximal x="0" y="0" z="0" diameter="15.6"/>
        <distal x="20" y="7" z="0" diameter="5"/>
    </segment>
    <segment id="11" name="Axon_2">
        <parent segment="10" fractionAlong="0.6"/>
        <distal x="12" y="9" z="0" diameter="5"/>
    </segment>
    ...

And the meaning of the attributes and the enclosed <proximal> and <distal> tags is as follows:

  • id: A distinct integer to name give the segment a name.

  • name: A more human-readable name to describe the segment with, but it is not actually used as part of the model; the id number is, instead.

  • Within the segment there are more tags:

    • Two tags called proximal and distal to define the extent of each frustum. Each has the following attributes:

      • x, y, z: The coordinates of the end’s centre in 3-d space, in microns.

      • diameter:The diameter of the frustum on that end, in microns.

    • For non-root segments, a <parent> tag specifying a the id of each’s parent segment. Said parent must have been declared previously in the <segment> sequence.

      • For more precision, the exact point on the parent segment to which this segment is attached, can also be specified. The point is given as the fractionAlong the parent segment’s length, from 0 (proximal “start” point) to 1 (distal “end” point).

        • The default value is 1 (attached to the parent’s distal end).

    • If the segment has a parent, defining its <proximal> end is optional; by default, it will be same as the point interpolated by fractionAlong, on the parent segment. the x, y, z as well as the diameter are interpolated between the two ends of the parent.

Standards in morphology coordinates

As expected, the 3-D coordinates for the proximal and distal ends of each segment must have an origin \((0,0,0)\); but what this is, is up to who traced the neurons. A common convention to place the origin near the middle of the soma.
Likewise, what the \(\overrightarrow x, ~ \overrightarrow y,~ \overrightarrow z ~\) directions represent is, unfortunately, up to the person who traced the neuron; there is no commonly agreed convention. Hence, you may have to tweak the coordinates to your favourite coordinate system:

  • Observe a discernible feature e.g. the axon;

  • then see which direction it extends to in the modeller’s coordinates;

  • see which direction it should be following in your coordinates, and you’ll understand how to map this axis.

  • If it’s important, repeat for the other two axes (e.g. with planarity, or asymmetrical shape) to find out how to transform the other axes.

Note: the origin used for the morphology is also used to (dis)place the cell <𝚒𝚗𝚜𝚝𝚊𝚗𝚌𝚎>s within a <𝚗𝚎𝚝𝚠𝚘𝚛𝚔>. See how this is done for example, on the Network tutorial and 𝙴𝚍𝚎𝚗𝙲𝚞𝚜𝚝𝚘𝚖𝚂𝚎𝚝𝚞𝚙 example.

Note: The NeuroML formulation is similar to very common SWC file format that describes morphologies. The additional features that NeuroML offers beyond SWC are:

  • The proximal end of a segment can be different from the parent’s distal; i.e. it can have different diameter and displacement in space.

  • The <segmentGroup> concept is much more powerful, in that the groups have names and descriptions (or metadata), and can contain and/or overlap with other groups.

  • The preferable discretisation of the cell can be provided, which is important for accurate and efficient simulation.

Specifying a location on a cell

Given the tree of <segment>s that outlines the shape of the neuron, a site on its surface can be designated as a NeuroML <segment> id, and fractionAlong the axial direction of said <segment>. This pair of segment and fractionAlong is used to:

  1. precisely attach synapses and stimulation probes to cells;

  2. record local action potentials with <EventOutputFile>s, and trajectories of quantities in <OutputFile>s — to specify fractionAlong together with the segment, use EDEN’s extended LEMS paths.

In all these cases:

  • If the segmentId is not specified, it is assumed 0 : typically the soma (but look out), or the entirety of a point neuron.

  • If the fractionAlong is not specified, it is assumed 0.5: the middle of the segment, doesn’t matter for point neurons. (or spatially detailed ones that comprise just one compartment, for that matter)

Aside: One could also ask for a specific depth into the cytosol and azimuth on the plane with the same fractionAlong, but NeuroML (and SNN simulators in general) hasn’t yet reached this level of spatial detail in the space of representable models (or there would be a way to model the soma in detail already).

Segment groups

The segments that build up the neuritic tree can be organised into groups that represent parts of the neuron. Groups are like mathematical sets in that they one can completely include another, or they can overlap at an intersection subset of segments. This also means that segments don’t have to be part of one group exclusively.

Segment groups thus demarcate interesting regions within a neuron, for looking at, or also relevant to our model.

In NeuroML, segment groups are defined within a <morphology>, as <segmentGroup> tags that include member segments, and/or other segment groups.

The <segmentGroup> tag has the following properties:

  • id: unique name for the group. Note: Unlike the ids of <segment>s, this can be any alphanumeral.

  • neuroLexId: A machine-readable NeuroLex keyword that identifies the type of the neuron’s region. For example, GO:0043025 stands for “soma”, GO:0043024 stands for “axon”, and GO:0097447 stands for “all the dendrites”.

    • The identifier sao864921383 has a special meaning for the implicit discretisation of the neuron, explained in a following section.

Each <segmentGroup> tag may also contain other tags, that list its properties:

  • <member segment="<seg id>"/>: Adds a single segment to this group.

  • <include segmentGroup="group id"/>: Adds a previously defined group to this group.

  • There also exist <𝚙𝚊𝚝𝚑> and <𝚜𝚞𝚋𝚃𝚛𝚎𝚎> tags to include segment spans in the standard, but EDEN doesn’t support them yet; they have not been used in practice. (Note: Segment groups can always be constructed programmatically to exactly fit the modeller’s intention.)

  • <property tag="<name>" value="<content>"> tags can also be used to add more information about the segment group. The tag="numberInternalDivisions" may be used to help with discretisation, when neuroLexId=sao864921383 “compartment”.

  • Finally, the <inhomogeneousParameter> tag can be used to associate a spatial scalar function with the group. This function can then be used to specify non-uniformity in ion channel distributions; hence the tag is explained in the relevant section.

The group "all"

If it is not otherwise defined by the modeller, the group "all" is defined as a shorthand that captures the entirety of the cell. It is, for example, useful to place mechanisms all over the cell, such as the membrane’s intrinsic electrical leak, or possibly an ion channel distribution that’s present throughout the cell’s membrane.

Simulation aspect: Discretisation into compartments

So far, we saw how the anatomical shape of the neuron is expressed in NeuroML. But as we saw at the start of the article, this shape is split (in mathematical language, discretised) into little bits, called compartments. The neuron is approximated by these compartments in the same way that an image can be approximated by pixels. There often are good reasons for the modeller to control how this discretisation happens; and NeuroML allows such control, albeit through more implicit than explicit rules.

Note: Although the exact discretisation method shouldn’t matter in theory, it is often important in practice. Naturally, a too big (or less often, too small) compartment may provide a poor reproduction of the ideal “real” thing, as we’ll see in a following exercise. But there are models where compartments (often containing certain sensitive biophysical mechanisms) have to have an exact size to reproduce the phenotype!

The implicit rules of NeuroML for discretisation are as follows:

  • <segmentGroup>s tagged with neuroLexId="sao864921383" follow “cable” rules;

  • Each <segment> that doesn’t belong in such a corresponds to one individual compartment.

    • This is how discretisation is always done, for example, in the GENESIS simulator.

Tip: If your modelling pipeline already offers its own discretisation (and it’s not as easy to convert as with NEURON), consider splitting your <segment>s (or equivalent) at compartment boundaries, and then joining the segment parts of each compartment into “cable” <segmentGroup>s (see below). If your modelled compartments include branching points within them, consult also the comp_area, comp_capacitance and comp_conductance_to_parent information that explain cell provides, to make equivalent NeuroML-supported compartments; look out for the influence of altered <segment>s on their surface area and electrical tortuosity.

The unbranched section directive

NeuroML offers an implicit rule for discretising spans of unbranched neurite cable. The rule applies when the span is gathered into a <segmentGroup> with the attribute neuroLexId="sao864921383". In that case, the length of all contained segments is added up and divided among a number of compartments, so that each compartment takes over an equal fraction of the total cable length. The number of compartments can be specified by adding inside the same <segmentGroup> the tag <property tag="numberInternalDivisions" value="<number>">.

  • If neuroLexId="sao864921383" and this tag is not specified, the whole group is assumed to map to one compartment (i.e. numberInternalDivisions = 1).

See also these tickets:

Tip

neuroLexId="sao864921383" corresponds to the <cable>s of the previous incarnation of NeuroML, NMLv1. It facilitates conversion of models that were originally written for the NEURON simulator, which shares the concept of continuous sections of neurite tube making up the neurons.

To prevent ambiguity, “cable” <segmentGroup>s are not allowed to overlap.

Refer also to NEURON’s user guide and forum on wise choices for numberInternalDivisions.
But also keep in mind that uniform spacing could be a suboptimal discretisation method if the neurite diameter, or parameters of active mechanisms, vary wildly. Take care when building models, and consider fine-tuning the discretisation around the important points.
As a final note: For some models and experiments, the discretisation can also be reduced to a smaller number of compartments (see this early paper and other advances since). This will affect the neuron’s dynamics (by more or by less) and drastically cut down the time it takes to run a simulation. Check if it’s applicable for your use case; but also be aware of the trade-offs, and the points where the decision will have to change.

Retrieving and using the discretisation

After specifying the discretisation, getting it back can also prove useful. Typical applications are:

EDEN can help with extracting the discretisation of neurons, through the auxiliary explain cell feature, as seen above. This feature produces a breakdown of cell types into the compartments that make them (which may have a many-to-many relationship with <segment>s in general), and the geometrical and physical properties of each compartment such as area and electrical capacitance of the membrane, and distance from the soma along the neurite path. It also produces additional compartment-related information, such as the list of compartments that belong to each segment group and a 3-D polygon mesh that represents the neuron’s morphology, and can be coloured per individual compartment.

Visualisation example

Now that we know what <segment>s and <segmentGroup>s are, let’s see the neuron that was painted per compartment before, painted here separately for each main region. Note that all dendrites happen to be classified as “basal” in this cell model.



# For each group: NML name, color, description
group_info = [
    ('somatic', (0.8, 0.7, 0.4), 'Soma'),
    ('axonal' , (0.8, 0.2, 0.2), 'Axon'     ),
    ('basal'  , (0.5, 0.1, 0.5), 'Dendrites'), ]

# First, set a color for comps that are in none of these groups.
comp_colors[:,:] = [[0.5,0.5,0.5]]
# Then, paint the compartments for each group.
for group_name, color, _ in group_info:
    comps_in_group = cell_info['segment_groups'][group_name]['comps']
    comp_colors[comps_in_group] = color

# Show the neuron:
plot = Groups_plot = ek3d.Plot(grid_visible=False, menu_visibility=False);
plot += ek3d.plot_neuron(cell_info, comp_colors)

# And add some annotations.
for i, (_, color, description) in enumerate(group_info):
    plot += k3d.text2d(description, (0,0.9-0.15*i), color=ek3d.RgbToInt(color),
                      is_html=True, label_box=False, size=2.5)

plot.show_html()

The biophysics of the cell

After specifying everything about the shape and discretisation of a neuron, it’s time to specify the biophysical properties that are present in (on, and around) it in our model. In NeuroML, these can be classified as:

In NeuroML, all these are defined within a <biophysicalProperties> tag. This contains the <membraneProperties>, <intracellularProperties> and <extracellularProperties> tags where each in turn contains some of the properties (as more appropriate location-wise).

Until a later revision of NeuroML includes them, other properties may still be model-able today as LEMS components that are designed to apply the desired effect.

All these are described for a cell type as distributions of parameters and mechanisms. The distributions apply all over the cell, or are bounded by a <segment> or <segmentGroup>, referring to the <morphology> that’s already defined for the <cell> type.

  • Each of the tags discussed in the following may have a respective segment or segmentGroup attribute; if neither is specified, it applies for the group all.

Caution

Make sure you don’t miss parts of the intended area when distributing over disparate segment groups; if you want a visual check, use the code above to verify.

Fun fact: The user guide of the simulator Arbor calls the process of placing mechanisms on neurons, “decoration” and specifically distributing attributes over the cell’s morphology, “painting”; this author approves of those metaphors.

Electrotonic properties

Some general electrical properties of the neuritic cable making up the neuron can be provided with the following tags:

  • <specificCapacitance value="(capacitance per area)"/>: The local specific capacitance of the membrane.

  • <resistivity value="(resistance times length)"/>: The local resistivity of the cytoplasm. (assuming current flows through the whole cross-section equally)

  • <initMembPotential value="(voltage value)" />: The initial value for the cross-membrane voltage \(V_m\) when the simulation starts.

  • <spikeThresh value="(voltage value)"/>: The threshold that \(V_m\) need to exceed, for an action potential to be registered in the model. This is also the point in time (plus possible delay) when event-based synapses are activated. (See also the comparison between “classic” and graded synapses.)

    • Note: <spikeThresh> is optional, but spikes will not be registered from the areas where it isn’t specified!

Except for <resistivity> which is located within <intracellularProperties>, the other properties are located within <membraneProperties>.

  • Important: These properties may be redefined by tags with overlapping areas. In this case, each spot is assigned the value of the last tag that applies to it, in sequence within the <biophysicalProprties> tag. This is similar to how EDEN’s <𝙴𝚍𝚎𝚗𝙲𝚞𝚜𝚝𝚘𝚖𝚂𝚎𝚝𝚞𝚙> extension works.

Note: Passive electrical leakage of the membrane is specified in NeuroML as a trans-membrane mechanism, typically a <channelDensity> of and <ionChannel> in the degenerate case (lacking any “<gate>”s).

Note

The dynamics of a spatially detailed cell are that of the passive neural cable, combined with the dynamics of the various mechanisms present on different parts of the cell.

Adding ion channels and other leaks to the membrane

More than simple physical properties, dynamic mechanisms may also be present over the neuron’s morphology. The most common type of trans-membrane mechanism is ion channels. These “open” and “close” for specific ion species (or sometimes, any form of current) and their state is influenced by electro-chemical factors, plus the state they were at right before. The constituent parts of the ion channel that participate in opening and closing often called “gates”.

Since a single ion channel is a quite small entity (compared to compartments of the neuron), populations of ion channels are usually modelled as distributions finely spread over areas.

For the same reason, the ion channels’ dynamics are usually modelled at the macroscopic (or smooth) level. In this level there are fractional gate variables that represent the average state of a gate population within the channel population, and there are continuous-timeODE‘s that govern the average rates of change between the “open” and “closed” states. (More faithful, stochastic models may still be implemented, but at this point in NeuroML they’ll be more practical to model as stochastic processes for the macroscopic gate variables.)

Designing ion channels

NeuroML describes ion channels through the <ionChannel> tag, and includes a multiplied-gate formulation of ion channel models for modellers to work with. (Do keep in mind that more formulations than those out of the box can be described, using LEMS).

All ion channel tags may have the following optional attributes:

  • type: Can be used in place of the tag’s name itself, to identify the particular type of ion channel; EDEN uses that or the tag name to use LEMS <ComponentType>s as ion channels.

  • species: The sort of substance (typically an ion) that passes through the channel. If missing, the ion channel is assumed to permit non-specific current; that is, current flows without really affecting the local concentrations of tracked ion species.

  • conductance: The conductance of a single ion channel. Must be set for a population of individual ion channels rather than a coarse density.

    • Note: The dynamics are the same either way, the “population” is still assumed to be lumped into one aggregate mechanism instance (per compartment) in NeuroML.

Unless specially modelled through LEMS, the <ionChannel> tag contains <gate>s of various types, each with its multiplicity and dynamics as described in the following:

Note: Constant leaks are also modelled in NeuroML as “degenerate” <ionChannel>s with no gates at all (thus a fixed conductance).

In NeuroML, an ion channel’s dynamics are modelled separately from the local reversal potential and the density of its instances that are spread over cells; said parameters are modelled as part of the particular channel distribution applied on cells.

Note that in the wild there may be more specific ion channel tag names such as <ionChannelPassive>, <ionChannelHH> and <ionChannelKS>; EDEN will parse and handle all of the contents regardless of stated type, so that potential hybrid models mixing up the gate types are also supported.

Refer also to the official NeuroML guide and the LEMS extension section of this guide on how to make your own ion channels in NeuroML.

In the following, the transition rates that govern the kinetics of gate variables are introduced, and in turn the various gate formulations that use these transition rates.

Transition rates

As mentioned above, the transition between states of individual gates can (in a big enough population) be macroscopically modelled as the fraction of gates that are in each state, and the rates of transition (as fraction per time, or \(T^{-1}\)) between states. These rates can be expressed in NeuroML as either “transition rate per existing fraction” functions (also seen as Hodgkin-Huxley \(\alpha, \beta\) “rate constants”). In chemical language, transition between two states can be expressed as the reversible reaction scheme:

\[\mathrm{Closed} ~ \xrightleftharpoons[\alpha]{\beta} ~ \mathrm{Open}\]

Another, equivalent way to define these rate functions is through the quantities “steady-state value to converge to” and “time-constant to converge to it”, which may written as \(x_\infty, \tau\). Often this is easier to get as a measure (say, from literature), or just more convenient to write down as a formula. These two are directly convertible to and from \(\alpha, \beta\) functions through the relations:

\[\tau = \frac{1}{\alpha + \beta},~ x_\infty = \frac{\alpha}{\alpha + \beta} \iff ~ \alpha = \frac{x_\infty}{\tau},~ \beta = \frac{1 - x_\infty}{\tau}\]

Keep in mind that \(\alpha, \beta\) (or \(\tau, x_\infty\) just the same) will (as pairs) vary against external factors, such as trans-membrane voltage. (Naturally, if the rates couldn’t change then ion channels would arrive and stay at a resting equilibrium, without oscillating.)

Rate formulae

NeuroML includes a small set of pre-defined rate formulae types as functions of membrane voltage only, in terms of: rate (for \(\alpha, \beta\) functions), steady state (for \(x_\infty\)) and time constant (for \(\tau\)). These are used as in <forwardRate>s, <reverseRate>s, <timeCourse>s and <steadyState>s of Hodgkin-Huxley-style gates, and forward, reverse and tauinf Transitions of Markov-style gates rates as applicable.
These pre-defined function types are:

  • 𝙷𝙷𝙴𝚡𝚙𝚅𝚊𝚛𝚒𝚊𝚋𝚕𝚎 and 𝙷𝙷𝙴𝚡𝚙𝚁𝚊𝚝𝚎, which are parameterisable exponential curves. The general formula is: \begin{gather*} f(V_m) = rate \cdot e^{\frac{V_m - midpoint}{scale}} \end{gather*}

    Their attributes are:

    • rate: The dependent-variable scaling factor for the quantity represented by the curve, from the dimensionless inner curve to the function’s dimensions. Is the value of the function when membrane voltage is midpoint.

    • midpoint: The reference point to center the curve around. The function has the value rate for this value of membrane voltage.

    • scale: The independent-variable scaling unit to evolve the curve by. Values scale apart on the independent variable are in a proportion of \(1\colon e\).

      • Notice that midpoint and rate are redundant for the exponential curve, set them in the way that most makes sense.

  • 𝙷𝙷𝙴𝚡𝚙𝙻𝚒𝚗𝚎𝚊𝚛𝚅𝚊𝚛𝚒𝚊𝚋𝚕𝚎 and 𝙷𝙷𝙴𝚡𝚙𝙻𝚒𝚗𝚎𝚊𝚛𝚁𝚊𝚝𝚎, which are parameterisable rectified linear curves with exponential taper to zero. The general formula is:

    \[\begin{gathered} f(V_m) = rate \cdot \frac{a} {1 - e^{-a}}, ~~ a = \frac{V_m - midpoint}{scale} \end{gathered}\]
    Their attributes are:

    • rate: The dependent-variable scaling factor for the quantity represented by the curve, from the dimensionless inner curve to the function’s dimensions. Is twice the value of the function when membrane voltage is midpoint.

    • midpoint: The reference point to center the curve around. The function has the value rate/2 for this value of membrane voltage.

    • scale: The independent-variable scaling unit to evolve the curve by. The extent of the “taper zone” from linear to zero is proportional to this value.

  • 𝙷𝙷𝚂𝚒𝚐𝚖𝚘𝚒𝚍𝚅𝚊𝚛𝚒𝚊𝚋𝚕𝚎 and 𝙷𝙷𝚂𝚒𝚐𝚖𝚘𝚒𝚍𝚁𝚊𝚝𝚎, which are parameterisable logisticcurves. The general formula is: \begin{gather} f(V_m) = rate \cdot \frac{1} {1 + e^{-a}}, ~~ a = \frac{V_m - midpoint}{scale} \end{gather}

    Their attributes are:

    • rate: The dependent-variable scaling factor for the quantity represented by the curve, from the dimensionless inner curve to the function’s dimensions. Is twice the value of the function when membrane voltage is midpoint.

    • midpoint: The reference point to center the curve around. The function has the value rate/2 for this value of membrane voltage.

    • scale: The independent-variable scaling unit to evolve the curve by. The extent of the “taper zone” between zero and maximum is proportional to this value.

  • 𝚏𝚒𝚡𝚎𝚍𝚃𝚒𝚖𝚎𝙲𝚘𝚞𝚛𝚜𝚎, which represents a fixed time constant to converge to the steady state by. The only attribute is tau, the value of the time constant (with units).

More then the pre-defined set, additional rate formulae can be supplied using LEMS.

Using these rate formulae, we can specify various types of ion channel gate dynamics.

General attributes of gates
In the pre-defined NeuroML formulation, <𝚒𝚘𝚗𝙲𝚑𝚊𝚗𝚗𝚎𝚕> models contain one or more gates. The fraction of ion channels that are open in an aggregate population is assumed to be the multiplicative product of the fractions of gates that are open, where the fraction of each gate that is open is the fraction of gate particles, multiplied over by the number or particle instances it takes to actually open a gate.
Every sort of <gate> in an <ionChannel> must specify the following attributes:

  • identifying name for the gate and its activation variable;

  • multiplicity as instances;

    • The gate’s activation is multiplied by itself this many times over, to scale the activation of the ion channels it belongs to (getting multiplied with the other gates in it as well).

  • and if the tag name is <gate>, an attribute that clarifies the specific dynamics type of the gate.

Classic Hodgkin-Huxley gate dynamics

The most common model for a gate assigns two possible states for it: either Open or Closed. In that case, one can consider \(x\), the fraction of gates of a certain type in an ion channel population that are open; it follows that \(1 - x\) of all such gates are closed. Hence, the state of that gate population is captured in just one number x, commonly called the “activation” or “gate variable”.

The transition rates between these two states can be expressed as \(\alpha, \beta\) functions, or alternatively \(\tau, x_\infty\) (or a mix of them!). The tag names for a gate under this formalism can be <𝚐𝚊𝚝𝚎𝙷𝙷𝚛𝚊𝚝𝚎𝚜>, <𝚐𝚊𝚝𝚎𝙷𝙷𝚝𝚊𝚞𝙸𝚗𝚏>, <𝚐𝚊𝚝𝚎𝙷𝙷𝚛𝚊𝚝𝚎𝚜𝚃𝚊𝚞>, <𝚐𝚊𝚝𝚎𝙷𝙷𝚛𝚊𝚝𝚎𝚜𝙸𝚗𝚏>, and even <𝚐𝚊𝚝𝚎𝙷𝙷𝚛𝚊𝚝𝚎𝚜𝚃𝚊𝚞𝙸𝚗𝚏> so that rates stands for \(\alpha, \beta\), tau stands for \(\tau\), inf stands for \(x_\infty\), and which functions are specified is contained in the tag name.

As explained above, each of the tags may contain the transition rate formulae as child tags:

  • \(\alpha\) as <forwardRate> and \(\beta\) as <reverseRate>, each with a type attribute referring to a per time-valued function;

  • \(\tau\) as <timeCourse>, with a type attribute referring to a time-valued function;

  • \(x_\infty\) as <steadyState>, with a type attribute referring to a dimensionless function.

If either \(\tau\) or \(x_\infty\) are specified, they are consulted to simulate the dynamics; otherwise, their values are derived from the \(\alpha, \beta\) formulae.
This means that the \(\alpha, \beta\) rates specified in the model may well be inconsistent with \(\tau\) and/or \(x_\infty\)! To prevent this, avoid the <gateHHratesTauInf> tag and prefer the <gateHHrates> and <gateHHtauInf> tags to describe the dynamics.
Fractional HH gate dynamics

Sometimes, gates are modelled as individual and independent sub-populations, where each “sub-gate” has Hodgkin-Huxley dynamics. Since the subpopulations are independent, they contribute additively (not multiplicatively!) to the total activation of the gate. The NeuroML tag to describe this is <𝚐𝚊𝚝𝚎𝙵𝚛𝚊𝚌𝚝𝚒𝚘𝚗𝚊𝚕>. It contains multiple <subGate>s, each with attributes id(entifier) and fractionalConductance (weight by which they contribute to the <gateFractional>’s activation). Each <subGate>is then modelled with the same attributes and semantics as a Hodgkin-Huxley-style <𝚐𝚊𝚝𝚎𝙷𝙷𝚝𝚊𝚞𝙸𝚗𝚏>.

Note: It may make more sense if the total of a <gateFractional>’s <subgate>s adds up to 1.

Multi-state Markov gate dynamics
Another way to model a gate population is to assume more than two possible (exclusive) states for each individual gate. In that case the state of a gate population can be described by the fractions of gates per state they’re in — but same as with Hodgkin-Huxley modelling, there is one degree of freedom less since all the fractions per state must add up to 1. The various pathways for gates to (statistically) move between states are then laid out as a general kinetic scheme (or continuous-time Markov process).
Tip: Non-Markov kinetic schemes can also be modelled just fine, by specifying non-stationary transition rates with LEMS.

The NeuroML tag to describe this is <𝚐𝚊𝚝𝚎𝙺𝚂>. It contains multiple <openState>s <closedState>s, each with an identifier name attribute. It also contains multiple transition paths between these states, as the tags:

  • <forwardTransition>, which contains a <rate> tag that has same attributes as the <forwardRate> in HH-style gate dynamics.

  • <reverseTransition>, which contains a <rate> tag with same attributes as <reverseRate>.

  • <TauInfTransition>, which expresses the rates in a different way (that would really make sense for HH dynamics). It contains the <timeCourse> and <steadyState> tags with the same semantics as in <gateHHtauInf>, but the true transition time constant and dynamic equilibrium between from and to also depend on all the other (direct and indirect) transition pathways as well. Use the \(\tau, x_\infty\) to \(\alpha, \beta\) conversion formulae to retrieve the forward and reverse rates that <TauInfTransition> stands for.

Each Transition must specify from which state it goes to which as attributes. For each <forwardTransition> between two states there should be a <reverseTransition> and vice versa, and multiple transitions between a pair of states are not allowed (sum then up instead).

  • Even if a transition goes in one direction only, this should be made explicit by specifying a <reverseTransition> with a <rate> of constant 0.

Examples

Take a look at the NeuroML description of the following gating kinetic scheme from Hirschberg et al. 1998, as shown in figure 10 of this paper:

\[\begin{split}\begin{alignedat}{4} & C_1 \rightleftharpoons && C_2 \rightleftharpoons && C_3 \rightleftharpoons && C_4 \\ & && && {\:\upharpoonleft\!\downharpoonright} && {\:\upharpoonleft\!\downharpoonright} \\ & && && O_1 && O_2 \end{alignedat}\end{split}\]

Compare and contrast with its description as a MOD file for NEURON.

Since the rate functions aren’t part of the pre-defined rate formulae, they have been described with LEMS. Note that all parameters are repeatedly defined for all rates, even though each rate actually cares for only one specific rate constant.

  • The interested reader may thus slim down the definitions of the LEMS rates, and even further, whittle the <ComponentType>s down to one calcium-dependent and one fixed transition rate, as seen with the second example below.

Another ion channel modelled with a <gateKS> is from Carter et al. 2012 figure 7a, with 12 separate states that arise from 6 mechanical states that in sequence open the ion channel, and a potential inactivation sub-state for each mechanical state. Note that all these rates follow one of just two base formulae, with only different parameter values changing — that’s much neater!

Exercise

Convert the Markov ion channel of from the MOD file description (or other cases of the kinetic model for NEURON, to a NeuroML description. Which description is easier to read and modify, in which use cases? Use the <q10Settings> tag to model the identical effect of temperature on all transitions.

Instantaneous gate dynamics

Finally, there are some ion channel gates that open or close instantaneously (or quick enough, in any case) in response to their surroundings. Hence they don’t carry an intrinsic state variable that affects their evolution moving forward; they are stateless.

The behaviour of such gates is described in NeuroML as a <𝚐𝚊𝚝𝚎𝙷𝙷𝙸𝚗𝚜𝚝𝚊𝚗𝚝𝚊𝚗𝚎𝚘𝚞𝚜> tag; this contains a <steadyState> tag with the type attribute that selects the dimensionless formula for the instantaneous activation fraction, and its corresponding parameters as attributes.

Gate rate scaling
As chemical reactions often do, an ion channel gate’s transitions may be accelerated or slowed down by the local temperature (or possibly other factors in the model). For convenience, NeuroML allows adding a <q10Settings> tag on any non-static type of <gate> (or <subGate>), which adds a temperature-dependence Q₁₀ factor to speed up the gate by (effectively scaling \(\tau\)).
The attributes of this tag may be either:

  • q10Factor (unitless, per 10°C) and experimentalTemp that the gate’s recorded dynamics refer to (i.e. the temperature at which the gate behaves the same as without the q10Settings mentioned).

  • or fixedQ10 (unitless scaling factor) for a reaction speed that actually stays the same within the relevant temperature zone of the model (or was calculated outside the model’s description, e.g. when generating a model file and having pre-calculated rates).

Conductance scaling

As we saw above, the open-close dynamics of ion channel gates may be affected by temperature and NeuroML offers a way to model that. Another effect of temperature that’s often encountered is that the conductance of certain ion channels may be altered as well. To capture this effect, NeuroML offers the <q10ConductanceScaling> tag, which can be inserted to <ionChannel>s and has attributes q10Factor and experimentalTemp, with the same meaning as for temperature-affected gate transition rates.

Refer also to how custom temperature-dependence models that may be affected by more parameters can be modelled with LEMS.

More dynamics models

Finally, there may be ion channel models (especially ones proposed for the first time) that don’t fully fall within one of the predefined formulations. If their dynamics allow, these can still be modelled as custom LEMS components, and also using EDEN’s extensions to NeuroML if necessary.

Ion channel distributions

After designing an (aggregate) ion channel mechanism, it’s time to place instances of this mechanism on cells. This is done in NeuroML through distribution tags that designate a type of ion channel, an extent over the cell, and parameters of the distribution, such as condDensity and erev(ersal). (Much like the tags for electrical biophysical properties do.)

There are many different distribution types in NeuroML (refer to the NeuroML reference and search the page for channelDensity and channelPopulation); what they all have in common is the following attributes:

  • ionChannel, which points to the ideal ion channel mechanism to replicate all over the place (or simply scale by conductanceDensity * compartment area for the usual macroscopic models).

  • ion (optionally) which designates the specific ion species to be pooled. May be non_specific when referring to an “<ionChannel>" lacking the ion property. Useful then there are e.g. two separate “pools” of calcium which activate certain mechanisms differently for some reason.

  • segment or segmentGroup (either optionally) which designates the part of the cell to add the ion channels to; all if neither attribute is present.

Note: for <channelPopulation>s which involve a specific number of ion channels, the per-individual conductance must have been specified for the <ionChannel>.

Caution

If two distributions with the same <ionChannel> name are defined for overlapping <segmentGroup>s, it is ambiguous whether the distributions should both be instantiated, or one should replace the other, in the zone of overlap. Make sure that distributions of the “same” ion channel mechanism do not overlap, to to avoid ambiguity!

Non-uniform ion channel distributions in NeuroML

The previous way to define distributions of mechanisms specifies, for each parameter, the same value all over the distribution; but some times we need more that that. Of course, one could always use a fine-grained set of tiny <segmentGroup>s and define a different distribution for each one of them for each non-uniform distribution, but this is cumbersome to write, read, store and use. NeuroML offers some provisions for non-uniformity through the <channelDensityNonUniform>, <channelDensityNonUniformNernst> and <channelDensityNonUniformGHK> tags for ion channel distributions.

Tip

For a fully general (albeit EDEN-specific) way to control the variability of all mechanisms on spatially-modelled cells, refer to the relevant usage example of the <EdenCustomSetup> extension.

Inhomogeneous value in the channelDensity

Instead of a segment or segmentGroup attribute, the aforementioned tags contain one or more <variableParameter> tags that define variability over segmentGroups, each through a relevant <inhomogeneousParameter> in the <segmentGroup> and <inhomogeneousValue> mathematical function of said parameter. The <variableParameter> tags have the following attributes:

  • segmentGroup: to distribute ion channels over, and search for <inhomogeneousParameter>s within;

  • parameter: The biophysical attribute to vary over the ion channel’s distribution; this can be condDensity for <channelDensityNonUniform> and <channelDensityNonUniformNernst>, and permeability for <channelDensityNonUniformGHK>.

And inside each <variableParameter> tag, there is a <inhomogeneousValue> tag with attributes:

  • inhomogeneousParameter: The id of the <inhomogeneousParameter> specified in the <segmentGroup>

  • value: A dimensionless LEMS expression that involves the <inhomogeneousParameter> by its variable name. Note: The numerical value of the expression stands for SI units, such as mhos per square meter, \({S}/{m^2}\) for condDensity and metres per second, \(m/s\) for permeability.

Inhomogeneous parameters in the morphology

To support mechanism distributions with non-uniform parameters, non-uniform dimensionless parameters can be defined over a <morphology> by inserting <inhomogeneousParameter> tags in it. These have the following attributes:

  • id: Must be unique within the <segmentGroup>.

  • variable: The actual variable name to be used in LEMS expressions referencing the <inhomogeneousParameter>.

  • metric: The name of the metric to give the variable values by, over the extent of the <segmentGroup>. The only permitted metric for now is Path Length from root. Note: Unlike most LEMS quantities, this variable is provided dimensionless and its value is in microns (if not made relative with normalizationEnd).

The following optional tags inside the <inhomogeneousParameter> apply transformations to the metric variable:

  • <proximal translationStart="<value>"/>: Offset the metric so that the minimum is value (typically 0).

  • <distal   normalizationEnd="<value>"/>: Divide the metric my its maximum and then multiply by value (typically 1), so that the variable is now a relative, dimensionless measure.

Adding concentration models near the membrane

Another type of local “mechanism” that can be present on (parts of) the neuron is concentrations of chemical species (typically ions) but near the cell’s membrane. The presence (and relevance) of such concentrations can be expressed in NeuroML through <species> tags inside <intracellularProperties> (i.e. along with axial <resistivity> tags).
Each <species> tag has the following attributes:

  • id: Name of the “distribution” of the concentration model. Useful when different concentration models are used for the same species.

  • segment or segmentGroup which designates where to add the concentration model to; all if neither attribute is present.

  • ion: Name of the substance whose concentration is being tracked (in theory, it shouldn’t be just ions though).

  • initialConcentration: The concentration of the chemical species near the inner side of the membrane, when the simulation starts.

  • initialExtConcentration: The concentration near the outer side of the membrane, when the simulation starts. This usually is kept steady, to model a well mixed buffer when neurons are isolated on a dish; but new <concentrationModels> can always be have more interesting dynamics.

The concentration model, i.e. relationship between incoming flux and local concentration over time for a chemical species, is described by NeuroML tags at the top level of the NeuroML files. NeuroML offers the following types of concentration models outside the box:

Same as with most NeuroML mechanisms, novel concentration models can also be modelled with LEMS.

Note

The total species influx arriving to the local concentration model, is the sum of all <ionChannel> distributions transporting the same chemical species.

Attaching probes and synapses on a spatially-modelled cell

As we saw in the Introduction section, mechanisms can also be individually attached to cells to complete the models, through <projection>s and <inputList>s in the <network>.

The exact spot to place these mechanisms on the cell is specified for:

  • <input>s in the <inputList>s, as <segment> and <fractionAlong>

  • <connection>s in the <projection>s, as {pre,post}Segment and {pre,post}FractionAlong attributes (one pair for each of the pre- and post-synaptic sites).

Recap

These are the things which make up a spatially-detailed biophysical neuron model: if you’ve managed to read this far, you now know enough to read models of neurons, as well as make your own.

Armed with this knowledge, look again at the .cell.nml files shown at the start. If you see anything that’s still not clear now, ask about it on the NeuroML discussion channels and tell the guide’s authors to include it in this article.

Example: Modelling and simulating a stretch of neural cable

As a first tutorial, we’ll implement in NeuroML and play with the simplest spatially-detailed neuron model: an unbranched section of neural cable (that could be part of the e.g. axon), with a spike that propagates along it following the Hodgkin-Huxley type, Na and K ion channel dynamics.

Modelling the cell’s morphology

The morphology of a straight cylinder is rather simple. We’ll define only one cylindrical <segment> and make use of neuroLexId="sao864921383" and numberInternalDivisions to have it chopped into finer compartments.

  • Note that, as explained above, the same morphology could be made by constructing individual segments and letting the default discretisation of 1 compartment per segment take hold.



cable_length_microns = 200
cable_diameter_microns = 1
cable_compartments = 50
cable_morph_file = f'''<neuroml>
<morphology id="Subdivided_Morph">
    <segment id ="0" name="Axon_0">
        <proximal x="0" y="0" z="0" diameter="1"/>
        <distal x="{cable_length_microns}" y="0" z="0" diameter="{cable_diameter_microns}"/>
    </segment>
    <segmentGroup id="axon" neuroLexId="sao864921383">
        <property tag="numberInternalDivisions" value="{cable_compartments}"/>
        <member segment="0"/>
    </segmentGroup>
</morphology>
</neuroml>'''
# print(cable_morph_file)
with open('Subdivided_Morph.nml', 'wt') as f: f.write(cable_morph_file)

Modelling the membrane mechanisms

Next, we’ll add the ion channels from the famous Hodgkin-Huxley model, as also seen on the official NeuroML guide. Note that more than by typing, NeuroML files can also be created with pynml as shown in the linked example, and any other suitable method.

We’ll make one .channel.nml file for each independent trans-membrane mechanism.

Sodium channel



%%writefile HH_Sodium.channel.nml
<neuroml>
<ionChannel id="na_chan" type="ionChannelHH" conductance="10pS" species="na">

    <gateHHrates id="m" instances="3">
        <q10Settings type="q10ExpTemp" q10Factor="3" experimentalTemp="6.3 degC"/>
        <forwardRate type="HHExpLinearRate" rate="1per_ms" midpoint="-40mV" scale="10mV"/>
        <reverseRate type="HHExpRate" rate="4per_ms" midpoint="-65mV" scale="-18mV"/>
    </gateHHrates>

    <gateHHrates id="h" instances="1">
        <q10Settings type="q10ExpTemp" q10Factor="3" experimentalTemp="6.3 degC"/>
        <forwardRate type="HHExpRate" rate="0.07per_ms" midpoint="-65mV" scale="-20mV"/>
        <reverseRate type="HHSigmoidRate" rate="1per_ms" midpoint="-35mV" scale="10mV"/>
    </gateHHrates>
</ionChannel>
</neuroml>

Writing HH_Sodium.channel.nml

Potassium channel



%%writefile HH_Potassium.channel.nml
<neuroml>
<ionChannel id="k_chan" type="ionChannelHH" conductance="10pS" species="k">

    <gateHHrates id="n" instances="4">
        <q10Settings type="q10ExpTemp" q10Factor="3" experimentalTemp="6.3 degC"/>
        <forwardRate type="HHExpLinearRate" rate="0.1per_ms" midpoint="-55mV" scale="10mV"/>
        <reverseRate type="HHExpRate" rate="0.125per_ms" midpoint="-65mV" scale="-80mV"/>
    </gateHHrates>

</ionChannel>
</neuroml>

Writing HH_Potassium.channel.nml

Passive membrane leak

Note that we only need the conductance attribute defined if we’re using distributions of a counted number of channels, such as <𝚌𝚑𝚊𝚗𝚗𝚎𝚕𝙿𝚘𝚙𝚞𝚕𝚊𝚝𝚒𝚘𝚗>.



%%writefile Passive_leak.channel.nml
<neuroml>
<ionChannelHH id="passiveChan">
            <notes>Leak conductance</notes>
</ionChannelHH>
</neuroml>

Writing Passive_leak.channel.nml

Modelling the biophysics distributions

Then, we’ll apply the ion channels over the whole cable, with erev, condDensity and other biophysical parameters again as per the example in the NeuroML guide.



%%writefile Subdivided_Cable.cell.nml
<neuroml>
<!-- Get the previous definitions -->
<include href="Subdivided_Morph.nml"/>
<include href="HH_Sodium.channel.nml"/>
<include href="HH_Potassium.channel.nml"/>
<include href="Passive_leak.channel.nml"/>

<!-- Set the biophysics -->
<biophysicalProperties id="bioPhys_Passive">
    <membraneProperties>
        <channelDensity id="na_all" ionChannel="na_chan"     ion="na"           condDensity="0.12   S_per_cm2" erev="+50.0 mV" />
        <channelDensity id="k_all"  ionChannel="k_chan"      ion="k"            condDensity="0.036  S_per_cm2" erev="-77.0 mV" />
        <channelDensity id="leak"   ionChannel="passiveChan" ion="non_specific" condDensity="0.0003 S_per_cm2" erev="-54.3 mV" />
        <spikeThresh value="0 mV"/>
        <specificCapacitance value="1.0 uF_per_cm2"/>
        <initMembPotential value="-65mV" />
    </membraneProperties>
    <intracellularProperties>
        <resistivity value="1 kohm_cm"/>
    </intracellularProperties>
</biophysicalProperties>

<!-- And compose into a cell type -->
<cell id="Subdivided_Cable" morphology="Subdivided_Morph" biophysicalProperties="bioPhys_Passive"/>
</neuroml>

Writing Subdivided_Cable.cell.nml

Modelling the experimental rig

We’ll put a current clamp on some \(3/10\) along the cable, to see how the induced spike will propagate along the cable towards both directions.



cable_pop_name='SingleAxon'
model_file = f'''<neuroml>
    <include href="Subdivided_Cable.cell.nml"/>

    <pulseGenerator id="pulseGen1" delay="5ms" duration="200ms" amplitude="0.1nA"/>
    <network id="Net">

        <population id="{cable_pop_name}" component="Subdivided_Cable" size="1" />
        <inputList id="stimInput_Subdivided_1" component="pulseGen1" population="{cable_pop_name}">
            <input id="0" target="../{cable_pop_name}[0]" segmentId="0" fractionAlong="0.3" destination="synapses"/>
        </inputList>

    </network>
</neuroml>'''
# print(model_file)
with open('Model_LonelyAxon.nml', 'wt') as f: f.write(model_file)

Recording over the cell’s full extent

We want to record the state of the cable throughout, over time. We’ll use EDEN’s explain_cell and GetLemsLocatorsForCell helper routines to tell us which spots on the neuron to record, so that all compartments are covered.

Note

Because official NeuroML allows recording on <segment>s but only on the middle of each <segment> (when there’s only one in our stylised morphology), we’ll have to use extended LEMS path syntax, which is supported by EDEN but not official NeuroML.



import eden_simulator
cells_info = eden_simulator.experimental.explain_cell('Model_LonelyAxon.nml')
cell_info = cells_info['Subdivided_Cable']



comp_locators = eden_simulator.experimental.GetLemsLocatorsForCell(cell_info)
cell_voltage_paths = [ f'{cable_pop_name}[0]/{loc}/v' for loc in comp_locators ]
tabline = '\n        '
sim_file = f'''<Lems>
<Include href="Model_LonelyAxon.nml"/>
<Simulation id="MySim" length="100 ms" step="50 us" target="Net">
    <OutputFile id="MyOutFile" fileName="results.gen.txt">
        {tabline.join([ f'<OutputColumn id="v_0_{loc}"  quantity= "{path}"/>' for loc, path in zip(comp_locators, cell_voltage_paths)])}
    </OutputFile>
</Simulation>
<Target component="MySim"/></Lems>'''
# print(sim_file)
with open('Sim_LonelyAxon.xml', 'wt') as f: f.write(sim_file)

Running the simulation and displaying results



results = eden_simulator.runEden('Sim_LonelyAxon.xml')



import numpy as np
rec_time_axis_sec = results['t']
neuron_waveforms = np.array([results[path] for path in cell_voltage_paths])
neuron_waveforms.shape



(50, 2001)

Since it’s a straight unbranched cable, we can intuitively assume a single spatial axis x along it, and plot a 2-D function over x and time.

When a neuron is more like a sprawling tree (or bush), we can still record along a certain linear section, but otherwise we may have to use a more skeuomorphic display.



from matplotlib import pyplot as plt
# help(plt.imshow)
im = plt.imshow(neuron_waveforms,
    extent=[ *rec_time_axis_sec[[0,-1]], cable_length_microns, 0 ],
           aspect='auto',interpolation='None')
cbar = plt.colorbar(im); cbar.set_label('Voltage (V)')
plt.xlabel('Time (sec)');plt.ylabel('Position (μm)')
plt.xlim((0, 0.02));
_images/intro_spatial_64_0.png
Finally, let’s visualise the “real” thing in space, painted with (of course) false colour that represents membrane voltage.
(And as a first step to that, reduce the amount of displayed data to a practical frame rate.)


from eden_simulator.display.animation import subsample_trajectories
_, anim_axis, sampled_time_axis_sec, (sampled_voltage,) = subsample_trajectories(rec_time_axis_sec, [neuron_waveforms.T])
print("Animation duration: %d frames, %.3f sec" % (len(anim_axis), anim_axis[-1]))
anim_axis

Animation duration: 1001 frames, 33.333 sec



array([0.00000000e+00, 3.33333333e-02, 6.66666667e-02, ...,
       3.32666667e+01, 3.33000000e+01, 3.33333333e+01], shape=(1001,))

Let’s see the painted neuron that you might have already visualised in your mind, while modelling it. We’ll use K3D for an animated 3-D display.



import k3d; from eden_simulator.display.spatial.k3d import Plot, plot_neuron
plot = Cable_plot = Plot()
plot += plot_neuron(cell_info, sampled_voltage, time_axis_sec=anim_axis, color_map='jet',);
# Add a nice text box to inform us about the time elapsed in the simulation.
k3d_label_dict = { str(real_time): f't = {sim_time*1000:.1f} ~ ms' for (real_time, sim_time) in zip(anim_axis, sampled_time_axis_sec)  }
plt_label = k3d.text2d(k3d_label_dict, (0.,0.)); plot += plt_label
# And display as a HTML inset.
from IPython.display import HTML, display
plot.snapshot_type = 'inline'; display(HTML(plot.get_snapshot()))

Exercise: When the cable is 200 μm long, we observe that the action potential propagates slower near its starting point, than away of it. (Enlarge the time axis or focus on a small part of it to see this clearly.)

  • Grow the cable to 1~2 mm (increase the numberInternalDivisions if necessary for an accurate result) and look again at membrane voltage over time and space. How does the picture change? Plot the subset, show it beside the original shorter cable so that the clamped points are aligned.

  • Why is AP propagation so sped up near the ends of the cable? Consider the system’s dynamics.

  • If an endless cable is to be modelled (to study small-scale propagation dynamics and such), how would you alter the model, to minimise the effect of the cable’s finite length? Look up how this is handled with other types of cables.

Exercise: Explore how the discretisation’s fineness affects the result. What happens when there are too few (or too many) compartments? What happens as the spatial resolution starts getting too low or too high? Compare and contrast with inappropriate values for the simulation’s time step.

  • Tip: Wrap the model-to-2D-plot procedure in a routine instead of running all the cells separately.

More resources on spatially-detailed modelling

Multi-compartment modelling in other simulators

The concept of spatially-detailed neurons that are modelled as coupled compartments is also supported by other neuro-simulators, and their platform-specific modelling formats. Readers who have learned (or are learning) to use another simulator, or want to compare and contrast the implementations of concepts and see more usage examples, may refer to:

  • Chapter 5 of the NEURON book (and the whole book and user’s guide for more details);

  • Arbor’s implementation of the concept, much in line with NEURON and NeuroML, as “cable” cells;

  • Chapters 2 and 5 of the Book of Genesis;

  • BRIAN’s multi-compartment capability — note that the compartment’s dynamics equations are all written explicitly (i.e. same as with point neurons), rather than being composed of the various distributed mechanisms.

Next steps

After learning how to model spatially-detailed neurons and their geometry-based biophysics, it’s time to model more neurons, with more mechanisms on them and more effects that control phenomena. Since each model and experiment typically introduces highly-specific effects, mechanisms that capture them can be written by the model authors, as LEMS extensions.

Following this tutorial, you may also try out the “Building a network” tutorial, which shows how to use and display a full-featured spatially-detailed cell. As mentioned previously, if a quantity needs to vary among neurons or across their extent for a desired model, EDEN permits full control over variability through <𝙴𝚍𝚎𝚗𝙲𝚞𝚜𝚝𝚘𝚖𝚂𝚎𝚝𝚞𝚙>.