Simulating a multi compartment OLM neuron

In this section we will model and simulate a multi-compartment Oriens-lacunosum moleculare (OLM) interneuron cell from the rodent hippocampal CA1 network model developed by Bezaire et al. ([BRB+16]). The complete network model can be seen here on GitHub, and here on Open Source Brain.

Membrane potential for neuron recorded from the simulation at the soma

Fig. 15 Membrane potential of the simulated OLM cell at the soma.

This plot, saved as olm_example_sim_seg0_soma0-v.png is generated using the following Python NeuroML script:

#!/usr/bin/env python3
"""
Multi-compartmental OLM cell example

File: olm-example.py

Copyright 2021 NeuroML contributors
Authors: Padraig Gleeson, Ankur Sinha
"""

from neuroml import (NeuroMLDocument, IncludeType, Population, PulseGenerator, ExplicitInput, Network, SegmentGroup, Member, Property, Include, Instance, Location)
from CellBuilder import (create_cell, add_segment, add_channel_density, set_init_memb_potential, set_resistivity, set_specific_capacitance, get_seg_group_by_id)
from pyneuroml import pynml
from pyneuroml.lems import LEMSSimulation
import numpy as np


def main():
    """Main function

    Include the NeuroML model into a LEMS simulation file, run it, plot some
    data.
    """
    # Simulation bits
    sim_id = "olm_example_sim"
    simulation = LEMSSimulation(sim_id=sim_id, duration=600, dt=0.01, simulation_seed=123)
    # Include the NeuroML model file
    simulation.include_neuroml2_file(create_olm_network())
    # Assign target for the simulation
    simulation.assign_simulation_target("single_olm_cell_network")

    # Recording information from the simulation
    simulation.create_output_file(id="output0", file_name=sim_id + ".dat")
    simulation.add_column_to_output_file("output0", column_id="pop0_0_v", quantity="pop0[0]/v")
    simulation.add_column_to_output_file("output0",
                                         column_id="pop0_0_v_Seg0_soma_0",
                                         quantity="pop0/0/olm/0/v")
    simulation.add_column_to_output_file("output0",
                                         column_id="pop0_0_v_Seg1_soma_0",
                                         quantity="pop0/0/olm/1/v")
    simulation.add_column_to_output_file("output0",
                                         column_id="pop0_0_v_Seg0_axon_0",
                                         quantity="pop0/0/olm/2/v")
    simulation.add_column_to_output_file("output0",
                                         column_id="pop0_0_v_Seg1_axon_0",
                                         quantity="pop0/0/olm/3/v")
    simulation.add_column_to_output_file("output0",
                                         column_id="pop0_0_v_Seg0_dend_0",
                                         quantity="pop0/0/olm/4/v")
    simulation.add_column_to_output_file("output0",
                                         column_id="pop0_0_v_Seg1_dend_0",
                                         quantity="pop0/0/olm/6/v")
    simulation.add_column_to_output_file("output0",
                                         column_id="pop0_0_v_Seg0_dend_1",
                                         quantity="pop0/0/olm/5/v")
    simulation.add_column_to_output_file("output0",
                                         column_id="pop0_0_v_Seg1_dend_1",
                                         quantity="pop0/0/olm/7/v")
    # Save LEMS simulation to file
    sim_file = simulation.save_to_file()

    # Run the simulation using the NEURON simulator
    pynml.run_lems_with_jneuroml_neuron(sim_file, max_memory="2G", nogui=True,
                                        plot=False, skip_run=False)
    # Plot the data
    plot_data(sim_id)


def plot_data(sim_id):
    """Plot the sim data.

    Load the data from the file and plot the graph for the membrane potential
    using the pynml generate_plot utility function.

    :sim_id: ID of simulaton

    """
    data_array = np.loadtxt(sim_id + ".dat")
    pynml.generate_plot([data_array[:, 0]], [data_array[:, 1]], "Membrane potential (soma seg 0)", show_plot_already=False, save_figure_to=sim_id + "_seg0_soma0-v.png", xaxis="time (s)", yaxis="membrane potential (V)")
    pynml.generate_plot([data_array[:, 0]], [data_array[:, 2]], "Membrane potential (soma seg 1)", show_plot_already=False, save_figure_to=sim_id + "_seg1_soma0-v.png", xaxis="time (s)", yaxis="membrane potential (V)")
    pynml.generate_plot([data_array[:, 0]], [data_array[:, 3]], "Membrane potential (axon seg 0)", show_plot_already=False, save_figure_to=sim_id + "_seg0_axon0-v.png", xaxis="time (s)", yaxis="membrane potential (V)")
    pynml.generate_plot([data_array[:, 0]], [data_array[:, 4]], "Membrane potential (axon seg 1)", show_plot_already=False, save_figure_to=sim_id + "_seg1_axon0-v.png", xaxis="time (s)", yaxis="membrane potential (V)")


def create_olm_network():
    """Create the network

    :returns: name of network nml file
    """
    net_doc = NeuroMLDocument(id="network",
                              notes="OLM cell network")
    net_doc_fn = "olm_example_net.nml"
    net_doc.includes.append(IncludeType(href=create_olm_cell()))
    # Create a population: convenient to create many cells of the same type
    pop = Population(id="pop0", notes="A population for our cell",
                     component="olm", size=1, type="populationList")
    pop.instances.append(Instance(id=1, location=Location(0., 0., 0.)))
    # Input
    pulsegen = PulseGenerator(id="pg_olm", notes="Simple pulse generator", delay="100ms", duration="100ms", amplitude="0.08nA")

    exp_input = ExplicitInput(target="pop0[0]", input="pg_olm")

    net = Network(id="single_olm_cell_network", note="A network with a single population")
    net_doc.pulse_generators.append(pulsegen)
    net.explicit_inputs.append(exp_input)
    net.populations.append(pop)
    net_doc.networks.append(net)

    pynml.write_neuroml2_file(nml2_doc=net_doc, nml2_file_name=net_doc_fn, validate=True)
    return net_doc_fn


