14 The Integral-equation Approach

14.5 Adaptation in Population Equations

Fig. 14.12: Population of adaptive neurons in response to a step stimulus. The results from quasi-renewal theory are compared to a direct simulation of 25000 neurons with parameters optimized for cortical cells. A. The solution of the integral equation Eqs. (14.86) and (14.96) for the population activity (thin black line) is overlaid on the population activity calculated by simulating 25 000 SRM neurons with parameters capturing properties of cortical cells (cf. Ch.11). The dotted line which gives the solution of the integral equation of non-adapting neurons (time-dependent renewal theory, Eq. (14.5)) indicates that adaptation makes a major contribution to the observed population activity. B. Comparing adapted activity AA_{\infty} from simulations (squares, error bars correspond to one standard deviation) to predictions from renewal (dotted line) and quasi-renewal (black line) theory. Modified from Naud and Gerstner (358).

The integral equations presented so far apply only to neurons whose state does not depend on the spiking history, except on the most recent spike. Cortical neurons, however, show pronounced adaptation, as discussed in Chapter 11. Since time-dependent renewal theory cannot account for adaptation, the predicted population activity does not match the one observed in simulations; see Fig. 14.12.

As indicated in Eq. (14.12), the intuitions that have led to the integral equation (14.5) of time-dependent renewal theory, are still applicable, but the interval distribution must take into account the past history.

Let PI,A(t|t^)P_{I,A}(t|\hat{t}) be, again, the probability of observing a spike at time tt given the last spike at time t^\hat{t}, the input-current and the activity history A(t)A(t) until time tt, then the integral equation for adapting neurons is Eq. (14.12) which we recall here for convenience

A(t)=-tPI,A(t|t^)A(t^)dt^.A(t)=\int_{-\infty}^{t}P_{I,A}(t|\hat{t})A(\hat{t})d\hat{t}. (14.86)

In this section, we develop the methods necessary to describe adapting neuron populations. First we present the systematic derivation of the integral equation Eq. (14.86). Then, in Section 14.5.2 we describe the event-based moment expansion which provides a framework for approximating PI,A(t|t^)P_{I,A}(t|\hat{t}).

14.5.1 Quasi-Renewal Theory (*)

Chapter 12 discussed the concept of a homogeneous population. We have seen that the activity of a homogeneous population of uncorrelated neurons can be seen as an average of the spike trains across the population, A(t)=S(t)NA(t)=\langle S(t)\rangle_{N} where the subscript NN denotes an average across the NN neurons. Because neurons are stochastic, all neurons have different spiking histories. Accordingly, an average over a population is equivalent to an average over all possible past spike trains,

A(t)=ρ(t|S)SA(t)=\langle\rho(t|S)\rangle_{S} (14.87)

where ρ(t|S)\rho(t|S) is the hazard (or instantaneous firing rate) at time tt given the previous spikes in the spike train S(t)S(t) and the average S\langle\cdot\rangle_{S} is over all spike train histories. In other words, A(t)A(t) can be calculated as the expected activity averaged over all possible spike train histories, whereas ρ(t|S)\rho(t|S) gives the stochastic intensity of generating a spike around time tt given the specific history summarized by the past spike train SS (358).

The spike train SS consists of a series of spikes t^1\hat{t}_{1}, t^2\hat{t}_{2}, … t^n\hat{t}_{n}. To average over all possible spike trains requires that we consider spike trains made of n=0n=0 to \infty spikes. Moreover, for a given number nn of spikes, all possible time sequences of the nn spikes must be considered. This average takes the form of a path integral (149) where the path here is the spike trains. Therefore, the average in Eq. (14.87) is

A(t)=n=0-t-t^1-t^n-1ρ(t|t^1,t^2,,t^n)P(t^1,t^2,)dt^ndt^1A(t)=\sum_{n=0}^{\infty}\int_{-\infty}^{t}\int_{-\infty}^{\hat{t}_{1}}...\int_% {-\infty}^{\hat{t}_{n-1}}\rho(t|\hat{t}_{1},\hat{t}_{2},...,\hat{t}_{n})P(\hat% {t}_{1},\hat{t}_{2},...){\text{d}}\hat{t}_{n}...{\text{d}}\hat{t}_{1} (14.88)

