# 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.

## 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. |