Equations objects are initialised with a string as follows:
eqs=Equations('''
dx/dt=(y-x)/tau + a : volt # differential equation
y=2*x : volt # equation
z=x # alias
a : volt/second # parameter
''')
It is possible to pass a string instead of an Equations object when initialising a neuron group. In that case, the string is implicitly converted to an Equations object. There are 4 different types of equations:
Right hand sides must be valid Python expressions, possibly including comments and multiline characters (\).
The units of all variables except aliases must be specified. Note that in first line, the units volt are meant for x, not dx/dt. The consistency of all units is checked with the method check_units(), which is automatically called when initialising a neuron group (through the method prepare()).
When an Equations object is finalised (through the method prepare(), automatically called the NeuronGroup initialiser), the names of variables defined by non-differential equations are replaced by their (string) values, so that differential equations are self-consistent. In the process, names of external variables are also modified to avoid conflicts (by adding a prefix).
The key idea for efficient simulations is to update synchronously the state variables of all identical neuron models. A neuron group is defined by the model equations, and optionally a threshold condition and a reset. For example for 100 neurons:
eqs=Equations('dv/dt=-v/tau : volt')
group=NeuronGroup(100,model=eqs,reset=0*mV,threshold=10*mV)
The model keyword also accepts strings (in that case it is converted to an Equations object), e.g.:
group=NeuronGroup(100,model='dv/dt=-v/tau : volt',reset=0*mV,threshold=10*mV)
The units of both the reset and threshold are checked for consistency with the equations. The code above defines a group of 100 integrate-and-fire neurons with threshold 10 mV and reset 0 mV. The second line defines an object named group which contains all the state variables, which can be accessed with the dot notation, i.e. group.v is a vector with the values of variable v for all of the 100 neurons. It is an array with units as defined in the equations (here, volt). By default, all state variables are initialised at value 0. It can be initialised by the user as in the following example:
group.v=linspace(0*mV,10*mV,100)
Here the values of v for all the neurons are evenly spaced between 0 mV and 10 mV (linspace is a NumPy function). The method group.rest() may also be used to set the resting point of the equations, but convergence is not always guaranteed.
Subgroups can be created with the slice operator:
subgroup1=group[0:50]
subgroup2=group[50:100]
Then subgroup2.v[i] equals group.v[50+i]. An alternative equivalent method is the following:
subgroup1=group.subgroup(50)
subgroup2=group.subgroup(50)
The parent group keeps track of the allocated subgroups. But note that the two methods are mutually exclusive, e.g. in the following example:
subgroup1=group[0:50]
subgroup2=group.subgroup(50)
both subgroups are actually identical.
Subgroups are useful when creating connections or monitoring the state variables or spikes. The best practice is to define groups as large as possible, then divide them in subgroups if necessary. Indeed, the larger the groups are, the faster the simulation runs. For example, for a network with a feedforward architecture, one should first define one group holding all the neurons in the network, then define the layers as subgroups of this big group.
For details, see the reference documentation for NeuronGroup.
More complex resets can be defined. The value of the reset keyword can be:
The simplest way to customise the reset is to define it as a Python statement, e.g.:
eqs='''
dv/dt=-v/tau : volt
dw/dt=-w/tau : volt
'''
group=NeuronGroup(100,model=eqs,reset="v=0*mV;w+=3*mV",threshold=10*mV)
The string must be a valid Python statement (possibly a multiline string). It can contain variables from the neuron group, units and any variable defined in the namespace (e.g. tau), as for equations. Be aware that if a variable in the namespace has the same name as a neuron group variable, then it masks the neuron variable. The way it works is that the code is evaluated with each neuron variable v replaced by v[spikes], where spikes is the array of indexes of the neurons that just spiked.
To define a specific reset, the generic method is define a function as follows:
def myreset(P,spikes):
P.v[spikes]=rand(len(spikes))*5*mV
group=NeuronGroup(100,model=eqs,reset=myreset,threshold=10*mV)
or faster:
def myreset(P,spikes):
P.v_[spikes]=rand(len(spikes))*5*mV
Every time step, the user-defined function is called with arguments P, the neuron group, and spikes, the list of indexes of the neurons that just spiked. The function above resets the neurons that just spiked to a random value.
It is possible to specify the reset variable explicitly:
group=NeuronGroup(100,model=eqs,reset=Reset(0*mV,state='w'),threshold=10*mV)
Here the variable w is reset.
The value of the reset can be given by another state variable:
group=NeuronGroup(100,model=eqs,reset=VariableReset(0*mV,state='v',resetvaluestate='w'),threshold=10*mV)
Here the value of the variable w is used to reset the variable v.
As for the reset, the threshold can be customised.
The simplest way to customise the threshold is to define it as a Python expression, e.g.:
eqs='''
dv/dt=-v/tau : volt
dw/dt=(v-w)/tau : volt
'''
group=NeuronGroup(100,model=eqs,reset=0*mV,threshold="v>=w")
The string must be an expression which evaluates to a boolean. It can contain variables from the neuron group, units and any variable defined in the namespace (e.g. tau), as for equations. Be aware that if a variable in the namespace has the same name as a neuron group variable, then it masks the neuron variable. The way it works is that the expression is evaluated with the neuron variables replaced by their vector values (values for all neurons), so that the expression returns a boolean vector.
The generic method to define a custom threshold condition is to pass a function of the state variables which returns a boolean (true if the threshold condition is met), for example:
eqs='''
dv/dt=-v/tau : volt
dw/dt=(v-w)/tau : volt
'''
group=NeuronGroup(100,model=eqs,reset=0*mV,threshold=lambda v,w:v>=w)
Here we used an anonymous function (lambda keyword) but of course a named function can also be used. In this example, spikes are generated when v is greater than w. Note that the arguments of the function must be the state variables with the same order as in the Equations string.
It is possible to specify the threshold variable explicitly:
group=NeuronGroup(100,model=eqs,reset=0*mV,threshold=Threshold(0*mV,state='w'))
Here the variable w is checked.
The same model as in the functional threshold example can be defined as follows:
group=NeuronGroup(100,model=eqs,reset=0*mV,threshold=\
VariableThreshold(state='v',threshold_state='w'))
For Hodgkin-Huxley models, one needs to determine the threshold empirically. Here the threshold should really be understood rather as the onset of the spikes (used to propagate the spikes to the other neurons), since there is no explicit reset. There is a Threshold subclass for this purpose:
group=NeuronGroup(100,model=eqs,threshold=EmpiricalThreshold(threshold=-20*mV,refractory=3*ms))
Spikes are triggered when the membrane potential reaches the value -20 mV, but only if it has not spiked in the last 3 ms (otherwise there would be spikes every time step during the action potential). The state keyword may be used to specify the state variable which should be checked for the threshold condition.
It is possible to generate spikes with a given probability rather than when a threshold condition is met, by using the class PoissonThreshold, as in the following example:
group=NeuronGroup(100,model='x : Hz',threshold=PoissonThreshold(state='x'))
x=linspace(0*Hz,10*Hz,100)
Here spikes are generated as Poisson processes with rates given by the variable x (the state keyword is optional: default = first variable defined). Note that x can change over time (inhomogeneous Poisson processes). The units of variable x must be Hertz.