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

Extending NeuroML with LEMS components

NeuroML is a rich language for describing neural models. However, the mechanisms in its standard library are not that extensive and to be fair, it would not be possible to anticipate all the possible mechanisms that could be proposed in neuroscience. Hence, NeuroML works together with another modelling language, LEMS. LEMS can describe arbitrary dynamical equations for mechanisms (arbitrary as in, any ODE’s plus instantaneous changes). In a sense, it is similar to the blockdiagram models seen in, e.g. MathWorks Simulink™.

Within NeuroML, LEMS may be used to describe new types of cells, synapses and much more (see below for the full list).

This article will introduce LEMS components, and how these can be used to create custom neural dynamics. The various NeuroML-related aspects of LEMS components are presented, and then an example is given on how to write custom dynamics for each sort of neural mechanism (neuron, synapse, &c.)

What is a LEMS component?

A LEMS component is a virtual physical entity (or mechanism) that interacts with its environment in ways described my mathematical expressions. Its internal state is captured by a number of scalar variables, and its interactions with the environment are in terms of scalar physical quantities or instantaneous events that it is affected by or determines. It may also include other physical quantities called parameters that remain fixed throughout a simulation, but still influence the state’s dynamics.

Side note: The scalar quantities can be used to implement continuous in space (also known as “density-type”) mechanisms, since all models are discretised to run on the computer; refer to the guide’s article on spatially detailed cells, and also the examples about modelling fields.

Within its Requirements (what comes in) and the Exposures (what comes out), a LEMS component behaves like a black box; this means that it does not otherwise interact with its surroundings. This property is helpful to keep track of what goes where in complex systems.

It is not necessary for a LEMS component to have state variables: a stateless component can determine and expose quantities that are a function of the component’s Requirements at each point in time. This the case for instantaneous ion channel gates, gap junctions and other “static” mechanisms, where the e.g. electrical conductance is an important parameter.

LEMS components are added to NeuroML files just like NeuroML-standard mechanisms. A tag named after their type is added inside the <neuroml> but not inside the <network>, and the tag’s attributes specify the parameters that determine its behaviour. Indeed, the NeuroML-standard mechanisms we’ve used in the previous articles also have LEMS descriptions.

What remains then, is to describe the types of the components.

What is a LEMS ComponentType?

Besides <Component>, the other core concept of LEMS is <ComponentType>. This describes the external interface and the internal structure of components, and the mathematical expressions that govern their behaviour.

Their usefulness comes from the frequent occasion where many types of mechanism can be captured by the same mathematical equations. A famous case is the Izhikevich neuron model where the same equations capture many and much different neural activity patterns: all that changes is the values of some parameters. By describing a ComponentType we can thus specify the essential dynamics once, instead of re-defining all the properties and behaviour for each component instance. ComponentTypes are also very useful when archiving or exchanging NeuroML models, since every user (or use case) may require a different set of parameters for the same general model type.

Thus, the ComponentType is defined once and then specialised, concrete components are defined as follows:

<ComponentType name="MyNewSynapseType">
    ...
    <Dynamics>
        ...
    </Dynamics>
    ...
</ComponentType>

...

<MyNewSynapseType this_parameter="60 ms" that_parameter="37 degC">

    or, equivalently:

<Component type="MyNewSynapseType"
    this_parameter="60 ms" that_parameter="37 degC">

Note that the components can indeed be further replicated, as in the case of multiple <connection>s in a <projection>, or <input>s in an <inputList>. In case that individual instances need to have different parameters within the same <inputList> or <projection>, this type of variability can be specified through EDEN’s <CustomSetup> extension.

Now let’s see how to write a new ComponentType.

The “components” of a ComponentType

Here is a reference of the LEMS properties that can be used to make NeuroML mechanisms (as far as EDEN is concerned).

• Attribute name; that’s used to refer to the ComponentType and hence must be unique for each.

• Attribute extends; This specifies the ComponentType as a ‘derived’ type that inherits properties from another, more abstract ComponentType. This way it’s not necessary to re-write the whole ComponentType from scratch. Refer to the NeuroML reference (check sidebar) for the standard LEMS components, and the NML component list below for a per-type analysis.

  Note: Components should extend a relevant base type of neuron, synapse, etc. to be considered by EDEN for said roles.

• Attribute description; that’s optional and up to taste. It is better than XML comments for explaining what the ComponentType does or represents. Note that the other LEMS tags may also accept a description attribute as well.

• Interface tags:

  • <Requirement>: This is a quantity that the component needs from its surroundings, to determine its own behaviour. For example, an ion channel may open or close depending on membrane potential.

    • Attributes are the name identifier (which mush match with an available <Exposure>), and the dimension (in the sense used in Physics) of said quantity.

  • <Exposure>: This is a quantity that the component determines based on its requirements and internal state. Information about this quantity may be thus exposed to related parts of the model. For example, an ion channel may determine and expose the electrical current flux through it.

    • Attributes are the name identifier (which must match with the <Requirement>s that dependent components have) and the dimension of the exposed quantity.

  • <EventPort>: This is a named “port” (or “letter box” for lack of a better word) that components can receive the occurrence of discrete events through, or alternatively emit discrete events In a sense, it is similar to a <Requirement> or <Exposure> for discrete events rather than quantities over time. For example, a post-synapse waits for the firing of the pre-synaptic neuron that’s a discrete event, and a (site of a) neuron may announce when a “spike” has happened on it. Attributes are the name identifier, and the direction which can be either in (the port receives events) or out (the port sends events).

• Static attributes tags:

  • <Parameter>: This is a quantity about the component, that remains the same for the component throughout the simulation. Its value must be specified as a XML attribute in every declared <Component> of this <ComponentType>. For example, different neurons may have different <Parameter> values, but the values stay the same for each cell. Attributes are name and dimension.

  • <Property>: This is same as the <Parameter> tag, except there can also be a defaultValue attribute that is used if the value is not specified in a declared <Component>. For example, the weight of synapses and input sources is 1 unless otherwise specified.

  • <Constant>: This is a quantity that’s just immutable. Its value is specific to the <ComponentType>; it may not be defined for each instance separately. They are typically used for natural constants such as \(\pi\), or units that are used in math expressions. But they can also be used to reduce the memory footprint of the component, so that more of them can fit in the computer during simulation. Attributes are name, dimension, and mandatory value.

