brian2tools documentation¶

Installation with Anaconda¶

Since brian2tools (and brian2 on which it depends) are not part of the main Anaconda distribution, you have to install it from the brian-team channel. To do so use:

conda install -c brian-team brian2tools

You can also permanently add the channel to your list of channels:

conda config --add channels brian-team

This has only to be done once. After that, you can install and update the brian2 packages as any other Anaconda package:

conda install brian2tools

conda install mayavi

Installation with pip¶

If you decide not to use Anaconda, you can install brian2tools from the Python package index: https://pypi.python.org/pypi/brian2tools

To do so, use the pip utility:

pip install brian2tools

You might want to add the --user flag, to install Brian 2 for the local user only, which means that you don’t need administrator privileges for the installation.

If you have an older version of pip, first update pip itself:

# On Linux/MacOsX:
pip install -U pip

# On Windows
python -m pip install -U pip

If you don’t have pip but you have the easy_install utility, you can use it to install pip:

easy_install pip

If you have neither pip nor easy_install, use the approach described here to install pip: https://pip.pypa.io/en/latest/installing.htm

Plotting recorded activity¶

We’ll use the following example (the CUBA example from Brian 2) as a demonstration.

from brian2 import *

Vt = -50*mV
Vr = -60*mV

eqs = '''dv/dt  = (ge+gi-(v + 49*mV))/(20*ms) : volt (unless refractory)
         dge/dt = -ge/(5*ms) : volt
         dgi/dt = -gi/(10*ms) : volt
      '''
P = NeuronGroup(4000, eqs, threshold='v>Vt', reset='v = Vr', refractory=5*ms,
                method='linear')
P.v = 'Vr + rand() * (Vt - Vr)'
P.ge = 0*mV
P.gi = 0*mV

we = (60*0.27/10)*mV # excitatory synaptic weight (voltage)
wi = (-20*4.5/10)*mV # inhibitory synaptic weight
Ce = Synapses(P[:3200], P, on_pre='ge += we')
Ci = Synapses(P[3200:], P, on_pre='gi += wi')
Ce.connect(p=0.02)
Ci.connect(p=0.02)

spike_mon = SpikeMonitor(P)
rate_mon = PopulationRateMonitor(P)
state_mon = StateMonitor(P, 'v', record=[0, 100, 1000])  # record three cells

run(1 * second)

We will also assume that brian2tools has been imported like this:

from brian2tools import *

Spikes¶

To plot a basic raster plot, you can call brian_plot with the SpikeMonitor as its argument:

brian_plot(spike_mon)

To have more control over the plot, or to plot spikes that are not stored in a SpikeMonitor, use plot_raster:

plot_raster(spike_mon.i, spike_mon.t, time_unit=second, marker=',', color='k')

Rates¶

Calling brian_plot with the PopulationRateMonitor will plot the rate smoothed with a Gaussian window with 1ms standard deviation.:

brian_plot(rate_mon)

To plot the rate with a different smoothing and/or to set other details of the plot use plot_raster:

plot_rate(rate_mon.t, rate_mon.smooth_rate(window='flat', width=10.1*ms),
          linewidth=3, color='gray')

State variables¶

Finally, calling brian_plot with the StateMonitor will plot the recorded voltage traces:

brian_plot(state_mon)

By indexing the StateMonitor, the plot can be restricted to a subset of the recorded neurons:

brian_plot(state_mon[1000])

Again, for more detailed control you can directly use the plot_state function. Here we also demonstrate the use of the returned Axes object to add a legend to the plot:

ax = plot_state(state_mon.t, state_mon.v.T, var_name='membrane potential', lw=2)
ax.legend(['neuron 0', 'neuron 100', 'neuron 1000'], frameon=False, loc='best')

Plotting synaptic connections and variables¶

For the following examples, we create synapses and synaptic weights according to “distances” (differences between the source and target indices):

from brian2 import *

group = NeuronGroup(100, 'dv/dt = -v / (10*ms) : volt',
                    threshold='v > -50*mV', reset='v = -60*mV')

synapses = Synapses(group, group, 'w : volt', on_pre='v += w')

# Connect to cells with indices no more than +/- 10 from the source index with
# a probability of 50% (but do not create self-connections)
synapses.connect(j='i+k for k in sample(-10, 10, p=0.5) if k != 0',
                 skip_if_invalid=True)  # ignore values outside of the limits
# Set synaptic weights depending on the distance (in terms of indices) between
# the source and target cell and add some randomness
synapses.w = '(exp(-(i - j)**2/10.) + 0.5 * rand())*mV'
# Set synaptic weights randomly
synapses.delay = '1*ms + 2*ms*rand()'

Connections¶

