Getting started

Reading model descriptions from XML files

NineML documents can contain abstraction layer models, user layer models (with references to abstraction layer models defined in other documents) or both.

To read a file containing only abstraction layer elements:

>>> import nineml, pprint
>>> doc = nineml.read("./BrunelIaF.xml")
>>> pprint(doc.items())
[('BrunelIaF', Dynamics(name='BrunelIaF')),
 ('current', Dimension(name='current', i=1)),
 ('resistance', Dimension(name='resistance', i=-2, m=1, t=-3, l=2)),
 ('time', Dimension(name='time', t=1)),
 ('voltage', Dimension(name='voltage', i=-1, m=1, t=-3, l=2)]

This gives us a Document instance, a dictionary-like object containing a Dynamics definition of an integrate-and-fire neuron model, together with the definitions of the physical dimensions of parameters and state variables used in the model.

Now for a file containing an entire user layer model (with references to other NineML documents containing the abstraction layer definitions):

>>> doc = nineml.read("./network/Brunel2000/AI.xml")
>>> pprint(doc.items())
[('All': Selection(name='All')),
 ('Exc': Population(name='Exc', number=4000, cell=nrn)),
 ('Excitation': Projection(name="Excitation", source=Population(name='Exc', number=4000, cell=nrn), destination=Selection(name='All'), connectivity=BaseComponent(name="RandomExc", componentclass="RandomFanIn"), response=BaseComponent(name="syn", componentclass="AlphaPSR")plasticity=BaseComponent(name="ExcitatoryPlasticity", componentclass="StaticConnection"), delay=Delay(value=1.5, unit=ms), with 2 port-connections)),
 ('Ext': Population(name='Ext', number=5000, cell=stim)),
 ('External': Projection(name="External", source=Population(name='Ext', number=5000, cell=stim), destination=Selection(name='All'), connectivity=BaseComponent(name="OneToOne", componentclass="OneToOne"), response=BaseComponent(name="syn", componentclass="AlphaPSR")plasticity=BaseComponent(name="ExternalPlasticity", componentclass="StaticConnection"), delay=Delay(value=1.5, unit=ms), with 2 port-connections)),
 ('Hz': Unit(name='Hz', dimension='per_time', power=0)),
 ('Inh': Population(name='Inh', number=1000, cell=nrn)),
 ('Inhibition': Projection(name="Inhibition", source=Population(name='Inh', number=1000, cell=nrn), destination=Selection(name='All'), connectivity=BaseComponent(name="RandomInh", componentclass="RandomFanIn"), response=BaseComponent(name="syn", componentclass="AlphaPSR")plasticity=BaseComponent(name="InhibitoryPlasticity", componentclass="StaticConnection"), delay=Delay(value=1.5, unit=ms), with 2 port-connections)),
 ('Mohm': Unit(name='Mohm', dimension='resistance', power=6)),
 ('current': Dimension(name='current', i=1)),
 ('mV': Unit(name='mV', dimension='voltage', power=-3)),
 ('ms': Unit(name='ms', dimension='time', power=-3)),
 ('nA': Unit(name='nA', dimension='current', power=-9)),
 ('per_time': Dimension(name='per_time', t=-1)),
 ('resistance': Dimension(name='resistance', i=-2, m=1, t=-3, l=2)),
 ('time': Dimension(name='time', t=1)),
 ('voltage': Dimension(name='voltage', i=-1, m=1, t=-3, l=2))]

Again we get a Document instance object containing all the NineML objects in the document. An alternative representation can be obtained by reading the file as a Network object:

>>> from nineml.user import Network
>>> net = doc.read("./network/Brunel2000/AI.xml").as_network('BrunelAI')
>>> print(net)
Network(name='BrunelAI')

This gives a much more structured representation. For example, all the Populations within the model are available through the populations attribute:

>>> pprint(list(net.populations))
[Population(name='Exc', number=4000, cell=nrn),
 Population(name='Ext', number=5000, cell=stim),
 Population(name='Inh', number=1000, cell=nrn)]

Introspecting NineML models

Introspecting abstraction layer models

Once we have loaded a model from an XML file we can begin to examine its structure.

>>> model = doc['BrunelIaF']
>>> model
Dynamics(name='BrunelIaF')

We can see a list of model parameters:

>>> pprint(list(model.parameters))
[Parameter(theta, dimension=voltage),
 Parameter(Vreset, dimension=voltage),
 Parameter(R, dimension=resistance),
 Parameter(tau_rp, dimension=time),
 Parameter(tau, dimension=time)]

a list of state variables:

>>> pprint(list(model.state_variables))
[StateVariable(V, dimension=voltage),
 StateVariable(t_rpend, dimension=time)]

and a list of the variables that are imported from/exposed to the outside world:

>>> pprint(list(model.ports))
[AnalogSendPort('V', dimension='Dimension(name='voltage', i=-1, m=1, t=-3, l=2)'),
 AnalogSendPort('t_rpend', dimension='Dimension(name='time', t=1)'),
 AnalogReducePort('Isyn', dimension='Dimension(name='current', i=1)', op='+'),
 EventSendPort('spikeOutput')]

Delving more deeply, we can examine the model’s regimes more closely:

>>> pprint(list(model.regimes))
[Regime(refractoryRegime),
 Regime(subthresholdRegime)]
>>> r_ref, r_sth = model.regimes

Looking first at the subthreshold regime, we can see the differential equations:

>>> list(r_sth.time_derivatives)
[TimeDerivative( dV/dt = (-V + R*Isyn)/tau )]

and the conditions under which the model will transition to the refractory regime:

>>> list(r_sth.transitions)
[OnCondition( V > theta )]
>>> tr_spike = next(r_sth.transitions)

The trigger for this transition is for the variable V to pass a threshold (parameter theta):

>>> tr_spike.trigger
Trigger('V > theta')

When the transition is initiated, the model will emit an output event (i.e. a spike) and discontinuously change the values of some of the state variables:

>>> tr_spike.output_events
[OutputEvent('spikeOutput')]
>>> tr_spike.state_assignments
[StateAssignment('t_rpend', 't + tau_rp'), StateAssignment('V', 'Vreset')]

Then it will move to the refractory regime:

>>> tr_spike.target_regime
Regime(refractoryRegime)

The refractory regime can be introspected in a similar way.

Introspecting user layer models

As shown above, once a complete network model has been loaded as a Network object, we can look at its neuron populations and the connections between these populations (“projections”):

>>> pprint(list(net.populations))
[Population(name='Exc', number=4000, cell=nrn),
 Population(name='Ext', number=5000, cell=stim),
 Population(name='Inh', number=1000, cell=nrn)]

>>> pprint(list(net.projections))
[Projection(name="Inhibition", pre=Population(name='Inh', size=2500, cell=nrn), post=Selection(name='All', Concatenate(Item(name='0'), Item(name='1'))), connectivity=Connectivity(rule=RandomFanIn, src_size=2500, dest_size=12500), response=DynamicsProperties(name="syn", component_class="Alpha")plasticity=DynamicsProperties(name="InhibitoryPlasticity", component_class="Static"), delay=1.5 * ms,event_port_connections=[EventPortConnection(sender=role:pre->spike_output, receiver=role:response->input_spike)], analog_port_connections=[AnalogPortConnection(sender=role:response->i_synaptic, receiver=role:post->i_synaptic), AnalogPortConnection(sender=role:plasticity->fixed_weight, receiver=role:response->weight)]),
 Projection(name="External", pre=Population(name='Ext', size=12500, cell=stim), post=Selection(name='All', Concatenate(Item(name='0'), Item(name='1'))), connectivity=Connectivity(rule=OneToOne, src_size=12500, dest_size=12500), response=DynamicsProperties(name="syn", component_class="Alpha")plasticity=DynamicsProperties(name="ExternalPlasticity", component_class="Static"), delay=1.5 * ms,event_port_connections=[EventPortConnection(sender=role:pre->spike_output, receiver=role:response->input_spike)], analog_port_connections=[AnalogPortConnection(sender=role:response->i_synaptic, receiver=role:post->i_synaptic), AnalogPortConnection(sender=role:plasticity->fixed_weight, receiver=role:response->weight)]),
 Projection(name="Excitation", pre=Population(name='Exc', size=10000, cell=nrn), post=Selection(name='All', Concatenate(Item(name='0'), Item(name='1'))), connectivity=Connectivity(rule=RandomFanIn, src_size=10000, dest_size=12500), response=DynamicsProperties(name="syn", component_class="Alpha")plasticity=DynamicsProperties(name="ExcitatoryPlasticity", component_class="Static"), delay=1.5 * ms,event_port_connections=[EventPortConnection(sender=role:pre->spike_output, receiver=role:response->input_spike)], analog_port_connections=[AnalogPortConnection(sender=role:response->i_synaptic, receiver=role:post->i_synaptic), AnalogPortConnection(sender=role:plasticity->fixed_weight, receiver=role:response->weight)])]

NineML also supports “selections”, groupings of neurons which span populations:

>>> pprint(list(net.selections))
[Selection(name='All', Concatenate(Item(name='0'), Item(name='1')))]

Note

in NineML version 1, the only type of selection is a concatenation of two or more populations. In future versions it will be possible to select and combine sub-populations.

Looking more closely at a population, we can see its name, the number of neurons it contains and the neuron model used (Component):

>>> p_exc = net.population('Exc')
>>> p_exc
Population(name='Exc', size=4000, cell=nrn)
>>> p_exc.size
4000
>>> p_exc.cell
DynamicsProperties(name="nrn", componentclass="BrunelIaF")

In the neuron model component we can see its abstraction layer definition (ComponentClass), it’s properties (parameter values), and the initial values of its state variables.

Note

the handling of initial values is likely to change in future versions of NineML.

>>> p_exc.cell.component_class
Dynamics(name='BrunelIaF')
>>> pprint(list(p_exc.cell.properties))
[Property(name=Vreset, value=10.0, unit=mV),
 Property(name=tau, value=20.0, unit=ms),
 Property(name=R, value=1.5, unit=Mohm),
 Property(name=tau_rp, value=2.0, unit=ms),
 Property(name=theta, value=20.0, unit=mV)]
>>> pprint(list(p_exc.cell.initial_values))
[Initial(name='t_rpend', value=0.0, unit=ms),
 Initial(name='V', value=0.0, unit=mV)]

Turning from a population to a projection:

>>> prj_inh = net.projection('Inhibition')
>>> prj_inh.pre
Population(name='Inh', number=1000, cell=nrn)
>>> prj_inh.post
Selection(name='All', Concatenate(Item(name='0'), Item(name='1')))
>>> prj_inh.response
DynamicsProperties(name="syn", componentclass="AlphaPSR")
>>> prj_inh.connectivity
DynamicsProperties(name="RandomInh", componentclass="RandomFanIn")
>>> prj_inh.plasticity
DynamicsProperties(name="InhibitoryPlasticity", componentclass="StaticConnection")
>>> prj_inh.delay
1.5 * ms
>>> pprint(list(prj_inh.port_connections))
[AnalogPortConnection(sender=role:response->i_synaptic, receiver=role:post->i_synaptic),
 AnalogPortConnection(sender=role:plasticity->fixed_weight, receiver=role:response->weight),
 EventPortConnection(sender=role:pre->spike_output, receiver=role:response->input_spike)]

Note that the pre and post attributes point to Populations or Projections, the connectivity rule, the post-synaptic response model and the synaptic plasticity model are all Components. The port_connections attribute indicates which ports in the different components should be connected together.

Writing model descriptions in Python

Writing abstraction layer models

subthreshold_regime = Regime(
    name="subthreshold_regime",
    time_derivatives=[
        "dV/dt = alpha*V*V + beta*V + zeta - U + Isyn / C_m",
        "dU/dt = a*(b*V - U)", ],

    transitions=[On("V > theta",
                    do=["V = c",
                        "U =  U+ d",
                        OutputEvent('spike')],
                    to='subthreshold_regime')]
)

ports = [AnalogSendPort("V", un.voltage),
         AnalogReducePort("Isyn", un.current, operator="+")]

parameters = [
    Parameter('theta', un.voltage),
    Parameter('a', un.per_time),
    Parameter('b', un.per_time),
    Parameter('c', un.voltage),
    Parameter('d', un.voltage / un.time),
    Parameter('C_m', un.capacitance),
    Parameter('alpha', un.dimensionless / (un.voltage * un.time)),
    Parameter('beta', un.per_time),
    Parameter('zeta', un.voltage / un.time)]

state_variables = [
    StateVariable('V', un.voltage),
    StateVariable('U', un.voltage / un.time)]

izhi = Dynamics(
    name="Izhikevich",
    parameters=parameters,
    state_variables=state_variables,
    regimes=[subthreshold_regime],
    analog_ports=ports)

Writing user layer models

# Meta-parameters
order = 1000       # scales the size of the network
Ne = 4 * order     # number of excitatory neurons
Ni = 1 * order     # number of inhibitory neurons
epsilon = 0.1      # connection probability
Ce = int(epsilon * Ne)  # number of excitatory synapses per neuron
Ci = int(epsilon * Ni)  # number of inhibitory synapses per neuron
Cext = Ce          # effective number of external synapses per neuron
delay = 1.5        # (ms) global delay for all neurons in the group
J = 0.1            # (mV) EPSP size
Jeff = 24.0 * J      # (nA) synaptic weight
Je = Jeff          # excitatory weights
Ji = -g * Je       # inhibitory weights
Jext = Je          # external weights
theta = 20.0       # firing thresholds
tau = 20.0         # membrane time constant
tau_syn = 0.1      # synapse time constant
# nu_thresh = theta / (Je * Ce * tau * exp(1.0) * tau_syn) # threshold rate
nu_thresh = theta / (J * Ce * tau)
nu_ext = eta * nu_thresh      # external rate per synapse
input_rate = 1000.0 * nu_ext * Cext   # mean input spiking rate

# Parameters
neuron_parameters = dict(tau=tau * ms,
                         v_threshold=theta * mV,
                         refractory_period=2.0 * ms,
                         v_reset=10.0 * mV,
                         R=1.5 * Mohm)  # units??
psr_parameters = dict(tau=tau_syn * ms)

# Initial Values
v_init = RandomDistributionProperties(
    "uniform_rest_to_threshold",
    ninemlcatalog.load("randomdistribution/Uniform",
                       'UniformDistribution'),
    {'minimum': (0.0, unitless),
     'maximum': (theta, unitless)})
neuron_initial_values = {"v": (v_init * mV),
                         "refractory_end": (0.0 * ms)}
synapse_initial_values = {"a": (0.0 * nA), "b": (0.0 * nA)}
tpoisson_init = RandomDistributionProperties(
    "exponential_beta",
    ninemlcatalog.load('randomdistribution/Exponential',
                       'ExponentialDistribution'),
    {"rate": (1000.0 / input_rate * unitless)})

# Dynamics components
celltype = DynamicsProperties(
    "nrn",
    ninemlcatalog.load('neuron/LeakyIntegrateAndFire',
                       'LeakyIntegrateAndFire'),
    neuron_parameters, initial_values=neuron_initial_values)
ext_stim = DynamicsProperties(
    "stim",
    ninemlcatalog.load('input/Poisson', 'Poisson'),
    dict(rate=(input_rate, Hz)),
    initial_values={"t_next": (tpoisson_init, ms)})
psr = DynamicsProperties(
    "syn",
    ninemlcatalog.load('postsynapticresponse/Alpha', 'Alpha'),
    psr_parameters,
    initial_values=synapse_initial_values)

# Connecion rules
one_to_one_class = ninemlcatalog.load(
    '/connectionrule/OneToOne', 'OneToOne')
random_fan_in_class = ninemlcatalog.load(
    '/connectionrule/RandomFanIn', 'RandomFanIn')

# Populations
exc_cells = Population("Exc", Ne, celltype, positions=None)
inh_cells = Population("Inh", Ni, celltype, positions=None)
external = Population("Ext", Ne + Ni, ext_stim, positions=None)

# Selections
all_cells = Selection(
    "All", Concatenate(exc_cells, inh_cells))

# Projections
input_prj = Projection(
    "External", external, all_cells,
    connectivity=ConnectionRuleProperties(
        "OneToOne", one_to_one_class),
    response=psr,
    plasticity=DynamicsProperties(
        "ExternalPlasticity",
        ninemlcatalog.load("plasticity/Static", 'Static'),
        properties={"weight": (Jext, nA)}),
    port_connections=[
        EventPortConnection(
            'pre', 'response', 'spike_output', 'spike'),
        AnalogPortConnection(
            "plasticity", "response", "fixed_weight", "weight"),
        AnalogPortConnection(
            "response", "destination", "i_synaptic", "i_synaptic")],
    delay=(delay, ms))

exc_prj = Projection(
    "Excitation", exc_cells, all_cells,
    connectivity=ConnectionRuleProperties(
        "RandomExc", random_fan_in_class, {"number": (Ce * unitless)}),
    response=psr,
    plasticity=DynamicsProperties(
        "ExcitatoryPlasticity",
        ninemlcatalog.load("plasticity/Static", 'Static'),
        properties={"weight": (Je, nA)}),
    port_connections=[
        EventPortConnection(
            'pre', 'response', 'spike_output', 'spike'),
        AnalogPortConnection(
            "plasticity", "response", "fixed_weight", "weight"),
        AnalogPortConnection(
            "response", "destination", "i_synaptic", "i_synaptic")],
    delay=(delay, ms))

inh_prj = Projection(
    "Inhibition", inh_cells, all_cells,
    connectivity=ConnectionRuleProperties(
        "RandomInh", random_fan_in_class, {"number": (Ci * unitless)}),
    response=psr,
    plasticity=DynamicsProperties(
        "InhibitoryPlasticity",
        ninemlcatalog.load("plasticity/Static", 'Static'),
        properties={"weight": (Ji, nA)}),
    port_connections=[
        EventPortConnection(
            'pre', 'response', 'spike_output', 'spike'),
        AnalogPortConnection(
            "plasticity", "response", "fixed_weight", "weight"),
        AnalogPortConnection(
            "response", "destination", "i_synaptic", "i_synaptic")],
    delay=(delay, ms))

# Save to document in NineML Catalog
network = Network(name if name else "BrunelNetwork")
network.add(exc_cells, inh_cells, external, all_cells, input_prj, exc_prj,
            inh_prj)