This page was generated from docs/intro_neuroml.ipynb. Interactive online version: . Download notebook.

NeuroML basics

The program that this user guide is about, EDEN, simulates the activity of spiking neural networks. The models of such networks and their neurons are described in NeuroML, which is a neural modelling language in the XML format. EDEN then works by reading NeuroML files and crunching the numbers the model described in these files.

This article will explain how models are written in NeuroML, and how to use EDEN to simulate them. The various parts of a NeuroML simulation are introduced one after the other, in the following sections.

Creating a simulation

First of all, to set up a NeuroML simulation and “run” it, it takes:

• a NeuroML network that contains all the things to simulate;

• a simulation that contains the parameters of how to simulate said things (i.e. for how long, what to observe and such).

These are created by writing down their descriptions as text-based, XML files, for various programs to read back.

Let’s now see how to write the smallest possible network and simulation.

In the following, the %%writefile line is a Jupyter magic so that when this code cell is run, the specified file will be (over)written with the contents as follows in the cell.

Tip

We’ll be changing the model files through the tutorial to show different features. When running the notebook, if you want to re-run a previous section, make sure to overwrite the %%writefile cells to that point before re-running the simulation.

%%writefile EmptyModel.nml
<!-- These properties on the top level <neuroml> tag are just to positively identify the file as NeuroML.  They are not mandatory -->
<neuroml xmlns="http://www.neuroml.org/schema/neuroml2"  xmlns:xs="http://www.w3.org/2001/XMLSchema" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:schemaLocation="http://www.neuroml.org/schema/neuroml2 https://raw.github.com/NeuroML/NeuroML2/development/Schemas/NeuroML2/NeuroML_v2.1.xsd">
<!-- Here comes the network -->
<!-- set the temperature to that you think your model is running under (depends on the species, whether it's cold blooded) -->
<network id="MyNetwork" type="networkWithTemperature" temperature="37degC" >
<!-- Nothing here yet -->
</network>
</neuroml>

Writing EmptyModel.nml

And the simulation file (see also the NeuroML guide) pointing to the EmptyModel, by convention there is one file for the <Simulation> and one for the <network>:

%%writefile LEMS_EmptySim.xml
<?xml version="1.0" encoding ="UTF-8" ?>
<Lems>
<!-- we will be <include>ing the nml file with MyNetwork -->
<include file="EmptyModel.nml" />
<!-- the Simulation is the main content, it needs a duration and a timestep.  It involves a certain "target" <network> -->
<Simulation id="MySim" length="10 s" step="10 us" target="MyNetwork" >
<!-- Nothing here yet -->
</Simulation>
<!-- Just in case there are multiple <Simulation>s described, specify the one to run here -->
<Target component="MySim"/>
</Lems>

Writing LEMS_EmptySim.xml

Now let’s run this completely minimal simulation.

EDEN can be used from Python with the eden_simulator package, just pip install it if you haven’t already in Quickstart:

from eden_simulator import runEden
%time runEden('LEMS_EmptySim.xml');

CPU times: user 54.4 ms, sys: 4.23 ms, total: 58.7 ms
Wall time: 652 ms

Now let’s see how we can make it a tad more interesting.

Adding a neuron and recording output

The first step to define a neural model is to define a neuron. Let’s add a neuron to the model and see how it behaves.

This is done by defining the neuron model as a cell type, and adding to the <network> a <population> with cells of said cell type (called the component for this <population>).

Instead of specifying the cell type from scratch, we’ll use the common linear integrate-and fire model as a base, and supply our own dynamical parameters for it. This cell type can be found in the standard NeuroML library as <iafCell> and has parameters C(apacitance of cell membrane), leakConductance (of membrane), leakReversal (potential), reset (potential right after firing) and thresh(old for firing).

We’ll set them as follows: C = 200 pF, leakConductance = 10 nmho, leakReversal = -70 mV, reset = -70 mV, thresh = -50 mV.

As we’ll see in an in-depth tutorial, the population can also be a populationList with specific 3-D positions for each cell, but this is not necessary for simulating the population.

%%writefile SingleCell_Model.nml
<!-- These properties on the top level <neuroml> tag are just to positively identify the file as NeuroML.  They are not mandatory -->
<neuroml xmlns="http://www.neuroml.org/schema/neuroml2"  xmlns:xs="http://www.w3.org/2001/XMLSchema" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:schemaLocation="http://www.neuroml.org/schema/neuroml2 https://raw.github.com/NeuroML/NeuroML2/development/Schemas/NeuroML2/NeuroML_v2.1.xsd">
<!-- Here comes the cell type ❗ -->
<iafCell id="MyFirstCellType" C="200 pF" leakConductance="10 nS" leakReversal="-70 mV" reset="-70mV" thresh="-50mV" />
<!-- Here comes the network -->
<network id="MyNetwork" type="networkWithTemperature" temperature="37degC" >
    <!-- Add a population with a single cell ❗ -->
    <population id="MyFirstPopulation" component="MyFirstCellType" size="1"/>
</network>
</neuroml>

Writing SingleCell_Model.nml

There are more point neuron types defined in base NeuroML (see for example, adExIaFCell and izhikevich2007Cell), and custom ones can be built as LEMS components in a following section.

More than point neurons, NeuroML and EDEN also support spatially modelled neurons with a detailed morphology and spatial compartments. Due to the added details and modelling options, these are covered in a following section.

In order to observe the cell, we’ll have to record one of its attributes over time.

This is done in the simulation file by adding an <OutputFile> with <OutputColumn>s to the <Simulation>.

• Each OutputColumn points to a quantity, which in the case of simplified neurons is described as population name[cell number in population starting from 0]/variable. (Refer to the path documentation for more possibilities.)

We see (click on “Dynamics” here) that v stands for the membrane voltage, let’s record that.

Hence we want to record the quantity MyFirstPopulation[0]/v.

Let’s write the new simulation file:

%%writefile LEMS_SingleCell_Sim.xml
<?xml version="1.0" encoding ="UTF-8" ?>
<Lems>
<!-- we will be <include>ing the nml file with the SingleCellModel -->
<include file="SingleCell_Model.nml" />
<!-- the Simulation is the main content, it needs a duration and a timestep.  It involves a certain "target" <network> -->
<Simulation id="MySim" length="1 s" step="10 us" target="MyNetwork" >
    <!-- Add an OutputFile ❗ -->
    <OutputFile id="MyFirstOutputFile" fileName="results.gen.txt">
        <OutputColumn id="my_first_vm" quantity="MyFirstPopulation[0]/v"/>
    </OutputFile>
</Simulation>
<!-- Just in case there are multiple <Simulation>s described, specify the one to run here -->
<Target component="MySim"/>
</Lems>

Writing LEMS_SingleCell_Sim.xml

And run the new simulation:

results = runEden('LEMS_SingleCell_Sim.xml')
results

{'MyFirstPopulation[0]/v': array([-0.07, -0.07, -0.07, ..., -0.07, -0.07, -0.07], shape=(100001,)),
 't': array([0.0000e+00, 1.0000e-05, 2.0000e-05, ..., 9.9998e-01, 9.9999e-01,
        1.0000e+00], shape=(100001,))}

The output of Eden is in the form of a dict. Its entries are t and the specified quantity of each OutputColumn, same as when using the general tool PyNeuroML.

The value of results['t'] stands for the list of recorded points in time for the time series, and every other entry value stands for the trajectory of the respective variable over the sampling points in time.

All values are are in SI units, that is seconds for t volts for v.

Let’s plot MyFirstPopulation[0]/v against t then:

from matplotlib import pyplot as plt; plt.figure(figsize=(4,2))
plt.plot(results['t'],results['MyFirstPopulation[0]/v'])
plt.xlabel('Time (sec)'); plt.ylabel('Membrane voltage (V)');

Apparently this type of neuron doesn’t fire intrinsically; it just stays polarised and quiet, producing this electrical flat line.

Other types of neurons fire or oscillate intrinsically, but we’ll need to apply some physiological stimulus on this one to see its response.

Adding stimuli

We typically study neurons in isolation, or at least as part of a small piece of tissue. But neurons normally react to electro-chemical stimuli supplied by other neurons, thus we’ll have to provide an artificial “equivalent” to study their behaviour. (We can also study the natural response of neurons to their natural stimuli in vivo, but this is much harder to understand and reason about.)

A basic way to play with cells is the current clamp, and in this virtual world we don’t even have to worry about injuring the cells. By injecting “positive” current, we can depolarise a cell and cause action potentials.

This is done by defining the input’s model as an input source type, and adding to the <network> an <inputList> with probes of said source type (called the component for this <inputList>).

Instead of specifying the input type from scratch (as we’ll be doing in a [following section]), we’ll use the common DC pulse current clamp (flat current injected over a time window) as a base, and supply our own dynamical parameters for it. This input type can be found in the standard NeuroML library as <pulseGenerator> and has parameters delay (before the pulse starts), duration (of the pulse), and amplitude.

We’ll set them as follows: delay = 100 msec, duration = 500 msec, amplitude = 210 pA.

Next, we’ll add an <inputList> to the <network>. Each <inputList> contains one or more copies (called <input> instances) of the same input source component, and is applied on a certain cell population.

Each <input> has a cell target which can be, for example, written as <population name>[<cell instance number>]. (See more about LEMS paths.) It also specifies a destination which is at the moment reserved to be the value synapses (as in, the input injects current to the cell just like synapses do). Finally, the alternative name inputW can also be used to apply a weight property on the amplitude of the specific <input> instance. (See how it’s done below).

As we’ll see in the following chapter on spatially detailed cells, the specific segment and fractionAlong (the segment) to apply the <input> on can also be optionally specified, otherwise the input will be placed on the “soma” of the cell.

Let’s see then, how to instantiate a current clamp and apply it on the cell in the following.

%%writefile SingleCell_Model.nml
<!-- These properties on the top level <neuroml> tag are just to positively identify the file as NeuroML.  They are not mandatory -->
<neuroml xmlns="http://www.neuroml.org/schema/neuroml2"  xmlns:xs="http://www.w3.org/2001/XMLSchema" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:schemaLocation="http://www.neuroml.org/schema/neuroml2 https://raw.github.com/NeuroML/NeuroML2/development/Schemas/NeuroML2/NeuroML_v2.1.xsd">
<!-- Here comes the cell type -->
<iafCell id="MyFirstCellType" C="200 pF" leakConductance="10 nS" leakReversal="-70 mV" reset="-70mV" thresh="-50mV" />
<!-- Here comes the input type ❗ -->
<pulseGenerator id="MyFirstInputSource" delay="100ms" duration="500ms" amplitude="0.21nA"/>
<!-- Here comes the network -->
<network id="MyNetwork" type="networkWithTemperature" temperature="37degC" >
    <!-- Add a population with a single cell -->
    <population id="MyFirstPopulation" component="MyFirstCellType" size="1"/>
    <!-- Apply an inputList on the population ❗ -->
    <inputList id="MyFirstInputList" population="MyFirstPopulation" component="MyFirstInputSource">
        <!-- With one stim placement on it -->
        <input id="0" target="MyFirstPopulation[0]" destination="synapses"/>
    </inputList>
</network>
</neuroml>

Overwriting SingleCell_Model.nml

And run the simulation with the updated model file:

results = runEden('LEMS_SingleCell_Sim.xml')
plt.plot(results['t'],results['MyFirstPopulation[0]/v'])
plt.xlabel('Time (sec)'); plt.ylabel('Membrane voltage (V)');

That’s more like it.

Change the amplitude of the input source (or equivalently ​weight​ with inputW), repeat with the new model files and see how the firing rate changes.

There are more input types defined in base NeuroML (see for example, ramp and sine current sources and voltage clamps), and custom ones can be built as LEMS components in a following section.

Adding more neurons and synapses

This section shows how to make models with more than one neuron. This way, one can explore in detail the mechanism of interaction between a pair of neurons, or set up whole populations of neurons in an inter-connected network and explore the group dynamics.

First of all, let’s extend the first population from out previous model to contain two neurons this time. This is quite easy to do, we can just set the size of the population to what we want.

Then, we’ll connect these two cells with a synapse.

This is done by defining the synapse’s model as an synapse type, and adding to the <network> a synaptic <projection> with synapses of said type (called the component for this <projection>).

Instead of specifying the synapse type from scratch (as we’ll be doing in a following section), we’ll use the common exponential-decay conductance model as a base, and supply our own dynamical parameters for it. This synapse type can be found in the standard NeuroML library as <expOneSynapse> and has parameters erev(ersal potential of the active synapse), gbase(peak post-synaptic conductance per spike), and decay (time constant).

We’ll set them as follows: erev = 0 mV, gbase = 10 nmho, decay = 5 ms for educational purposes.

Next, we’ll add a <projection> to the <network>. Each <projection> contains one or more copies (called <connection>s) of the same synapse type, and is applied from a certain cell presynapticPopulation to a certain postsynapticPopulation.

Each <connection> has a distinct serial id and its pre and post site. Each site is located using the CellId (of the pre or post population). The alternative names connectionW and connectionWD can also be used to apply a weight property on the amplitude of the specific <connection> for the former, and weight as well as propagation delay for the latter.

As we’ll see in the following section on spatially detailed cells, the specific SegmentId and FractionAlong (the segment) to apply the <synapse> between can also be optionally specified for the pre and post sites, otherwise the sites will be located on the “soma” of each corresponding cell.

Fun fact: One can even model a self-synapse that way.

Let’s see then, how to instantiate a synapse and connect our two cells with it in the following. We’ll use a huge delay of 10 msec to make its effect visible.

%%writefile MultiCell_Model.nml
<neuroml xmlns="http://www.neuroml.org/schema/neuroml2"  xmlns:xs="http://www.w3.org/2001/XMLSchema" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:schemaLocation="http://www.neuroml.org/schema/neuroml2 https://raw.github.com/NeuroML/NeuroML2/development/Schemas/NeuroML2/NeuroML_v2.1.xsd">
<iafCell id="MyFirstCellType" C="200 pF" leakConductance="10 nS" leakReversal="-70 mV" reset="-70mV" thresh="-50mV" />
<pulseGenerator id="MyFirstInputSource" delay="100ms" duration="500ms" amplitude="0.21nA"/>

<expOneSynapse id="MyFirstSynapseType" gbase="10nS" erev="0V" tauDecay="5ms"/>
<network id="MyNetwork" type="networkWithTemperature" temperature="37degC" >

    <!-- Add a population with TWO cells ❗ -->
    <population id="MyFirstPopulation" component="MyFirstCellType" size="2"/>

    <!-- Apply an inputList on the population -->
    <inputList id="MyFirstInputList" population="MyFirstPopulation" component="MyFirstInputSource">
        <!-- With one stim placement on it -->
        <input id="0" target="MyFirstPopulation[0]" destination="synapses"/>
    </inputList>

    <!-- Apply a synaptic projection ❗ -->
    <projection id="MyFirstApProjection" presynapticPopulation="MyFirstPopulation" postsynapticPopulation="MyFirstPopulation" synapse="MyFirstSynapseType">
        <!-- With one connection on it -->
        <connectionWD id="0" preCellId="0" postCellId="1" weight="1" delay="10 ms"/>
    </projection>

</network>
</neuroml>

Writing MultiCell_Model.nml

There are more synapse types defined in base NeuroML, (see for example, alphaCurrentSynapse and expThreeSynapse), and custom ones can be built as LEMS components in a following section.

In order to observe the two cells, we’ll have to record the new cell as well. This is done by adding one more <OutputColumn>s to the <Simulation>, as follows:

%%writefile LEMS_MultiCell_Sim.xml
<?xml version="1.0" encoding ="UTF-8" ?>
<Lems>
<include file="MultiCell_Model.nml" />
<Simulation id="MySim" length="1 s" step="100 us" target="MyNetwork" >
    <OutputFile id="MyFirstOutputFile" fileName="results.gen.txt">
        <!-- Add TWO outputColumns, one for each cell ❗ -->
        <OutputColumn id="my_first_vm"  quantity="MyFirstPopulation[0]/v"/>
        <OutputColumn id="my_second_vm" quantity="MyFirstPopulation[1]/v"/>
    </OutputFile>
</Simulation>
<Target component="MySim"/>
</Lems>

Writing LEMS_MultiCell_Sim.xml

Let’s run this simulation:

results = runEden('LEMS_MultiCell_Sim.xml')
plt.figure(figsize=(10,5))
plt.plot(results['t'],results['MyFirstPopulation[0]/v'], label='My first cell')
plt.plot(results['t'],results['MyFirstPopulation[1]/v'], label='My second cell')
plt.xlabel('Time (sec)'); plt.ylabel('Membrane voltage (V)'); plt.legend(); plt.xlim((0, .8))

(0.0, 0.8)

Adding more cell types and spike sources

There are many SNN models that are formed by one population of similar cells, but there also exist heterogeneous networks comprising multiple populations of different cells, connected by various synaptic projections. Moreover, some of these populations may receive action potentials from extra-model spike sources.

In this section, we’ll introduce another point-neuron type and its corresponding population, see how to access different cell populations, and add some external spike input to our network as well.

Adding a new cell type and population

Adding another point-neuron type and population is straightforward: We’ll introduce the new neuron type with a tag outside the <network> definition, and add a second population inside the <network>, just like we did with our first cell type and population.

This time, we’ll use the common and very expressive Izhikevich cell model (a simple non-linear system that captures many waveform types) as a base, and supply our own dynamical parameters for it.

This input type can be found in the standard NeuroML library as <izhikevich2007Cell> (an update of the original <izhikevichCell> model that used mostly dimensionless variables) and has parameters C(apacitance), a (rate of adaptation), b (coupling of adaptation to voltage), c (after-spike reset voltage), d (after-spike adaptation increment), k(onstant of voltage growth rate), v0 (initial membrane voltage), vpeak (max voltage before reset), vr (resting membrane potential), and vt (threshold potential).

After consulting the Izhikevich model gallery repo on OpenSourceBrain that shows different activity types with their corresponding parameters, we select the “Regular Spiking” activity type and modify it a bit to show interaction between neurons. We’ll set them as follows:

C = 100 pF, a = 0.03 Hz⁻¹, b = -2 nS, c = -50 mV, d = 100 pA, k = 0.7 μS/V, vr = -60 mV, vpeak = 35 mV, vt = -40 mV, v0 = vr

%%writefile MultiCell_Model.nml
<neuroml xmlns="http://www.neuroml.org/schema/neuroml2"  xmlns:xs="http://www.w3.org/2001/XMLSchema" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:schemaLocation="http://www.neuroml.org/schema/neuroml2 https://raw.github.com/NeuroML/NeuroML2/development/Schemas/NeuroML2/NeuroML_v2.1.xsd">
<iafCell id="MyFirstCellType" C="200 pF" leakConductance="10 nS" leakReversal="-70 mV" reset="-70mV" thresh="-50mV" />

<!-- Add a new a new cell type ❗ -->
<izhikevich2007Cell id="MySecondCellType" v0 = "-60mV" C="100 pF" k = "0.7 nS_per_mV"
                    vr = "-60 mV" vt = "-40 mV" vpeak = "35 mV"
                    a = "0.03 per_ms" b = "-2 nS" c = "-50 mV" d = "100 pA">
    <notes>Regular spiking cell</notes>
</izhikevich2007Cell>

<pulseGenerator id="MyFirstInputSource" delay="100ms" duration="500ms" amplitude="0.21nA"/>
<expOneSynapse id="MyFirstSynapseType" gbase="10nS" erev="0V" tauDecay="5ms"/>

<network id="MyNetwork" type="networkWithTemperature" temperature="37degC" >
    <population id="MyFirstPopulation"  component="MyFirstCellType"  size="2"/>

    <!-- Add another population with two cells ❗ -->
    <population id="MySecondPopulation" component="MySecondCellType" size="2"/>

    <!-- Apply an inputList on the first population -->
    <inputList id="MyFirstInputList" population="MyFirstPopulation" component="MyFirstInputSource">
        <!-- With one stim placement on it -->
        <input id="0" target="MyFirstPopulation[0]" destination="synapses"/>
    </inputList>

     <!-- Apply an inputList on the SECOND population ❗ -->
    <inputList id="MySecondInputList" population="MySecondPopulation" component="MyFirstInputSource">
        <!-- With two stim placements on it, one for each cell -->
        <inputW id="0" target="0" destination="synapses" weight="0.4"/>
        <inputW id="1" target="1" destination="synapses" weight="0.3"/>
    </inputList>

    <!-- Apply the AP synaptic projection -->
    <projection id="MyFirstApProjection" presynapticPopulation="MyFirstPopulation" postsynapticPopulation="MyFirstPopulation" synapse="MyFirstSynapseType">
        <!-- With one connection on it -->
        <connectionWD id="0" preCellId="0" postCellId="1" weight="1" delay="10 ms"/>
    </projection>

</network>
</neuroml>

Overwriting MultiCell_Model.nml

And let’s add more one more <OutputColumn>s to the <Simulation> to keep track of the new cells.

Setting the randomisation seed

Also, because we’ll be using simulated random input in the following, we can set another parameter of the <Simulation> to “fix” the random numbers that will be pulled out of the hat and hence get the same result for a given model, every time.

• This is the optional seed that should be an integer number. If seed is not specified, then EDEN will pick the time and date as the seed number, making all the random numbers change every time.

• Important: Due to multiple technological reasons, different computing setups come with different numerical perturbations. Being (usually) chaotic systems, neural networks may show different activity from the same initial state, because of these perturbations; though a good model’s qualitative function should not be affected by this. Read more about the limitations of numerical reproducibility, in a frequently answered question.

Caution

Please do not rely on the seed to make your model work, random variables should be random. A working model should be working for at least most of the seeds, always check before publishing.

%%writefile LEMS_MultiCell_Sim.xml
<?xml version="1.0" encoding ="UTF-8" ?>
<Lems>
<include file="MultiCell_Model.nml" />
<!-- Set the SEED to freeze the random numbers ❗ -->
<Simulation id="MySim" length="1 s" step="100 us" seed="20190109" target="MyNetwork" >
    <OutputFile id="MyFirstOutputFile" fileName="results.gen.txt">
        <OutputColumn id="my_first_vm"  quantity= "MyFirstPopulation[0]/v"/>
        <OutputColumn id="my_second_vm" quantity= "MyFirstPopulation[1]/v"/>
        <!-- Add TWO more outputColumns, one for each new cell ❗ -->
        <OutputColumn id="my_third_vm"  quantity="MySecondPopulation[0]/v"/>
        <OutputColumn id="my_fourth_vm" quantity="MySecondPopulation[1]/v"/>
    </OutputFile>
</Simulation>
<Target component="MySim"/>
</Lems>

Overwriting LEMS_MultiCell_Sim.xml

And see our new cells behave. They exhibit an actual AP shape, whereas the LIF cells reset right on the threshold-to-spike instead of the spike’s peak.

def show_multicell_results(results):
    from matplotlib import pyplot as plt;
    plt.figure(figsize=(10,5))
    plt.plot(results['t'],results[ 'MyFirstPopulation[0]/v'], label='LIF #0, v')
    plt.plot(results['t'],results[ 'MyFirstPopulation[1]/v'], label='LIF #1, v')
    plt.plot(results['t'],results['MySecondPopulation[0]/v'], label='Izh #0, v', linestyle='--')
    plt.plot(results['t'],results['MySecondPopulation[1]/v'], label='Izh #1, v', linestyle='--')
    plt.xlabel('Time (sec)'); plt.ylabel('Membrane voltage (V)'); plt.legend(); plt.xlim((0, .8))

results = runEden('LEMS_MultiCell_Sim.xml')
show_multicell_results(results)

Adding abstract spike sources and applying them to the network

As we noted previously, our virtual neurons may need to be stimulated with realistic input; this way that we can watch and play with them at a “natural-like” state, instead of crushing loneliness on a glass dish. Of course, these stimuli are normally provided by adjacent neurons, which motivates us to model those neurons as well as well as the stimuli that these neurons receive… At some point we end up out of scope and/or technical resources, and will have to replace external neuron input with a rougher abstraction. One such abstraction is entities that emit “firing events” independently of what the network does, according to a mathematical model. These firing events can then activate synapses on target cells in our neural network.

Let’s see how this works in NeuroML, by defining two different spike source types, adding populations of said types in the network, and delivering their potentials through projections.

We will add two spike source types from the NeuroML library: a deterministic <spikeArray> and a stochastic <SpikeSourcePoisson>:

• For the former, we will specify the exact occurrence time of every single spike by adding <spike> child tags inside it.

• For the latter, we’ll supply parameters as usual in the previous parts, namely: start = 200 ms, duration = 300 ms, rate = 20 Hz.

Next, we’ll instantiate them with <population> tags as if they were cells. (It so happens.)

Finally, we’ll add <projection>s with synapse connections from those not-cell populations, to real cells.

Note: we cannot specify spike sources as post-synaptic targets since they are not neurons, and hence we can’t just attach biophysical synapses on them. Neither can we connect them with graded synapses for that matter.

%%writefile MultiCell_Model.nml
<neuroml xmlns="http://www.neuroml.org/schema/neuroml2"  xmlns:xs="http://www.w3.org/2001/XMLSchema" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:schemaLocation="http://www.neuroml.org/schema/neuroml2 https://raw.github.com/NeuroML/NeuroML2/development/Schemas/NeuroML2/NeuroML_v2.1.xsd">
<iafCell id="MyFirstCellType" C="200 pF" leakConductance="10 nS" leakReversal="-70 mV" reset="-70mV" thresh="-50mV" />

<izhikevich2007Cell id="MySecondCellType" v0 = "-60mV" C="100 pF" k = "0.7 nS_per_mV"
                    vr = "-60 mV" vt = "-40 mV" vpeak = "35 mV"
                    a = "0.03 per_ms" b = "-2 nS" c = "-50 mV" d = "100 pA">
    <notes>Regular spiking cell</notes>
</izhikevich2007Cell>
<pulseGenerator id="MyFirstInputSource" delay="100ms" duration="500ms" amplitude="0.21nA"/>

<!-- Add new spike sources ❗ -->
<spikeArray id="MySpikeList">
    <spike id="0" time="210 ms"/>
    <spike id="1" time="350 ms"/>
    <spike id="2" time="450 ms"/>
</spikeArray>
<SpikeSourcePoisson id="MyRandomSpikeSource" start="200ms" duration="300ms" rate="15Hz"/>

<expOneSynapse id="MyFirstSynapseType" gbase="10nS" erev="0V" tauDecay="5ms"/>

<network id="MyNetwork" type="networkWithTemperature" temperature="37degC" >
    <population id="MyFirstPopulation"  component="MyFirstCellType"  size="2"/>
    <population id="MySecondPopulation" component="MySecondCellType" size="2"/>

    <!-- Add "populations" of SPIKE SOURCES beside the cells ❗ -->
    <population id="Spikes_Counted" component="MySpikeList"         size="1" />
    <population id="Spikes_Random"  component="MyRandomSpikeSource" size="1" />

    <!-- Apply an inputList on the first population -->
    <inputList id="MyFirstInputList" population="MyFirstPopulation" component="MyFirstInputSource">
        <!-- With one stim placement on it -->
        <input id="0" target="MyFirstPopulation[0]" destination="synapses"/>
    </inputList>
     <!-- Apply an inputList on the second population -->
    <inputList id="MySecondInputList" population="MySecondPopulation" component="MyFirstInputSource">
        <!-- With two stim placements on it, one for each cell -->
        <inputW id="0" target="0" destination="synapses" weight="0.4"/>
        <inputW id="1" target="1" destination="synapses" weight="0.3"/>
    </inputList>

    <!-- Apply the AP synaptic projection -->
    <projection id="MyFirstApProjection" presynapticPopulation="MyFirstPopulation" postsynapticPopulation="MyFirstPopulation" synapse="MyFirstSynapseType">
        <!-- With one connection on it -->
        <connectionWD id="0" preCellId="0" postCellId="1" weight="1" delay="10 ms"/>
    </projection>

    <!-- And, finally, project synapses from spike sources to cells ❗ -->
    <projection id="MyFirstSourceProjection" presynapticPopulation="Spikes_Counted" postsynapticPopulation="MySecondPopulation" synapse="MyFirstSynapseType">
        <connection id="0" preCellId="0" postCellId="0"/>
    </projection>
    <projection id="MySecondSourceProjection" presynapticPopulation="Spikes_Random" postsynapticPopulation="MyFirstPopulation" synapse="MyFirstSynapseType">
        <connection id="0" preCellId="0" postCellId="1"/>
    </projection>
</network>
</neuroml>

Overwriting MultiCell_Model.nml

There is this interesting detail about spike sources: although functionally they are inputs and are classed as such in the NeuroML specification, model-wise they act like a lesser form of neurons and in fact they are instantiated in the network as such.

There are more spike sources defined in base NeuroML, (see for example, the regular spikeGenerator and the uniformly-distributed spikeGeneratorRandom). Custom ones can be built as LEMS components in a following section.

If the desired spike series are too complicated to describe with a base NeuroML tag pr with a LEMS spike source, they can also be generated outside EDEN and streamed into the simulation, using the IO extension. The same applies for other input, in the form of biophysical time series.

Let’s see what happens.

results = runEden('LEMS_MultiCell_Sim.xml')
show_multicell_results(results)

The spikes are applied to one cell of each type cells, and the exact times of occurrence are seen as very abrupt increases (i.e. depolarisations) in membrane potential.

Observe that the events from the <spikeArray> trigger the synapses at the precise times specified in the <spike>s, whereas the spikes from the <SpikeSourcePoisson> are stochastically distributed.

Remove the seed mentioned above and re-run the cell, to see how the <spikeArray> will remain the same and the <SpikeSourcePoisson> will vary.

Adding graded interactions

In the previous sections, we added a synaptic connection that follows the traditional “action potential” model: if the pre-synaptic neuron passes a threshold, then an invariant quantity of neurotransmitters is released, the neurotransmitters slowly cross the synaptic gap and the post-synaptic potential is applied to the post-neuron. Information travels one way (pre 🡒 post) only, exclusively at the instants that the pre-neuron fires.

(See also Figure 3.1 from the Neuronal Dynamics book.)

This is not the only way that neurons are seen to interact in the wild. Transmitter release may happen gradually with membrane potential (it’s not all or nothing), and neurons often also interact through graded potentials. These are ways that neurons engage in continuous interaction, and we can also model that in NeuroML.

In this section, we will add graded interactions between cells in the network, and look closely at their effects.

Caution

Use graded synapses only with neurons that have realistic action potential shapes (not LIF, for example), lest you end up with wrong conclusions.

We will add two different graded synaptic component types from the NeuroML library: an ohmic <gapJunction> model, and a dummy <silentSynapse> model that applies no effect (for one-way interactions).

As usual, we’ll define the synapse types inside the <neuroml> and outside the <network>, but this time we’ll different projection tags to state that unlike the previous ‘classical’ synapses, these are bi-directional continuous connections. The new tags used are <electricalProjection> for symmetric connections and <continuousProjection> for not necessarily symmetric connections. These contain different tags, <electricalConnectionInstance> and <continuousConnectionInstance> respectively.

The difference with <projection> and <connection> is that the synapse is not an attribute of the projection but individually specified for each connection instance:

• <electricalConnectionInstance>s have a synapse type parameter. This synapse gets instantiated twice, once for each side of the connection, and the two halves are symmetrically attached to the (only nominally) pre and post synaptic sites.

• <continuousConnectionInstance>s have two different synapse type parameters, preComponent and postComponent. In this case, the pre component may be different from the post component, each is instantiated and attached on its respective side of the synapse.

So then, the tags are slightly different from the classic-type <connection> and the {pre,post}synapticComponents are specified individually for each synapse, that’s because the style changed during the development of NeuroML. EDEN may be rather permissive on mixing and matching from either style. Just remember that you may have to use a slightly differently style for action and graded potentials when making your NeuroML.

In order to specify weight, <electricalConnectionInstance> and <continuousConnectionInstance> can be replaced with the similar <electricalConnectionInstanceW> and <continuousConnectionInstanceW> tags.

• There also are two equivalent but a little bit shorter tags <electricalConnection> and <continuousConnection>; but tags like <electricalConnectionW> and <continuousConnectionW> are not in the standard and will be just ignored by EDEN, watch out.

Note that these continuous connections can’t have a delay associated with them, though that would help make neural mass models. If you are interested in continuous interactions with pure lag, do state your requirements in a feature request for NeuroML.

%%writefile MultiCell_Model.nml
<neuroml xmlns="http://www.neuroml.org/schema/neuroml2"  xmlns:xs="http://www.w3.org/2001/XMLSchema" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:schemaLocation="http://www.neuroml.org/schema/neuroml2 https://raw.github.com/NeuroML/NeuroML2/development/Schemas/NeuroML2/NeuroML_v2.1.xsd">
<iafCell id="MyFirstCellType" C="200 pF" leakConductance="10 nS" leakReversal="-70 mV" reset="-70mV" thresh="-50mV" />

<!-- Add a new a new cell type -->
<izhikevich2007Cell id="MySecondCellType" v0 = "-60mV" C="100 pF" k = "0.7 nS_per_mV"
                    vr = "-60 mV" vt = "-40 mV" vpeak = "35 mV"
                    a = "0.03 per_ms" b = "-2 nS" c = "-50 mV" d = "100 pA">
    <notes>Regular spiking cell</notes>
</izhikevich2007Cell>

<pulseGenerator id="MyFirstInputSource" delay="100ms" duration="500ms" amplitude="0.21nA"/>
<expOneSynapse id="MyFirstSynapseType" gbase="10nS" erev="0V" tauDecay="5ms"/>

<!-- Add a new synapse type, a gap junction ❗ -->
<gapJunction id="Syn_gapJunction" conductance="2000pS"/> <!-- With a pretty huge conductance, to demonstrate the effect -->

<!-- Add a new synapse type, a dummy graded synapse ❗ -->
<silentSynapse id="Syn_SilentSynapse"/> <!-- No parameters, no effect, no worries -->

<network id="MyNetwork" type="networkWithTemperature" temperature="37degC" >
    <population id="MyFirstPopulation"  component="MyFirstCellType"  size="2"/>
    <population id="MySecondPopulation" component="MySecondCellType" size="2"/>

    <!-- Apply an inputList on the first population -->
    <inputList id="MyFirstInputList" population="MyFirstPopulation" component="MyFirstInputSource">
        <!-- With one stim placement on it -->
        <input id="0" target="MyFirstPopulation[0]" destination="synapses"/>
    </inputList>
     <!-- Apply an inputList on the second population -->
    <inputList id="MySecondInputList" population="MySecondPopulation" component="MyFirstInputSource">
        <!-- With two stim placements on it, one for each cell -->
        <inputW id="0" target="0" destination="synapses" weight="0.4"/>
        <inputW id="1" target="1" destination="synapses" weight="0.3"/>
    </inputList>

    <!-- Apply the AP synaptic projection -->
    <projection id="MyFirstApProjection" presynapticPopulation="MyFirstPopulation" postsynapticPopulation="MyFirstPopulation" synapse="MyFirstSynapseType">
        <!-- With one connection on it -->
        <connectionWD id="0" preCellId="0" postCellId="1" weight="1" delay="10 ms"/>
    </projection>

    <!-- Apply the new graded synaptic projections ❗ -->
    <electricalProjection id="MyGapProjection" presynapticPopulation="MySecondPopulation" postsynapticPopulation="MySecondPopulation">
        <electricalConnectionInstanceW id="0" preCell="0" postCell="1" preSegment="0" postSegment="0" synapse="Syn_gapJunction" weight="1"/> <!-- seg 0 is the only one for point neurons -->
    </electricalProjection>
    <continuousProjection id="MyUnsymmetricGap" presynapticPopulation="MySecondPopulation" postsynapticPopulation="MyFirstPopulation">
        <continuousConnection id="0" preCell="1" postCell="0" preComponent="Syn_SilentSynapse" postComponent="Syn_gapJunction"/>
    </continuousProjection>

</network>
</neuroml>

Overwriting MultiCell_Model.nml

The standard NeuroML library is a bit sparse in continuous synapses: more than the component types we just used, it includes one detailed chemical receptor model that has been used for modelling the pyloric rhythm of the crustacean stomatogastric ganglion. But this is only because there are no other typical, commonly agreed interaction models yet. To implement your preferred interaction model, see the LEMS tutorial and the implementation of said STG model, and EDEN’s multi-flux capability for more options.

Let’s run this simulation:

results = runEden('LEMS_MultiCell_Sim.xml')
show_multicell_results(results)

Observe how the membrane voltage each Izhikevich cell is sharply affected whenever the other one fires; the one-way graded potential applied to the first LIF cell has an even more pronounced effect. The steady leak due to different resting potentials for the two cells is also visible at the start (before the current clamps activate), feel free to change vr to same in both cells if it’s a concern.

NeuroML file structure conventions

In the previous examples, we used two files per model. This is not the only option; we could also have split the parts of the simulation differently, as long as the file for each part is <Include>d by the file we pointed the simulator to in runEden, or by a another file <Include>d by a file that’s ultimately <Include>d by the simulation file, and so forth. Or (when using EDEN) we could just roll all model parts side by side with the <Simulation> in a single file if that’s more convenient. (By the way, if the same file is <Include>d by different files, that’s no problem.)

As a neural modelling project gets more detailed, and its parts get re-mixed in different ways (also between projects), there are some NeuroML file namingconventions that may help manage the model contents:

• each <cell> type is defined in its own file with extension .cell.nml;

• each synapse type used is defined in its own file with extension .synapse.nml;

• the <network> to be simulated is placed in its own file with extension .net.nml;

• the <simulation> information are placed in their own file with a name like .xml.

• There is no common convention for <input>s; they are usually added to the .net.nml files, but they could also be added to NeuroML files in some other way that makes sense.

Of course, if there are overriding reasons the strcuture can be different in practice, as mentioned previously.

For example, if a set of different model parts is also meaningful as a named collection, there could be an additional .nml file that <Includes> the files for said parts.

The following code cells show how the parts comprising the last model can also be laid out in a more modular file structure.

%%writefile MultiCell.net.nml
<neuroml xmlns="http://www.neuroml.org/schema/neuroml2"  xmlns:xs="http://www.w3.org/2001/XMLSchema" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:schemaLocation="http://www.neuroml.org/schema/neuroml2 https://raw.github.com/NeuroML/NeuroML2/development/Schemas/NeuroML2/NeuroML_v2.1.xsd">

<!-- Include all the re-usable components ❗ -->
<include file="MultiCell_Iaf.cell.nml" />
<include file="MultiCell_Izh.cell.nml" />
<include file="SimpleChemical.synapse.nml" />
<include file="SimpleGapJunction.synapse.nml" />
<include file="Graded_Dummy.synapse.nml" />

<!-- Leave the input sources here i guess, or not ... -->
<pulseGenerator id="MyFirstInputSource" delay="100ms" duration="500ms" amplitude="0.21nA" />

<!-- And leave here just the network ❗ (and auxiliaries) -->
<network id="MyNetwork" type="networkWithTemperature" temperature="37degC" >
    <population id="MyFirstPopulation"  component="MyFirstCellType"  size="2"/>
    <population id="MySecondPopulation" component="MySecondCellType" size="2"/>

    <inputList id="MyFirstInputList" population="MyFirstPopulation" component="MyFirstInputSource">
        <input id="0" target="MyFirstPopulation[0]" destination="synapses"/>
    </inputList>
    <inputList id="MySecondInputList" population="MySecondPopulation" component="MyFirstInputSource">
        <inputW id="0" target="0" destination="synapses" weight="0.4"/>
        <inputW id="1" target="1" destination="synapses" weight="0.3"/>
    </inputList>

    <projection id="MyFirstApProjection" presynapticPopulation="MyFirstPopulation" postsynapticPopulation="MyFirstPopulation" synapse="MyFirstSynapseType">
        <connectionWD id="0" preCellId="0" postCellId="1" weight="1" delay="10 ms"/>
    </projection>

    <electricalProjection id="MyGapProjection" presynapticPopulation="MySecondPopulation" postsynapticPopulation="MySecondPopulation">
        <electricalConnectionInstanceW id="0" preCell="0" postCell="1" preSegment="0" postSegment="0" synapse="Syn_gapJunction" weight="1"/> <!-- seg 0 is the only one for point neurons -->
    </electricalProjection>
    <continuousProjection id="MyUnsymmetricGap" presynapticPopulation="MySecondPopulation" postsynapticPopulation="MyFirstPopulation">
        <continuousConnection id="0" preCell="1" postCell="0" preComponent="Syn_SilentSynapse" postComponent="Syn_gapJunction"/>
    </continuousProjection>
</network>
</neuroml>

Writing MultiCell.net.nml

%%writefile MultiCell_Iaf.cell.nml
<neuroml>
<iafCell id="MyFirstCellType" C="200 pF" leakConductance="10 nS" leakReversal="-70 mV" reset="-70mV" thresh="-50mV" />
</neuroml>

Writing MultiCell_Iaf.cell.nml

%%writefile MultiCell_Izh.cell.nml
<neuroml>
<izhikevich2007Cell id="MySecondCellType" v0 = "-60mV" C="100 pF" k = "0.7 nS_per_mV"
                    vr = "-60 mV" vt = "-40 mV" vpeak = "35 mV"
                    a = "0.03 per_ms" b = "-2 nS" c = "-50 mV" d = "100 pA">
    <notes>Regular spiking cell</notes>
</izhikevich2007Cell>
</neuroml>

Writing MultiCell_Izh.cell.nml

%%writefile SimpleChemical.synapse.nml
<neuroml>
<expOneSynapse id="MyFirstSynapseType" gbase="10nS" erev="0V" tauDecay="5ms"/>
</neuroml>

Writing SimpleChemical.synapse.nml

%%writefile SimpleGapJunction.synapse.nml
<neuroml>
<gapJunction id="Syn_gapJunction" conductance="2000pS"/> <!-- With a pretty huge conductance, to demonstrate the effect -->
</neuroml>

Writing SimpleGapJunction.synapse.nml

%%writefile Graded_Dummy.synapse.nml
<neuroml>
<silentSynapse id="Syn_SilentSynapse"/> <!-- No parameters, no effect, no worries -->
</neuroml>

Writing Graded_Dummy.synapse.nml

Make another LEMS file that points to the new .net.nml:

%%writefile LEMS_MultiCell_Modular_Sim.xml
<?xml version="1.0" encoding ="UTF-8" ?>
<Lems>
<include file="MultiCell.net.nml" />
<!-- Set the seed to freeze the random numbers -->
<Simulation id="MySim" length="1 s" step="100 us" seed="20190109" target="MyNetwork" >
    <OutputFile id="MyFirstOutputFile" fileName="results.gen.txt">
        <OutputColumn id="my_first_vm"  quantity= "MyFirstPopulation[0]/v"/>
        <OutputColumn id="my_second_vm" quantity= "MyFirstPopulation[1]/v"/>
        <OutputColumn id="my_third_vm"  quantity="MySecondPopulation[0]/v"/>
        <OutputColumn id="my_fourth_vm" quantity="MySecondPopulation[1]/v"/>
    </OutputFile>
</Simulation>
<Target component="MySim"/>
</Lems>

Writing LEMS_MultiCell_Modular_Sim.xml

And re-run the modular version of the model.

results = runEden('LEMS_MultiCell_Modular_Sim.xml')
show_multicell_results(results)

---

This concludes the introduction to the basics of NeuroML. You may proceed to lean more about spatially detailed cells and making your own mechanisms with LEMS.