A call of brian_plot with a Synapses object will plot all connections, plotting either the matrix as an image, the connections as a scatter plot, or a 2-dimensional histogram (using matplotlib’s hexbin function). The decision which type of plot to use is based on some heuristics applied to the number of synapses and might possibly change in future versions:

brian_plot(synapses)

As explained above, for a large connection matrix this would instead use an approach based on a hexagonal 2D histogram:

big_group = NeuronGroup(10000, '')
many_synapses = Synapses(big_group, big_group)
many_synapses.connect(j='i+k for k in range(-2000, 2000) if rand() < exp(-(k/1000.)**2)',
                      skip_if_invalid=True)
brian_plot(many_synapses)

Under the hood brian_plot calls plot_synapses which can also be used directly for more control:

plot_synapses(synapses.i, synapses.j, plot_type='scatter', color='gray', marker='s')

Synaptic variables (weights, delays, etc.)¶

Synaptic variables such as synaptic weights or delays can also be plotted with brian_plot:

subplot(1, 2, 1)
brian_plot(synapses.w)
subplot(1, 2, 2)
brian_plot(synapses.delay)
tight_layout()

Again, using plot_synapses provides more control. The following code snippet also calls the add_background_pattern function to make the distinction between white color values and the background clearer:

ax = plot_synapses(synapses.i, synapses.j, synapses.w, var_name='synaptic weights',
                   plot_type='scatter', cmap='hot')
add_background_pattern(ax)
ax.set_title('Recurrent connections')

Multiple synapses per source-target pair¶

In Brian, source-target pairs can be connected by more than a single synapse. In this case you cannot plot synaptic state variables (because it is ill-defined what to plot) but you can still plot connections which will show how many synapses exists. For example, if we make the same connect from above a second time, the new synapses will be added to the existing ones so some source-target pairs are now connected by two synapses:

synapses.connect(j='i+k for k in sample(-10, 10, p=0.5) if k != 0',
                 skip_if_invalid=True)

Calling brian_plot or plot_synapses will now show the number of synapses between each pair of neurons:

brian_plot(synapses)

Plotting neuronal morphologies¶

In the following, we’ll use a reconstruction from the Destexhe lab (a neocortical pyramidal neuron from the cat brain [1]) that we load into Brian:

from brian2 import *

morpho = Morphology.from_file('51-2a.CNG.swc')

Dendograms¶

Calling brian_plot with a Morphology will plot a dendogram:

brian_plot(morpho)

The plot_dendrogram function does the same thing, but in contrast to the other plot functions it does not allow any customization at the moment, so there is no benefit over using brian_plot.

Morphologies in 2D or 3D¶

In addition to the dendogram which only plots the general structure but not the actual morphology of the neuron in space, you can plot the morphology using plot_morphology. For a 3D morphology, this will plot the morphology in 3D using the Mayavi package

plot_morphology(morpho)

For artificially created morphologies (where one might only use coordinates in 2D) or to get a quick view of a morphology, you can also plot it in 2D (this will be done automatically if the coordinates are 2D only):

plot_morphology(morpho, plot_3d=False)

Both 2D and 3D morphology plots can be further customized, e.g. they can show the width of the compartments and do not use the default alternation between blue and red for each section:

plot_morphology(morpho, plot_3d=True, show_compartments=True,
                show_diameter=True, colors=('darkblue',))

[1]	Available at http://neuromorpho.org/neuron_info.jsp?neuron_name=51-2a

Working example¶

As a demonstration, we use a simple unconnected Integrate & Fire neuron model with refractoriness and given initial values.

from brian2 import *
import brian2tools.nmlexport

set_device('neuroml2', filename="nml2model.xml")

n = 100
duration = 1*second
tau = 10*ms

eqs = '''
dv/dt = (v0 - v) / tau : volt (unless refractory)
v0 : volt
'''
group = NeuronGroup(n, eqs, threshold='v > 10*mV', reset='v = 0*mV',
                    refractory=5*ms, method='linear')
group.v = 0*mV
group.v0 = '20*mV * i / (N-1)'

rec_idx = [2, 63]
statemonitor = StateMonitor(group, 'v', record=rec_idx)
spikemonitor = SpikeMonitor(group, record=rec_idx)

run(duration)

The use of the exporter requires only a few changes to an existing Brian 2 script. In addition to the standard brian2 import at the beginning of your script, you need to import the brian2tools.nmlexport package. You can then set a “device” called neuroml2 which will generate NeuroML2/LEMS code instead of executing your model. You will also have to specify a keyword argument filename with the desired name of the output file.

The above code will result in a file nml2model.xml and an additional file LEMSUnitsConstants.xml with units definitions in form of constants (necessary due to the way units are handled in LEMS equations).

The file nml2model.xml will look like this:

<Lems>
  <Include file="NeuroML2CoreTypes.xml"/>
  <Include file="Simulation.xml"/>
  <Include file="LEMSUnitsConstants.xml"/>
  <ComponentType extends="baseCell" name="neuron1">
    <Property dimension="voltage" name="v0"/>
    <Property dimension="time" name="tau"/>
    <EventPort direction="out" name="spike"/>
    <Exposure dimension="voltage" name="v"/>
    <Dynamics>
      <StateVariable dimension="voltage" exposure="v" name="v"/>
      <OnStart>
        <StateAssignment value="0" variable="v"/>
      </OnStart>
      <Regime name="refractory">
        <StateVariable dimension="time" name="lastspike"/>
        <OnEntry>
          <StateAssignment value="t" variable="lastspike"/>
        </OnEntry>
        <OnCondition test="t .gt. ( lastspike + 5.*ms )">
          <Transition regime="integrating"/>
        </OnCondition>
      </Regime>
      <Regime initial="true" name="integrating">
        <TimeDerivative value="(v0 - v) / tau" variable="v"/>
        <OnCondition test="v .gt. (10 * mV)">
          <EventOut port="spike"/>
          <StateAssignment value="0*mV" variable="v"/>
          <Transition regime="refractory"/>
        </OnCondition>
      </Regime>
    </Dynamics>
  </ComponentType>
  <ComponentType extends="basePopulation" name="neuron1Multi">
    <Parameter dimension="time" name="tau_p"/>
    <Parameter dimension="none" name="N"/>
    <Constant dimension="voltage" name="mVconst" symbol="mVconst" value="1mV"/>
    <Structure>
      <MultiInstantiate componentType="neuron1" number="N">
        <Assign property="v0" value="20*mVconst * index /  ( N-1 ) "/>
        <Assign property="tau" value="tau_p"/>
      </MultiInstantiate>
    </Structure>
  </ComponentType>
  <network id="neuron1MultiNet">
    <Component N="100" id="neuron1Multipop" tau_p="10. ms" type="neuron1Multi"/>
  </network>
  <Simulation id="sim1" length="1s" step="0.1 ms" target="neuron1MultiNet">
    <Display id="disp0" timeScale="1ms" title="v" xmax="1000" xmin="0" ymax="11" ymin="0">
      <Line id="line3" quantity="neuron1Multipop[3]/v" scale="1mV" timeScale="1ms"/>
      <Line id="line64" quantity="neuron1Multipop[64]/v" scale="1mV" timeScale="1ms"/>
    </Display>
    <OutputFile fileName="recording_nml2model.dat" id="of0">
      <OutputColumn id="3" quantity="neuron1Multipop[3]/v"/>
      <OutputColumn id="64" quantity="neuron1Multipop[64]/v"/>
    </OutputFile>
    <EventOutputFile fileName="recording_nml2model.spikes" format="TIME_ID" id="eof">
      <EventSelection eventPort="spike" id="line3" select="neuron1Multipop[3]"/>
      <EventSelection eventPort="spike" id="line64" select="neuron1Multipop[64]"/>
    </EventOutputFile>
  </Simulation>
  <Target component="sim1"/>
</Lems>

The exporting device creates a new ComponentType for each cell definition implemented as a Brian 2 NeuronGroup. Later that particular ComponentType is bundled with the initial value assignment into a a new ComponentType (here called neuron1Multi) by MultiInstantiate and eventually a network (neuron1MultiNet) is created out of a defined Component (neuron1Multipop).

Note that the integration method does not matter for the NeuroML export, as NeuroML/LEMS only describes the model not how it is numerically integrated.

To validate the output, you can use the tool jNeuroML. Make sure that jnml has access to the NeuroML2CoreTypes folder by setting the JNML_HOME environment variable.

With jnml installed you can run the simulation as follows:

jnml nml2model.xml

Supported Features¶

Currently, the NeuroML2 export is restricted to simple neural models and only supports the following classes (and a single run statement per script):

• NeuronGroup - The definition of a neuronal model. Mechanisms like threshold, reset and refractoriness are taken into account. Moreover, you may set the initial values of the model parameters (like v0 above).
• StateMonitor - If your script uses a StateMonitor to record variables, each recorded variable is transformed into to a Line tag of the Display in the NeuroML2 simulation and an OutputFile tag is added to the model. The name of the output file is recording_<<filename>>.dat.
• SpikeMonitor - A SpikeMonitor is transformed into an EventOutputFile tag, storing the spikes to a file named recording_<<filename>>.spikes.

Limitations¶

As stated above, the NeuroML2 export is currently quite limited. In particular, none of the following Brian 2 features are supported:

• Synapses
• Network input (PoissonGroup, SpikeGeneratorGroup, etc.)
• Multicompartmental neurons (SpatialNeuronGroup)
• Non-standard simulation protocols (multiple runs, store/restore mechanism, etc.).

NMLExporter¶

