## Abstract

We demonstrate the mathematical equivalence of two commonly used forms of firing rate model equations for neural networks. In addition, we show that what is commonly interpreted as the firing rate in one form of model may be better interpreted as a low-pass-filtered firing rate, and we point out a conductance-based firing rate model.

At least since the pioneering work of Wilson & Cowan (1972), it has been common to study neural circuit behavior using rate equations—equations that specify neural activities simply in terms of their rates of firing action potentials, as opposed to spiking models, in which the actual emissions of action potentials, or spikes, are modeled. Rate models can be derived as approximations to spiking models in a variety of ways (Wilson & Cowan, 1972; Mattia & Del Giudice, 2002; Shriki, Hansel, & Sompolinsky, 2003; Ermentrout, 1994; La Camera, Rauch, Luscher, Senn, & Fusi, 2004; Aviel and Gerstner, 2006; Ostojic & Brunel, 2011; reviewed in Ermentrout & Terman, 2010; Gerstner & Kistler, 2002; and Dayan & Abbott, 2001).

Two forms of rate model most commonly used to model neural circuits are the following, which we will refer to as the v-equation and r-equation respectively:
1
2
Here, v and r are each vectors representing neural activity, with each element representing the activity of one neuron in the modeled circuit. v is commonly thought of as representing voltage, while r is commonly thought of as representing firing rate (probability of spiking per unit time). f(x) is a nonlinear input-output function that acts element-by-element on the elements of x, that is, it has ith element (f(x))i = f(xi) for some nonlinear function of one variable f. f typically takes such forms as an exponential, a power law, or a sigmoid function, and f(vi) is typically regarded as a static nonlinearity converting the voltage of the ith cell vi to the cell's instantaneous firing rate. W is the matrix of synaptic weights between the neurons in the modeled circuit. and I are the vectors of external inputs to the neurons in the v or r networks, respectively, which may be time dependent. In the appendix, we illustrate a simple heuristic derivation of the v-equation, starting from the biophysical equation for the voltages v. Along the way, we also point to a conductance-based version of the rate equation.

When developing a rate model of a network, it can be unclear which form of equation to use or whether it makes a difference. Here we demonstrate that the choice between equations 1 and 2 makes no difference: the two models are mathematically equivalent, and so will display the same set of behaviors. It has been noted previously (Beer, 2006) that when I is constant and W is invertible, the two equations are equivalent under the relationship v = Wr + I, . We generalize this result to demonstrate the equivalence of the two equations when W is not invertible and inputs may be time dependent.

The v-equation is defined when we specify the input across time, , and the initial condition v(0); we will call the combination of these and equation 1 a v-model. The r-equation is defined when we specify I(t) and r(0); we will call the combination of these and equation 2 an r-model. We will show that any v-model can be mapped to an r-model and any r-model can be mapped to a v-model such that the solutions to equations 1 and 2 satisfy v = Wr + I.

As we will see, the inputs in equivalent models are related by , or . That is, I is a low-pass-filtered version of . Note that there is an equivalence class of I, parameterized by I(0), that all correspond to the same under this equivalence. We assume that the equivalence class has been specified, that is, has been specified (if I has been specified, can be found as ). Then a v-model is defined by specifying v(0), while an r-model is defined by specifying the set {r(0), I(0)}. If W is D × D, then v(0) is D-dimensional, while {r(0), I(0)} is 2D-dimensional, so we can guess that the map from r to v takes a D-dimensional space of r-models to a single v-model, and conversely the map from v to r takes a single v-model back to a D-dimensional space of r-models, and we will show that this is true.

We first show that if r evolves according to the r-equation, then Wr + I evolves according to the v-equation. Setting v = Wr + I, we find:
3
4
5
Therefore, if v evolves according to the v-equation and r evolves according to the r-equation and v(0) = Wr(0) + I(0), then, since the v-equation propagates Wr + I forward in time, v = Wr + I at all times t > 0. We will thus have established the desired equivalence if we can solve v(0) = Wr(0) + I(0) for any v-model, specified by v(0), or for any r-model, specified by {r(0), I(0)}.

Note that, as expected, a D-dimensional space of r-models converges on the same v-model. Since {r(0), I(0)} forms a 2D-dimensional space, which is constrained by the D-dimensional equation v(0) = Wr(0) + I(0), the D-dimensional subspace of r-models {r(0), I(0)} that satisfy this equation all converge on the same v-model.

To go from an r-model to a v-model is straightforward: we simply set v(0) = Wr(0) + I(0). To go from a v-model to an r-model, we first define some useful notation:1

• •.

is the null space of W, that is, the subspace of all vectors that W maps to 0. PN is the projection operator into .

• •.