def create_olm_cell():
    """Create the complete cell.

    :returns: cell object
    """
    nml_cell_doc = NeuroMLDocument(id="oml_cell")
    cell = create_cell("olm")
    nml_cell_file = cell.id + ".cell.nml"

    # Add two soma segments
    diam = 10.0
    soma_0 = add_segment(cell,
                         prox=[0.0, 0.0, 0.0, diam],
                         dist=[0.0, 10., 0.0, diam],
                         name="Seg0_soma_0",
                         group="soma_0")

    soma_1 = add_segment(cell,
                         prox=None,
                         dist=[0.0, 10. + 10., 0.0, diam],
                         name="Seg1_soma_0",
                         parent=soma_0,
                         group="soma_0")

    # Add axon segments
    diam = 1.5
    axon_0 = add_segment(cell,
                         prox=[0.0, 0.0, 0.0, diam],
                         dist=[0.0, -75, 0.0, diam],
                         name="Seg0_axon_0",
                         parent=soma_0,
                         fraction_along=0.0,
                         group="axon_0")
    axon_1 = add_segment(cell,
                         prox=None,
                         dist=[0.0, -150, 0.0, diam],
                         name="Seg1_axon_0",
                         parent=axon_0,
                         group="axon_0")

    # Add 2 dendrite segments

    diam = 3.0
    dend_0_0 = add_segment(cell,
                           prox=[0.0, 20, 0.0, diam],
                           dist=[100, 120, 0.0, diam],
                           name="Seg0_dend_0",
                           parent=soma_1,
                           fraction_along=1,
                           group="dend_0")

    dend_1_0 = add_segment(cell,
                           prox=None,
                           dist=[177, 197, 0.0, diam],
                           name="Seg1_dend_0",
                           parent=dend_0_0,
                           fraction_along=1,
                           group="dend_0")

    dend_0_1 = add_segment(cell,
                           prox=[0.0, 20, 0.0, diam],
                           dist=[-100, 120, 0.0, diam],
                           name="Seg0_dend_1",
                           parent=soma_1,
                           fraction_along=1,
                           group="dend_1")
    dend_1_1 = add_segment(cell,
                           prox=None,
                           dist=[-177, 197, 0.0, diam],
                           name="Seg1_dend_1",
                           parent=dend_0_1,
                           fraction_along=1,
                           group="dend_1")

    # XXX: For segment groups to be correctly mapped to sections in NEURON,
    # they must include the correct neurolex ID
    for section_name in ["soma_0", "axon_0", "dend_0", "dend_1"]:
        section_group = get_seg_group_by_id(section_name, cell)
        section_group.neuro_lex_id = 'sao864921383'

    den_seg_group = get_seg_group_by_id("dendrite_group", cell)
    den_seg_group.includes.append(Include(segment_groups="dend_0"))
    den_seg_group.includes.append(Include(segment_groups="dend_1"))
    den_seg_group.properties.append(Property(tag="color", value="0.8 0 0"))

    ax_seg_group = get_seg_group_by_id("axon_group", cell)
    ax_seg_group.includes.append(Include(segment_groups="axon_0"))
    ax_seg_group.properties.append(Property(tag="color", value="0 0.8 0"))

    soma_seg_group = get_seg_group_by_id("soma_group", cell)
    soma_seg_group.includes.append(Include(segment_groups="soma_0"))

    soma_seg_group.properties.append(Property(tag="color", value="0 0 0.8"))

    # Other cell properties
    set_init_memb_potential(cell, "-67mV")
    set_resistivity(cell, "0.15 kohm_cm")
    set_specific_capacitance(cell, "1.3 uF_per_cm2")

    # channels
    # leak
    add_channel_density(cell, nml_cell_doc,
                        cd_id="leak_all",
                        cond_density="0.01 mS_per_cm2",
                        ion_channel="leak_chan",
                        ion_chan_def_file="olm-example/leak_chan.channel.nml",
                        erev="-67mV",
                        ion="non_specific")
    # HCNolm_soma
    add_channel_density(cell, nml_cell_doc,
                        cd_id="HCNolm_soma",
                        cond_density="0.5 mS_per_cm2",
                        ion_channel="HCNolm",
                        ion_chan_def_file="olm-example/HCNolm.channel.nml",
                        erev="-32.9mV",
                        ion="h",
                        group="soma_group")
    # Kdrfast_soma
    add_channel_density(cell, nml_cell_doc,
                        cd_id="Kdrfast_soma",
                        cond_density="73.37 mS_per_cm2",
                        ion_channel="Kdrfast",
                        ion_chan_def_file="olm-example/Kdrfast.channel.nml",
                        erev="-77mV",
                        ion="k",
                        group="soma_group")
    # Kdrfast_dendrite
    add_channel_density(cell, nml_cell_doc,
                        cd_id="Kdrfast_dendrite",
                        cond_density="105.8 mS_per_cm2",
                        ion_channel="Kdrfast",
                        ion_chan_def_file="olm-example/Kdrfast.channel.nml",
                        erev="-77mV",
                        ion="k",
                        group="dendrite_group")
    # Kdrfast_axon
    add_channel_density(cell, nml_cell_doc,
                        cd_id="Kdrfast_axon",
                        cond_density="117.392 mS_per_cm2",
                        ion_channel="Kdrfast",
                        ion_chan_def_file="olm-example/Kdrfast.channel.nml",
                        erev="-77mV",
                        ion="k",
                        group="axon_group")
    # KvAolm_soma
    add_channel_density(cell, nml_cell_doc,
                        cd_id="KvAolm_soma",
                        cond_density="4.95 mS_per_cm2",
                        ion_channel="KvAolm",
                        ion_chan_def_file="olm-example/KvAolm.channel.nml",
                        erev="-77mV",
                        ion="k",
                        group="soma_group")
    # KvAolm_dendrite
    add_channel_density(cell, nml_cell_doc,
                        cd_id="KvAolm_dendrite",
                        cond_density="2.8 mS_per_cm2",
                        ion_channel="KvAolm",
                        ion_chan_def_file="olm-example/KvAolm.channel.nml",
                        erev="-77mV",
                        ion="k",
                        group="dendrite_group")
    # Nav_soma
    add_channel_density(cell, nml_cell_doc,
                        cd_id="Nav_soma",
                        cond_density="10.7 mS_per_cm2",
                        ion_channel="Nav",
                        ion_chan_def_file="olm-example/Nav.channel.nml",
                        erev="50mV",
                        ion="na",
                        group="soma_group")
    # Nav_dendrite
    add_channel_density(cell, nml_cell_doc,
                        cd_id="Nav_dendrite",
                        cond_density="23.4 mS_per_cm2",
                        ion_channel="Nav",
                        ion_chan_def_file="olm-example/Nav.channel.nml",
                        erev="50mV",
                        ion="na",
                        group="dendrite_group")
    # Nav_axon
    add_channel_density(cell, nml_cell_doc,
                        cd_id="Nav_axon",
                        cond_density="17.12 mS_per_cm2",
                        ion_channel="Nav",
                        ion_chan_def_file="olm-example/Nav.channel.nml",
                        erev="50mV",
                        ion="na",
                        group="axon_group")

    nml_cell_doc.cells.append(cell)
    pynml.write_neuroml2_file(nml_cell_doc, nml_cell_file, True, True)
    return nml_cell_file


if __name__ == "__main__":
    main()

As we will see, we repeat the same operations in the script while adding segments and ion-channels to our model, so we also write some helper functions to make it easier for ourselves:

#!/usr/bin/env python3
"""
Utility functions to help build cells in NeuroML

File: CellBuilder.py

Copyright 2021 NeuroML contributors
Author: Ankur Sinha <sanjay DOT ankur AT gmail DOT com>
"""


from typing import List
from neuroml import (Cell, Morphology, MembraneProperties, IntracellularProperties, BiophysicalProperties, Segment, SegmentGroup, Point3DWithDiam, SegmentParent, Member, InitMembPotential, Resistivity, SpecificCapacitance, NeuroMLDocument, IncludeType, ChannelDensity) # type: ignore  # noqa
from pyneuroml.pynml import print_function  # type: ignore


neuro_lex_ids = {
    'axon': "GO:0030424",
    'dend': "GO:0030425",
    'soma': "GO:0043025",
}


def create_cell(cell_id):
    # type: (str) -> Cell
    """Create a NeuroML Cell.

    Initialises the cell with these properties assigning IDs where applicable:
    - Morphology: "morphology"
    - BiophysicalProperties: "biophys"
    - MembraneProperties
    - IntracellularProperties
    - SegmentGroups: "all", "soma_group", "dendrite_group", "axon_group" which
      can be used to include all, soma, dendrite, and axon segments
      respectively.

    Note that since this cell does not currently include a segment in its
    morphology, it is *not* a valid NeuroML construct. Use the `add_segment`
    function to add segments. `add_segment` will also populate the default
    segment groups this creates.

    :param cell_id: id of the cell
    :type cell_id: str
    :returns: created cell object of type neuroml.Cell

    """
    cell = Cell(id=cell_id)
    cell.morphology = Morphology(id='morphology')
    membrane_properties = MembraneProperties()
    intracellular_properties = IntracellularProperties()

    cell.biophysical_properties = BiophysicalProperties(
        id="biophys", intracellular_properties=intracellular_properties,
        membrane_properties=membrane_properties)

    seg_group_all = SegmentGroup(id='all')
    seg_group_soma = SegmentGroup(id='soma_group',
                                  neuro_lex_id=neuro_lex_ids["soma"],
                                  notes="Default soma segment group for the cell")
    seg_group_axon = SegmentGroup(id='axon_group',
                                  neuro_lex_id=neuro_lex_ids["axon"],
                                  notes="Default axon segment group for the cell")
    seg_group_dend = SegmentGroup(id='dendrite_group',
                                  neuro_lex_id=neuro_lex_ids["dend"],
                                  notes="Default dendrite segment group for the cell")
    cell.morphology.segment_groups.append(seg_group_all)
    cell.morphology.segment_groups.append(seg_group_soma)
    cell.morphology.segment_groups.append(seg_group_axon)
    cell.morphology.segment_groups.append(seg_group_dend)

    return cell


def add_segment(cell, prox, dist, name=None, parent=None, fraction_along=1.0, group=None):
    # type: (Cell, List[float], List[float], str, SegmentParent, float, SegmentGroup) -> Segment
    """Add a segment to the cell.

    Convention: use axon_, soma_, dend_ prefixes for axon, soma, and dendrite
    segments respectivey. This will allow this function to add the correct
    neurolex IDs to the group.

    The created segment is also added to the default segment groups that were
    created by the `create_cell` function: "all", "dendrite_group",
    "soma_group", "axon_group" if the convention is followed.

    :param cell: cell to add segment to
    :type cell: Cell
    :param prox: proximal segment information
    :type prox: list with 4 float entries: [x, y, z, diameter]
    :param dist: dist segment information
    :type dist: list with 4 float entries: [x, y, z, diameter]
    :param name: name of segment
    :type name: str
    :param parent: parent segment
    :type parent: SegmentParent
    :param fraction_along: where the new segment is connected to the parent (0: distal point, 1: proximal point)
    :type fraction_along: float
    :param group: segment group to add the segment to
    :type group: SegmentGroup
    :returns: the created segment

    """
    try:
        if prox:
            p = Point3DWithDiam(x=prox[0], y=prox[1], z=prox[2], diameter=prox[3])
        else:
            p = None
    except IndexError as e:
        print_function("{}: prox must be a list of 4 elements".format(e))
    try:
        d = Point3DWithDiam(x=dist[0], y=dist[1], z=dist[2], diameter=dist[3])
    except IndexError as e:
        print_function("{}: dist must be a list of 4 elements".format(e))

    segid = len(cell.morphology.segments)

    sp = SegmentParent(segments=parent.id, fraction_along=fraction_along) if parent else None
    segment = Segment(id=segid, proximal=p, distal=d, parent=sp)

    if name:
        segment.name = name

    if group:
        seg_group = None
        seg_group = get_seg_group_by_id(group, cell)
        seg_group_all = get_seg_group_by_id("all", cell)
        seg_group_default = None
        neuro_lex_id = None

        if "axon_" in group:
            neuro_lex_id = neuro_lex_ids["axon"]  # See http://amigo.geneontology.org/amigo/term/GO:0030424
            seg_group_default = get_seg_group_by_id("axon_group", cell)
        if "soma_" in group:
            neuro_lex_id = neuro_lex_ids["soma"]
            seg_group_default = get_seg_group_by_id("soma_group", cell)
        if "dend_" in group:
            neuro_lex_id = neuro_lex_ids["dend"]
            seg_group_default = get_seg_group_by_id("dendrite_group", cell)

        if seg_group is None:
            seg_group = SegmentGroup(id=group, neuro_lex_id=neuro_lex_id)
            cell.morphology.segment_groups.append(seg_group)

        seg_group.members.append(Member(segments=segment.id))
        # Ideally, these higher level segment groups should just include other
        # segment groups using Include, which would result in smaller NML
        # files. However, because these default segment groups are defined
        # first, they are printed in the NML file before the new segments and their
        # groups. The jnml validator does not like this.
        # TODO: clarify if the order of definition is important, or if the jnml
        # validator needs to be updated to manage this use case.
        if seg_group_default:
            seg_group_default.members.append(Member(segments=segment.id))

    seg_group_all.members.append(Member(segments=segment.id))

    cell.morphology.segments.append(segment)
    return segment