We emphasize that in Eq. (14.88) the spikes are ordered counting backward in time; therefore t^1\hat{t}_{1} is the most recent spike so that t>t^1>t^2>>t^nt>\hat{t}_{1}>\hat{t}_{2}>\dots>\hat{t}_{n}.

Each spike train can occur with a probability entirely defined by the hazard function

P(t^1,t^2,t^n)=[t^iSnρ(t^i|Sn)]e--tρ(x|Sn)dx.P(\hat{t}_{1},\hat{t}_{2},...\hat{t}_{n})=\left[\prod_{\hat{t}_{i}\in S_{n}}% \rho(\hat{t}_{i}|S_{n})\right]e^{-\int_{-\infty}^{t}\rho(x|S_{n})dx}. (14.89)

The equality follows from the product of two probabilities: the probabilities of observing the spikes at the times t^\hat{t} given the rest of the history (the hazard function) and the probability of not spiking between each of the spikes (the survivor function). Again, it is a generalization of the quantities seen in Chapters 7 and 9. Writing the complete path integral, we have

A(t)=ρ(t)e--tρ(x)dx+n=1-t-t^1-t^n-1[e-t^1tρ(x|t^1,,t^n)dxρ(t|t^1,,t^n)]×ρ(t^1|t^2,,t^n)ρ(t^2|t^3,,t^n)ρ(t^n)e-t^2t^1ρ(x|t^2,,t^n)dx---t^nρ(x)dxdt^ndt^2dt^1.\begin{split}A(t)=&\rho(t)e^{-\int_{-\infty}^{t}\rho(x)dx}+\sum_{n=1}^{\infty}% \int_{-\infty}^{t}\int_{-\infty}^{\hat{t}_{1}}...\int_{-\infty}^{\hat{t}_{n-1}% }\left[e^{-\int_{\hat{t}_{1}}^{t}\rho(x|\hat{t}_{1},...,\hat{t}_{n})dx}\rho(t|% \hat{t}_{1},...,\hat{t}_{n})\right]\\ &\times\rho(\hat{t}_{1}|\hat{t}_{2},...,\hat{t}_{n})\rho(\hat{t}_{2}|\hat{t}_{% 3},...,\hat{t}_{n})...\rho(\hat{t}_{n})e^{-\int_{\hat{t}_{2}}^{\hat{t}_{1}}% \rho(x|\hat{t}_{2},...,\hat{t}_{n})dx-...-\int_{-\infty}^{\hat{t}_{n}}\rho(x)% dx}{\text{d}}\hat{t}_{n}...{\text{d}}\hat{t}_{2}{\text{d}}\hat{t}_{1}.\end{split} (14.90)

The case with zero spikes in the past (the first term on the right-hand-side of Eq. (14.90)) can be neglected, because we can always formally assign a firing time at t^=-\hat{t}=-\infty to neurons that have not fired in the recent past. We now focus on the factor enclosed in square brackets on the right-hand side of Eq. (14.90). This factor resembles PI(t|t^)P_{I}(t|\hat{t}). In fact, if, for a moment, we made a renewal assumption, the most recent spike were the only one that mattered. This would imply that we could set ρ(t|t^1,,t^n)=ρ(t|t^)\rho(t|\hat{t}_{1},\dots,\hat{t}_{n})=\rho(t|\hat{t}) inside the square brackets (where t^=t^1\hat{t}=\hat{t}_{1} is the most recent spike) and shift the factor enclosed by square brackets in front of the n-1n-1 integral over t^2,t^3t^n\hat{t}_{2},\hat{t}_{3}\dots\hat{t}_{n}. In this case, the factor in square brackets would become exactly PI(t|t^)P_{I}(t|\hat{t}). The remaining integrals over the n-1n-1 variables can be recognized as the average of ρ(t^1|S)\rho(\hat{t}_{1}|S) over the possible histories, that is, A(t^1)A(\hat{t}_{1}). Therefore the renewal assumption leads back to Eq. (14.5), as it should.