• The <Dynamics> tag (or tags), that describes the actual behaviour of the <Component>. Its structure will be described in the following sub-section.

The components of <Dynamics>

These LEMS properties are used inside the <Dynamics> tag of a <ComponentType> to define its behaviour:

“Dynamic” variables

• <StateVariable>: This is a self-standing quantity inside the component that may change during the simulation. Its value must be assigned <OnStart>, see below.

  • Required attributes are name and dimension. Optionally the exposure attribute can be used to have it be a named <Exposure>.

• <DerivedVariable>: This is a dependent quantity whose value is a function of other quantities. It can’t depend directly on indirectly on itself, of course.

  • Required attributes are name and either:

    • value, which is the mathematical expression that yields the <DerivedVariable>’s value;

    • or select, which may be the sum of the <Exposure>s on the mechanism’s <Children>. This will be explained further through EDEN’s Multiflux extension.

• <ConditionalDerivedVariable>: This is a <DerivedVariable> which may be evaluate as different expressions, depending on whether certain mathematical conditions apply.

  Note: It is notably used for ion channel gate equations when the formula can evaluate as 0/0 for certain Vm. In these cases, the formula is fixed by instead taking the limit and linearisation, or just visually substituting with the curve’s value, near the problematic point.

  • Required attributes are name and dimension.

  • It also must have one or more <Case> tags, that each specify as attributes the applicable condition and the associated math expression yielding the value.

  • The last case should not have a condition; it is the default <Case>, and its value will apply when all other <Case>s don’t.

  • Note: The character < is not allowed in XML attributes (though EDEN’s parser may still take it). Consider instead of <, >, <=, >=, the alternative FORTRAN spellings .lt., .gt., .le., .ge..

  • The conditions are evaluated in the provided order: in case multiple conditions apply, the first stated <Case> takes precedence. But the default <Case> comes last (and should be stated last). The algorithm goes like: if <first case> ... else if <second case> ... else <default case>.

  • Due to how the original LEMS interpreter was made, absent a default <Case>, a ConditionalVariable slipping all cases will be 0 but please don’t rely on this!

• <DerivedParameter>: This is a <DerivedVariable> which depends only on fixed <Parameter>s, <Property>s and <Constants>, and hence its value remains the same for each component.

  • Required attributes are name, dimension and value.

Dynamic behaviour

