We consider a population of nonlinear integrate-and-fire models in the diffusion limit. We start with the continuity equation (13.7)

In the stationary state, the membrane potential density does not depend on time, $p(u,t)=p_{0}(u)$, so that the left-hand side of (13.38) vanishes. Eq. (13.38) therefore simplifies to

Therefore the flux takes a constant value, except at two values: at the numerical threshold $\theta_{\rm reset}$, at which the membrane potential is reset; and at the potential $u_{r}$, to which the reset occurs. The flux vanishes for $u>\theta_{\rm reset}$ and for $u<u_{r}$. For $u_{r}<u<\theta_{\rm reset}$ the constant value $J(u,t)=c>0$ still needs to be determined.

A neuron in the population under consideration is driven by the action potentials emitted by presynaptic neurons from the same or other populations. For stochastic spike arrivals at rate $\nu_{k}$, where each spike causes a jump of the voltage by an amount $w_{k}$, the contributions to the flux have been determined in Eqs. (13.10) and (13.11). In the diffusion limit the flux according to Eq. (13.21) is

In the stationary state, we have $p(u,t)=p_{0}(u)$ and $J(u,t)=c$ for $u_{r}<u<\theta_{\rm reset}$. (For $u<u_{r}$ we have $J(u,t)=0$.) Hence Eq. (13.40) simplifies to a first-order differential equation

The differential equation (13.41) can be integrated numerically starting at the upper bound $u=\theta_{\rm reset}$ with initial condition $p_{0}(\theta_{\rm reset})=0$ and ${\text{d}}p_{0}/{\text{d}}u|_{\theta_{\rm reset}}=-1=-2c\tau_{m}/\sigma^{2}$. When the integration passes $u=u_{r}$ the constant switches from c to zero. The integration is stopped at a lower bound $u_{\rm low}$ when $p_{0}$ has approached zero. The exact value of the lower bound is of little importance.

At the end of the integration, the surface under the voltage distribution $\int_{u_{\rm low}}^{\theta_{\rm reset}}p_{0}(u){\text{d}}u=1$ which determines the constant $c$ and therefore the firing rate $A_{0}$ in the stationary state.

So far we restricted the discussion to a constant input $I_{0}$. We now add a small periodic perturbation

where $\omega=2\pi/T$; the parameter $T$ denotes the period of the modulation. We expect the periodic drive to lead to a small periodic change in the population activity

Our aim is to calculate the complex number $\hat{G}$ that characterizes the linear response or ‘gain’ at frequency $\omega$

The result for a population of exponential integrate-and-fire neurons is shown in Fig. 13.9.

Once we have found the gain $\hat{G}$ for arbitrary frequencies $\omega$, an inverse Fourier transform of $\hat{G}$ gives the linear filter $G(s)$. The linear response of the population activity $A(t)=A_{0}+\Delta A(t)$ to an arbitrary input current $I(t)=I_{0}+\Delta I(t)$ is

A linear response is a valid description if the change is small:

In our case of periodic stimulation, the amplitude of the input current is scaled by $\epsilon$. We are free to choose $\epsilon$ small enough to fulfill condition (13.47).

The small periodic drive at frequency $\omega$ leads to a small periodic change in the density of membrane potentials

which has a phase lead ($\phi_{p}>0$) or lag ($\phi_{p}<0$) with respect to the periodic drive. We assume that the modulation amplitude in the input current $\epsilon$ is small so that $p_{1}(u)\ll p_{0}(u)$ for most values of $u$. We say that the change is at most of ‘order $\epsilon$’.

Similarly to the membrane potential density, the flux $J(u,t)$ will also exhibit a small perturbation of order $\epsilon$. For the exponential integrate-and-fire model of Eq. (13.42) with $u_{\rm rest}=0$ the flux is (cf. Eq. (13.40))

where we have introduced on the right-hand side a linear operator $Q$ that comprises all the terms inside the square brackets. The stationary state, that we have analyzed in the previous subsection under the assumption of constant input $I(t)=I_{0}$ has a flux $J_{0}(u)=Q_{0}(u)\,p_{0}(u)$. In the presence of the periodic perturbation, the flux is

with

where $Q_{1}(u,t)=\epsilon\,R\,\cos(\omega\,t)/\tau_{m}$ is the change of the operator $Q$ to order $\epsilon$. Note that the flux through the threshold $\theta_{\rm reset}$ gives the periodic modulation of the population activity.

We insert our variables $A,p,J$ into the differential equation (13.38) and find that the stationary terms with subscript zero cancel each other, as it should be, and only the terms with subscript ‘1’ survive. To simplify the notation, it is convenient to switch from real to complex numbers and include the phase into the definition of $A_{1},p_{1}(u),J_{1}(u)$, e.g. $\hat{A}_{1}=A_{1}\,\exp(i\phi_{A})$; the hat indicates the complex number. If we take the Fourier transform over time, Eq. (13.38) becomes