In quasi-renewal theory, instead of assuming ρ(t|S)=ρ(t|t^1)\rho(t|S)=\rho(t|\hat{t}_{1}) we assume that

ρ(t|S)ρ(t|t^1,S)S,\langle\rho(t|S)\rangle\approx\langle\rho(t|\hat{t}_{1},S^{\prime})\rangle_{S^% {\prime}}, (14.91)

where SS^{\prime} is made of all the spikes in SS but t^1\hat{t}_{1}. This assumption is reasonable given the following two observations in real neurons. First, the strong effect of the most recent spike needs to be considered explicitly. Second, the rest of the spiking history only introduces a self-inhibition that is similar for all neurons in the population and that depends only on the expected distribution of spikes in the past (358). The approximation (14.91) is not appropriate for intrinsically bursting neurons, but should apply well to other cell types (fast-spiking, non-fast-spiking, adapting, delayed, …). If we insert Eq. (14.91) into the terms inside the square brackets of Eq. (14.90), we can define the interspike interval distribution

PI,A(t|t^)=ρ(t|t^,S)Sexp(-t^tρ(x|t^,S)Sdx),P_{I,A}(t|\hat{t})=\langle\rho(t|\hat{t},S^{\prime})\rangle_{S^{\prime}}\exp% \left(-\int_{\hat{t}}^{t}\langle\rho(x|\hat{t},S^{\prime})\rangle_{S^{\prime}}% dx\right), (14.92)

such that Eq. (14.90) becomes Eq. (14.86).

14.5.2 Event-based Moment Expansion (*)

The development in the previous subsection is general and does not rely on a specific neuron model nor on a specific noise model. In this section, we now introduce approximation methods that are valid for the broad class of Generalized Linear Models of spiking neurons with exponential escape noise, e.g., the Spike Response Model (see Chapter 9). The aim is to find theoretical expressions for PI,A(t|t^)P_{I,A}(t|\hat{t}) according to quasi-renewal theory (Eq. (14.92)). In particular we need to evaluate the average of the likelihood over the past history ρ(t|t^1,S)S\langle\rho(t|\hat{t}_{1},S^{\prime})\rangle_{S^{\prime}}.

We recall the SRM model with exponential escape noise (Chapter 9). The instantaneous stochastic intensity at time tt, given the spike history is

ρ(t|S)=ρ¯exp(h(t)+t^kSη(t-t^k))=ρ¯eh(t)+[ηS](t).\rho(t|S)=\overline{\rho}\exp\left(h(t)+\sum_{\hat{t}_{k}\in S}\eta(t-\hat{t}_% {k})\right)=\overline{\rho}e^{h(t)+[\eta\ast S](t)}. (14.93)

Here, h(t)=0κ(s)I(t-s)dsh(t)=\int_{0}^{\infty}\kappa(s)I(t-s)\,{\text{d}}s is the input potential, η(s)\eta(s) describes the spike-after potential introduced by each spike and u(t)=h(t)+t^kSη(t-t^k)u(t)=h(t)+\sum_{\hat{t}_{k}\in S}\eta(t-\hat{t}_{k}) is the total membrane potential. The standard formulation of exponential escape noise, ρ(t)=f(u(t)-ϑ)=ρ0exp[β(u(t)-ϑ)]\rho(t)=f(u(t)-\vartheta)=\rho_{0}\exp[\beta\,(u(t)-\vartheta)], contains two extra parameters (ϑ\vartheta for the threshold and 1/β1/\beta for the noise level), but we can rescale the voltage and rate units so as to include the parameters β\beta and ϑ\vartheta into the definition of uu and ρ¯\overline{\rho}, respectively. The time course of η\eta and κ\kappa and the parameter ρ¯\bar{\rho} can be optimized so as to describe the firing behavior of cortical neurons; cf. Chapter 11.

In the model defined by Eq. (14.93), all the history dependence is contained in a factor, eηSe^{\eta\ast S}, which can be factorized into the contribution from the last spike and that of all previous spikes,