• <TimeDerivative>: This is the rate of change over time that applies for a <StateVariable>.

  • Required attributes are the variable the rate is for, and the value expression for the rate (its dimension should be (the variable's dimension)/time).

  • It can be specified once per <StateVariable>. If you want the rate to be Conditional, then make a <ConditionalDerivedVariable> for the rate and use that as the value. Note: It is not necessary for every StateVariable to have a rate; some variables may be updated only in specific occasions through <OnCondition> or <OnEvent> tags.

• <OnStart>: This tag contains the necessary actions to take when the <Component> first starts to exist. That is, assign an initial value for each of its <StateVariables>.

  • It contains one or more <StateAssignment> tags:

    • Required attributes for each are the variable and the value expression to be assigned to it.

  • Try to not depend on the sequence these are evaluated in.

• <OnCondition>: This tag contains the necessary actions to take when a mathematical condition is met. For example, when a LIF <iafCell>’s v reaches the firing threshold, then the cell emits a spike and v’s value is reset.

  • The necessary attribute is the condition to test. It may contain the following tags:

  • <StateAssignment>: Set a state variable to the value of a math expression, same as with <OnStart>.

  • <EventOut>: Emit a spike on the specified output port.

• <OnEvent>: This tag is similar to <OnCondition>, but instead of the test expression it waits on an <EventPort> (with direction in) named by the port attribute. When a spike arrives on that port, the enclosed actions are taken as with OnCondition. For example, when a spike reaches a synapse’s in port, its neurotransmitter cascade is started or stimulated further.

How can LEMS components be used in NeuroML models?

As mentioned at the start, NeuroML has some models for all types of mechanisms, but that’s hardly ever enough. For our modelling, we may prefer a slightly different model; it could be a model from a recent paper, or a brand new proposed set of dynamics equations. We can make such models in NeuroML by describing them in the LEMS dynamics language.

In practice, this means describing individual mechanisms as LEMS components, and using these components as options for model parts in our NeuroML model; the NeuroML formulation remains in charge of the high-level structure of the model.

This section will present the needed LEMS context, and the method to model in LEMS each type of mechanism that EDEN supports.

For even more <Dynamics> options and other capabilities to model custom interactions, refer to the “EDEN extensions” section of this user guide.

Overview of LEMS mechanisms

EDEN allows LEMS components for all these types of mechanisms:

• “Point” neurons without detailed anatomy

• Synapses, on their pre- and post-synaptic parts

  • AP based, graded, or mixed interaction

• Input sources of the experimental rig such as current, voltage, and other clamps and electrodes

• Ion channels and their properties:

  • Custom whole-gate dynamics

    • Or gate transition rates (expressed as \(\alpha\), \(\beta\) rates or, \(\tau, \infty\) time to stabilise and steady state) for HH, Markov and instantaneous gates

  • Custom scaling factor for ion channel conductance

  • Entire ion channels with fully-custom dynamics

• Ion pool models

Let’s see how to work for each one of them.

Commonly available quantities

Beside the quantities that are provided for each type, there are some common, intrinsically exposed quantities that are always available:

• time: The virtual time that the simulation has been running for, starting at 0. Also aliased as t.

• temperature: It’s the same across the simulated <network>.

Neuron-local quantities

For mechanisms that are located on a neuron, its biophysical properties nearby are also available:

• v: Membrane voltage

• caConc: Calcium ion concentration (named ca) for the mechanisms that care about it. Note: There is also a way to specify an alternate Ca²⁺ species ca2 which does not affect mechanisms.

Custom point neurons

Base <ComponentType>s for this type to extend include <𝚋𝚊𝚜𝚎𝙲𝚎𝚕𝚕>, <𝚋𝚊𝚜𝚎𝚂𝚙𝚒𝚔𝚒𝚗𝚐𝙲𝚎𝚕𝚕>, <𝚋𝚊𝚜𝚎𝙲𝚎𝚕𝚕𝙼𝚎𝚖𝚋𝙿𝚘𝚝>, <𝚋𝚊𝚜𝚎𝙲𝚎𝚕𝚕𝙼𝚎𝚖𝚋𝙿𝚘𝚝𝙲𝚊𝚙>, <𝚋𝚊𝚜𝚎𝙲𝚎𝚕𝚕𝙼𝚎𝚖𝚋𝙿𝚘𝚝𝙳𝙻>, <𝚋𝚊𝚜𝚎𝙸𝚊𝚏>, <𝚋𝚊𝚜𝚎𝙸𝚊𝚏𝙲𝚊𝚙𝙲𝚎𝚕𝚕>, <𝚋𝚊𝚜𝚎𝙲𝚎𝚕𝚕>, and <𝚋𝚊𝚜𝚎𝙿𝚢𝙽𝙽𝙲𝚎𝚕𝚕>, <𝚋𝚊𝚜𝚎𝙿𝚢𝙽𝙽𝙸𝚊𝙵𝙲𝚎𝚕𝚕>, <𝚋𝚊𝚜𝚎𝙿𝚢𝙽𝙽𝙸𝚊𝙵𝙲𝚘𝚗𝚍𝙲𝚎𝚕𝚕>.

The LEMS interface for point neurons is as follows:

Available quantities to have as<Requirement>s:

• iSyn: The total inward current flow, from attached synapses and (unlike its name would suggest) input sources.

• ISyn: Similar, but used for dimensionless current flowing into cells with dimensionless state variables (e.g. FitzHugh-Nagumo or the original Izhikevich formulation)

Optional<Exposure>:

• v: The membrane voltage of the cell, needed by synapses and certain input sources.

Optional<EventPort>:

• spike (direction: out): The one and only (for now) spike output of cells.

Note that neurons does not receive events directly hence they don’t have an inward <EventPort>; incoming spikes are mediated by synapses instead.

If the model requires that spikes change the internal state of a neuron directly, see EDEN’s <WritableRequirement> extension.

Example: The Quadratic Integrate-and-Fire neuron model

As we saw in the NeuroML tutorial, NeuroML includes in its standard component library linear integrate-and fire neurons (as it does AdEx neurons and some classic models).

But we may prefer a higher-quality model for our experiments. Here we’ll see how to make a custom neuron type, using the Quadratic Integrate-and-Fire model as example.

Its dynamics equations are:

\[\begin{split}\begin{gathered} C_m~\frac{dV_m}{dt} = k~(V_m - V_{\rm rest})~(V_m - V_{\rm crit}) + i_{\rm syn}\\ \\ V_m > V​_{\rm peak}~~⇒~~V_m ~\colon= ​​​​V_{\rm reset}​​ \end{gathered}\end{split}\]

And this is how it’s written in LEMS:

%%writefile Custom_Cell_Model.nml
<neuroml>

<!-- Let's make a new ComponentType ❗ -->
<ComponentType name="QifCell"
    extends="baseCellMembPotCap"
    description="Integrate-and-fire cell with quadratic Vm dynamics">

    <!-- Reminder: The following are inherited from "baseCellMembPotCap" -->
    <!-- <Exposure name="v" dimension="voltage" description="Membrane potential"/> -->
    <!-- <EventPort name="spike" direction="out" description="Spike event"/> -->
    <!-- <Parameter name="C" dimension="capacitance" description="Total capacitance of the cell membrane"/> -->
    <!-- <Exposure name="iMemb" dimension="current" description="Total current crossing the cell membrane"/> -->
    <!-- And iSyn is provided and added as requirement by EDEN -->
    <!-- <Requirement name="iSyn" dimension="current" description="Total current due to synaptic inputs"/> -->

    <!-- Now let's add our own stuff -->
    <Parameter name="v_rest" dimension="voltage" description="Resting potential"/>
    <Parameter name="v_crit" dimension="voltage" description="Critical potential for initiating spike"/>
    <Parameter name="v_peak" dimension="voltage" description="Threshold potential to register the spike and reset"/>
    <Parameter name="v_reset" dimension="voltage" description="Reset (repolarisation) potential"/>
    <Parameter name="v2_factor" dimension="conductance_per_voltage" description="Intrinsic rate scaling factor"/>

    <Dynamics>

        <!-- It is a one-dimensional model, it has one state variable -->
        <StateVariable name="v" dimension="voltage" exposure="v"/>

        <!-- Extract iMemb from the rate equation to clear things up,
            and satisfy the inherited Exposure as well -->
        <DerivedVariable name="iMemb" dimension="current" exposure="iMemb"
          value="v2_factor * (v-v_rest) * (v-v_crit) + iSyn"/>

        <!-- Here comes the rate, it's reasonable -->
        <TimeDerivative variable="v" value="iMemb / C"/>

        <OnStart> <!-- Start at resting state -->
            <StateAssignment variable="v" value="v_rest"/>
        </OnStart>

        <OnCondition test="v > v_peak"> <!-- Send spike and reset -->
            <StateAssignment variable="v" value="v_reset"/>
            <EventOut port="spike"/>
        </OnCondition>

    </Dynamics>
</ComponentType>

<!-- Here comes the concrete cell type with parameters ❗ -->
<QifCell id="MyFirstQif" C="200 pF" v2_factor="0.7 nS_per_mV"
    v_rest="-60 mV" v_crit="-30 mV" v_peak="+30 mV" v_reset="-70mV" />

<!-- Add a current clamp to rouse the cell -->
<pulseGenerator id="A_DC_Clamp" delay="100ms" duration="500ms" amplitude="0.21nA"/>
<!-- Here comes the network -->
<network id="Net" type="networkWithTemperature" temperature="37degC" >
    <!-- Add a population with a single cell -->
    <population id="Pop" component="MyFirstQif" size="1"/>
    <!-- Apply a stim on said cell -->
    <inputList id="Inps" population="Pop" component="A_DC_Clamp">
        <input id="0" target="Pop[0]" destination="synapses"/></inputList>
</network>
</neuroml>

Writing Custom_Cell_Model.nml

And add a <Simulation> file, don’t forget to record the cell’s vm:

%%writefile LEMS_Custom_Cell_Sim.xml
<Lems>
<include file="Custom_Cell_Model.nml" /> <!-- Remember to check the filename -->
<Simulation id="Sim" length="0.7 s" step="0.1 ms" target="Net" >
    <OutputFile id="MyFirstOutputFile" fileName="results.gen.txt">
        <OutputColumn id="vm" quantity="Pop[0]/v"/>
    </OutputFile>
</Simulation>
<Target component="Sim"/>
</Lems>

Writing LEMS_Custom_Cell_Sim.xml

Look at the action potentials that this cell generates:

from eden_simulator import runEden; from matplotlib import pyplot as plt
results = runEden('LEMS_Custom_Cell_Sim.xml')
plt.plot(results['t'],results['Pop[0]/v'])
plt.xlabel('Time (sec)'); plt.ylabel('Membrane voltage (V)');

/home/docs/checkouts/readthedocs.org/user_builds/eden-simulator/envs/latest/lib/python3.10/site-packages/eden_simulator/embeden.py:3: UserWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html. The pkg_resources package is slated for removal as early as 2025-11-30. Refrain from using this package or pin to Setuptools<81.
  import pkg_resources

Custom synapses

Note that every synapse components is attached to a particular neuron, either on the pre- or post-synaptic side. Hence bi-directional synapses are made up of two components.

Base <ComponentType>s for this type to extend include <𝚋𝚊𝚜𝚎𝚂𝚢𝚗𝚊𝚙𝚜𝚎>, <𝚋𝚊𝚜𝚎𝚅𝚘𝚕𝚝𝚊𝚐𝚎𝙳𝚎𝚙𝚂𝚢𝚗𝚊𝚙𝚜𝚎>, <𝚋𝚊𝚜𝚎𝙲𝚞𝚛𝚛𝚎𝚗𝚝𝙱𝚊𝚜𝚎𝚍𝚂𝚢𝚗𝚊𝚙𝚜𝚎>, <𝚋𝚊𝚜𝚎𝙲𝚘𝚗𝚍𝚞𝚌𝚝𝚊𝚗𝚌𝚎𝙱𝚊𝚜𝚎𝚍𝚂𝚢𝚗𝚊𝚙𝚜𝚎>, <𝚋𝚊𝚜𝚎𝙲𝚘𝚗𝚍𝚞𝚌𝚝𝚊𝚗𝚌𝚎𝙱𝚊𝚜𝚎𝚍𝚂𝚢𝚗𝚊𝚙𝚜𝚎𝚃𝚠𝚘>, <𝚋𝚊𝚜𝚎𝙿𝚢𝚗𝚗𝚂𝚢𝚗𝚊𝚙𝚜𝚎>, and <𝚋𝚊𝚜𝚎𝙶𝚛𝚊𝚍𝚎𝚍𝚂𝚢𝚗𝚊𝚙𝚜𝚎>, <𝚐𝚊𝚙𝙹𝚞𝚗𝚌𝚝𝚒𝚘𝚗>.

The LEMS interface for synapses is as follows:

Available quantities to have as<Requirement>s:

• All the intra-cell quantities, plus:

• vpeer (voltage): The membrane potential of the non-local cell (other side of the synaptic connection). Available only if exposed by the peer cell.

Required<Exposure>s:

• i: The electrical current that this synapse injects into the cell. Must be of the same dimension as the cell’s iSyn. May also be capital I as per the convention for dimensionless current.

Optional<Exposure>s:

• g: The electrical conductance of the mechanism, basically \(\frac{d i}{d V_m}\). May help with simulation accuracy when available.

• Other flows that a synapse may inject to a cell using the Multiflux extension.

Optional<EventPort>:

• in (direction: in): The incoming spikes that (re)activate the synapse.

The following examples show how to make a custom chemical synapse, gap junction and more.

Example: A time-varying chemical synapse

It is known that sometimes, the effect of chemical synapses is modulated by the concurrent presence of other chemical signals. For example, the spike-induced increment in neurotransmitter release could be more or less depending on external conditions.

Here we’ll show an exponential-decay post-synaptic current synapse, whose weight is a function of elapsed simulatedtime. The generated current decays by a time constant of \(\tau_{\rm syn}\), and it is incremented by \(I_{\rm base} \cdot sin(k t)​​\) on every spike.

Therefore, its dynamics equations are:

\[\begin{split}\begin{gathered} \frac{dI_{\rm syn}}{dt} = -\frac{I_{\rm syn}}{\tau_{\rm syn}}\\ \\ spike~in ~~\Rightarrow~~ I_{\rm syn} ~\colon= I_{\rm syn} + I_{\rm base} ~ sin( 2 \pi {time \over T_{weight}})​​ \end{gathered}\end{split}\]

And this is how it’s written in LEMS:

%%writefile Custom_ChemSyn_Model.nml
<neuroml>

<!-- Let's make a new ComponentType ❗ -->
<ComponentType name="SinSyn"
    extends="baseCurrentBasedSynapse"
    description="Time-varying exponential current synapse">

    <!-- Reminder: The following are inherited from "baseCurrentBasedSynapse" -->
    <!-- <Exposure name="i" dimension="current" description="Generated current"/> -->
    <!-- <EventPort name="in" direction="in" description="Synapse trigger"/> -->

    <!-- Now let's add our own stuff -->
    <Parameter name="i_base" dimension="current" description="Base peak-current level (mefore modulation)"/>
    <Parameter name="tau_syn" dimension="time" description="Time constant of exponential decay"/>
    <Parameter name="T_weight" dimension="time" description="Period of time-varying weight"/>

    <!-- Note: we may need some unit constants because physical values are not directly allowed in LEMS expressions -->
    <Constant name="pA" dimension="current" value="1 pA"/>
    <Constant name="pi" dimension="none" value="3.14159"/> <!-- Add the dimensionless constant π -->

    <!-- Add dependence on "time" ❗ -->
    <Requirement name="time" dimension="time" description="Time elapsed during the simulation"/>

    <Dynamics>

        <!-- It is a one-dimensional model, it has one state variable -->
        <StateVariable name="i" dimension="current" exposure="i"/>
        <!-- And a rate of change -->
        <TimeDerivative variable="i" value="- i / tau_syn"/>

        <OnStart> <!-- Start at idle, no generated current -->
            <StateAssignment variable="i" value="0 * pA"/>
        </OnStart>

        <OnEvent port="in"> <!-- Increase current flow -->
            <StateAssignment variable="i" value="i + i_base * sin(2*pi*time/T_weight)"/>
        </OnEvent>

    </Dynamics>
</ComponentType>

<!-- Here comes the concrete synapse type with parameters ❗ -->
<SinSyn id="MyCustomChemSynapse" i_base="3 nA" tau_syn="1 msec" T_weight="450 msec" />

<!-- Add a simple cell, and spike source to use our new synapse with -->
<iafCell id="MyFirstCellType" C="200 pF" leakConductance="10 nS" leakReversal="-70 mV" reset="-70mV" thresh="-50mV" />
<spikeGenerator id="spikeGenRegular" period="70 ms"/>
<!-- Here comes the network -->
<network id="Net" type="networkWithTemperature" temperature="37degC" >
    <!-- Add a single cell and spike source -->
    <population id="Pop" component="MyFirstCellType" size="1"/>
    <population id="Spi" component="spikeGenRegular" size="1" />
    <!-- Connect them with a synapse -->
    <projection id="Pro" presynapticPopulation="Spi" postsynapticPopulation="Pop" synapse="MyCustomChemSynapse">
        <connection id="0" preCellId="0" postCellId="0"/> </projection>
</network>
</neuroml>

Writing Custom_ChemSyn_Model.nml

And let’s see what happens:

%%writefile LEMS_Custom_ChemSyn_Sim.xml
<Lems>
<include file="Custom_ChemSyn_Model.nml" /> <!-- Remember to check the filename -->
<Simulation id="Sim" length="1 s" step="0.1 ms" target="Net" >
    <OutputFile id="MyFirstOutputFile" fileName="results.gen.txt">
        <OutputColumn id="vm" quantity="Pop[0]/v"/>
    </OutputFile>
</Simulation>
<Target component="Sim"/>
</Lems>

Writing LEMS_Custom_ChemSyn_Sim.xml

from eden_simulator import runEden; from matplotlib import pyplot as plt
results = runEden('LEMS_Custom_ChemSyn_Sim.xml')
plt.plot(results['t'],results['Pop[0]/v'])
plt.xlabel('Time (sec)'); plt.ylabel('Membrane voltage (V)');

The same technique may also be used for factors other than time, and for other mechanism types as well.

Example: A rectifying gap junction

It is also known that the electrical gap junctions coupling cells aren’t always ohmic; sometimes the magnitude of current flow depends on polarity! Importantly, this means that the current flow function is also not symmetric.

Here we’ll model a non-linear gap junction just like that. We’ll roughly follow the current vs. voltage curve of a diode with internal resistance, based on the HHExpLinearRate formula.

The conductance equation will then be of the form: \begin{gather*} I = I_0~ (\frac {x} { 1 - e^{x} } - 1) ~, ~~~ x = {{V} \over {V_0}} \end{gather*}

for \(I\) the current from the anode to the cathode, \(V\) the voltage between cathode and anode respectively.

And this is how it’s written in LEMS:

%%writefile Custom_GradSyn_Model.nml
<neuroml>

<!-- Let's make a new ComponentType ❗ -->
<ComponentType name="GapDiode_Cathode"
    extends="baseGradedSynapse"
    description="Non-symmetrical gap junction, cathode half">

    <!-- Reminder: The following are inherited from "baseGradedSynapse" -->
    <!-- <Exposure name="i" dimension="current" description="Generated current"/> -->
    <!-- And vpeer is provided and added as requirement by EDEN -->
    <!-- <Requirement name="vpeer" dimension="voltage" description="Potential of peer cell"/> -->

    <!-- Now let's add our own stuff -->
    <Parameter name="I_rev_max" dimension="current" description="Asymptotic max reverse current"/>
    <Parameter name="forward_conductance" dimension="conductance" description="Marginal (ie small-signal) conductance when fwd polarised"/>

    <Dynamics>
        <!-- It is a *static* model, it has no state variables ❗ -->
        <!-- Use a derived variable to simplify the formula -->
        <DerivedVariable name="x" dimension="none" value="(vpeer - v)/(I_rev_max/forward_conductance)"/>

        <ConditionalDerivedVariable name="i" exposure="i" dimension="current">

            <!-- Special case for near the singularity, linearise ❗ -->
            <Case condition="abs(x) < 0.0001" value="x * I_rev_max / 2"/>

            <Case  value="I_rev_max * ( (x / (1 - exp(-x))) - 1 )"/> <!-- Usual case, but beware of exp(x) = 1 ❗ -->
        </ConditionalDerivedVariable>

    </Dynamics>
</ComponentType>
<!-- And one more for the other side ❗ -->
<ComponentType name="GapDiode_Anode" extends="baseGradedSynapse" description="Non-symmetrical gap junction, anode half">
    <Parameter name="I_rev_max" dimension="current" description="Asymptotic max reverse current"/>
    <Parameter name="forward_conductance" dimension="conductance" description="Marginal (ie small-signal) conductance when fwd polarised"/>
    <Dynamics>

        <!-- Now things get interesting. FLIP x to sample the same asymmetric curve,
            and then FLIP the result for the opposite current into this side ❗  -->
        <DerivedVariable name="x" dimension="none" value="(v - vpeer)/(I_rev_max/forward_conductance)"/>
        <ConditionalDerivedVariable name="i" exposure="i" dimension="current">=
            <Case condition="abs(x) < 0.05" value=" - I_rev_max * x"/>
            <Case  value=" - I_rev_max * ( (x / (1 - exp(-x))) - 1 )"/> <!-- Usual case -->
        </ConditionalDerivedVariable>

    </Dynamics>
</ComponentType>

<!-- Here comes the concrete synapse type with parameters ❗ -->
<GapDiode_Cathode id="MyCathode" forward_conductance="10 nS" I_rev_max="0.001 nA"/>
<GapDiode_Anode   id="MyAnode"   forward_conductance="10 nS" I_rev_max="0.001 nA"/>

<!-- Add a simple cell type, and a pair of pulse sources -->
<iafCell id="MyFirstCellType" C="200 pF" leakConductance="10 nS" leakReversal="-70 mV" reset="-70mV" thresh="-50mV" />
<pulseGenerator id="Inp0" delay="150ms" duration="100ms" amplitude="0.1nA"/>
<pulseGenerator id="Inp1" delay="750ms" duration="100ms" amplitude="0.1nA"/>
<!-- Here comes the network -->
<network id="Net" type="networkWithTemperature" temperature="37degC" >
    <!-- Add a pair of cells -->
    <population id="Pop" component="MyFirstCellType" size="2"/>
    <!-- Connect them with a synapse -->
    <continuousProjection id="Pro" presynapticPopulation="Pop" postsynapticPopulation="Pop">
        <continuousConnection id="0" preCell="0" postCell="1" preComponent="MyAnode" postComponent="MyCathode"/>
    </continuousProjection>
    <!-- And stimulate them -->
    <inputList id="Inpli0" population="Pop" component="Inp0"><input id="0" target="0" destination="synapses"/></inputList>
    <inputList id="Inpli1" population="Pop" component="Inp1"><input id="0" target="1" destination="synapses"/></inputList>
</network>
</neuroml>

Writing Custom_GradSyn_Model.nml

And let’s see what happens:

%%writefile LEMS_Custom_GradSyn_Sim.xml
<Lems>
<include file="Custom_GradSyn_Model.nml" /> <!-- Remember to check the filename -->
<Simulation id="Sim" length="1 s" step="0.1 ms" target="Net" >
    <OutputFile id="MyFirstOutputFile" fileName="results.gen.txt">
        <OutputColumn id="vm" quantity="Pop[0]/v"/>
        <OutputColumn id="v2" quantity="Pop[1]/v"/>
    </OutputFile>
</Simulation>
<Target component="Sim"/>
</Lems>

Writing LEMS_Custom_GradSyn_Sim.xml

from eden_simulator import runEden; from matplotlib import pyplot as plt
results = runEden('LEMS_Custom_GradSyn_Sim.xml')
plt.plot(results['t'],results['Pop[0]/v'], label='My first cell')
plt.plot(results['t'],results['Pop[1]/v'], label='My second cell')
plt.xlabel('Time (sec)'); plt.ylabel('Membrane voltage (V)'); plt.legend();

When current is injected to the first cell, a significant part flows to the second cell; but when the second cell is stimulated, barely any flows “back” to the first cell. Observe that the potential peak is in the former case is less than in the latter, since the current is divided among both cells.

Note that a single ComponentType could also be used here: the \(V-I\) curve can be flipped by means of an additional “polarity” Parameter. And in fact, in this model the I_rev_max parameter happens to affect x linearly in both places where a flip is needed: flip the sign on MyCathode and see what happens vs. MyAnode !

More synapse examples

For further LEMS examples of on complicated synapse models, check out an auxiliary spike counter implemented as a synapse, and highly-retailed models of multi-path plasticity (Ebner et al. 2019), and different interacting receptor-transmitter pairs (Schwartz et al. 2012).

Custom input sources

There are two types of input sources: neuron-attached probes, and self-standing spike generators.

Base <ComponentType>s for this type to extend include for the former <𝚋𝚊𝚜𝚎𝙿𝚘𝚒𝚗𝚝𝙲𝚞𝚛𝚛𝚎𝚗𝚝>, <𝚋𝚊𝚜𝚎𝚅𝚘𝚕𝚝𝚊𝚐𝚎𝙳𝚎𝚙𝙿𝚘𝚒𝚗𝚝𝙲𝚞𝚛𝚛𝚎𝚗𝚝>, <𝚋𝚊𝚜𝚎𝚅𝚘𝚕𝚝𝚊𝚐𝚎𝙳𝚎𝚙𝙿𝚘𝚒𝚗𝚝𝙲𝚞𝚛𝚛𝚎𝚗𝚝𝚂𝚙𝚒𝚔𝚒𝚗𝚐>, <𝚋𝚊𝚜𝚎𝙿𝚘𝚒𝚗𝚝𝙲𝚞𝚛𝚛𝚎𝚗𝚝𝙳𝙻>, <𝚋𝚊𝚜𝚎𝚅𝚘𝚕𝚝𝚊𝚐𝚎𝙳𝚎𝚙𝙿𝚘𝚒𝚗𝚝𝙲𝚞𝚛𝚛𝚎𝚗𝚝𝙳𝙻>, and for the latter <𝚋𝚊𝚜𝚎𝚂𝚙𝚒𝚔𝚎𝚂𝚘𝚞𝚛𝚌𝚎>.

Note that if some neuron-attached input sources are implied to emit spikes in the LEMS definitions, they actually don’t: the existence of ports is just a technicality of the LEMS underpinnings, outside the NeuroML context. Therefore <baseVoltageDepPointCurrentSpiking> is equivalent to <baseVoltageDepPointCurrent>.

The LEMS interface for probes is as follows:

Available quantities to have as<Requirement>s:

• All the intra-cell quantities, plus:

Suggested<Exposure>s:

• i: The current that this probe injects into the cell.

• g: The electrical conductance of the mechanism, basically \(\frac{d i}{d V_m}\). May help with simulation accuracy when available.

• I: Similar, but used for dimensionless current flowing into cells with dimensionless state variable.)

• Other flows that a probe may inject to a cell using the Multiflux extension.

For spike sources, requirements and exposure are as for point neurons; refer to the same section on making custom spike sources.

Example: Ornstein-Uhlenbeck noise

An popular electrical signal to stimulate cells with is Ornstein-Uhlenbeck noise which is thought to better approximate neural activity than, say, white noise. This is a stochastic process that evolves as a random walk. A very interesting feature of random walks is that due to their fractal nature, the change of its state non-linear on the simulation’s timestep. Hence a <TimeDerivative> won’t cut it, regardless of (first order accurate or better) ODE solver. We need to update the value with an explicit formula on every timestep!

Besides showing how to make a custom stimulus probe, this model also shows two commonly requested tricks:

• Random numbers following the normal distribution instead of the uniform that random() provides;

• Discrete updates that happen on every timestep. Note: This is outside the LEMS specification and does not necessarily behave the same on all simulators! But this is a proposed feature for LEMS and it may get codified in the future.

%%writefile OU_Input_Model.nml
<neuroml>
<!-- Let's make a new ComponentType ❗ -->
<ComponentType name="Noise_OU" extends="basePointCurrent">
    <Parameter name="mean" dimension="current" description="Average current to revert to"/>
    <Parameter name="sigma" dimension="current" description="Standard deviation of current"/>
    <Parameter name="tau" dimension="time" description="Correlation time constant"/>
    <Parameter name="delta_t" dimension="time" description="Sampling period for noise"/> <!-- Should be exactly  a <Simulation> step ❗ -->
    <Dynamics>
        <StateVariable name="i" exposure="i" dimension="current"/>
        <DerivedVariable name="randn" exposure="randn" dimension="none" value="sqrt(-2*ln(random(1)))*cos(2*3.14159265359*random(1))" description="A normally random variable for every moment ❗"/>
        <OnCondition test="1+1==2"> <StateAssignment variable="i" value="i + (mean - i)*delta_t/tau + sigma*randn*sqrt(2*delta_t/tau)"/> </OnCondition>
        <OnStart> <StateAssignment variable="i" value="mean"/> <!-- or to taste --> </OnStart>
    </Dynamics>
</ComponentType>

<!-- Add a simple cell type, and this input source -->
<iafCell id="MyFirstCellType" C="200 pF" leakConductance="10 nS" leakReversal="-70 mV" reset="-70mV" thresh="-50mV" />
<Noise_OU id="MyNoise" mean="0.12 nA" sigma="0.02 nA" tau="300 msec" delta_t="0.1 ms"/>

<network id="Net" type="networkWithTemperature" temperature="37degC" >
    <!-- Add a cell -->
    <population id="Pop" component="MyFirstCellType" size="1"/>
    <!-- And stimulate it -->
    <inputList id="Inpli0" population="Pop" component="MyNoise"><input id="0" target="0" destination="synapses"/></inputList>
</network>
</neuroml>

Writing OU_Input_Model.nml

%%writefile LEMS_OU_Input_Sim.xml
<Lems>
<include file="OU_Input_Model.nml" /> <!-- Remember to check the filename -->
<Simulation id="Sim" length="1 s" step="0.1 ms" target="Net" seed="20190109" > <!-- seed used -->
    <OutputFile id="MyFirstOutputFile" fileName="results.gen.txt">
        <OutputColumn id="vm" quantity="Pop[0]/v"/>
        <OutputColumn id="I" quantity="Inpli0[0]/i"/>
    </OutputFile>
</Simulation>
<Target component="Sim"/>
</Lems>

Writing LEMS_OU_Input_Sim.xml

And let’s see what happens.

Note: We’ll use extended LEMS paths to access the source’s current, so we can show it side by side with voltage.

from eden_simulator import runEden; from matplotlib import pyplot as plt
results = runEden('LEMS_OU_Input_Sim.xml')
def halflines(ax,x, y, color, label):
    ax.set_ylabel(label, color=color)
    ax.plot(x,y, color=color)
    ax.tick_params(axis='y', labelcolor=color)
col1 = 'tab:blue'
halflines(plt.gca(), results['t'],results['Pop[0]/v']*1000, color='tab:blue', label='Membrane voltage (mV)')
ax2 = plt.gca().twinx()
halflines(plt.gca(), results['t'],results['Inpli0[0]/i']*1e9, color='tab:orange', label='OU noise input (nA)')
plt.xlabel('Time (sec)');

Observe the delay between swings of input current and current. The membrane becomes less polarised than normally, due to the “positive” input current.

Set mean input to 0 or increase sigma and see what happens to the plots: autoscale can be misleading !

For an arbitrary sampling rate (and more commonly agreed behaviour), it just takes keeping track of the time of last update. Refer to the example from the OpenSourceBrain Stochasticity Showcase.

Custom ion channels, gates and rates

In biophysically modelled neurons, NeuroML offers a quite rich variety of modelling options. The summary is that ion channels are made up of gates, which open and close according to transition rates. Each of these parts can be modelled by a component from the core NeuroML library, or by a component described in LEMS. Details and examples on how to model each part follow:

Custom gate transition rates

The most common customisation is to just select different \(\alpha, ~\beta\) or \(\tau, ~x_\infty\) rates for the gate variables or ion channels, under either the Hodgkin-Huxley or Markov formulation.

Base <ComponentType>s for this type to extend include <𝚋𝚊𝚜𝚎𝚅𝚘𝚕𝚝𝚊𝚐𝚎𝙳𝚎𝚙𝚁𝚊𝚝𝚎>, <𝚋𝚊𝚜𝚎𝚅𝚘𝚕𝚝𝚊𝚐𝚎𝙲𝚘𝚗𝚌𝙳𝚎𝚙𝚁𝚊𝚝𝚎>, <𝚋𝚊𝚜𝚎𝚅𝚘𝚕𝚝𝚊𝚐𝚎𝙳𝚎𝚙𝚃𝚒𝚖𝚎>, <𝚋𝚊𝚜𝚎𝚅𝚘𝚕𝚝𝚊𝚐𝚎𝙲𝚘𝚗𝚌𝙳𝚎𝚙𝚃𝚒𝚖𝚎>, <𝚋𝚊𝚜𝚎𝚅𝚘𝚕𝚝𝚊𝚐𝚎𝙳𝚎𝚙𝚅𝚊𝚛𝚒𝚊𝚋𝚕𝚎>, and <𝚋𝚊𝚜𝚎𝚅𝚘𝚕𝚝𝚊𝚐𝚎𝙲𝚘𝚗𝚌𝙳𝚎𝚙𝚅𝚊𝚛𝚒𝚊𝚋𝚕𝚎>. Different <Requirement>s are inherited, depending on the base <componentType>.

The LEMS interface for gate rates is as follows:

Available quantities to have as<Requirement>s:

• All the intra-cell quantities.

Required<Exposure>s: One of either, depending on formula type:

• r (dimension: per_time): for forwardRate and reverseRate

• t (dimension: time): for timeCourse

• x (dimension: none): for steadyState

Example: Custom HH rates for a calcium channel

An example from the NeuroML-DB on how to write your own \(\alpha, ~\beta\) rates for your gate variables can be seen here (source).

Example: Custom timeCourse for a sodium channel

An example using both a core NeuroML steadyState and a custom timeCourse can be seen here (source).

Example: A Ca2⁺-dependent potassium channel

An example using calcium-dependent dynamics can be seen here (source).

Both timeCourse and steadyState are customised into components that extend <𝚋𝚊𝚜𝚎𝚅𝚘𝚕𝚝𝚊𝚐𝚎𝙲𝚘𝚗𝚌𝙳𝚎𝚙𝚃𝚒𝚖𝚎> and <𝚋𝚊𝚜𝚎𝚅𝚘𝚕𝚝𝚊𝚐𝚎𝙲𝚘𝚗𝚌𝙳𝚎𝚙𝚅𝚊𝚛𝚒𝚊𝚋𝚕𝚎> and expose t and x respectively (though a <𝚏𝚒𝚡𝚎𝚍𝚃𝚒𝚖𝚎𝙲𝚘𝚞𝚛𝚜𝚎> could have been used in this case).

Note that in all cases, v and caConc is normalised into dimensionless quantities(of implicit milliVolt and kiloMolar units), so that they can be used in f̵u̵n̵k̵y̵empirical formulae which wouldn’t make sense with dimensional quantities.

Custom gates

One may prefer a fully-custom model for ion channel gates, when customising its rates is not enough. The <componentType> to extend is <𝚋𝚊𝚜𝚎𝙶𝚊𝚝𝚎>.

The LEMS interface for gates is as follows:

Available quantities to have as<Requirement>s:

• All the intra-cell quantities.

Required<Exposure>s:

• q (dimension: none): The ‘state’ of the gate population, by convention the fraction of gates that are open.

• fcond (dimension: none): The activation of the gate (its multiplicative factor on channel conductance). May be, for example, q to the power of the gate’s multiplicity. Think of it as “fraction of conductivity”.

Custom conductance scaling for ion channels

More than ion channel gates, the conductance of an ion channel may be affected by other biophysical factors. Along with the temperature-dependent <q10ConductanceScaling> one may also specify custom LEMS components, that may also depend on other factors. Base <ComponentType>s for this type to extend include <𝚋𝚊𝚜𝚎𝙲𝚘𝚗𝚍𝚞𝚌𝚝𝚊𝚗𝚌𝚎𝚂𝚌𝚊𝚕𝚒𝚗𝚐> and <𝚋𝚊𝚜𝚎𝙲𝚘𝚗𝚍𝚞𝚌𝚝𝚊𝚗𝚌𝚎𝚂𝚌𝚊𝚕𝚒𝚗𝚐𝙲𝚊𝙳𝚎𝚙𝚎𝚗𝚍𝚎𝚗𝚝> .

The LEMS interface for conductance scaling is as follows:

Available quantities to have as<Requirement>s:

• All the intra-cell quantities.

Required<Exposure>s:

• factor (dimension: none): The number that this scaling effect multiplies the ion channel’s base conductance by.

Example: Ca2⁺-dependent base conductance

An example of calcium-dependent baseConductance for an ion channel can be seen here (source).

Note that in this case it is inserted in the ion channel as a <baseConductanceScalingCaDependent type="cond_scaling_kc_fast"/>, the way EDEN detects it is through the type attribute of the tag (or the tag name that should match the ComponentType if type), and the fact that it extends baseConductanceScalingCaDependent.

Custom ion channel models

Finally, one may prefer a fully-custom model for the whole ion channel. The <componentType> to extend is <𝚋𝚊𝚜𝚎𝙸𝚘𝚗𝙲𝚑𝚊𝚗𝚗𝚎𝚕>.

The LEMS interface for gates is as follows:

Available quantities to have as<Requirement>s:

• All the intra-cell quantities.

Required<Exposure>s:

• fopen (dimension: none): The fraction of ion channels that are open, after accounting for conductance scaling.

• g (dimension: conductance): The electrical conductance of the ion channel, in absolute terms.

Custom ion concentration models

Ion concentration models determine the dynamical relation between local ion concentration and the ion fluxes that feed it. For custom dynamics, the base <ComponentType>to extend is <𝚌𝚘𝚗𝚌𝚎𝚗𝚝𝚛𝚊𝚝𝚒𝚘𝚗𝙼𝚘𝚍𝚎𝚕>.

The LEMS interface for concentration models is as follows:

Available quantities to have as<Requirement>s:

• All the intra-cell quantities, plus:

• surfaceArea: The total surface area of the compartment’s membrane.

• initialConcentration: The initial value for concentration, as specified with <species> tags inside <intracellularProperties>.

• initialExtConcentration: The same for the outer side of the membrane.

• iCa, iCa2: the total per-compartment inflow in terms of current for the two independent calcium ion species that NeuroML supports. Perhaps the other ion flows should also be exposed like this? What about species without a charge?

Required<Exposure>s:

• concentration: The current value of internal concentration.

• extConcentration: Likewise for the outer side.

Example: An extended concentration model from Hay et al. 2011

An alternative concentration model from Hay et al. 2011 can be seen here (source). It is better than the similar <𝚍𝚎𝚌𝚊𝚢𝚒𝚗𝚐𝙿𝚘𝚘𝚕𝙲𝚘𝚗𝚌𝚎𝚗𝚝𝚛𝚊𝚝𝚒𝚘𝚗𝙼𝚘𝚍𝚎𝚕> from core NeuroML, in that the latter assumes the compartment to be a sphere and calculates the near-membrane volume accordingly. Hence the core NeuroML model is better suited for neurons modelled as spheres, whereas the <concentrationModelHayEtAl> model is better suited to anatomically detailed cells.

Note: Do take into account what this model does before using it, and customise to fit your own model!

Example: A concentration model for an independent calcium pool

Some models have not one, but two separate species of calcium ion that have different effects and dynamics (see also <cell2CaPools> and friends.) A concentration model for the second species can be specified as in here (source).

---

This article listed the sorts of ways for EDEN users to specify whole new mechanisms in NeuroML.

Refer also to the official NeuroML user guide:

• On how NeuroML and LEMS intersect

• More on making new components

• The core NeuroML component type library

• The user guide for pure LEMS (also beyond what NeuroML covers)

Further than the NeuroML specification, EDEN also offers some extensions that give even more powers to LEMS components.