def set_init_memb_potential(cell, v, group="all"):
    # type: (Cell, str, str) -> None
    """Set the initial membrane potential of the cell.

    :param cell: cell to modify
    :type cell: Cell
    :param v: value to set for membrane potential with units
    :type v: str
    :param group: id of segment group to modify
    :type group: str
    """
    cell.biophysical_properties.membrane_properties.init_memb_potentials = \
        [InitMembPotential(value=v, segment_groups=group)]


def set_resistivity(cell, resistivity, group="all"):
    # type: (Cell, str, str) -> None
    """Set the resistivity of the cell

    :param cell: cell to modfify
    :param resistivity: value resistivity to set with units
    :type resistivity: str
    :param group: segment group to modify
    :type group: str
    """
    cell.biophysical_properties.intracellular_properties.resistivities = [Resistivity(value=resistivity, segment_groups=group)]


def set_specific_capacitance(cell, spec_cap, group="all"):
    # type: (Cell, str, str) -> None
    """Set the specific capacitance for the cell.

    :param cell: cell to set specific capacitance for
    :type cell: Cell
    :param spec_cap: value of specific capacitance with units
    :type spec_cap: str
    :param group: segment group to modify
    :type group: str
    """
    cell.biophysical_properties.membrane_properties.specific_capacitances.append(SpecificCapacitance(value=spec_cap, segment_groups=group))


def add_channel_density(cell, nml_cell_doc, cd_id, cond_density, ion_channel, ion_chan_def_file="", erev="0.0 mV", ion="non_specific", group="all"):
    # type: (Cell, NeuroMLDocument, str, str, str, str, str, str, str) -> None
    """Add channel density.

    :param cell: cell to be modified
    :type cell: Cell
    :param nml_cell_doc: cell NeuroML document to which channel density is to be added
    :type nml_cell_doc: NeuroMLDocument
    :param cd_id: id for channel density
    :type cd_id: str
    :param cond_density: value of conductance density with units
    :type cond_density: str
    :param ion_channel: name of ion channel
    :type ion_channel: str
    :param ion_chan_def_file: path to NeuroML2 file defining the ion channel, if empty, it assumes the channel is defined in the same file
    :type ion_chan_def_file: str
    :param erev: value of reversal potential with units
    :type erev: str
    :param ion: name of ion
    :type ion: str
    :param group: segment groups to add to
    :type group: str
    """
    cd = ChannelDensity(id=cd_id,
                        segment_groups=group,
                        ion=ion,
                        ion_channel=ion_channel,
                        erev=erev,
                        cond_density=cond_density)

    cell.biophysical_properties.membrane_properties.channel_densities.append(cd)

    if len(ion_chan_def_file) > 0:
        if IncludeType(ion_chan_def_file) not in nml_cell_doc.includes:
            nml_cell_doc.includes.append(IncludeType(ion_chan_def_file))


def get_seg_group_by_id(sg_id, cell):
    # type (str, Cell) -> SegmentGroup
    """Return the SegmentGroup object for the specified segment group id.

    :param sg_id: id of segment group to find
    :type sg_id: str
    :param cell: cell to look for segment group in
    :type cell: Cell
    :returns: SegmentGroup object of specified ID or None if not found

    """
    if not sg_id or not cell:
        print_function("Please specify both a segment group id and a Cell")
        return None

    for sg in cell.morphology.segment_groups:
        if sg.id == sg_id:
            return sg

    return None


def get_seg_group_by_id_substring(sg_id_substring, cell):
    # type (str, Cell) -> list
    """Return segment groups that include the specified substring.

    :param sg_id_substring: substring to match
    :type sg_id_substring: str
    :param cell: cell to look for segment group in
    :type cell: Cell
    :returns: list of SegmentGroups whose IDs include the given substring

    """
    sg_group_list = []
    if not sg_id_substring or not cell:
        print_function("Please specify both a segment group id and a Cell")
        return None
    for sg in cell.morphology.segment_groups:
        if sg_id_substring in sg.id:
            sg_group_list.append(sg)
    return sg_group_list

These helper functions will be included in the Python NeuroML libraries in their next release. Currently, we import them into our Python script at the top:

from CellBuilder import (create_cell, add_segment, add_channel_density, set_init_memb_potential, set_resistivity, set_specific_capacitance, get_seg_group_by_id)

Declaring the model in NeuroML

Similar to previous examples, we will first declare the model, visualise it, and then simulate it. The OLM model is slightly more complex than the HH neuron model we had worked with in the previous tutorial since it includes multiple compartments. However, where we had declared the ion-channels ourselves in the previous example, here will will not do so. We will include channels that have been pre-defined in NeuroML to demonstrate how components defined in NeuroML can be easily re-used in models.

We will follow the same method as before. We will first define the cell, create a network with one instance of the cell, and then simulate it to record and plot the membrane potential from different segments.

Declaring the cell

To keep our Python script modularised, we start constructing our cell in a separate function.

def create_olm_cell():
    """Create the complete cell.

    :returns: cell object
    """
    nml_cell_doc = NeuroMLDocument(id="oml_cell")
    cell = create_cell("olm")
    nml_cell_file = cell.id + ".cell.nml"

    # Add two soma segments
    diam = 10.0
    soma_0 = add_segment(cell,
                         prox=[0.0, 0.0, 0.0, diam],
                         dist=[0.0, 10., 0.0, diam],
                         name="Seg0_soma_0",
                         group="soma_0")

    soma_1 = add_segment(cell,
                         prox=None,
                         dist=[0.0, 10. + 10., 0.0, diam],
                         name="Seg1_soma_0",
                         parent=soma_0,
                         group="soma_0")

    # Add axon segments
    diam = 1.5
    axon_0 = add_segment(cell,
                         prox=[0.0, 0.0, 0.0, diam],
                         dist=[0.0, -75, 0.0, diam],
                         name="Seg0_axon_0",
                         parent=soma_0,
                         fraction_along=0.0,
                         group="axon_0")
    axon_1 = add_segment(cell,
                         prox=None,
                         dist=[0.0, -150, 0.0, diam],
                         name="Seg1_axon_0",
                         parent=axon_0,
                         group="axon_0")

    # Add 2 dendrite segments

    diam = 3.0
    dend_0_0 = add_segment(cell,
                           prox=[0.0, 20, 0.0, diam],
                           dist=[100, 120, 0.0, diam],
                           name="Seg0_dend_0",
                           parent=soma_1,
                           fraction_along=1,
                           group="dend_0")

    dend_1_0 = add_segment(cell,
                           prox=None,
                           dist=[177, 197, 0.0, diam],
                           name="Seg1_dend_0",
                           parent=dend_0_0,
                           fraction_along=1,
                           group="dend_0")

    dend_0_1 = add_segment(cell,
                           prox=[0.0, 20, 0.0, diam],
                           dist=[-100, 120, 0.0, diam],
                           name="Seg0_dend_1",
                           parent=soma_1,
                           fraction_along=1,
                           group="dend_1")
    dend_1_1 = add_segment(cell,
                           prox=None,
                           dist=[-177, 197, 0.0, diam],
                           name="Seg1_dend_1",
                           parent=dend_0_1,
                           fraction_along=1,
                           group="dend_1")

    # XXX: For segment groups to be correctly mapped to sections in NEURON,
    # they must include the correct neurolex ID
    for section_name in ["soma_0", "axon_0", "dend_0", "dend_1"]:
        section_group = get_seg_group_by_id(section_name, cell)
        section_group.neuro_lex_id = 'sao864921383'

    den_seg_group = get_seg_group_by_id("dendrite_group", cell)
    den_seg_group.includes.append(Include(segment_groups="dend_0"))
    den_seg_group.includes.append(Include(segment_groups="dend_1"))
    den_seg_group.properties.append(Property(tag="color", value="0.8 0 0"))

    ax_seg_group = get_seg_group_by_id("axon_group", cell)
    ax_seg_group.includes.append(Include(segment_groups="axon_0"))
    ax_seg_group.properties.append(Property(tag="color", value="0 0.8 0"))

    soma_seg_group = get_seg_group_by_id("soma_group", cell)
    soma_seg_group.includes.append(Include(segment_groups="soma_0"))

    soma_seg_group.properties.append(Property(tag="color", value="0 0 0.8"))

    # Other cell properties
    set_init_memb_potential(cell, "-67mV")
    set_resistivity(cell, "0.15 kohm_cm")
    set_specific_capacitance(cell, "1.3 uF_per_cm2")

    # channels
    # leak
    add_channel_density(cell, nml_cell_doc,
                        cd_id="leak_all",
                        cond_density="0.01 mS_per_cm2",
                        ion_channel="leak_chan",
                        ion_chan_def_file="olm-example/leak_chan.channel.nml",
                        erev="-67mV",
                        ion="non_specific")
    # HCNolm_soma
    add_channel_density(cell, nml_cell_doc,
                        cd_id="HCNolm_soma",
                        cond_density="0.5 mS_per_cm2",
                        ion_channel="HCNolm",
                        ion_chan_def_file="olm-example/HCNolm.channel.nml",
                        erev="-32.9mV",
                        ion="h",
                        group="soma_group")
    # Kdrfast_soma
    add_channel_density(cell, nml_cell_doc,
                        cd_id="Kdrfast_soma",
                        cond_density="73.37 mS_per_cm2",
                        ion_channel="Kdrfast",
                        ion_chan_def_file="olm-example/Kdrfast.channel.nml",
                        erev="-77mV",
                        ion="k",
                        group="soma_group")
    # Kdrfast_dendrite
    add_channel_density(cell, nml_cell_doc,
                        cd_id="Kdrfast_dendrite",
                        cond_density="105.8 mS_per_cm2",
                        ion_channel="Kdrfast",
                        ion_chan_def_file="olm-example/Kdrfast.channel.nml",
                        erev="-77mV",
                        ion="k",
                        group="dendrite_group")
    # Kdrfast_axon
    add_channel_density(cell, nml_cell_doc,
                        cd_id="Kdrfast_axon",
                        cond_density="117.392 mS_per_cm2",
                        ion_channel="Kdrfast",
                        ion_chan_def_file="olm-example/Kdrfast.channel.nml",
                        erev="-77mV",
                        ion="k",
                        group="axon_group")
    # KvAolm_soma
    add_channel_density(cell, nml_cell_doc,
                        cd_id="KvAolm_soma",
                        cond_density="4.95 mS_per_cm2",
                        ion_channel="KvAolm",
                        ion_chan_def_file="olm-example/KvAolm.channel.nml",
                        erev="-77mV",
                        ion="k",
                        group="soma_group")
    # KvAolm_dendrite
    add_channel_density(cell, nml_cell_doc,
                        cd_id="KvAolm_dendrite",
                        cond_density="2.8 mS_per_cm2",
                        ion_channel="KvAolm",
                        ion_chan_def_file="olm-example/KvAolm.channel.nml",
                        erev="-77mV",
                        ion="k",
                        group="dendrite_group")
    # Nav_soma
    add_channel_density(cell, nml_cell_doc,
                        cd_id="Nav_soma",
                        cond_density="10.7 mS_per_cm2",
                        ion_channel="Nav",
                        ion_chan_def_file="olm-example/Nav.channel.nml",
                        erev="50mV",
                        ion="na",
                        group="soma_group")
    # Nav_dendrite
    add_channel_density(cell, nml_cell_doc,
                        cd_id="Nav_dendrite",
                        cond_density="23.4 mS_per_cm2",
                        ion_channel="Nav",
                        ion_chan_def_file="olm-example/Nav.channel.nml",
                        erev="50mV",
                        ion="na",
                        group="dendrite_group")
    # Nav_axon
    add_channel_density(cell, nml_cell_doc,
                        cd_id="Nav_axon",
                        cond_density="17.12 mS_per_cm2",
                        ion_channel="Nav",
                        ion_chan_def_file="olm-example/Nav.channel.nml",
                        erev="50mV",
                        ion="na",
                        group="axon_group")

    nml_cell_doc.cells.append(cell)
    pynml.write_neuroml2_file(nml_cell_doc, nml_cell_file, True, True)
    return nml_cell_file