ρ(t|I,t^,S)S=ρ¯eh(t)+η(t-t^)e[ηS](t)S.\langle\rho(t|I,\hat{t},S^{\prime})\rangle_{S^{\prime}}=\overline{\rho}e^{h(t)% +\eta(t-\hat{t})}\langle e^{[\eta\ast S^{\prime}](t)}\rangle_{S^{\prime}}. (14.94)

In order to evaluate the expected values on the right-hand side of Eq. (14.92), we need to calculate the quantity MS[η]=eηSSM_{S}[\eta]=\langle e^{\eta\ast S}\rangle_{S}.

MSM_{S} is called a moment generating functional because functional derivatives with respect to η(t)\eta(t) and evaluated at η(t)=0\eta(t)=0 yields the moments of SS: δMSδη[η=0]=S(t)\frac{\delta M_{S}}{\delta\eta}[\eta=0]=\langle S(t)\rangle, δ2MSδη2[η=0]=S(t)S(t)\frac{\delta^{2}M_{S}}{\delta\eta^{2}}[\eta=0]=\langle S(t)S(t^{\prime})\rangle, …   . Explicit formulas for the moment generating functional are known (529). One of the expansion schemes unfolds in terms of the correlation functions gk(t1,t2,,tk)g_{k}(t_{1},t_{2},...,t_{k}) and provides a useful framework for the approximation of our functional at hand.

We have already seen two examples of such correlation functions in various chapters of the book. The first term is simply the population activity g1(t1)=S(t1)=A(t1)g_{1}(t_{1})=\langle S(t_{1})\rangle=A(t_{1}), i.e., the expectation value or ‘mean’ of the spike count in a short interval Δt0\Delta t\to 0. The second term in the expansion scheme, is the second-order correlation function which we have encountered in Chapter 7 in the context of autocorrelation and the noise spectrum. There, the quantity Cii0(s)C_{ii}^{0}(s) is the stationary version of the slightly more general time-dependent correlation term g2(t1,t2)=S(t1)S(t2)-S(t1)S(t2)g_{2}(t_{1},t_{2})=\langle S(t_{1})S(t_{2})\rangle-\langle S(t_{1})\rangle% \langle S(t_{2})\rangle. Higher orders would follow the same pattern(529), but they will play no role in the following.

Using the moment expansion, the expected value in Eq. (14.94) becomes

e[ηS](t)S=exp(m=11m!-t^(eη(t-s1)-1)(eη(t-sm)-1)gm(s1,,sm)ds1dsm).\langle e^{[\eta\ast S^{\prime}](t)}\rangle_{S^{\prime}}=\exp\left(\sum_{m=1}^% {\infty}\frac{1}{m!}\int_{-\infty}^{\hat{t}}\left(e^{\eta(t-s_{1})}-1\right)..% .\left(e^{\eta(t-s_{m})}-1\right)g_{m}(s_{1},...,s_{m})ds_{1}...ds_{m}\right). (14.95)

We call the expansion the event-based moment expansion. As a physical rule of thumb, the contribution of terms m2m\geq 2 in Eq. (14.95) decreases rapidly with increasing mm, whenever the events (i.e., spike times at times tt and tt^{\prime}) are weakly coupled. For the purpose of Eq. (14.92), we can effectively ignore all second- and higher-order correlations and only keep the term with m=1m=1. This gives

ρ(t|t^,S)S=ρ¯eh(t)+η(t-t^)+-t^(eη(t-z)-1)A(z)dz.\langle\rho(t|\hat{t},S^{\prime})\rangle_{S^{\prime}}=\overline{\rho}e^{h(t)+% \eta(t-\hat{t})+\int_{-\infty}^{\hat{t}}(e^{\eta(t-z)}-1)A(z)dz}. (14.96)

Eq. (14.96) can be used as an excellent approximation in the survivor function as well as in formula for the interval distribution PI,A(t|t^)P_{I,A}(t|\hat{t}) in Eq. (14.92). Used within the integral equation (14.12) it gives an implicit description of the population activity for infinite an infinite population of adapting neurons.