Introduction to Brian part 3: Simulations¶
If you haven’t yet read parts 1 and 2 on Neurons and Synapses, go read them first.
This tutorial is about managing the slightly more complicated tasks that crop up in research problems, rather than the toy examples we’ve been looking at so far. So we cover things like inputting sensory data, modelling experimental conditions, etc.
As before we start by importing the Brian package and setting up matplotlib for IPython:
This tutorial is a static non-editable version. You can launch an interactive, editable version without installing any local files using the Binder service (although note that at some times this may be slow or fail to open):
Alternatively, you can download a copy of the notebook file
to use locally:
See the tutorial overview page for more details.
Let’s start by looking at a very common task: doing multiple runs of a simulation with some parameter that changes. Let’s start off with something very simple, how does the firing rate of a leaky integrate-and-fire neuron driven by Poisson spiking neurons change depending on its membrane time constant? Let’s set that up.
Now if you’re running the notebook, you’ll see that this was a little slow to run. The reason is that for each loop, you’re recreating the objects from scratch. We can improve that by setting up the network just once. We store a copy of the state of the network before the loop, and restore it at the beginning of each iteration.
That’s a very simple example of using store and restore, but you can use it in much more complicated situations. For example, you might want to run a long training run, and then run multiple test runs afterwards. Simply put a store after the long training run, and a restore before each testing run.
You can also see that the output curve is very noisy and doesn’t
increase monotonically like we’d expect. The noise is coming from the
fact that we run the Poisson group afresh each time. If we only wanted
to see the effect of the time constant, we could make sure that the
spikes were the same each time (although note that really, you ought to
do multiple runs and take an average). We do this by running just the
Poisson group once, recording its spikes, and then creating a new
SpikeGeneratorGroup that will output those recorded spikes each
You can see that now there is much less noise and it increases monotonically because the input spikes are the same each time, meaning we’re seeing the effect of the time constant, not the random spikes.
Note that in the code above, we created
Network objects. The reason
is that in the loop, if we just called
run it would try to simulate
all the objects, including the Poisson neurons
P, and we only want
to run that once. We use
Network to specify explicitly which objects
we want to include.
The techniques we’ve looked at so far are the conceptually most simple way to do multiple runs, but not always the most efficient. Since there’s only a single output neuron in the model above, we can simply duplicate that output neuron and make the time constant a parameter of the group.
WARNING "tau" is an internal variable of group "neurongroup_2", but also exists in the run namespace with the value 10. * msecond. The internal variable will be used. [brian2.groups.group.Group.resolve.resolution_conflict]
You can see that this is much faster again! It’s a little bit more complicated conceptually, and it’s not always possible to do this trick, but it can be much more efficient if it’s possible.
Let’s finish with this example by having a quick look at how the mean and standard deviation of the interspike intervals depends on the time constant.
Notice that we used the
spike_trains() method of
This is a dictionary with keys being the indices of the neurons and
values being the array of spike times for that neuron.
Changing things during a run¶
Imagine an experiment where you inject current into a neuron, and change the amplitude randomly every 10 ms. Let’s see if we can model that using a Hodgkin-Huxley type neuron.
In the code above, we used a loop over multiple runs to achieve this.
That’s fine, but it’s not the most efficient way to do it because each
time we call
run we have to do a lot of initialisation work that
slows everything down. It also won’t work as well with the more
efficient standalone mode of Brian. Here’s another way.
We’ve replaced the loop that had multiple
run calls with a
run_regularly. This makes the specified block of code run every
run_regularly lets you run code specific to a
NeuronGroup, but sometimes you might need more flexibility.
For this, you can use
network_operation which lets you run arbitrary
Python code (but won’t work with the standalone mode).
Now let’s extend this example to run on multiple neurons, each with a different capacitance to see how that affects the behaviour of the cell.
So that runs, but something looks wrong! The injected currents look like they’re different for all the different neurons! Let’s check:
Sure enough, it’s different each time. But why? We wrote
group.run_regularly('I = rand()*50*nA', dt=10*ms) which seems like
it should give the same value of I for each neuron. But, like threshold
and reset statements,
run_regularly code is interpreted as being run
separately for each neuron, and because I is a parameter, it can be
different for each neuron. We can fix this by making I into a shared
variable, meaning it has the same value for each neuron.
Ahh, that’s more like it!
Now let’s think about a neuron being driven by a sinusoidal input. Let’s go back to a leaky integrate-and-fire to simplify the equations a bit.
So far, so good and the sort of thing we saw in the first tutorial. Now,
what if that input current were something we had recorded and saved in a
file? In that case, we can use
TimedArray. Let’s start by
reproducing the picture above but using
Note that for the example where we put the
sin function directly in
the equations, we had to use the
method='euler' argument because the
exact integrator wouldn’t work here (try it!). However,
is considered to be constant over its time step and so the linear
integrator can be used. This means you won’t get the same behaviour from
these two methods for two reasons. Firstly, the numerical integration
euler give slightly different results.
sin is not constant over a timestep whereas
Now just to show that
TimedArray works for arbitrary currents, let’s
make a weird “recorded” current and run it on that.
Finally, let’s finish on an example that actually reads in some data from a file. See if you can work out how this example works.