Let us walk through this function:

    nml_cell_doc = NeuroMLDocument(id="oml_cell")
    cell = create_cell("olm")
    nml_cell_file = cell.id + ".cell.nml"

First, we create a new NeuroML document that we will use to save this cell. Then, we proceed to create a new cell using the Cell NeuroML component type, and define the name of the file we will use to store the cell in NeuroML. To create the cell, we use the create_cell utility function:

def create_cell(cell_id):
    cell = Cell(id=cell_id)
    cell.morphology = Morphology(id='morphology')
    membrane_properties = MembraneProperties()
    intracellular_properties = IntracellularProperties()

    cell.biophysical_properties = BiophysicalProperties(
        id="biophys", intracellular_properties=intracellular_properties,
        membrane_properties=membrane_properties)

    seg_group_all = SegmentGroup(id='all')
    seg_group_soma = SegmentGroup(id='soma_group',
                                  neuro_lex_id=neuro_lex_ids["soma"],
                                  notes="Default soma segment group for the cell")
    seg_group_axon = SegmentGroup(id='axon_group',
                                  neuro_lex_id=neuro_lex_ids["axon"],
                                  notes="Default axon segment group for the cell")
    seg_group_dend = SegmentGroup(id='dendrite_group',
                                  neuro_lex_id=neuro_lex_ids["dend"],
                                  notes="Default dendrite segment group for the cell")
    cell.morphology.segment_groups.append(seg_group_all)
    cell.morphology.segment_groups.append(seg_group_soma)
    cell.morphology.segment_groups.append(seg_group_axon)
    cell.morphology.segment_groups.append(seg_group_dend)

    return cell

We now know that a Cell component has two children: morphology, and biophysical properties. The function simply creates the new Cell component for us and initialises the morphology and biophysical properties. Additionally, it creates some default segment groups for us that we use to organise our our segments in later.

We now have an empty cell. Since we are building a multi-compartmental cell, we now proceed to define the detailed morphology of the cell. We do this by adding segments and grouping them in to segment groups which are both children elements of the morphology of the cell. This is done using the add_segment utility function as required:

    soma_0 = add_segment(cell,
                         prox=[0.0, 0.0, 0.0, diam],
                         dist=[0.0, 10., 0.0, diam],
                         name="Seg0_soma_0",
                         group="soma_0")

The utility function takes the dimensions of the segment—it’s proximal and distal co-ordinates and the diameter to create a segment of the provided name. Additionally, since segments need to be contiguous, it adds the segment to a parent segment. Finally, it places the segment into the specified segment group and the default groups that we also have and adds the segment to the cell’s morphology.

def add_segment(cell, prox, dist, name=None, parent=None, fraction_along=1.0, group=None):
    try:
        if prox:
            p = Point3DWithDiam(x=prox[0], y=prox[1], z=prox[2], diameter=prox[3])
        else:
            p = None
    except IndexError as e:
        print_function("{}: prox must be a list of 4 elements".format(e))
    try:
        d = Point3DWithDiam(x=dist[0], y=dist[1], z=dist[2], diameter=dist[3])
    except IndexError as e:
        print_function("{}: dist must be a list of 4 elements".format(e))

    segid = len(cell.morphology.segments)

    sp = SegmentParent(segments=parent.id, fraction_along=fraction_along) if parent else None
    segment = Segment(id=segid, proximal=p, distal=d, parent=sp)

    if name:
        segment.name = name

    if group:
        seg_group = None
        seg_group = get_seg_group_by_id(group, cell)
        seg_group_all = get_seg_group_by_id("all", cell)
        seg_group_default = None
        neuro_lex_id = None

        if "axon_" in group:
            neuro_lex_id = neuro_lex_ids["axon"]  # See http://amigo.geneontology.org/amigo/term/GO:0030424
            seg_group_default = get_seg_group_by_id("axon_group", cell)
        if "soma_" in group:
            neuro_lex_id = neuro_lex_ids["soma"]
            seg_group_default = get_seg_group_by_id("soma_group", cell)
        if "dend_" in group:
            neuro_lex_id = neuro_lex_ids["dend"]
            seg_group_default = get_seg_group_by_id("dendrite_group", cell)

        if seg_group is None:
            seg_group = SegmentGroup(id=group, neuro_lex_id=neuro_lex_id)
            cell.morphology.segment_groups.append(seg_group)

        seg_group.members.append(Member(segments=segment.id))
        if seg_group_default:
            seg_group_default.members.append(Member(segments=segment.id))

    seg_group_all.members.append(Member(segments=segment.id))

    cell.morphology.segments.append(segment)
    return segment

We call the same function multiple times to add soma, dendritic, and axonal segments to our cell. Note how the segments connect to each other to form the contiguous cell morphology.

    # Add two soma segments
    diam = 10.0
    soma_0 = add_segment(cell,
                         prox=[0.0, 0.0, 0.0, diam],
                         dist=[0.0, 10., 0.0, diam],
                         name="Seg0_soma_0",
                         group="soma_0")

    soma_1 = add_segment(cell,
                         prox=None,
                         dist=[0.0, 10. + 10., 0.0, diam],
                         name="Seg1_soma_0",
                         parent=soma_0,
                         group="soma_0")

    # Add axon segments
    diam = 1.5
    axon_0 = add_segment(cell,
                         prox=[0.0, 0.0, 0.0, diam],
                         dist=[0.0, -75, 0.0, diam],
                         name="Seg0_axon_0",
                         parent=soma_0,
                         fraction_along=0.0,
                         group="axon_0")
    axon_1 = add_segment(cell,
                         prox=None,
                         dist=[0.0, -150, 0.0, diam],
                         name="Seg1_axon_0",
                         parent=axon_0,
                         group="axon_0")

    # Add 2 dendrite segments

    diam = 3.0
    dend_0_0 = add_segment(cell,
                           prox=[0.0, 20, 0.0, diam],
                           dist=[100, 120, 0.0, diam],
                           name="Seg0_dend_0",
                           parent=soma_1,
                           fraction_along=1,
                           group="dend_0")

    dend_1_0 = add_segment(cell,
                           prox=None,
                           dist=[177, 197, 0.0, diam],
                           name="Seg1_dend_0",
                           parent=dend_0_0,
                           fraction_along=1,
                           group="dend_0")

    dend_0_1 = add_segment(cell,
                           prox=[0.0, 20, 0.0, diam],
                           dist=[-100, 120, 0.0, diam],
                           name="Seg0_dend_1",
                           parent=soma_1,
                           fraction_along=1,
                           group="dend_1")
    dend_1_1 = add_segment(cell,
                           prox=None,
                           dist=[-177, 197, 0.0, diam],
                           name="Seg1_dend_1",
                           parent=dend_0_1,
                           fraction_along=1,
                           group="dend_1")