is the subspace perpendicular to . This is the subspace spanned by the rows of W. PN is the projection operator into .

• •.

is the range of W, that is, the subspace of vectors that can be written Wx for some x. This is the subspace spanned by the columns of W. PR is the projection operator into .

• •.

is the subspace perpendicular to , also called the left null space. PR is the projection operator into .

For any vector x, we define xNPNx, xNPNx, xRPRx, xRPRx. We rely on the fact that x = xN + xN = xR + xR.

Given a v-model, the equation v(0) = Wr(0) + I(0) has a solution if and only if , which is true if and only if vR(0) − IR(0) = 0,2 so we must choose
9
Letting DR be the dimension of and DN the dimension of , the fundamental theorem of linear algebra states that DR + DN = D. So IR(0) has dimension DN. This leaves unspecified IR(0), which has dimension DR.
To solve for rN(0), we note that the equation v = Wr + I can equivalently be written v = WrN + I (because WrN = 0, so Wr = WrN). That is, knowledge of v specifies only rN. We define W−1 to be the Moore-Penrose pseudo-inverse of W. This is the matrix that gives the one-to-one mapping of into that inverts the one-to-one mapping of to induced by W, and that maps all vectors in to 0.3 The pseudo-inverse has the property that W−1W = PN while WW−1 = PR. Then we can solve for rN(0) as
10
This is a DR-dimensional equation for the 2DR-dimensional set of unknowns {rN(0), IR(0)}, so it determines DR of these parameters and leaves DR free. For example, it could be solved by freely choosing IR(0) and then setting rN(0) = W−1(vR(0) − IR(0)), or by freely choosing rN(0) and then setting IR(0) = vR(0) − WrN(0).

Equations 10 and 9 together ensure the equality v(0) = Wr(0) + I(0). Applying W to both sides of equation 10 yields vR(0) = WrN(0) + IR(0) = Wr(0) + IR(0). This states that the equality holds within the range of W; orthogonal to the range of W, we have PRWr = 0 and vR(0) = IR(0). Together, these yield v(0) = Wr(0) + I(0).

Finally, we can freely choose rN(0), which has no effect on the equation v(0) = Wr(0) + I(0). rN(0) has DN dimensions, so we have freely chosen DR + DN = D dimensions in finding an r-model that is equivalent to the v-model. That is, we have found a D-dimensional subspace of such r-models—those that satisfy v(0) = Wr(0) + I(0).

To summarize, we have established the equivalence between r-models and v-models. For each fixed choice of W, τ, and , an r-model is specified by {r(0), I(0)} and equation 2, while a v-model is specified by v(0) and equation 1. The equivalence is established by setting v(0) = Wr(0) + I(0), which yields a D-dimensional subspace of equivalent r-models for a given v-model. Under this equivalence, v obeys equation 1, r obeys equation 2, and the two are related at all times by v = Wr + I, with . To go from an r-model to its equivalent v-model, we simply set v(0) = Wr(0) + I(0). To go from a v-model to one of its equivalent r-models, we set IR(0) = vR(0), freely choose rN(0), and freely choose {rN(0), IR(0)} from the DR-dimensional subspace of such choices that satisfy rN(0) = W−1(vR(0) − IR(0)), where W−1 is the pseudoinverse of W.

Finally, note that equation 2 can be written . That is, if we regard v as a voltage and f(v) as a firing rate, as suggested by the “derivation” in the appendix, then r is a low-pass-filtered version of the firing rate, just as I is a low-pass-filtered version of the input .

## Appendix:  Simple “Derivation” of the v-Equation

As an example of an unsophisticated and heuristic derivation of these equations (more sophisticated derivations can be found in the references in the main text), the v-equation can be “derived” as follows. We start with the equation for the membrane voltage of the ith neuron:
A.1
where Ci is the capacitance of the ith neuron and gij is the jth conductance onto the neuron, with reversal potential Eij. We assume that the gij’s are composed of an intrinsic conductance, gLi, with reversal potential ELi; extrinsic input with reversal potential ; and within-network synaptic conductances, with representing input from neuron j with reversal potential . Dividing by ∑kgik and defining τi(t) = Ci/∑kgik gives
A.2
We now make a number of further simplifying assumptions. We assume that is proportional to the firing rate rj of neuron j, with proportionality constant : . This ignores synaptic time courses, among other things. We assume that rj is given by the static nonlinearity rj = f(vj) (see Miller & Troyer, 2002; Hansel & van Vreeswijk, 2002; Priebe, Mechler, Carandini, & Ferster, 2004, for such a relationship between firing rate and voltage averaged over a few tens of milliseconds). We assume synapses are either excitatory with reversal potential EE or inhibitory with reversal potential EI, and linearly transform the units of voltage so that EE = 1 and EI = −1. We define . This is now a synaptic weight that is positive for excitatory synapses and negative for inhibitory synapses. We define and define . This yields the conductance-based rate equation,
A.3
with τi(t) = Ci/(gi + ∑k|Wik|f(vk)).
Finally, we assume that the total conductance, represented by the denominator in the last term of equation A.3, can be taken to be constant, for example, if gLi is much larger than synaptic and external conductances or if inputs tend to be push-pull, with withdrawal of some inputs compensating for addition of others. We absorb the constant denominator into the definitions of and Wij and note that this also implies that τi is constant, to arrive finally at the v-equation:
A.4

