More on equations¶
Equations may contain external variables. When an
Equations object is initialised,
a dictionary is built with the values of all external variables. These values are taken
from the namespace where the
Equations object was defined. It is possible to go one or
several levels up in the namespaces by specifying the keyword
The value of these parameters can in general be changed during the simulation and the
modifications are taken into account, except in two situations: when the equations are
frozen (see below) or when the integration is exact (linear equations). In those cases,
the values of the parameters are the ones at initialisation time.
Alternatively, the string defining the equations can be evaluated within a given namespace by providing keywords at initialisation time, e.g.:
eqs=Equations('dx/dt=-x/tau : volt',tau=10*ms)
In that case, the values of all external variables are taken from the specified dictionary (given by the keyword arguments), even if variables with the same name exist in the namespace where the string was defined. The two methods for passing the values of external variables are mutually exclusive, that is, either all external variables are explicitly specified with keywords (if not, they are left unspecified even if there are variables with the same names in the namespace where the string was defined), or all values are taken from the calling namespace.
More can be done with keyword arguments. If the value is a string, then the name of the variable is replaced, e.g.:
eqs=Equations('dx/dt=-x/tau : volt',tau=10*ms,x='Vm')
changes the variable name x to Vm. This is useful for writing functions which return equations where the variable name is provided by the user.
Finally, if the value is
None then the name of the variable is replaced by a unique
eqs=Equations('dx/dt=-x/tau : volt',tau=10*ms,x=None)
This is useful to avoid conflicts in the names of hidden variables.
Equations can be combined using the sum operator. For example:
eqs=Equations('dx/dt=(y-x)/tau : volt') eqs+=Equations('dy/dt=-y/tau: volt')
Note that some variables may be undefined when defining the first equation. No error is
raised when variables are undefined and absent from the calling namespace.
Equations objects are added, the consistency is checked. For example it is not
possible to add two
Equations objects which define the same variable.
Which variable is the membrane potential?¶
Several objects, such as
can be initialised without specifying which variable is the membrane potential, in which
case it is assumed that it is the first variable.
Internally, the variables of an
Equations object are reorderered so that the first one
is most likely to be the membrane potential (using
The first variable is, with decreasing priority :
- the first defined variable.
The currently available integration methods are:
- Exact integration when the equations are linear.
- Euler integration (explicit, first order).
- Runge-Kutta integration (explicit, second order).
- Exponential Euler integration (implicit, first order).
The method is selected when a
NeuronGroup is initialized.
If the equations are linear, exact integration is automatically selected.
Otherwise, Euler integration is selected by default, unless the keyword
implicit=True is passed, which selects the exponential Euler method. A second-order method
can be selected using the keyword
order=2 (explicit Runge-Kutta method, midpoint estimation).
It is possible to override this behaviour with the
method keyword when initialising
NeuronGroup. Possible values are
If the differential equations are linear, then the update phase
X(t)->X(t+dt) can be calculated exactly with a matrix product.
First, the equations are examined to determine whether they are linear
with the method
islinear() and the function
is_affine() (this is currently done using dynamic typing).
Second, the matrix M and the vector B such that dX/dt=M(X-B) are calculated with
Third, the matrix A such that X(t+dt)=A*(X(t)-B)+B is calculated at initialisation
of a specific state updater object,
as A=expm(M*dt), where expm is the matrix exponential.
Important remark: since the update matrix and vector are precalculated, the values of all external variables in the equations are frozen at initialisation. If external variables are modified after initialisation, those modifications are not taken into account during the simulation.
Inexact exact integration: If the equation cannot be put into the form dX/dt=M(X-B), for example if the equation is dX/dt=MX+A where M is not invertible, then the equations are not integrated exactly, but using a system equivalent to Euler integration but with dt 100 times smaller than specified. Updates are of the form X(t+dt)=A*X(t)+C where the matrix A and vector C are computed by applying Euler integration 100 times to the differential equations.
The Euler is a first order explicit integration method. It is the default one for nonlinear equations. It is simply implemented as X(t+dt)=X(t)+f(X)*dt.
Exponential Euler integration¶
The exponential Euler method is used for Hodgkin-Huxley type equations, are which stiff. Equations of that type are conditionally linear, that is, the differential equation for each variable is linear in that variable (i.e., linear if all other variables are considered constant). The idea is thus to solve the differential equation for each variable over one time step, assuming that all other variables are constant over that time step. The numerical scheme is still first order, but it is more stable than the forward Euler method. Each equation can be written as dx/dt=a*x+b, where a and b depend on the other variables and thus change after each time step. The values of a and b are obtained during the update phase by calculating a*x+b for x=0 and x=1 (note that these values are different for every neuron, thus we calculate vectors A and B). Then x(t+dt) is calculated in the same way as for the exact integration method above.
Stochastic differential equations¶
Noise is introduced in differential equations with the keyword
xi, which means normalised gaussian noise (the derivative of the Brownian term).
Currently, this is implemented simply by adding a normal random number to the variable
at the end of the integration step (independently for each neuron).
The unit of white noise is non-trivial, it is
second**(-.5). Thus, a typical stochastic
sigma is in the same units as
x. We note the following two facts:
The noise term is independent between neurons. Thus, one cannot use this method to analyse the response to frozen noise (where all neurons receive the same input noise). One would need to use an external variable representing the input, updated by a user-defined operation.
The noise term is independent between equations. This can however be solved by the following trick:
dx/dt=-x/tau+sigmax*u/tau**.5 : volt dy/dt=-y/tau+sigmay*u/tau**.5 : volt u=xi : second**(-.5)
It is not possible to modulate the noise term with a variable (e.g.
v*xi/tau**.5). One reason is that, with multiplicative noise, there is
an ambiguity between the Ito and the Stratonovich interpretation.
Unfortunately, this limitation also applies to parameters, i.e.,
sigma*xi/tau**.5 is not possible if sigma is a parameter, as in the
eqs=Equations('''dx/dt=-x/tau + sigma*xi/tau**.5 : volt sigma : volt''') group = NeuronGroup(1, eqs, threshold='x>vt')
However, the problem can usually be solved by some rewriting:
eqs=Equations('''dy/dt=-y/tau + xi/tau**.5 : 1 x=sigma*y : volt sigma : volt''') group = NeuronGroup(1, eqs, threshold = 'x>vt')
The time variable
t can be directly inserted into an equation string.
It is replaced at run time by the current value of the time variable for the relevant
neuron group, and also appears as a state variable of the neuron group.
External variables can be frozen by passing the keyword
False) at initialization of a
Then when the string defining the equations are compiled into Python functions
the external variables are replaced by their float values (units are discarded).
This can result in a significant speed-up.
TODO: more on the implementation.
State updates can be compiled into Python code objects by passing the keyword
compile=True at initialization of a a
Note that this is different from the method
which compiles the equation for every variable into a Python function
(not the whole state update process).
compile keyword is set, the method
exponential_euler_code() is called. It generates
a string containing the Python code for the update of all state variables (one time step),
then compiles it into Python code object. That compiled object is then called at every time step.
All external variables are frozen in the process (regardless of the value of the
This results in a significant speed-up (although the exponential Euler code is not
quite optimised yet). Note that only Python code is generated, thus a
C compiler is not required.
Working with equations¶
eqs=Equations(''' dx/dt=(y-x)/(10*ms) : volt dy/dt=-z/(5*ms) : volt z=2*(x+y) : volt ''')
Applying an equation¶
The value of z can be calculated using the
The second argument is a dictionary containing the values of all dependent variables
(here the result is
The right-hand side of differential equations can also be calculated in the same way:
Note in the second case that only the values of the dynamic variables should be passed.
Calculating a fixed point¶
A fixed point of the equations can be calculated as follows:
where the optional keywords give the initial point (zero if not provided).
Internally, the function
optimize.fsolve from the Scipy package is used to
find a zero of the set of differential equations (thus, convergence is not
guaranteed; in that case, the initial values are returned).
A dictionary with the values of the dynamic variables at the fixed point is returned.
If the equations were previously frozen, then the units disappear from the equations and unit consistency problems may arise.
Equationsobjects need to be “prepared” before use, as follows:
This is automatically called by the NeuronGroup initialiser.
|||Note that this approach raises an issue when dX/dt=B. We currently (temporarily) solve this problem by adding a small diagonal matrix to M to make it invertible.|