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 slong 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 equantions 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 comapred 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 desctibed 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)
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 <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; theid
number is, instead.Within the segment there are more tags:
Two tags called
proximal
anddistal
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 frusum on that end, in microns.
For non-root segments, a
<parent>
tag specifying a theid
of each’s parentsegment
. 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, from0
(proximal
“start” point) to1
(distal
“end” point).The default value is
1
(attached to the parent’sdistal
end).
If the segment has a
parent
, defining its<proximal>
end is optional; by default, it will be same as the point interpolated byfractionAlong
, on theparent
segment. thex
,y
,z
as well as thediameter
are interpolated between the two ends of the parent.
Standards in morphology coordinates¶
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.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 (eg 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 asegment
can be different from theparent
’sdistal
; 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 fractionlong
the axial direction of said <segment>
. This pair of segment
and fractionAlong
is used to:
precisely attach synapses and stimulation probes to cells;
record local action potentials with
<EventOutputFile>
s, and trajectories of quantities in<OutputFile>
s — to specifyfractionAlong
together with thesegment
, use EDEN’s extended LEMS paths.
In all these cases:
If the
segmentId
is not specified, it is assumed0
: typically the soma (but look out), or the entirety of a point neuron.If the
fractionAlong
is not specified, it is assumed0.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 theid
s of<segment>
s, this can be any alphanumeral.neuroLexId
: A machine-readable NeuroLex keyword that identifie 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 infromation about the segment group. Thetag="numberInternalDivisions"
may be used to help with discretisation, whenneuroLexId=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 withneuroLexId="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.
numberInternalDivisions
.Retrieving and using the discretisation¶
After specifying the discretisation, getting it back can also prove useful. Typical applications are:
displaying a changing quantity across the neuron;
visually inspecting the discretisation’s fineness while modelling;
controlling non-uniform biophysical properties at the finest level;
better reconciling the equations between domains of the simulation, when the surroundings of the neuron matter; see also the LFP example.
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 aditional 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:
Electrical properties of the neuritic cable, such as surface capacitance, axial resistance, and voltage threshold to register an action potential (to trigger post-synapses);
Distributions of trans-membrane mechanisms, such as ion channels, passive leaks, and concentration models (and pumps) of substances (such as ions) around the membrane;
Individual (usually trans-membrane) instances of mechanisms, such as synapses and stimulation probes.
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
orsegmentGroup
attribute; if neither is specified, it applies for the groupall
.
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:
Classic Hodgkin-Huxley gate models:
<gateHHrates>
,<gateHHtauInf>
,<gateHHratesTau>
,<gateHHratesInf>
,<gateHHratesTauInf>
HH sub-gate models:
<gateFractional>
Multi-state gates with general kinetics:
<gateKS>
Stateless gates:
<gateInstantaneous>
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 poetential 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:
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 direcly convertible to and from \(\alpha, \beta\) functions through the relations:
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
type
s 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
Transition
s of
Markov-style gates rates as applicable.𝙷𝙷𝙴𝚡𝚙𝚅𝚊𝚛𝚒𝚊𝚋𝚕𝚎 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 ismidpoint
.midpoint
: The reference point to center the curve around. The function has the valuerate
for this value of membrane voltage.scale
: The independent-variable scaling unit to evolve the curve by. Valuesscale
apart on the independent variable are in a proportion of \(1\colon e\).Notice that
midpoint
andrate
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 ismidpoint
.midpoint
: The reference point to center the curve around. The function has the valuerate/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 ismidpoint
.midpoint
: The reference point to center the curve around. The function has the valuerate/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¶
<gate>
in an <ionChannel>
must specify the following attributes:id
entifying 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 dynamicstype
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 altenatively \(\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 atype
attribute referring to a per time-valued function;\(\tau\) as
<timeCourse>
, with atype
attribute referring to a time-valued function;\(x_\infty\) as
<steadyState>
, with atype
attribute referring to a dimensionless function.
<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¶
The NeuroML tag to describe this is <𝚐𝚊𝚝𝚎𝙺𝚂>. It contains multiple <openState>
s <closedState>
s, each with an id
entifier name attibute. It also contains multiple transition paths between these state
s, 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 betweenfrom
andto
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 constant0
.
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:
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¶
<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\)).q10Factor
(unitless, per 10°C) andexperimentalTemp
that the gate’s recorded dynamics refer to (i.e. the temperature at which the gate behaves the same as without theq10Settings
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 insterted to <ionChannel>
s and has attributess 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 channnel 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 ther place (or simply scale byconductanceDensity * compartment area
for the usual macroscopic models).ion
(optionally) which designates the specific ion species to be pooled. May benon_specific
when referring to an “<ionChannel>"
lacking theion
property. Useful then there are e.g. two separate “pools” of calcium which activate certain mechanisms differently for some reason.segment
orsegmentGroup
(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 segmentGroup
s, 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 becondDensity
for<channelDensityNonUniform>
and<channelDensityNonUniformNernst>
, andpermeability
for<channelDensityNonUniformGHK>
.
And inside each <variableParameter>
tag, there is a <inhomogeneousValue>
tag with attributes:
inhomogeneousParameter
: Theid
of the<inhomogeneousParameter>
specified in the<segmentGroup>
value
: A dimensionless LEMS expression that involves the<inhomogeneousParameter>
by itsvariable
name. Note: The numerical value of the expression stands for SI units, such as mhos per square meter, \({S}/{m^2}\) forcondDensity
and metres per second, \(m/s\) forpermeability
.
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 isPath Length from root
. Note: Unlike most LEMS quantities, this variable is provided dimensionless and its value is in microns (if not made relative withnormalizationEnd
).
The following optional tags inside the <inhomogeneousParameter>
apply transformations to the metric
variable:
<proximal translationStart="<value>"/>
: Offset the metric so that the minimum isvalue
(typically0
).<distal normalizationEnd="<value>"/>
: Divide the metric my its maximum and then multiply byvalue
(typically1
), so that the variable is now a relative, dimensionless measure.
Adding concentration models near the membrane¶
<species>
tags inside <intracellularProperties>
(i.e. along with axial <resistivity>
tags).<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
orsegmentGroup
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:
<𝚍𝚎𝚌𝚊𝚢𝚒𝚗𝚐𝙿𝚘𝚘𝚕𝙲𝚘𝚗𝚌𝚎𝚗𝚝𝚛𝚊𝚝𝚒𝚘𝚗𝙼𝚘𝚍𝚎𝚕>, for modelling a thin shell of a certain
shellThickness
near the compartment’s membrane, and adecayConstant
of time over which the concentration converges torestingConc
;<𝚏𝚒𝚡𝚎𝚍𝙵𝚊𝚌𝚝𝚘𝚛𝙲𝚘𝚗𝚌𝚎𝚗𝚝𝚛𝚊𝚝𝚒𝚘𝚗𝙼𝚘𝚍𝚎𝚕>, for similar dynamics, but with a fixed linear relationship between the incoming current density and the resulting rate of change of the local concentration (instead of considering the Faraday constant and
shellVolume
, as the previous model does).
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 thepre
- andpost
-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));
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])
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()))