Next, we add extra information to our segments and organise them so that they can be correctly exported to the NEURON format for simulation later.

    for section_name in ["soma_0", "axon_0", "dend_0", "dend_1"]:
        section_group = get_seg_group_by_id(section_name, cell)
        section_group.neuro_lex_id = 'sao864921383'

    den_seg_group = get_seg_group_by_id("dendrite_group", cell)
    den_seg_group.includes.append(Include(segment_groups="dend_0"))
    den_seg_group.includes.append(Include(segment_groups="dend_1"))
    den_seg_group.properties.append(Property(tag="color", value="0.8 0 0"))

    ax_seg_group = get_seg_group_by_id("axon_group", cell)
    ax_seg_group.includes.append(Include(segment_groups="axon_0"))
    ax_seg_group.properties.append(Property(tag="color", value="0 0.8 0"))

    soma_seg_group = get_seg_group_by_id("soma_group", cell)
    soma_seg_group.includes.append(Include(segment_groups="soma_0"))

    soma_seg_group.properties.append(Property(tag="color", value="0 0 0.8"))

We have now completed adding the morphological information to our cell. Next, we proceed to our biophysical properties, which are split into two:

We also use a few simple helper functions defined in the CellBuilder.py module to add these to our cell:

    # Other cell properties
    set_init_memb_potential(cell, "-67mV")
    set_resistivity(cell, "0.15 kohm_cm")
    set_specific_capacitance(cell, "1.3 uF_per_cm2")

    # channels
    # leak
    add_channel_density(cell, nml_cell_doc,
                        cd_id="leak_all",
                        cond_density="0.01 mS_per_cm2",
                        ion_channel="leak_chan",
                        ion_chan_def_file="olm-example/leak_chan.channel.nml",
                        erev="-67mV",
                        ion="non_specific")
    # HCNolm_soma
    add_channel_density(cell, nml_cell_doc,
                        cd_id="HCNolm_soma",
                        cond_density="0.5 mS_per_cm2",
                        ion_channel="HCNolm",
                        ion_chan_def_file="olm-example/HCNolm.channel.nml",
                        erev="-32.9mV",
                        ion="h",
                        group="soma_group")
    # Kdrfast_soma
    add_channel_density(cell, nml_cell_doc,
                        cd_id="Kdrfast_soma",
                        cond_density="73.37 mS_per_cm2",
                        ion_channel="Kdrfast",
                        ion_chan_def_file="olm-example/Kdrfast.channel.nml",
                        erev="-77mV",
                        ion="k",
                        group="soma_group")
    # Kdrfast_dendrite
    add_channel_density(cell, nml_cell_doc,
                        cd_id="Kdrfast_dendrite",
                        cond_density="105.8 mS_per_cm2",
                        ion_channel="Kdrfast",
                        ion_chan_def_file="olm-example/Kdrfast.channel.nml",
                        erev="-77mV",
                        ion="k",
                        group="dendrite_group")
    # Kdrfast_axon
    add_channel_density(cell, nml_cell_doc,
                        cd_id="Kdrfast_axon",
                        cond_density="117.392 mS_per_cm2",
                        ion_channel="Kdrfast",
                        ion_chan_def_file="olm-example/Kdrfast.channel.nml",
                        erev="-77mV",
                        ion="k",
                        group="axon_group")
    # KvAolm_soma
    add_channel_density(cell, nml_cell_doc,
                        cd_id="KvAolm_soma",
                        cond_density="4.95 mS_per_cm2",
                        ion_channel="KvAolm",
                        ion_chan_def_file="olm-example/KvAolm.channel.nml",
                        erev="-77mV",
                        ion="k",
                        group="soma_group")
    # KvAolm_dendrite
    add_channel_density(cell, nml_cell_doc,
                        cd_id="KvAolm_dendrite",
                        cond_density="2.8 mS_per_cm2",
                        ion_channel="KvAolm",
                        ion_chan_def_file="olm-example/KvAolm.channel.nml",
                        erev="-77mV",
                        ion="k",
                        group="dendrite_group")
    # Nav_soma
    add_channel_density(cell, nml_cell_doc,
                        cd_id="Nav_soma",
                        cond_density="10.7 mS_per_cm2",
                        ion_channel="Nav",
                        ion_chan_def_file="olm-example/Nav.channel.nml",
                        erev="50mV",
                        ion="na",
                        group="soma_group")
    # Nav_dendrite
    add_channel_density(cell, nml_cell_doc,
                        cd_id="Nav_dendrite",
                        cond_density="23.4 mS_per_cm2",
                        ion_channel="Nav",
                        ion_chan_def_file="olm-example/Nav.channel.nml",
                        erev="50mV",
                        ion="na",
                        group="dendrite_group")
    # Nav_axon
    add_channel_density(cell, nml_cell_doc,
                        cd_id="Nav_axon",
                        cond_density="17.12 mS_per_cm2",
                        ion_channel="Nav",
                        ion_chan_def_file="olm-example/Nav.channel.nml",
                        erev="50mV",
                        ion="na",
                        group="axon_group")

    nml_cell_doc.cells.append(cell)

This completes the definition of our cell. We write it to a NeuroML file, and validate it.

    pynml.write_neuroml2_file(nml_cell_doc, nml_cell_file, True, True)

The resulting NeuroML file is:

