Custom per-element variability¶
As we saw in the NeuroML introduction tutorial, models comprise a <network>
of neuron populations connected by synapse projections, and some experimental rig around them to interact with through.
In practice, a problem that emerges often when modelling is that NeuroML defines populations of identical neurons, and projections with identical components. The state of each instance may of course vary through the simulation, but the parameters and initial state are the same for all instances of a cell (or synapse or input type, or generally LEMS component instance).
However, this is not always sufficient to make a realistic model, as neurons can show considerable physiological variability in nature. And even when neurons and synapses are of the same type, a small amount of variability is desired in models to make them more robust.
For these reasons, modellers often require that each neuron or synapse may have its parameters slightly different than the rest; and what’s more, the kind of variability, or anatomy-informed) is often unique to that model part, or at least project or research laboratory. Typical examples of variability models include mathematical formulae, anatomical data or random sampling of statistical distributions. What’s then needed is a way to apply any kind of variability to the model directly. And in the general case, this includes assigning an explicit set of values, one by one, to a collection of mechanisms.
EDEN supports models with arbitrary per-instance variability, in a way that allows for any modelling pipeline by design. It does this through the <EdenCustomSetup>
extension, and a special setup language for supplying and designating the variability data.
Tip
This is a simulator-specific interim data format. It may be replaced with equivalent NeuroML capabilities as they appear.
The first step in using custom variability is to add the tag <EdenCustomSetup filename="<file name>" />
inside the usual <Simulation>
tag.
The file under file name
then, contains the customisation information as follows.
The EdenCustomSetup
file specification¶
The setup language is (to this point) essentially a sequence of set
statements, which set the values of quantities (i.e. <Property>
s ) over a set of mechanisms.
Note that these statements are considered in strict order: a set
statement may supersede values assigned by a previous statement to the same mechanism!
In order to specify a collection of mechanisms instead of just a single one, the LEMS path machinery is used in a kind of re-mixed form. It is broken down into:
the part that locates an entity in the
<network>
, which in turn is broken down into:the “collection” of entities in the network, that is a neuron
<population>
, synaptic<projection>
, or<inputList>
,and the “identifiers” of entities within said collection. A
set
statement may name multiple identifiers, to share the same quantity value or get each one value from a list ofmulti
ple, as we’ll see below.For spatially-detailed
<cell>
populations, there is also a list of different cell locations that comes after the list of cell identifiers. The mechanisms in these locations may in turn take the same value, or different ones.
and the path within the entity that locates the quantity in terms of the sub-mechanism that it belongs to, and its identifier name within.
For the extended set of LEMS paths that EDEN understands (and a hint as to how they may be appled here), refer to EDEN’s LEMS paths reference.
General structure¶
<EdenCustomSetup>
files are written as human-readable text lines; an equivalent binary format for higher data density is also planned.
The contents of the file are arranged in lines; each line contains tokens separated with space and tab characters, and other “white-space” characters.
The types of lines are as follows:
set
statements, for:cell
populationssynapse
projectionsinput
lists
values
lines, which follow theset
statements that specifymulti
ple values;empty and “comment” lines to aid readability for humans.
The structure of these lines is described below.
set cell
statement¶
A set cell
statement line is a sequence of tokens as follows:
set cell <population> <cell list> <location list> <attribute> <value>
The tokens have the following meanings, in order:
population
: The name of the neuron<population>
whose cells’ quantities are customised.cell list
: Either a comma-separated list of cell instanceid
’s in strictly ascending order, orall
to designate the whole population.location list
: Can be either:a comma-separated list of segment locators on cell;
all
to apply all over the cell;or the name of a “cable” segment group, to sample values along a track of sampling points.
For point neurons, options are the equivalent
all
or0
(the implicit singular “segment”id
, as per LEMS paths).attribute
: The part of the LEMS path within the neuron, that follows afterpopulation[id]/segment/
value
: The single quantity value ,ormulti
to specify more, in thevalues
lines that follow. Examples for howvalue
s look like:100 pA
to specify an amount of electrical current;multi mV
to specify one or more voltage values, in one or morevalues
lines that follow;3.14159
to specify a dimensionless, pure number;multi
to specify one or more dimensionless numbers, or <𝚅𝚊𝚛𝚒𝚊𝚋𝚕𝚎𝚁𝚎𝚚𝚞𝚒𝚛𝚎𝚖𝚎𝚗𝚝> paths.
If the location list
is “all” and multi
ple values are specified, one values
line follows, with one value per cell instance in the cell list
.
If the cell list
includes only one cell, the location list
contains multiple segment locators and multi
ple values are specified, one values
line follows, with one value per segment locator in the location list
.
When both the cell list
contains multiple cells and the location list
contains multiple segment locators, there are different ways to specify multiplicity. The multi
value specifier has to be replaced by:
multi_location
, in which case the same values are applied to all cells, and onevalues
line follows with the value for each segment locator;multi_cell
, in which case the same value is applied to all segment locators per cell, and onevalues
line follows with the value for each cell;multi_multi
, in which case a different value is applied to each cell on each segment locator. As manyvalues
lines as the number of mentioned cells follow, each line with the value for each segment locator on the listed cell.
Cable mode¶
When the location list is a “cable” <segmentGroup>
, a special cable interpolation mode is engaged, wherein the mechanism value is sampled across the cable’s extent for each compartment, based on a piecewise-linear function that’s provided in the values
lines that follow.
The value
token must be cable
or cable_multi
. At least two values
lines follow:
The first
values
line is special; it contains the sampling points for the values that follow. These sampling points represent the relative fractional length along the cable, starting from the proximal end and up to the distal end, and thus may have values from \(0\) to \(1\).If
value
isvalue_multi
, onevalues
line for each cell specified in thecell list
follows, each line with the piecewise-linear function’s values at the aplready specified sampling points.
Tip
This mode is targeted toward porting models written for NEURON, when values for DENSITY
mechanism quantities are assigned over a neuron’s section
with any sort of variability. The values instantiated on the NEURON model can be extracted with:
access <section>; for (x) print x, <property>(x)
as seen in chapter 5 of the NEURON book, section “Working with range variables”.
set synapse
statement¶
A set synapse
statement line is a sequence of tokens as follows:
set cell <projection> <synapse list> <site> <attribute> <value>
The tokens have the following meanings, in order:
projection
: The name of the synaptic<projection>
whose components’ quantities are customised.synapse list
: Either a comma-separated list of<connection>
id
’s in strictly ascending order, orall
to designate the whole projection.site
: either orpre
orpost
. Since classical event-based<projection>
s have a synaptic component only on the post-synaptic site,post
is the only allowed site for them.attribute
: The part of the LEMS path within the synaptic component, that follows afterprojection[id]/post/
.
Examples for value
s are the same as with cells.
When multi
ple values are specified, one values
line follows, with one value per mechanism instance in the specified list.
set input
statement¶
A set input
statement line is a sequence of tokens as follows:
set cell <list name> <input list> <attribute> <quantity with units>
The tokens have the following meanings, in order:
projection
: The name of the<inputList>
whose probes’ quantities are customised.input list
: Either a comma-separated list of<input>
id
’s in strictly ascending order, orall
to designate the whole<inputList>
.attribute
: The part of the LEMS path within the input source component, that follows afterinputlist[id]/
.
Examples for value
s are the same as with cells.
When multi
ple values are specified, one values
line follows, with one value per mechanism instance in the specified list.
Comment lines and inline comments¶
In the text-based format, there may be some space between, and on, the statement lines to write comments on, to aid people (or perhaps an automated process):
Empty lines may be present anywhere in the text; they are ignored, but may not split other lines.
Within each line, the part starting from the first token that starts with a
#
character, is also ignored.
Setting VariableRequirement
s¶
Along with numerical values for regular quantities, set
statements can also be used to assign the targets of <𝚅𝚊𝚛𝚒𝚊𝚋𝚕𝚎𝚁𝚎𝚚𝚞𝚒𝚛𝚎𝚖𝚎𝚗𝚝>s; in fact, this is the only way for them to take a value (which they must!). In this case, the value to specify is a LEMS path instead of a number. There is no need to specify units, as the model’s quantities are already automatically managed by EDEN.
Note: When a VariableRequirement
set of values is interpolated over an unbranched section using cable mode with sample points, nearest-neighbour interpolation is used instead of the linear interpolation that’s used for quantities. (If there is one sampling point per compartment, the effect is the same.)
Example: A network with variability over neurons, synapses and probes¶
To show how the above specification works, let’s make a 2-D connected grid of neurons and randomise parameters throughout:
import numpy as np; from numpy.linalg import norm; from matplotlib import pyplot as plt
rng = np.random.default_rng(seed=43264326)
offsets = [(1,0), (-1,1), (0,1), (1,1)] # sequential-first neighborhood
# for bidirectional: offsets = [(x,y) for y in [-1, 0, +1] for x in [-1, +1]]+[(0,y) for y in [-1, +1]]
# Make a 2D grid and then linearise it
grid_step_size_microns = 50
X_steps = 15; x_space = range(X_steps); X_hw = (X_steps-1)*grid_step_size_microns/2
Y_steps = 12; y_space = range(Y_steps); Y_hw = (Y_steps-1)*grid_step_size_microns/2
xmat, ymat = np.meshgrid(x_space, y_space, indexing='ij')
neuron_points_steps = np.stack((xmat,ymat),axis=-1).reshape((-1,2)); neurons = len(neuron_points_steps)
def xy_to_index(x,y): return x*Y_steps + y
# and in spatial, not abstract coordinates
neuron_points_microns = (np.array(neuron_points_steps) + .5 - (X_steps/2, Y_steps/2)) * grid_step_size_microns
pts = neuron_points_microns # shorthand
# Now gather the possible synapses under Neumann adjacency (ie 8-way)
syn_list = [(xy_to_index(x,y), xy_to_index(xx,yy))
for (x,y) in neuron_points_steps
for (dx, dy) in offsets for (xx, yy) in [(x+dx, y+dy)]
if 0 <= xx < X_steps and 0 <= yy < Y_steps
]
# Randomly flip the directions of some connections (could have both directions but it would be hard to show that)
syn_list = np.array(syn_list)
which_to_flip = rng.uniform(size=len(syn_list)) < 0.5
syn_list[which_to_flip] = [(y, x) for (x,y) in syn_list[which_to_flip]]
# And set up some parameters
input_stim_idxs = np.where( 0
| ( norm((pts - (+3*X_hw/14,+Y_hw/2.75))/(1), axis=-1) < Y_hw*.5/2.75)
| ( norm((pts - (-3*X_hw/14,+Y_hw/2.75))/(1), axis=-1) < Y_hw*.5/2.75)
| (
( norm((pts - (0,-Y_hw/2.75))/(1,.55), axis=-1) < X_hw/2)
& (pts[:,1] < -50)
)
)[0]
input_stim_delay_msec = 20; input_stim_duration_msec = 200+50*rng.uniform(size=len(input_stim_idxs));
neuron_Vrest_mV = ( -70
+ 30* (norm(pts - (0,0), axis=-1) < Y_hw*.95)
+ 10*rng.normal(size=neurons))
neuron_Vrest_mV[input_stim_idxs]=-70 # Don't differentiate the
syn_tau_msec = 0.1+10*rng.uniform(size=len(syn_list)); max_syn_tau_msec = np.max(syn_tau_msec)
# And display
fig, ax = plt.subplots(figsize=(16,10)); ax.set_aspect('equal')
for idx,(i,j) in enumerate(syn_list):
plt.arrow(*pts[i], *pts[j]-pts[i], color=(0.9*(syn_tau_msec[idx]/max_syn_tau_msec),0,0),
head_width=8, width = 1,lw = 0, length_includes_head=True)
plt.scatter(*neuron_points_microns[input_stim_idxs].T, color='xkcd:sky', s=225, lw=0, )
plt.scatter(*neuron_points_microns.T, c=neuron_Vrest_mV, s=100, lw=0)
ax.set_facecolor('#f0f4ff')
There’s a network with variability all over the place.
Now let’s generate the model file using f-strings and list comprehensions:
tabline = '\n ' # annoyingly enough, \ may not exist at all in a f string, not even in brackets. Fixed in Python 3.12+
model_file = f'''
<neuroml>
<izhikevich2007Cell id="IzhCell" v0 = "-60mV" C="100 pF" k = "0.7 nS_per_mV"
vr = "-60 mV" vt = "-30 mV" vpeak = "35 mV"
a = "0.03 per_ms" b = "-2 nS" c = "-50 mV" d = "100 pA"/>
<pulseGenerator id="DcProbe" delay="{input_stim_delay_msec} ms" duration="0 ms" amplitude="0.5nA"/>
<expOneSynapse id="ExpCondSynapse" gbase="8.1nS" erev="20mV" tauDecay="5ms"/>
<network id="Net" type="networkWithTemperature" temperature="37degC" >
<population id="Pop" component="IzhCell">
{tabline.join([ f'<instance id="{i}"><location x="{x}" y="{y}" z="{0}"/></instance>' for i, (x,y) in enumerate(pts)])}
</population>
<inputList id="Inp" population="Pop" component="DcProbe">
{tabline.join([ f'<inputW id="{i}" target="Pop[{idx}]" destination="synapses" weight="1"/>' for i, idx in enumerate(input_stim_idxs)])}
</inputList>
<projection id="GridProjection" presynapticPopulation="Pop" postsynapticPopulation="Pop" synapse="ExpCondSynapse">
{tabline.join([ f'<connectionWD id="{ii}" preCellId="{pre}" postCellId="{post}" weight="1" delay="5 ms"/>' for ii, (pre,post) in enumerate(syn_list)])}
</projection>
</network>
<Simulation id="MySim" length="1 s" step="100 us" target="Net"> <!-- no need for a seed yet, add one if you want ❗ -->
<OutputFile id="MyOutFile" fileName="results.gen.txt">
{tabline.join([ f'<OutputColumn id="v_{i}" quantity= "Pop[{i}]/v"/>' for i, _ in enumerate(pts)])}
</OutputFile>
<!-- Use EdenCustomSetup ❗ -->
<EdenCustomSetup filename="Example_CustomSetup.txt"/>
</Simulation>
<Target component="MySim"/></neuroml>
'''
# print(model_file)
with open('Example_CustomSetup.nml', 'wt') as f: f.write(model_file)
And the Example_CustomSetup.txt
that goes with it:
setup_file = ''
setup_file += 'set cell Pop all all vr multi mV\n'
setup_file += 'values '+ ' '.join(f'{x}' for x in neuron_Vrest_mV)+'\n'
setup_file += 'set synapse GridProjection all post tauDecay multi msec\n'
setup_file += 'values '+ ' '.join(f'{x}' for x in syn_tau_msec)+'\n'
setup_file += 'set input Inp all duration multi msec\n'
setup_file += 'values '+ ' '.join(f'{x}' for x in input_stim_duration_msec)+'\n'
# print(setup_file)
with open('Example_CustomSetup.txt', 'wt') as f: f.write(setup_file)
And run the simulation:
from eden_simulator import runEden
results = runEden('Example_CustomSetup.nml')
neuron_waveforms = np.array([ results[f"Pop[{i}]/v"] for i, _ in enumerate(pts) ])
# plt.imshow(neuron_waveforms, interpolation='none', aspect='auto');plt.colorbar()#plt.clim((-0.07, -0.04))
# plt.plot(neuron_waveforms[input_stim_idxs].T);
Let’s see the neuron activity in pretty 4-D (surprisingly more compact than matplotlib animations):
from eden_simulator.display.animation import subsample_trajectories
samples, anim_axis, sampled_time_axis, [sampled_soma_voltage] = subsample_trajectories(
results['t'], [neuron_waveforms.T * 1000], animation_speed=0.0096, animation_frames_per_second=48)
!pip install -q k3d
import k3d; from eden_simulator.display.spatial.k3d import Plot
points_plot = plot = Plot(background_color=0xf0f4ff, camera_auto_fit=False) # to override camera orientation
k3d_anim_dict = { str(real_time): x.astype(np.float32)
for (real_time, x) in zip(anim_axis, sampled_soma_voltage) }
k3d_label_dict = { str(real): f't = {sim:.3f} ~ s' for (real, sim) in zip(anim_axis, sampled_time_axis) }
plt_points = k3d.points(positions=np.c_[pts,[0]*len(pts)].astype(np.float32),attribute=k3d_anim_dict,
color_range=[-80, -20],point_size=30,color_map=k3d.matplotlib_color_maps.Hot); plot += plt_points
plot.camera = plot.get_auto_camera(pitch=30, yaw=10)[:6]+[0,1,0] # set 'y' to up ! vecs are pos, tgt, up
plt_label = k3d.text2d(k3d_label_dict, (0.,0.), label_box=False); plot += plt_label # add 2d elements AFTER setting auto camera
plot.fps = 48;
from IPython.display import display, HTML, IFrame
plot.snapshot_type = 'inline'; display(HTML(plot.get_snapshot()))