## Acknowledgments

This work was supported by R01-EY11001 from the National Eye Institute and by the Gatsby Charitable Foundation through the Gatsby Initiative in Brain Circuitry at Columbia University.

## Notes

1

If W is normal, the eigenvectors are orthogonal, so the null space is precisely the space orthogonal to the range: PN = PR and PN = PR. However, if W is nonnormal, then vectors orthogonal to the null space can be mapped into the null space; the range always has the dimension of the full space minus the dimension of the null space, but it need not be orthogonal to the null space.

2
Note that the condition , meaning that v = Wr + I can be solved, is true for all time if it is true in the initial condition. We compute:
6
7
Applying PR to equation 7 and noting that PRW = 0, we find
8
If , then vR(0) − IR(0) = 0, and hence vRIR = 0 at all subsequent times so at all subsequent times. Note also that for any initial conditions, the condition is true asymptotically as t → ∞.
3

If the singular value decomposition of a matrix M is M = USV†, where S is the diagonal matrix of singular values and U and v are unitary matrices, then its pseudo-inverse is , where is the pseudoinverse of S, obtained by inverting all nonzero singular values in S.

## References

Aviel
,
Y.
, &
Gerstner
,
W.
(
2006
).
From spiking neurons to rate models: A cascade model as an approximation to spiking neuron models with refractoriness
.
Phys. Rev. E
,
73
,
051908
.
Beer
,
R. D.
(
2006
).
Parameter space structure of continuous-time recurrent neural networks
.
Neural Comput.
,
18
,
3009
3051
.
Dayan
,
P.
, &
Abbott
,
L. F.
(
2001
).
Theoretical neuroscience
.
Cambridge, MA
:
MIT Press
.
Ermentrout
,
B.
(
1994
).
Reduction of conductance based models with slow synapses to neural nets
.
Neural Comput.
,
6
,
679
695
.
Ermentrout
,
G. B.
, &
Terman
,
D. H.
(
2010
).
Mathematical foundations of neuroscience
.
New York
:
Springer
.
Gerstner
,
W.
, &
Kistler
,
W.
(
2002
).
Spiking neuron models
.
Cambridge
:
Cambridge University Press
.
Hansel
,
D.
, &
van Vreeswijk
,
C.
(
2002
).
How noise contributes to contrast invariance of orientation tuning in cat visual cortex
.
J. Neurosci.
,
22
,
5118
5128
.
La Camera
,
G.
,
Rauch
,
A.
,
Luscher
,
H. R.
,
Senn
,
W.
, &
Fusi
,
S.
(
2004
).
Minimal models of adapted neuronal response to in vivo–like input currents
.
Neural Comput.
,
16
,
2101
2124
.
Mattia
,
M.
, &
Del Giudice
,
P.
(
2002
).
Population dynamics of interacting spiking neurons
.
Phys. Rev. E
,
66
,
051917
.
Miller
,
K. D.
, &
Troyer
,
T. W.
(
2002
).
Neural noise can explain expansive, power-law nonlinearities in neural response functions
.
J. Neurophysiol.
,
87
,
653
659
.
Ostojic
,
S.
, &
Brunel
,
N.
(
2011
).
From spiking neuron models to linear-nonlinear models
.
PLoS Comput. Biol.
,
7
,
e1001056
.
Priebe
,
N.
,
Mechler
,
F.
,
Carandini
,
M.
, &
Ferster
,
D.
(
2004
).
The contribution of spike threshold to the dichotomy of cortical simple and complex cells
.
Nat. Neurosci.
,
7
(
10
),
1113
1122
.
Shriki
,
O.
, &
Hansel
,
D.
, &
Sompolinsky
,
H.
(
2003
).
Rate models for conductance-based cortical neuronal networks
.
Neural Comput.
,
15
,
1809
1841
.
Wilson
,
H. R.
, &
Cowan
,
J. D.
(
1972
).
Excitatory and inhibitory interactions in localized populations of model neurons
.
Biol. Cybern.
,
12
,
1
24
.