<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.2.xsd" id="oml_cell">
    <include href="olm-example/leak_chan.channel.nml"/>
    <include href="olm-example/HCNolm.channel.nml"/>
    <include href="olm-example/Kdrfast.channel.nml"/>
    <include href="olm-example/KvAolm.channel.nml"/>
    <include href="olm-example/Nav.channel.nml"/>
    <cell id="olm">
        <morphology id="morphology">
            <segment id="0" name="Seg0_soma_0">
                <proximal x="0.000000e+00" y="0.000000e+00" z="0.000000e+00" diameter="10.0"/>
                <distal x="0.000000e+00" y="1.000000e+01" z="0.000000e+00" diameter="10.0"/>
            </segment>
            <segment id="1" name="Seg1_soma_0">
                <parent segment="0"/>
                <distal x="0.000000e+00" y="2.000000e+01" z="0.000000e+00" diameter="10.0"/>
            </segment>
            <segment id="2" name="Seg0_axon_0">
                <parent segment="0" fractionAlong="0.0"/>
                <proximal x="0.000000e+00" y="0.000000e+00" z="0.000000e+00" diameter="1.5"/>
                <distal x="0.000000e+00" y="-7.500000e+01" z="0.000000e+00" diameter="1.5"/>
            </segment>
            <segment id="3" name="Seg1_axon_0">
                <parent segment="2"/>
                <distal x="0.000000e+00" y="-1.500000e+02" z="0.000000e+00" diameter="1.5"/>
            </segment>
            <segment id="4" name="Seg0_dend_0">
                <parent segment="1"/>
                <proximal x="0.000000e+00" y="2.000000e+01" z="0.000000e+00" diameter="3.0"/>
                <distal x="1.000000e+02" y="1.200000e+02" z="0.000000e+00" diameter="3.0"/>
            </segment>
            <segment id="5" name="Seg1_dend_0">
                <parent segment="4"/>
                <distal x="1.770000e+02" y="1.970000e+02" z="0.000000e+00" diameter="3.0"/>
            </segment>
            <segment id="6" name="Seg0_dend_1">
                <parent segment="1"/>
                <proximal x="0.000000e+00" y="2.000000e+01" z="0.000000e+00" diameter="3.0"/>
                <distal x="-1.000000e+02" y="1.200000e+02" z="0.000000e+00" diameter="3.0"/>
            </segment>
            <segment id="7" name="Seg1_dend_1">
                <parent segment="6"/>
                <distal x="-1.770000e+02" y="1.970000e+02" z="0.000000e+00" diameter="3.0"/>
            </segment>
            <segmentGroup id="all">
                <member segment="0"/>
                <member segment="1"/>
                <member segment="2"/>
                <member segment="3"/>
                <member segment="4"/>
                <member segment="5"/>
                <member segment="6"/>
                <member segment="7"/>
            </segmentGroup>
            <segmentGroup neuroLexId="GO:0043025" id="soma_group">
                <notes>Default soma segment group for the cell</notes>
                <property tag="color" value="0 0 0.8"/>
                <member segment="0"/>
                <member segment="1"/>
                <include segmentGroup="soma_0"/>
            </segmentGroup>
            <segmentGroup neuroLexId="GO:0030424" id="axon_group">
                <notes>Default axon segment group for the cell</notes>
                <property tag="color" value="0 0.8 0"/>
                <member segment="2"/>
                <member segment="3"/>
                <include segmentGroup="axon_0"/>
            </segmentGroup>
            <segmentGroup neuroLexId="GO:0030425" id="dendrite_group">
                <notes>Default dendrite segment group for the cell</notes>
                <property tag="color" value="0.8 0 0"/>
                <member segment="4"/>
                <member segment="5"/>
                <member segment="6"/>
                <member segment="7"/>
                <include segmentGroup="dend_0"/>
                <include segmentGroup="dend_1"/>
            </segmentGroup>
            <segmentGroup neuroLexId="sao864921383" id="soma_0">
                <member segment="0"/>
                <member segment="1"/>
            </segmentGroup>
            <segmentGroup neuroLexId="sao864921383" id="axon_0">
                <member segment="2"/>
                <member segment="3"/>
            </segmentGroup>
            <segmentGroup neuroLexId="sao864921383" id="dend_0">
                <member segment="4"/>
                <member segment="5"/>
            </segmentGroup>
            <segmentGroup neuroLexId="sao864921383" id="dend_1">
                <member segment="6"/>
                <member segment="7"/>
            </segmentGroup>
        </morphology>
        <biophysicalProperties id="biophys">
            <membraneProperties>
                <channelDensity id="leak_all" ionChannel="leak_chan" condDensity="0.01 mS_per_cm2" erev="-67mV" ion="non_specific"/>
                <channelDensity id="HCNolm_soma" ionChannel="HCNolm" condDensity="0.5 mS_per_cm2" erev="-32.9mV" segmentGroup="soma_group" ion="h"/>
                <channelDensity id="Kdrfast_soma" ionChannel="Kdrfast" condDensity="73.37 mS_per_cm2" erev="-77mV" segmentGroup="soma_group" ion="k"/>
                <channelDensity id="Kdrfast_dendrite" ionChannel="Kdrfast" condDensity="105.8 mS_per_cm2" erev="-77mV" segmentGroup="dendrite_group" ion="k"/>
                <channelDensity id="Kdrfast_axon" ionChannel="Kdrfast" condDensity="117.392 mS_per_cm2" erev="-77mV" segmentGroup="axon_group" ion="k"/>
                <channelDensity id="KvAolm_soma" ionChannel="KvAolm" condDensity="4.95 mS_per_cm2" erev="-77mV" segmentGroup="soma_group" ion="k"/>
                <channelDensity id="KvAolm_dendrite" ionChannel="KvAolm" condDensity="2.8 mS_per_cm2" erev="-77mV" segmentGroup="dendrite_group" ion="k"/>
                <channelDensity id="Nav_soma" ionChannel="Nav" condDensity="10.7 mS_per_cm2" erev="50mV" segmentGroup="soma_group" ion="na"/>
                <channelDensity id="Nav_dendrite" ionChannel="Nav" condDensity="23.4 mS_per_cm2" erev="50mV" segmentGroup="dendrite_group" ion="na"/>
                <channelDensity id="Nav_axon" ionChannel="Nav" condDensity="17.12 mS_per_cm2" erev="50mV" segmentGroup="axon_group" ion="na"/>
                <specificCapacitance value="1.3 uF_per_cm2"/>
                <initMembPotential value="-67mV"/>
            </membraneProperties>
            <intracellularProperties>
                <resistivity value="0.15 kohm_cm"/>
            </intracellularProperties>
        </biophysicalProperties>
    </cell>
</neuroml>

We can now already inspect our cell using the NeuroML tools:

pynml -png olm.cell.png
...
pyNeuroML >>> Writing to: /home/asinha/Documents/02_Code/00_mine/2020-OSB/NeuroML-Documentation/source/Userdocs/NML2_examples/olm.cell.png

This gives us a figure of the morphology of our cell:

Figure showing the morphology of the OLM cell generated from the NeuroML definition.

Fig. 16 Figure showing the morphology of the OLM cell generated from the NeuroML definition.

Declaring the network

We now use our cell in a network. Similar to our previous example, we are going to only create a network with one cell, and an explicit input to the cell:

def create_olm_network():
    """Create the network

    :returns: name of network nml file
    """
    net_doc = NeuroMLDocument(id="network",
                              notes="OLM cell network")
    net_doc_fn = "olm_example_net.nml"
    net_doc.includes.append(IncludeType(href=create_olm_cell()))
    # Create a population: convenient to create many cells of the same type
    pop = Population(id="pop0", notes="A population for our cell",
                     component="olm", size=1, type="populationList")
    pop.instances.append(Instance(id=1, location=Location(0., 0., 0.)))
    # Input
    pulsegen = PulseGenerator(id="pg_olm", notes="Simple pulse generator", delay="100ms", duration="100ms", amplitude="0.08nA")

    exp_input = ExplicitInput(target="pop0[0]", input="pg_olm")

    net = Network(id="single_olm_cell_network", note="A network with a single population")
    net_doc.pulse_generators.append(pulsegen)
    net.explicit_inputs.append(exp_input)
    net.populations.append(pop)
    net_doc.networks.append(net)

    pynml.write_neuroml2_file(nml2_doc=net_doc, nml2_file_name=net_doc_fn, validate=True)
    return net_doc_fn

We start in the same way, by creating a new NeuroML document and including our cell file into it. We then create a population comprising of a single cell. We create a pulse generator as an explicit input, which targets our population. Note that as the schema documentation for ExplicitInput notes, any current source (any component that extends basePointCurrent) can be used as an ExplicitInput.

We add all of these to the network and save (and validate) our network file. The NeuroML file generated is below:

<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.2.xsd" id="network">
    <notes>OLM cell network</notes>
    <include href="olm.cell.nml"/>
    <pulseGenerator id="pg_olm" delay="100ms" duration="100ms" amplitude="0.08nA">
        <notes>Simple pulse generator</notes>
    </pulseGenerator>
    <network id="single_olm_cell_network">
        <population id="pop0" component="olm" size="1" type="populationList">
            <notes>A population for our cell</notes>
            <instance id="1">
                <location x="0." y="0." z="0."/>
            </instance>
        </population>
        <explicitInput target="pop0[0]" input="pg_olm"/>
    </network>
</neuroml>

The generated NeuroML model

Before we look at simulating the model, we can inspect our model to check for correctness. All our NeuroML files were validated when they were created already, so we do not need to run this step again. However, if required, this can be easily done:

pynml -validate olm*nml

Next, we can visualise our model using the information noted in the visualising NeuroML models page (including the -v verbose option for more information on the cell):

