More on equations

The Equations class is a central part of Brian, since models are generally specified with an Equations object. Here we explain advanced aspects of this class.

External variables

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 level (default=0). 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 name, e.g.:

eqs=Equations('dx/dt=-x/tau : volt',tau=10*ms,x=None)

This is useful to avoid conflicts in the names of hidden variables.

Issues

  • There can be problems if a variable with the same name as the variable of a differential equation exists in the namespace where the Equations object was defined.

Combining equations

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. When two 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 Threshold or Reset objects 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 Equations.get_Vm()). The first variable is, with decreasing priority :

  • v
  • V
  • vm
  • Vm
  • the first defined variable.

Numerical integration

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 a NeuronGroup. Possible values are linear, nonlinear, Euler, RK, exponential_Euler.

Exact integration

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 the function get_linear_equations() [1]. Third, the matrix A such that X(t+dt)=A*(X(t)-B)+B is calculated at initialisation of a specific state updater object, LinearStateUpdater, 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.

Euler integration

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 equation reads:

dx/dt=-x/tau+sigma*xi/tau**.5

where 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)
    

Important notice

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 following example:

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')

Non-autonomous equations

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.

Freezing

External variables can be frozen by passing the keyword freeze=True (default = False) at initialization of a NeuronGroup object. Then when the string defining the equations are compiled into Python functions (method compile_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.

Compilation

State updates can be compiled into Python code objects by passing the keyword compile=True at initialization of a a NeuronGroup. Note that this is different from the method compile_functions(), which compiles the equation for every variable into a Python function (not the whole state update process).

When the compile keyword is set, the method forward_euler_code() or 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 freeze keyword). 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

Equations object can also be used outside simulations. In the following, we suppose that an Equations object is defined as follows:

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 apply() method:

z=eqs.apply('z',dict(x=3*mV,y=5*mV))

The second argument is a dictionary containing the values of all dependent variables (here the result is 8*mV). The right-hand side of differential equations can also be calculated in the same way:

x=eqs.apply('x',dict(x=2*mV,y=3*mV))
y=eqs.apply('y',dict(x=2*mV,y=3*mV))

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:

fp=eqs.fixedpoint(x=2*mV,y=3*mV)

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.

Issues

  • If the equations were previously frozen, then the units disappear from the equations and unit consistency problems may arise.

  • Equations objects need to be “prepared” before use, as follows:

    eqs.prepare()
    

    This is automatically called by the NeuronGroup initialiser.

Footnotes

[1]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.