# 13.5 Networks of Nonlinear Integrate-and-Fire Neurons

In Chapter 5 it was shown that nonlinear integrate-and-fire neurons such as the exponential integrate-and-fire model provide an excellent approximation to the dynamics of single neurons, much better than the standard leaky integrate-and-fire model. This section indicates how the Fokker-Planck approach can be used to analyze networks of such nonlinear integrate-and-fire models.

We determine, for arbitrary nonlinear integrate-and-fire models driven by a diffusive input with constant mean $\mu=R\,I_{0}$ and noise $\sigma^{2}$, the distribution of membrane potentials $p_{0}(u)$ as well as the linear response of the population activity

 $A(t)=A_{0}+A_{1}(t)=A_{0}+\int_{0}^{\infty}G(s)\,I_{1}(t-s)\,{\text{d}}s,$ (13.37)

to a drive $\mu(t)=R\,[I_{0}+I_{1}(t)]$ (434). As explained in Section 13.4, the knowledge of $p_{0}$ is sufficient to predict, for coupled populations of integrate-and-fire neurons, the activity $A_{0}$ in the stationary state of asynchronous firing. Moreover, we will see that the filter $G(s)$ contains all the information need to analyze the stability of the stationary state of asynchronous firing. In this Section we focus on one-dimensional nonlinear integrate-and-fire models, and add adaptation later on in Section 13.6.

# 13.5.1 Steady state population activity

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

 $\qquad{\partial\over\partial t}p(u,t)=-{\partial\over\partial u}J(u,t)+A(t)\,% \delta(u-u_{r})-A(t)\,\delta(u-\theta_{\rm reset})\,.$ (13.38)

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

 ${\partial\over\partial u}J(u,t)=A(t)\,\delta(u-u_{r})-A(t)\,\delta(u-\theta_{% \rm reset})\,.$ (13.39)

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. For $u_{r} 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

 $J(u,t)={1\over\tau_{m}}\,[f(u)+\mu(t)-{1\over 2}\sigma^{2}(t){\partial\over% \partial u}]\,p(u,t)\,.$ (13.40)

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

 ${{\text{d}}p_{0}(u)\over{\text{d}}u}={2\tau_{m}\over\sigma^{2}}\left[{f(u)+\mu% \over\tau_{m}}p_{0}(u)-c\right]\,.$ (13.41)

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.

Example: Exponential integrate-and-fire model

The exponential integrate-and-fire model (156) as defined in Eq. (5.6) of Chapter 5 is characterized by the differential equation

 $\tau{{\text{d}}\over{\text{d}}t}u=-(u-u_{\rm rest})+\Delta_{T}\,\exp\left({u-% \vartheta_{rh}\over\Delta_{T}}\right)+\mu\,,$ (13.42)

where $\mu$ is the total drive in units of the membrane potential. The result of the numerical integration of Eqs. (13.39) and (13.41) depends critically on the relative value of the reset potential, the driving potential $\mu$ with respect to the reset $u_{r}$, and the rheobase firing threshold $\vartheta_{rh}$, as well as on the level $\sigma$ of diffusive noise. Two representative scenarios are shown in Fig. 13.8A. Case (i) has very little noise ($\sigma=2$mV) and the total drive ($\mu=-45$mV) is above the rheobase firing threshold ($\vartheta_{rh}=-53$mV). In this case, the membrane potential density a few mV below the reset potential ($u_{r}=-60$mV) is negligible. Because of the strong drive, the model neuron model is in the super-threshold regime (cf. Section 8.3 in Ch. 8) and fires regularly with a rate of about 44Hz.

Case (ii) is different, because it has a larger amount of noise ($\sigma=6$mV) and a negligible drive of $\mu=u_{r}=-60$mV. Therefore, the distribution of membrane potentials is broad, with a maximum at $u_{r}$ and firing is noise driven and occurs at a low rate of 5.6Hz.

For both noise levels, the frequency-current curve $A_{0}=\nu=g_{\sigma}(I_{0})$ is shown in Fig. 13.8B, as a function of the total the drive $\mu=R\,I$. We emphasize that for very strong super-threshold drive, the population firing rate $A_{0}$ of noisy neurons is lower than that of noise-free neurons.

# 13.5.2 Response to modulated input (*)

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

 $I(t)=I_{0}+\epsilon\,\cos(\omega\,t),$ (13.43)

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

 $A(t)=A_{0}+A_{1}(\omega)\,\cos(\omega\,t+\phi_{A}(\omega))\,.$ (13.44)

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

 $\hat{G}(\omega)={A_{1}(\omega)\over\epsilon}\,e^{i\phi_{A}(\omega)}\,.$ (13.45)

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

 $\Delta A(t)=\int_{0}^{\infty}G(s)\,\Delta I(t-s)\,{\text{d}}s\,.$ (13.46)

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

 $\Delta A(t)\ll A(t)\,.$ (13.47)

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

 $p(u,t)=p_{0}(u)+p_{1}(u)\,\cos(\omega\,t+\phi_{p}(u))$ (13.48)

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

 $J(u,t)=\,\left[{-u+R\,I(t)\over\tau_{m}}+{\Delta_{T}\over\tau_{m}}\exp\left({u% -\vartheta_{rh}\over\Delta_{T}}\right)-{\sigma^{2}(t)\over 2\tau_{m}}{\partial% \over\partial u}\right]\,p(u,t)=Q(u,t)\,p(u,t)\,,$ (13.49)

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

 $J(u,t)=J_{0}(u)+J_{1}(u)\cos(\omega t+\phi_{J}(u))$ (13.50)

with

 $J_{1}(u)\cos(\omega t+\phi_{J}(u))=Q_{0}(u)\,p_{1}(u)\cos(\omega t+\phi_{p}(u)% )+Q_{1}(u,t)\,p_{0}(u)+0_{\rm order}(\epsilon^{2}),$ (13.51)

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

 $\qquad-{\partial\over\partial u}\hat{J}_{1}(u)=i\omega\hat{p}_{1}(u)+\hat{A}_{% 1}\,[\delta(u-\theta_{\rm reset})-\delta(u-u_{r})]\,.$ (13.52)

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

 ${\partial\over\partial u}\hat{p}_{1}(u)={2\over\sigma^{2}}\,\left[-u+R\,I_{0}+% \Delta_{T}\exp\left({u-\vartheta_{rh}\over\Delta_{T}}\right)\right]\hat{p}_{1}% (u)+{2R\,\epsilon\over\sigma^{2}}\,p_{0}(u)-{2\tau\over\sigma^{2}}\,\hat{J}_{1% }(u).$ (13.53)

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

 $0=J(u_{\rm low})=a_{1}\,J_{1}^{\rm free}(u_{\rm low})+a_{2}\,J_{1}^{\epsilon}(% u_{\rm low})\,.$ (13.54)

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)

 $\hat{G}(\omega)={a_{1}\over\epsilon}=-{a_{2}\over\epsilon}\,{J_{1}^{\epsilon}(% u_{\rm low})\over J_{1}^{\rm free}(u_{\rm low})}\,.$ (13.55)

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