pynml-summary olm_example_net.nml -v
*******************************************************
* NeuroMLDocument: network
*
*  ComponentType: ['Bezaire_HCNolm_tau', 'Bezaire_Kdrfast_betaq', 'Bezaire_KvAolm_taub', 'Bezaire_Nav_alphah']
*  IonChannel: ['HCNolm', 'Kdrfast', 'KvAolm', 'Nav', 'leak_chan']
*  PulseGenerator: ['pg_olm']
*
*  Cell: olm
*    <Segment|0|Seg0_soma_0>
*      Parent segment: None (root segment)
*      (0.0, 0.0, 0.0), diam 10.0um -> (0.0, 10.0, 0.0), diam 10.0um; seg length: 10.0 um
*      Surface area: 314.1592653589793 um2, volume: 785.3981633974482 um3
*    <Segment|1|Seg1_soma_0>
*      Parent segment: 0
*      None -> (0.0, 20.0, 0.0), diam 10.0um; seg length: 10.0 um
*      Surface area: 314.1592653589793 um2, volume: 785.3981633974482 um3
*    <Segment|2|Seg0_axon_0>
*      Parent segment: 0
*      (0.0, 0.0, 0.0), diam 1.5um -> (0.0, -75.0, 0.0), diam 1.5um; seg length: 75.0 um
*      Surface area: 353.4291735288517 um2, volume: 132.53594007331938 um3
*    <Segment|3|Seg1_axon_0>
*      Parent segment: 2
*      None -> (0.0, -150.0, 0.0), diam 1.5um; seg length: 75.0 um
*      Surface area: 353.4291735288517 um2, volume: 132.53594007331938 um3
*    <Segment|4|Seg0_dend_0>
*      Parent segment: 1
*      (0.0, 20.0, 0.0), diam 3.0um -> (100.0, 120.0, 0.0), diam 3.0um; seg length: 141.4213562373095 um
*      Surface area: 1332.8648814475098 um2, volume: 999.6486610856323 um3
*    <Segment|5|Seg1_dend_0>
*      Parent segment: 4
*      None -> (177.0, 197.0, 0.0), diam 3.0um; seg length: 108.89444430272832 um
*      Surface area: 1026.3059587145826 um2, volume: 769.7294690359369 um3
*    <Segment|6|Seg0_dend_1>
*      Parent segment: 1
*      (0.0, 20.0, 0.0), diam 3.0um -> (-100.0, 120.0, 0.0), diam 3.0um; seg length: 141.4213562373095 um
*      Surface area: 1332.8648814475098 um2, volume: 999.6486610856323 um3
*    <Segment|7|Seg1_dend_1>
*      Parent segment: 6
*      None -> (-177.0, 197.0, 0.0), diam 3.0um; seg length: 108.89444430272832 um
*      Surface area: 1026.3059587145826 um2, volume: 769.7294690359369 um3
*    Total length of 8 segments: 670.6316010800756 um; total area: 6053.518558099847 um2
*
*    SegmentGroup: all, 8 member(s),    0 included group(s);    contains 8 segments in total
*    SegmentGroup: soma_group,  2 member(s),    1 included group(s);   contains 2 segments in total
*    SegmentGroup: axon_group,  2 member(s),    1 included group(s);   contains 2 segments in total
*    SegmentGroup: dendrite_group,      4 member(s),    2 included group(s);    contains 4 segments in total
*    SegmentGroup: soma_0,      2 member(s),    0 included group(s);   contains 2 segments in total
*    SegmentGroup: axon_0,      2 member(s),    0 included group(s);   contains 2 segments in total
*    SegmentGroup: dend_0,      2 member(s),    0 included group(s);   contains 2 segments in total
*    SegmentGroup: dend_1,      2 member(s),    0 included group(s);   contains 2 segments in total
*
*    Channel density: leak_all on all;  conductance of 0.01 mS_per_cm2 through ion chan leak_chan with ion non_specific, erev: -67mV
*      Channel is on <Segment|0|Seg0_soma_0>,   total conductance: 0.1 S_per_m2 x 3.1415926535897934e-10 m2 = 3.1415926535897936e-11 S (31.41592653589794 pS)
*      Channel is on <Segment|1|Seg1_soma_0>,   total conductance: 0.1 S_per_m2 x 3.1415926535897934e-10 m2 = 3.1415926535897936e-11 S (31.41592653589794 pS)
*      Channel is on <Segment|2|Seg0_axon_0>,   total conductance: 0.1 S_per_m2 x 3.534291735288517e-10 m2 = 3.534291735288518e-11 S (35.34291735288518 pS)
*      Channel is on <Segment|3|Seg1_axon_0>,   total conductance: 0.1 S_per_m2 x 3.534291735288517e-10 m2 = 3.534291735288518e-11 S (35.34291735288518 pS)
*      Channel is on <Segment|4|Seg0_dend_0>,   total conductance: 0.1 S_per_m2 x 1.3328648814475097e-09 m2 = 1.3328648814475097e-10 S (133.28648814475096 pS)
*      Channel is on <Segment|5|Seg1_dend_0>,   total conductance: 0.1 S_per_m2 x 1.0263059587145826e-09 m2 = 1.0263059587145826e-10 S (102.63059587145825 pS)
*      Channel is on <Segment|6|Seg0_dend_1>,   total conductance: 0.1 S_per_m2 x 1.3328648814475097e-09 m2 = 1.3328648814475097e-10 S (133.28648814475096 pS)
*      Channel is on <Segment|7|Seg1_dend_1>,   total conductance: 0.1 S_per_m2 x 1.0263059587145826e-09 m2 = 1.0263059587145826e-10 S (102.63059587145825 pS)
*    Channel density: HCNolm_soma on soma_group;        conductance of 0.5 mS_per_cm2 through ion chan HCNolm with ion h, erev: -32.9mV
*      Channel is on <Segment|0|Seg0_soma_0>,   total conductance: 5.0 S_per_m2 x 3.1415926535897934e-10 m2 = 1.5707963267948968e-09 S (1570.796326794897 pS)
*      Channel is on <Segment|1|Seg1_soma_0>,   total conductance: 5.0 S_per_m2 x 3.1415926535897934e-10 m2 = 1.5707963267948968e-09 S (1570.796326794897 pS)
*    Channel density: Kdrfast_soma on soma_group;       conductance of 73.37 mS_per_cm2 through ion chan Kdrfast with ion k, erev: -77mV
*      Channel is on <Segment|0|Seg0_soma_0>,   total conductance: 733.7 S_per_m2 x 3.1415926535897934e-10 m2 = 2.3049865299388314e-07 S (230498.65299388315 pS)
*      Channel is on <Segment|1|Seg1_soma_0>,   total conductance: 733.7 S_per_m2 x 3.1415926535897934e-10 m2 = 2.3049865299388314e-07 S (230498.65299388315 pS)
*    Channel density: Kdrfast_dendrite on dendrite_group;       conductance of 105.8 mS_per_cm2 through ion chan Kdrfast with ion k, erev: -77mV
*      Channel is on <Segment|4|Seg0_dend_0>,   total conductance: 1058.0 S_per_m2 x 1.3328648814475097e-09 m2 = 1.4101710445714652e-06 S (1410171.0445714653 pS)
*      Channel is on <Segment|5|Seg1_dend_0>,   total conductance: 1058.0 S_per_m2 x 1.0263059587145826e-09 m2 = 1.0858317043200284e-06 S (1085831.7043200284 pS)
*      Channel is on <Segment|6|Seg0_dend_1>,   total conductance: 1058.0 S_per_m2 x 1.3328648814475097e-09 m2 = 1.4101710445714652e-06 S (1410171.0445714653 pS)
*      Channel is on <Segment|7|Seg1_dend_1>,   total conductance: 1058.0 S_per_m2 x 1.0263059587145826e-09 m2 = 1.0858317043200284e-06 S (1085831.7043200284 pS)
*    Channel density: Kdrfast_axon on axon_group;       conductance of 117.392 mS_per_cm2 through ion chan Kdrfast with ion k, erev: -77mV
*      Channel is on <Segment|2|Seg0_axon_0>,   total conductance: 1173.92 S_per_m2 x 3.534291735288517e-10 m2 = 4.1489757538898964e-07 S (414897.57538898964 pS)
*      Channel is on <Segment|3|Seg1_axon_0>,   total conductance: 1173.92 S_per_m2 x 3.534291735288517e-10 m2 = 4.1489757538898964e-07 S (414897.57538898964 pS)
*    Channel density: KvAolm_soma on soma_group;        conductance of 4.95 mS_per_cm2 through ion chan KvAolm with ion k, erev: -77mV
*      Channel is on <Segment|0|Seg0_soma_0>,   total conductance: 49.5 S_per_m2 x 3.1415926535897934e-10 m2 = 1.5550883635269477e-08 S (15550.883635269476 pS)
*      Channel is on <Segment|1|Seg1_soma_0>,   total conductance: 49.5 S_per_m2 x 3.1415926535897934e-10 m2 = 1.5550883635269477e-08 S (15550.883635269476 pS)
*    Channel density: KvAolm_dendrite on dendrite_group;        conductance of 2.8 mS_per_cm2 through ion chan KvAolm with ion k, erev: -77mV
*      Channel is on <Segment|4|Seg0_dend_0>,   total conductance: 28.0 S_per_m2 x 1.3328648814475097e-09 m2 = 3.7320216680530273e-08 S (37320.21668053028 pS)
*      Channel is on <Segment|5|Seg1_dend_0>,   total conductance: 28.0 S_per_m2 x 1.0263059587145826e-09 m2 = 2.8736566844008313e-08 S (28736.566844008314 pS)
*      Channel is on <Segment|6|Seg0_dend_1>,   total conductance: 28.0 S_per_m2 x 1.3328648814475097e-09 m2 = 3.7320216680530273e-08 S (37320.21668053028 pS)
*      Channel is on <Segment|7|Seg1_dend_1>,   total conductance: 28.0 S_per_m2 x 1.0263059587145826e-09 m2 = 2.8736566844008313e-08 S (28736.566844008314 pS)
*    Channel density: Nav_soma on soma_group;   conductance of 10.7 mS_per_cm2 through ion chan Nav with ion na, erev: 50mV
*      Channel is on <Segment|0|Seg0_soma_0>,   total conductance: 107.0 S_per_m2 x 3.1415926535897934e-10 m2 = 3.361504139341079e-08 S (33615.04139341079 pS)
*      Channel is on <Segment|1|Seg1_soma_0>,   total conductance: 107.0 S_per_m2 x 3.1415926535897934e-10 m2 = 3.361504139341079e-08 S (33615.04139341079 pS)
*    Channel density: Nav_dendrite on dendrite_group;   conductance of 23.4 mS_per_cm2 through ion chan Nav with ion na, erev: 50mV
*      Channel is on <Segment|4|Seg0_dend_0>,   total conductance: 234.0 S_per_m2 x 1.3328648814475097e-09 m2 = 3.118903822587173e-07 S (311890.3822587173 pS)
*      Channel is on <Segment|5|Seg1_dend_0>,   total conductance: 234.0 S_per_m2 x 1.0263059587145826e-09 m2 = 2.401555943392123e-07 S (240155.59433921232 pS)
*      Channel is on <Segment|6|Seg0_dend_1>,   total conductance: 234.0 S_per_m2 x 1.3328648814475097e-09 m2 = 3.118903822587173e-07 S (311890.3822587173 pS)
*      Channel is on <Segment|7|Seg1_dend_1>,   total conductance: 234.0 S_per_m2 x 1.0263059587145826e-09 m2 = 2.401555943392123e-07 S (240155.59433921232 pS)
*    Channel density: Nav_axon on axon_group;   conductance of 17.12 mS_per_cm2 through ion chan Nav with ion na, erev: 50mV
*      Channel is on <Segment|2|Seg0_axon_0>,   total conductance: 171.20000000000002 S_per_m2 x 3.534291735288517e-10 m2 = 6.050707450813942e-08 S (60507.07450813943 pS)
*      Channel is on <Segment|3|Seg1_axon_0>,   total conductance: 171.20000000000002 S_per_m2 x 3.534291735288517e-10 m2 = 6.050707450813942e-08 S (60507.07450813943 pS)
*
*    Specific capacitance on all: 1.3 uF_per_cm2
*      Capacitance of <Segment|0|Seg0_soma_0>,  total capacitance: 0.013000000000000001 F_per_m2 x 3.1415926535897934e-10 m2 = 4.084070449666732e-12 F (4.084070449666732 pF)
*      Capacitance of <Segment|1|Seg1_soma_0>,  total capacitance: 0.013000000000000001 F_per_m2 x 3.1415926535897934e-10 m2 = 4.084070449666732e-12 F (4.084070449666732 pF)
*      Capacitance of <Segment|2|Seg0_axon_0>,  total capacitance: 0.013000000000000001 F_per_m2 x 3.534291735288517e-10 m2 = 4.594579255875073e-12 F (4.594579255875073 pF)
*      Capacitance of <Segment|3|Seg1_axon_0>,  total capacitance: 0.013000000000000001 F_per_m2 x 3.534291735288517e-10 m2 = 4.594579255875073e-12 F (4.594579255875073 pF)
*      Capacitance of <Segment|4|Seg0_dend_0>,  total capacitance: 0.013000000000000001 F_per_m2 x 1.3328648814475097e-09 m2 = 1.732724345881763e-11 F (17.32724345881763 pF)
*      Capacitance of <Segment|5|Seg1_dend_0>,  total capacitance: 0.013000000000000001 F_per_m2 x 1.0263059587145826e-09 m2 = 1.3341977463289574e-11 F (13.341977463289574 pF)
*      Capacitance of <Segment|6|Seg0_dend_1>,  total capacitance: 0.013000000000000001 F_per_m2 x 1.3328648814475097e-09 m2 = 1.732724345881763e-11 F (17.32724345881763 pF)
*      Capacitance of <Segment|7|Seg1_dend_1>,  total capacitance: 0.013000000000000001 F_per_m2 x 1.0263059587145826e-09 m2 = 1.3341977463289574e-11 F (13.341977463289574 pF)
*
*  Network: single_olm_cell_network
*
*   1 cells in 1 populations
*     Population: pop0 with 1 components of type olm
*       Locations: [(0, 0, 0), ...]
*
*   0 connections in 0 projections
*
*   0 inputs in 0 input lists
*
*   1 explicit inputs (outside of input lists)
*     Explicit Input of type pg_olm to pop0(cell 0), destination: unspecified
*
*******************************************************