Before we proceed further, let us have a closer look at Eqs. (13.51) and (13.52). We highlight three aspects (434). The first observation is that we have, quite arbitrarily, normalized the membrane potential density to an integral of unity. We could have chosen, with equal rights, a normalization to the total number $N$ of neurons in the population. Then the flux on the left-hand side of Eq. (13.52) would be enlarged by a factor of $N$, but so would the membrane potential density and the population activity, which appear on the right-hand side. We could also multiply both sides of the equation by any other (complex) number. As a consequence, we can, for example, try to solve Eq. (13.52) for a population activity modulation $\hat{A}_{1}=1$ and take care of the correct normalization only later.

The second observation is that the flux $J_{1}$ in Eq. (13.51) can be quite naturally separated into two components. The first contribution to the flux is proportional to the perturbation $p_{1}$ of the membrane potential density. The second component is caused by the direct action of the external current $Q_{1}(u,t)=\epsilon\,R\,\cos(\omega\,t)/\tau_{m}$. We will exploit this separability of the flux in the following.

The third observation is that, for the case of exponential integrate-and-fire neurons, the explicit expression for the operators $Q_{0}$ and $Q_{1}$ in can be inserted into Eq. (13.51) which yields a first-order differential equation

We therefore have two first-order differential equations (13.52) and (13.53) which are coupled to each other. In order to solve the two equations, we now exploit our second observation and split the flux into two components. We drop the hats on $J_{1},p_{1},A_{1}$ to lighten the notation.

(i) A ‘free’ component $J_{1}^{\rm free}(u)$ describes the flux that would occur in the absence of external drive ($\epsilon=0$), but in the presence of a periodically modulated population activity $A_{1}$. Intuitively, the distribution of the membrane potential density $p_{1}^{\rm free}(u)$ and the flux $J_{1}^{\rm free}$ must exhibit some periodic ‘breathing pattern’ to enable the periodic modulation of the flux through the threshold of strength $\hat{A}_{1}$. We find the ‘free’ component by integrating Eq. (13.53) with $\epsilon=0$ in parallel with Eq. (13.52) with parameter $A_{1}=1$ (i.e. the periodic flow through threshold is of unit strength). The integration starts at the initial condition $p_{1}(\theta_{\rm reset})=0$ and $J_{1}^{\rm free}(\theta_{\rm reset})=A_{1}=1$ and continues toward decreasing voltage values. Integration stops at a lower bound $u_{\rm low}$ which we place at an arbitrarily large negative value.

(ii) A ‘driven’ component $J_{1}^{\rm\epsilon}$ accounts for the fact that the periodic modulation is in fact caused by the input. We impose for the ‘driven’ component that the modulation of the flux through threshold vanishes: $A_{1}=0$. We can therefore find the ‘driven’ component by integrating Eq. (13.53) with a finite $\epsilon>0$ in parallel with Eq. (13.52) with parameter $A_{1}=0$ starting at the initial condition $p_{1}(\theta_{\rm reset})=0$ and $J_{1}^{\rm\epsilon}(\theta_{\rm reset})=A_{1}=0$ and integrating toward decreasing voltage values. As before, integration stops at a lower bound $u_{\rm low}$ which we can shift to arbitrarily large negative values.

Finally we add the two solutions together: Any combination of parameters $a_{1},a_{2}$ for the total flux $a_{1}\,J_{1}^{\rm free}(u)+a_{2}\,J_{1}^{\epsilon}(u)$ and the total density $a_{1}\,p_{1}^{\rm free}(u)+a_{2}\,p_{1}^{\epsilon}(u)$ will solve the pair of differential equations (13.53) and (13.52). Which one is the right combination?

The total flux must vanish at large negative values. Therefore we require a boundary condition

We recall that $J_{1}^{\epsilon}(u)$ is proportional to the drive $\epsilon$. Furthermore, the initial condition was chosen such that $J_{1}^{\rm free}(\theta_{\rm reset})$ yields a population response $A_{1}=1$ so that, in the combined solution, the factor $a_{1}$ is the population response!

We are interested in the ’gain factor’ at stimulation frequency $\omega$. We started this section by applying a periodic stimulus of strength $\epsilon$ and observing, at the same periodic frequency a modulation of strength $A_{1}$. The gain factor is defined as $\hat{G}(\omega)=A_{1}/\epsilon$. With the above arguments the gain factor is (434)

The gain factor $\hat{G}(\omega)$ has an amplitude $|\hat{G}(\omega)|$ and a phase $\phi_{G}$. Amplitude and phase of the gain factor of a population of exponential integrate-and-fire neurons are plotted in Fig. 13.9A.

The same numerical integration scheme that starts at the numerical threshold $\theta_{\rm reset}$ and integrates downward with appropriate initial conditions can also be applied to other, linear as well as nonlinear, neuron models. Furthermore, with the same scheme it is also possible to calculate the gain factor $G_{\sigma}$ in response to a modulation of the variance $\sigma^{2}(t)$ (Fig. 13.9B), or the response $G_{\vartheta}$ to periodic modulation of model parameters such as the rheobase threshold $\vartheta_{rh}$ (434). For earlier results on $G_{\sigma}$ see also (78; 298; 477).