The whole process of building NeuroML model starts with calling the create_lems_model method. It selects crucial Brian 2 objects to further parse and pass them to respective methods.

if network is None:
    net = Network(collect(level=1))
else:
    net = network

if not constants_file:
    self._model.add(lems.Include(LEMS_CONSTANTS_XML))
else:
    self._model.add(lems.Include(constants_file))
includes = set(includes)
for incl in INCLUDES:
    includes.add(incl)
neuron_groups  = [o for o in net.objects if type(o) is NeuronGroup]
state_monitors = [o for o in net.objects if type(o) is StateMonitor]
spike_monitors = [o for o in net.objects if type(o) is SpikeMonitor]

for o in net.objects:
    if type(o) not in [NeuronGroup, StateMonitor, SpikeMonitor,
                       Thresholder, Resetter, StateUpdater]:
        logger.warn("""{} export functionality
                       is not implemented yet.""".format(type(o).__name__))

# Thresholder, Resetter, StateUpdater are not interesting from our perspective
if len(netinputs)>0:
    includes.add(LEMS_INPUTS)
for incl in includes:
    self.add_include(incl)
# First step is to add individual neuron deifinitions and initialize
# them by MultiInstantiate
for e, obj in enumerate(neuron_groups):
    self.add_neurongroup(obj, e, namespace, initializers)

# DOM structure of the whole model is constructed below
self._dommodel = self._model.export_to_dom()
# input support - currently only Poisson Inputs
for e, obj in enumerate(netinputs):
    self.add_input(obj, counter=e)
# A population should be created in *make_multiinstantiate*
# so we can add it to our DOM structure.
if self._population:
    self._extend_dommodel(self._population)
# if some State or Spike Monitors occur we support them by
# Simulation tag
self._model_namespace['simulname'] = "sim1"
self._simulation = NeuroMLSimulation(self._model_namespace['simulname'],
                                     self._model_namespace['networkname'])
for e, obj in enumerate(state_monitors):
    self.add_statemonitor(obj, filename=recordingsname, outputfile=True)
for e, obj in enumerate(spike_monitors):
    self.add_spikemonitor(obj, filename=recordingsname)

simulation = self._simulation.build()
self._extend_dommodel(simulation)
target = NeuroMLTarget(self._model_namespace['simulname'])
target = target.build()
self._extend_dommodel(target)

LEMSDevice¶

LEMSDevice allows you to take advantage of Brian 2’s code generation mechanism. It makes usage of the module easier, as it means for user that they just need to import brian2tools.nmlexport and set the device neuroml2 like this:

import brian2lems.nmlexport

set_device('neuroml2', filename="ifcgmtest.xml")

In the class init a flag self.build_on_run was set to True which means that exporter starts working immediately after encountering the run statement.

def __init__(self):
    super(LEMSDevice, self).__init__()
    self.runs = []
    self.assignments = []
    self.build_on_run = True
    self.build_options = None
    self.has_been_run = False

First of all method network_run is called which gathers of necessary variables from the script or function namespaces and passes it to build method. In build we select all needed variables to separate dictionaries, create a name of the recording files and eventually build the exporter.

initializers = {}
for descriptions, duration, namespace, assignments in self.runs:
    for assignment in assignments:
        if not assignment[2] in initializers:
            initializers[assignment[2]] = assignment[-1]
if len(self.runs) > 1:
    raise NotImplementedError("Currently only single run is supported.")
if len(filename.split("."))!=1:
    filename_ = 'recording_' + filename.split(".")[0]
else:
    filename_ = 'recording_' + filename
exporter = NMLExporter()
exporter.create_lems_model(self.network, namespace=namespace,
                                         initializers=initializers,
                                         recordingsname=filename_)
exporter.export_to_file(filename)

if lems_const_save:
    with open(os.path.join(nmlcdpath, LEMS_CONSTANTS_XML), 'r') as f:
        with open(LEMS_CONSTANTS_XML, 'w') as fout:
            fout.write(f.read())

Other modules¶

If you want to know more about other scripts included in package ( lemsrendering, supporting, cgmhelper), please read their docstrings or comments included in the code.

TODO¶

• synapses support;

First attempt to make synapses export work was made during GSOC period. The problem with that feature is related to the fact that NeuroML and brian2 internal synapses implementation differs substantially. For instance, in NeuroML there are no predefined rules for connections, but user needs to explicitly define a synapse. Moreover, in Brian 2, for efficiency reasons, postsynaptic potentials are normally modeled in the post-synaptic cell (for linearly summating synapses, this is equivalent but much more efficient), whereas in NeuroML they are modeled as part of the synapse (simulation speed is not an issue here).

• network input support;

Although there are some classes supporting PoissonInput in the supporting.py, full functionality of input is still not provided, as it is stongly linked with above synapses problems.