We can check the connectivity graph of the model:

pynml -graph 10 olm_example_net.nml

which will give us this figure:

Level 10 network graph generated by pynml

Fig. 17 Level 10 network graph generated by pynml

Simulating the model

Now that we have declared and inspected our network model and all its components, we can proceed to simulate it. We do this in the main function:

def main():
    """Main function

    Include the NeuroML model into a LEMS simulation file, run it, plot some
    data.
    """
    # Simulation bits
    sim_id = "olm_example_sim"
    simulation = LEMSSimulation(sim_id=sim_id, duration=600, dt=0.01, simulation_seed=123)
    # Include the NeuroML model file
    simulation.include_neuroml2_file(create_olm_network())
    # Assign target for the simulation
    simulation.assign_simulation_target("single_olm_cell_network")

    # Recording information from the simulation
    simulation.create_output_file(id="output0", file_name=sim_id + ".dat")
    simulation.add_column_to_output_file("output0", column_id="pop0_0_v", quantity="pop0[0]/v")
    simulation.add_column_to_output_file("output0",
                                         column_id="pop0_0_v_Seg0_soma_0",
                                         quantity="pop0/0/olm/0/v")
    simulation.add_column_to_output_file("output0",
                                         column_id="pop0_0_v_Seg1_soma_0",
                                         quantity="pop0/0/olm/1/v")
    simulation.add_column_to_output_file("output0",
                                         column_id="pop0_0_v_Seg0_axon_0",
                                         quantity="pop0/0/olm/2/v")
    simulation.add_column_to_output_file("output0",
                                         column_id="pop0_0_v_Seg1_axon_0",
                                         quantity="pop0/0/olm/3/v")
    simulation.add_column_to_output_file("output0",
                                         column_id="pop0_0_v_Seg0_dend_0",
                                         quantity="pop0/0/olm/4/v")
    simulation.add_column_to_output_file("output0",
                                         column_id="pop0_0_v_Seg1_dend_0",
                                         quantity="pop0/0/olm/6/v")
    simulation.add_column_to_output_file("output0",
                                         column_id="pop0_0_v_Seg0_dend_1",
                                         quantity="pop0/0/olm/5/v")
    simulation.add_column_to_output_file("output0",
                                         column_id="pop0_0_v_Seg1_dend_1",
                                         quantity="pop0/0/olm/7/v")
    # Save LEMS simulation to file
    sim_file = simulation.save_to_file()

    # Run the simulation using the NEURON simulator
    pynml.run_lems_with_jneuroml_neuron(sim_file, max_memory="2G", nogui=True,
                                        plot=False, skip_run=False)
    # Plot the data
    plot_data(sim_id)

Here we first create a LEMSSimulation instance and include our network NeuroML file in it. We must inform LEMS what the target of the simulation is. In our case, it’s the id of our network, single_olm_cell_network:

    sim_id = "olm_example_sim"
    simulation = LEMSSimulation(sim_id=sim_id, duration=600, dt=0.01, simulation_seed=123)
    # Include the NeuroML model file
    simulation.include_neuroml2_file(create_olm_network())
    # Assign target for the simulation
    simulation.assign_simulation_target("single_olm_cell_network")

We also want to record some information, so we create an output file first with an id of output0:

    simulation.create_output_file(id="output0", file_name=sim_id + ".dat")

Now, we can record any quantity that is exposed by NeuroML (any exposure). Here, for example, we add columns for the membrane potentials v of the different segments of the cell.

    simulation.add_column_to_output_file("output0",
                                         column_id="pop0_0_v_Seg0_soma_0",
                                         quantity="pop0/0/olm/0/v")
    simulation.add_column_to_output_file("output0",
                                         column_id="pop0_0_v_Seg1_soma_0",
                                         quantity="pop0/0/olm/1/v")
    simulation.add_column_to_output_file("output0",
                                         column_id="pop0_0_v_Seg0_axon_0",
                                         quantity="pop0/0/olm/2/v")
    simulation.add_column_to_output_file("output0",
                                         column_id="pop0_0_v_Seg1_axon_0",
                                         quantity="pop0/0/olm/3/v")
    simulation.add_column_to_output_file("output0",
                                         column_id="pop0_0_v_Seg0_dend_0",
                                         quantity="pop0/0/olm/4/v")
    simulation.add_column_to_output_file("output0",
                                         column_id="pop0_0_v_Seg1_dend_0",
                                         quantity="pop0/0/olm/6/v")
    simulation.add_column_to_output_file("output0",
                                         column_id="pop0_0_v_Seg0_dend_1",
                                         quantity="pop0/0/olm/5/v")
    simulation.add_column_to_output_file("output0",
                                         column_id="pop0_0_v_Seg1_dend_1",
                                         quantity="pop0/0/olm/7/v")

The path required to point to the quantity (exposure) to be recorded needs to be correctly provided. Here, where we use a population list that includes an instance of the cell, it is: population_id/instance_id/cell component type/segment id/exposure. (See tickets 15 and 16)

We then save the LEMS simulation file, and run our simulation with the NEURON simulator (since the default jNeuroML simulator can only simulate single compartment cells).

    # Run the simulation using the NEURON simulator
    pynml.run_lems_with_jneuroml_neuron(sim_file, max_memory="2G", nogui=True,
                                        plot=False, skip_run=False)

Plotting the recorded variables

To plot the variables that we recorded, we write a simple function that reads the data and uses the generate_plot utility function which generates the membrane potential graphs for different segments.

def plot_data(sim_id):
    """Plot the sim data.

    Load the data from the file and plot the graph for the membrane potential
    using the pynml generate_plot utility function.

    :sim_id: ID of simulaton

    """
    data_array = np.loadtxt(sim_id + ".dat")
    pynml.generate_plot([data_array[:, 0]], [data_array[:, 1]], "Membrane potential (soma seg 0)", show_plot_already=False, save_figure_to=sim_id + "_seg0_soma0-v.png", xaxis="time (s)", yaxis="membrane potential (V)")
    pynml.generate_plot([data_array[:, 0]], [data_array[:, 2]], "Membrane potential (soma seg 1)", show_plot_already=False, save_figure_to=sim_id + "_seg1_soma0-v.png", xaxis="time (s)", yaxis="membrane potential (V)")
    pynml.generate_plot([data_array[:, 0]], [data_array[:, 3]], "Membrane potential (axon seg 0)", show_plot_already=False, save_figure_to=sim_id + "_seg0_axon0-v.png", xaxis="time (s)", yaxis="membrane potential (V)")
    pynml.generate_plot([data_array[:, 0]], [data_array[:, 4]], "Membrane potential (axon seg 1)", show_plot_already=False, save_figure_to=sim_id + "_seg1_axon0-v.png", xaxis="time (s)", yaxis="membrane potential (V)")

This concludes this example. Here we have seen how to create, simulate, record, and visualise a multi-compartment neuron. In the next section, you will find an interactive notebook where you can play with this example.