Perfect spike detection via time reversal

TL;DR

This paper introduces a novel method for perfect retrospective spike detection in neuronal simulations, improving accuracy by propagating thresholds with time-inverted dynamics, applicable to linear neuron models.

Contribution

It develops and benchmarks a new geometric method for exact spike detection that surpasses traditional approximate methods in accuracy and efficiency.

Findings

The method guarantees perfect spike detection within linear neuron models.

It is faster than root-finding based perfect detection methods.

Missed spikes are extremely rare with existing approximate methods.

Abstract

Spiking neuronal networks are usually simulated with three main simulation schemes: the classical time-driven and event-driven schemes, and the more recent hybrid scheme. All three schemes evolve the state of a neuron through a series of checkpoints: equally spaced in the first scheme and determined neuron-wise by spike events in the latter two. The time-driven and the hybrid scheme determine whether the membrane potential of a neuron crosses a threshold at the end of of the time interval between consecutive checkpoints. Threshold crossing can, however, occur within the interval even if this test is negative. Spikes can therefore be missed. The present work derives, implements, and benchmarks a method for perfect retrospective spike detection. This method can be applied to neuron models with affine or linear subthreshold dynamics. The idea behind the method is to propagate the threshold with a time-inverted dynamics, testing whether the threshold crosses the neuron state to be evolved, rather than vice versa. Algebraically this translates into a set of inequalities necessary and sufficient for threshold crossing. This test is slower than the imperfect one, but faster than an alternative perfect tests based on bisection or root-finding methods. Comparison confirms earlier results that the imperfect test rarely misses spikes (less than a fraction $1/10^8$ of missed spikes) in biologically relevant settings. This study offers an alternative geometric point of view on neuronal dynamics.

Figures

Equations

$$V(t+h)⩾\theta ⟹\text{threshold crossing}$$

$$V(\mathbf{s})=\theta .$$

$$\begin{array}{c}{\mathbf{k}}^{⊺}=(0^{⊺},1),\kappa =\theta ,\\ {\mathbf{k}}^{⊺}\mathbf{s}=\kappa \text{on threshold,}{\mathbf{k}}^{⊺}\mathbf{s}<\kappa \text{below threshold.}\end{array}$$

$$\mathbf{s}\mapsto \text{M}\mathbf{s}+\mathbf{m},$$

$${{\mathbf{k}}^{\prime }}^{⊺}={\mathbf{k}}^{⊺}{\text{M}}^{-1},{\kappa }^{\prime }=\kappa +{\mathbf{k}}^{⊺}{\text{M}}^{-1}\mathbf{m}$$

$$\dot{\mathbf{s}}(t)=\text{A}\mathbf{s}(t)+\mathbf{q}.$$

$$\bm{s}(t+h)=\mathrm{e}^{h\textbf{{{A}}}}\bm{s}(t)+\bigl{(}\mathrm{e}^{h\textbf{{{A}}}}-1\bigr{)}\textbf{{{A}}}^{-1}\bm{q},\quad\bm{s}(t)=\bm{s}_{0}.$$

$$\text{A}=\left(\begin{array}{cc}\text{B}& \mathbf{d}\\ {\mathbf{c}}^{⊺}& \alpha \end{array}\right),\mathbf{q}=\left(\begin{array}{c}\mathbf{r}\\ \beta \end{array}\right),$$

$$V\bigl{[}\mathrm{e}^{t\textbf{{{A}}}}\bm{s}(t)+\bigl{(}\mathrm{e}^{t\textbf{{{A}}}}-1\bigr{)}\textbf{{{A}}}^{-1}\bm{q}\bigr{]}=\theta,\qquad t\in]0,h],$$

$$\mathbf{x}\mapsto (\mathbf{I},V)=(\mathbf{x},\theta ).$$

$$E\colon(\bm{x},t)\mapsto\mathrm{e}^{-t\textbf{{{A}}}}\begin{pmatrix}\bm{x}\\ \theta\end{pmatrix}+\bigl{(}\mathrm{e}^{-t\textbf{{{A}}}}-1\bigr{)}\textbf{{{A}}}^{-1}\begin{pmatrix}\bm{r}\\ \beta\end{pmatrix},$$

$$\begin{split}\bm{k}^{\intercal}_{t}&=\begin{pmatrix}\bm{0}^{\intercal}&1\end{pmatrix}\mathrm{e}^{t\textbf{{{A}}}},\\ \kappa_{t}&=\theta+\begin{pmatrix}\bm{0}^{\intercal}&1\end{pmatrix}\bigl{(}1-\mathrm{e}^{t\textbf{{{A}}}}\bigr{)}\textbf{{{A}}}^{-1}\begin{pmatrix}\bm{r}\\ \beta\end{pmatrix}.\end{split}$$

$$V=\theta ,$$

$$\displaystyle\begin{pmatrix}\bm{0}^{\intercal}&1\end{pmatrix}\mathrm{e}^{h\textbf{{{A}}}}\begin{pmatrix}\bm{I}\\ V\end{pmatrix}=\theta+\begin{pmatrix}\bm{0}^{\intercal}&1\end{pmatrix}\bigl{(}1-\mathrm{e}^{h\textbf{{{A}}}}\bigr{)}\textbf{{{A}}}^{-1}\begin{pmatrix}\bm{r}\\ \beta\end{pmatrix}.$$

$$\mathrm{T}E(\mathbf{x},t)=({\text{de}}_{\mathbf{x}}E,{\text{de}}_{t}E)=-{\mathrm{e}}^{-t\text{A}}\left(\begin{array}{cc}\text{I}& \text{B}\mathbf{x}+\mathbf{d}\theta +\mathbf{r}\\ 0^{⊺}& {\mathbf{c}}^{⊺}\mathbf{x}+\alpha \theta +\mu \end{array}\right).$$

$$\det \mathrm{T}E(\mathbf{x},t)=0⟺{\mathbf{c}}^{⊺}\mathbf{x}+\alpha \theta +\beta =0$$

$${\mathbf{c}}^{⊺}\mathbf{x}+\alpha \theta +\beta =0,0⩽t⩽h.$$

$$\begin{array}{c}{x}_1=y_1,\dots ,x_{N-2}=y_{N-2},\\ x_{N-1}=-\frac{c_1{y}_1+\cdots +c_{N-2}{y}_{N-2}+\beta +\alpha \theta }{c_{N-1}}.\end{array}$$

$${\mathbf{c}}^{⊺}{\text{de}}_{\mathbf{y}}\mathbf{x}=0^{⊺},$$

$$\bm{\varGamma}\colon(\bm{y},t)\mapsto\mathrm{e}^{-t\textbf{{{A}}}}\begin{pmatrix}\bm{x}(\bm{y})\\ \theta\end{pmatrix}+\bigl{(}\mathrm{e}^{-t\textbf{{{A}}}}-1\bigr{)}\textbf{{{A}}}^{-1}\begin{pmatrix}\bm{r}\\ \beta\end{pmatrix},\quad 0<t<h.$$

$$\{\mathbf{s}|\chi (\mathbf{s})=0,\mathbf{s}\in D_{\chi }\}=\{\Gamma (\mathbf{y},t)|\mathbf{y}\in {\mathbf{R}}^{N-2},0⩽t⩽h\}.$$

$${\mathbf{k}}^{⊺}(\text{M}\mathbf{s}+\mathbf{m})<\kappa \iff ({\mathbf{k}}^{⊺}\text{M})\mathbf{s}<(\kappa -{\mathbf{k}}^{⊺}\mathbf{m})$$

$${\mathbf{s}}_0=(\mathbf{I},V)\in \text{no-spike region}\iff$$

$$V<\theta$$

$$\displaystyle\!\begin{pmatrix}\bm{0}^{\intercal}&1\end{pmatrix}\mathrm{e}^{h\textbf{{{A}}}}\begin{pmatrix}\bm{I}\\ V\end{pmatrix}<\theta+\begin{pmatrix}\bm{0}^{\intercal}&1\end{pmatrix}\bigl{(}1-\mathrm{e}^{h\textbf{{{A}}}}\bigr{)}\textbf{{{A}}}^{-1}\begin{pmatrix}\bm{r}\\ \beta\end{pmatrix}$$

$$\chi (\mathbf{I},V)<0,(\mathbf{I},V)\in D_{\chi }.$$

$$\mathbf{s}=E(\mathbf{x},t)=E({\mathbf{x}}^{\prime },t^{\prime }),t\ne t^{\prime }.$$

$$V[\lambda {\mathbf{s}}_1(t)+(1-\lambda ){\mathbf{s}}_2(t)]<\theta ,0<\lambda <1,$$

$$\begin{array}{rl}\dot{I}(t)& =-\frac{1}{{\tau }_{\text{s}}}I(t),\\ \dot{V}(t)& =\frac{1}{C}I(t)-\frac{1}{\tau }V(t)+\frac{1}{C}{I}_{\text{e}},\end{array}$$

$$\frac{\mathrm{d}}{\mathrm{d}t}\left(\begin{array}{c}I\\ V\end{array}\right)$$

$$=\left(\begin{array}{c}-\frac{1}{{\tau }_{\text{s}}}\end{array}\right),$$

$${\mathbf{c}}^{⊺}$$

$$\exp(-t\textbf{{{A}}})=\begin{pmatrix}\mathrm{e}^{\frac{t}{\tau_{\text{s}}}}&0\\ \frac{\bigl{(}\mathrm{e}^{\frac{t}{\tau}}-\mathrm{e}^{\frac{t}{\tau_{\text{s}}}}\bigr{)}\tau\tau_{\text{s}}}{C(\tau-\tau_{\text{s}})}&\mathrm{e}^{\frac{t}{\tau}}\end{pmatrix},$$

$${\mathbf{k}}_t^{⊺}$$

Full text

1

Spiking neuronal networks are usually simulated with three main simulation schemes: the classical time-driven and event-driven schemes, and the more recent hybrid scheme. All three schemes evolve the state of a neuron through a series of checkpoints: equally spaced in the first scheme and determined neuron-wise by spike events in the latter two. The time-driven and the hybrid scheme determine whether the membrane potential of a neuron crosses a threshold at the end of of the time interval between consecutive checkpoints. Threshold crossing can, however, occur within the interval even if this test is negative. Spikes can therefore be missed.

The present work derives, implements, and benchmarks a method for perfect retrospective spike detection. This method can be applied to neuron models with affine or linear subthreshold dynamics. The idea behind the method is to propagate the threshold with a time-inverted dynamics, testing whether the threshold crosses the neuron state to be evolved, rather than vice versa. Algebraically this translates into a set of inequalities necessary and sufficient for threshold crossing. This test is slower than the imperfect one, but faster than an alternative perfect tests based on bisection or root-finding methods. Comparison confirms earlier results that the imperfect test rarely misses spikes (less than a fraction $1/10^{8}$ of missed spikes) in biologically relevant settings.

This study offers an alternative geometric point of view on neuronal dynamics.

\helveticabold

2 Keywords:

State-space analysis, NEST, time-driven, event-driven, simulation, LIF neuron, differential geometry

3 Introduction

In the last decade, considerable work has been devoted to improve the accuracy of simulators that are capable of efficiently simulating large networks of spiking neurons (Hansel et al., 1998; Mattia and Del Giudice, 2000; Shelley and Tao, 2001; Dehaene and Changeux, 2005; Morrison et al., 2007; Brette, 2007; D’Haene et al., 2009; van Elburg and van Ooyen, 2009; Zheng et al., 2009; Hanuschkin et al., 2010). The field is driven by the ideal of combining the capability to cope with the high-frequency of synaptic events arriving at a neuron in nature with a mathematically accurate implementation of the threshold process a wide class of neuron models is based on.

Two classical schemes to simulate neuronal networks are the time-driven and the event-driven schemes (Fujimoto, 2000; Zeigler et al., 2000; Ferscha, 1996). Both schemes describe the state of the neurons by a set of variables and the action potentials as events that mediate the interaction between them.

In a time-driven scheme, the state of a neuron is updated on a time grid defined by the simulation step (for a review see Morrison and Diesmann, 2008). After all neurons are updated, their membrane potential is checked for threshold crossings. If the membrane potential of a neuron is above the threshold at this checkpoint, a spike is delivered to all neurons it is connected to. Subsequently, a new iteration step begins. The step size stipulates how frequently occurrences of threshold crossings are inspected during the simulation. The choice of the step size a trade-off between spike accuracy and the speed of the simulation Morrison et al. (2007). Such grid-constrained simulations force each spike event to a position on the equidistant temporal grid spanned by the step size and therefore induce artificial synchronization of the network dynamics (Hansel et al., 1998; Shelley and Tao, 2001; Morrison et al., 2007; Brette, 2007; van Elburg and van Ooyen, 2009; Hanuschkin et al., 2010).

In an event-driven scheme, the state of a neuron is updated only when it receives a spike. A central queue of events is maintained and each spike is inserted into this queue with its own time stamp. Upon update a neuron predicts when its next spike will occur in the absence of further input. This preliminary event is inserted into the queue and confirmed if it becomes due or removed when invalidated by further input. Efficient and elegant predicition methods have been developed for classes of neuron models without invertable dynamics (Brette, 2007; van Elburg and van Ooyen, 2009; D’Haene et al., 2009; D’Haene and Schrauwen, 2010; Ferscha, 1996). However, maintaining a central queue in a distributed simulation is challenging and may compromise the time performance of the simulator (Hanuschkin et al., 2010, for a detailed review see).

A hybrid scheme circumvents the shortcomings of both these schemes by embedding a locally event-driven algorithm for each neuron into a globally time-driven scheme (Morrison et al., 2007). The arrival of a spike at a neuron introduces one additional update and check point. The dynamics of a given neuron is then propagated from incoming spike to incoming spike and eventually to the end point of the global timestep. If the membrane potential of a neuron is above the threshold at a local or global check point, the precise point of threshold crossing is determined in continuous time, and a spike is emitted. Next to their location on the time grid, in this scheme spike events carry a floating point offset. Thus, in contrast to an event-driven scheme, the hybrid scheme does not predict future spike times but identifies threshold crossings only retrospectively. Hanuschkin et al. (2010) demonstrate that the latter scheme is equally accurate as the former at lower computational costs.

The hybrid scheme still has a loophole, however: spikes can be missed. The reason is that, as in the time-driven scheme, crossing of the threshold voltage is tested by inspecting whether the membrane potential of the neuron is above threshold at the end of a checkpoint (see 1 for an illustration of the scenario).

Nevertheless, the membrane potential $V$, evolved by the subthreshold dynamics, can be below threshold $\theta$ at two consecutive checkpoints $t$ and $t+h$ while a double, quadruple, etc. threshold-crossing occurred in between. The first crossing constitutes a missed spike. Symbolically,

[TABLE]

so the test $V(t+h)\geqslant\theta$ is a sufficient but not a necessary condition for the occurrence of a threshold crossing during the time interval $]t,t+h]$. For future reference we call this test “standard test”.

The conceptual question therefore remains whether a globally time-driven scheme can formulated such that it detects every threshold crossing. Although Hanuschkin et al. (2010) argue that the loss of spikes of the standard test is not of practical relevance in natural parameter regimes, the availability of a method perfect by construction would free the researcher from inquietude and costly controls when faced with previously unexplored neuron models or network architectures.

In this work we propose a new spike-detection method, which we term “lossless method” or “lossless test” because it is a necessary and sufficient condition for threshold crossing to occur in a given time interval. The method is based on state-space analysis and works with any neuron model with affine or linear subthreshold dynamics. It consists of a system of inequalities – some linear, some non-linear in the state-space variables – that together determine whether the initial state of a neuron will or will not reach threshold within the time interval until the next checkpoint. The lossless method replaces the standard test (1) in the time-driven and the hybrid scheme. Alone, it does not solve the problem of artificial synchronization the time-driven scheme suffers from. Hence the method is most meaningful within the hybrid scheme. Thanks to its perfect spike detection the lossless method can in fact be used to benchmark the hybrid scheme based on the standard test; Hanuschkin et al. (2010) use the method of D’Haene et al. (2009) for this purpose.

In 4 we present the idea behind the lossless method for a general neuron model with an affine or linear dynamics, and develop its mathematical construction. Parts of this construction must be addressed on a case-by-case basis; therefore in 5 we provide a concrete implementation of the lossless method for the leaky integrate-and-fire model with exponential synaptic currents (Fourcaud and Brunel, 2002a), within the hybrid scheme. The method can be algorithmically expressed in different ways. We explore two alternative cascades of inequalities and assess their costs in terms of time-to-completion relative to each other and to the hybrid scheme based on the standard test (1). For the latter scheme we also assess the number of missed spikes in commonly considered network regimes. The hybrid scheme based on the lossless method delivers the desired exact implementation of the mathematical definition of the neuron model without any further approximation up to floating point precision.

Preliminary results have been published in abstract form (Kunkel et al., 2011; Krishnan et al., 2016). The technology described in the present article will be made available with one of next major releases of the open-source simulation software NEST. The conceptual and algorithmic work described here is a module in the long-term collaborative project to provide the technology for neural systems simulations (Gewaltig and Diesmann, 2007).

4 A time-reversed state-space analysis

4.1 Idea: moving a surface backwards instead of a point forward

Let us summarize the problem mentioned in the previous section. We assume that a neuron’s state evolves according to three different dynamics:

(a) an integrable subthreshold dynamic as long as the neuron’s membrane potential is below threshold and there are no changes in input currents;

(b) discrete jumps in the subthreshold dynamics at predetermined times, corresponding to incoming spikes or to sudden changes in external currents; these can be formally incorporated into the subthreshold dynamics (a) via delta functions; and

(c) a “spike”, i.e. an instantaneous jump of the membrane potential from threshold to a reset value, as soon as the potential reaches the threshold value. The jump may be followed by a refractory period in which the membrane potential remains constant at the reset value. Then the integrable dynamics takes place again.

The advantage of the integrable dynamics is that the state of the neuron at a time $t+h$ can be analytically determined by that at time $t$; here $h$ can be negative or positive. The evolution can thus be calculated in discrete time steps, in particular in between times at which jumps (b) occur. The spike part of the dynamics, however, forces us to check whether the membrane potential $V$ reached a threshold value $\theta$ within the timestep interval $]t,t+h]$. We call this event threshold crossing (by “crossing” we also mean tangency). A sufficient condition for threshold crossing is that the membrane potential be above threshold at the end $t+h$ of the time step: by continuity, it must have assumed the threshold value at some time in the interval $]t,t+h]$. But this condition is not necessary: during the time step the potential may touch or surpass the threshold value and then go below it again, as in 1, an even number of times. Its value is then below threshold at both ends of the time step, $t$ and $t+h$. A test that only relies on the sufficient condition $V(t+h)\geqslant\theta$ – the standard test (1) – can therefore miss some spikes, leading to an incorrect dynamics. We need a test based on a necessary and sufficient condition.

A necessary and sufficient condition for threshold crossing is that the trajectory of the state during $]t,t+h]$ intersect the hypersurface “$\text{membrane potential}=\text{threshold value}$”. Translated into analytic geometry this means finding the solutions of a system of parametric equations – one representing the threshold hypersurface, the other the trajectory – and to test whether its solution set is empty (no threshold crossing) or not (threshold crossing). This idea is illustrated in 2A for a two-dimensional state space. This system is usually transcendental and its solutions have to be found numerically. Unfortunately, numerical solutions typically rely on bisection algorithms (Press et al., 2007, ch. 9), involving an increasingly finer timestepping of the dynamics. This nullifies the advantage of having an integrable dynamics with coarse-grained time steps.

The problem is that the test of intersection between trajectory and threshold tells us not only whether a threshold crossing occurs, but also the time at which it does. The test’s high computational cost partly comes from delivering this additional information. This problem is avoided if we formulate a different geometric test that tells us whether a threshold crossing occurs but does not deliver the crossing time. Mathematics often offer non-constructive proofs: “there exists a solution to the problem, but we do not know what the solution is”. It turns out that this ignorance is exactly what we need in our problem.

Instead of evolving the state of the neuron forwards in time, tracing a trajectory in state-space, and checking if and when it crosses the threshold hypersurface, we can evolve the threshold hypersurface backwards in time, sweeping a hypervolume, and check if and when it “crosses” the initial neuron state, which is a point. In other words we are testing whether a point belongs to a particular state-space region. The test for the intersection of a 1-dimensional curve with an $(N-1)$-dimensional surface is replaced by the test for the membership of a point in an $N$-dimensional volume. This idea is illustrated in 2B for a two-dimensional state-space.

Mathematically the latter test translates into a system of inequalities that the initial state at $t$ must satisfy if it does not cross threshold within $]t,t+h]$. Even if transcendental functions appear in this system, we do not need to find their roots: we only need to test whether the inequalities are satisfied by simply inserting the value of the state-space variables in the functions.

The equations corresponding to these inequalities represent the piecewise-differentiable boundary between the set of states that will cross the threshold within the timestep $h$ – which we call spike region – and the set of those that won’t – which we call no-spike region. Finding these inequalities in explicit form is the most important point of this method, and can be achieved with this heuristic procedure:

1. I.

Find the hypervolume swept by the threshold, in parametric form. This is done by representing the backpropagation of the threshold in state-space as a map between two manifolds: the product manifold threshold${}\times{}$time, and the state-space manifold. 2. II.

Find the boundary of the hypervolume, in parametric form. This is done by determining the placement of the images of the boundaries of the product manifold and, most important, of the images of the critical points of the map. The latter are the points at which the map becomes singular, thus mapping hypervolume elements into hypersurface elements (Spivak 1999, ch. 2; Choquet-Bruhat et al. 1996, § II.A.1). 3. III.

Transform the equations of the boundary above from parametric to implicit form. If this is not possible, the boundary can still be approximated by a triangularization, as finely as we please. 4. IV.

Finally, exclude the parts of the boundary that lie in the interior of the hypervolume (so technically not a boundary). Some guidelines to achieve this can be given based on the the convex properties of the hypervolume.

In the following analysis we mathematically develop the first two steps in a general way for any neuron model with an affine dynamic. Let us point out that they can be generalized to other kinds of dynamics. The procedure for the last two steps depends on the particular neuron model, so we can only give general guidelines. 5 provides a concrete example of all steps. An advantage of this procedure is that it can be geometrically explained and needs to be carried out only once for any given neuron model.

The advantage of this novel test is that its computational cost is distributed differently from the approach where the trajectory of the initial point is followed. In hand-waving terms, in the trajectory-based test the crossing time is immediately available from the parametrization of the orbit; the threshold-crossing condition is represented by computationally costly inequalities, because the crossing time, when it exists, can be read from them immediately. In the novel test the crossing time is a more complicated function of the hypervolume parametrization; the threshold-crossing condition in this case is represented by inequalities that are computationally less costly, because they do not explicitly deliver the crossing time. Knowledge of the latter requires an additional computational cost. The novel test thus allows us to avoid this extra cost when we check for the existence of a threshold-crossing; we only need to pay it if the crossing does occur.

From a mathematical point of view the two approaches are equivalent. The state-space technique is just a useful reformulation of the conventional trajectory view that is obtained by suitable mathematical manipulations. The conceptual device of propagating the threshold backwards in time is useful because it performs those manipulations in a transparent manner and gives them an intuitive dynamical meaning.

4.2 Mathematical preliminaries

The final equations to be obtained, eq. (23), can be derived by concepts from vector analysis, Cartesian geometry, and functional analysis; but the derivation is lengthy. To shorten it we use concepts and terminology from affine spaces (Coxeter, 1969; Rockafellar, 1972; Nomizu and Sasaki, 1994; Artin, 1955; Porta Mana, 2011), and differential manifolds (Choquet-Bruhat et al., 1996; Burke, 1987, 1995; Simon et al., 1992; Nomizu and Sasaki, 1994; Marsden and Ratiu, 2007; Bossavit, 1991, 2002; Schouten, 1989; de Rham, 1984; Dodson and Poston, 1991; Ramanan, 2005).

The state-space $S$ of a neuron has a natural vector-space structure (an affine-space structure would also suffice), inherited from the physical quantities that define it: the membrane potential $V$ and other $N-1$ physical quantities $\bm{I}$ whose exact number and definition depend on the specific neuron model (e.g., $\bm{I}$ could represent currents or additional voltages, for example of different compartments). For our purpose it is useful to consider $S$ as an $N$-dimensional differential manifold: its points $\{\bm{s}\}$ are the neuron states, and the quantities $(\bm{I},V)$ are coordinates $V\colon S\to\mathbf{R}$ and $\bm{I}\colon S\to\mathbf{R}^{N-1}$ e.g., the membrane potential of a state $\bm{s}$ is $V(\bm{s})$. These coordinates respect the vector structure of the state-space, i.e. $V(\bm{s}_{1}+\bm{s}_{2})=V(\bm{s}_{1})+V(\bm{s}_{2})$ and likewise for $\bm{I}$.

Every hyperplane in the state-space is defined by an affine equation $\bm{k}^{\intercal}\bm{s}=-\kappa$, the covector $\bm{k}^{\intercal}$ being the normal to the hyperplane, and $-\kappa$ being the affine term. The inequality $\bm{k}^{\intercal}\bm{s}>-\kappa$ defines one of the two half-spaces delimited by the hyperplane. The threshold hyperplane is especially important: it is the set of states $\bm{s}$ whose membrane potential has the threshold value:

[TABLE]

Its equation $\bm{k}^{\intercal}\bm{s}=\kappa$ in coordinates $(\bm{I},V)$ has coefficients

[TABLE]

An affine transformation of the state-space onto itself,

[TABLE]

where M is a linear transformation and $\bm{m}$ a state, maps each hypersurface and half-space $\bm{k}^{\intercal}\bm{s}\geqslant\kappa$ to a hypersurface and half-space ${\bm{k}^{\prime}}^{\intercal}\bm{s}\geqslant\kappa^{\prime}$ with

[TABLE]

(the transformation of the normal $\bm{k}^{\intercal}$ shows why it is a covector rather than a vector).

We now show that the integrable part of the neuron dyamics within a finite time step $h$ is an affine transform. Similarly, the integrable part of of the neuron dynamics we consider an affine evolution determined by the equation

[TABLE]

In a time interval $h$, this dynamics propagates an initial state $\bm{s}_{0}$ at time $t$ into the final state

[TABLE]

This, for each $h$, is an affine transformation of the form (4). In coordinates $(\bm{I},V)$ the linear operator A and vector $\bm{q}$ have the block form

[TABLE]

where B is an $(N-1,N-1)$ matrix, $\alpha$ and $\beta$ numbers, and the dimensionalities of the rectangular matrices $\bm{c}^{\intercal}$, $\bm{d}$ follow accordingly.

4.3 Derivation of the threshold-crossing condition

Let us mathematically summarize the first threshold-crossing test discussed in 4.1. We said that the state evolution (7) can be efficiently used in a time-step scheme in numerical simulations, but we need to test whether a threshold crossing occurred at some time $t\in]t,t+h]$. A necessary and sufficient condition would be the existence of solutions of the trascendental equation in $t$

[TABLE]

which corresponds to the intersection of the trajectory (7) and the threshold hyperplane (2); but it is a costly condition to test.

We now develop the second kind of threshold-crossing test discussed in 4.1, according to the steps I–IV. Steps I and II are performed in full generality for an affine dynamics. Steps III and IV have to be solved on a case-by-case basis, so their analysis below is only a guideline; a concrete example on how to perform them is given in 5.3.3 and 5.3.4.

4.3.1 Hypervolume in parametric form

Consider the threshold hyperplane $V(\bm{s})=\theta$ as an $(N-1)$-dimensional manifold with coordinates $\bm{x}$. It is embedded in state-space via the map

[TABLE]

Each state $(\bm{x},\theta)$ on the threshold, when propagated backwards in time for an interval $h$, traces a curve in state-space (the yellow lines of 3).

The union of these curves is an $N$-dimensional product manifold, called the extrusion (Bossavit, 2003) of the threshold hyperplane. We can use coordinates $(\bm{x},t)\in(\mathbf{R}^{N-1}\times[0,h])$ on this manifold. Its mapping into state-space is given, with the help of eq. (7), by

[TABLE]

where $t>0$ is the direction of the past. This map is analytic, but generally not an embedding because it can have self-intersections. We will see the significance of this in 4.4. It is not an immersion either because it can have singular points; these will be especially important for us because they constitute part of the boundary between spike and no-spike regions. See 3 for a two-dimensional example.

For fixed $t$, the map $\bm{x}\mapsto E(\bm{x},t)$ is affine, and its image is a hyperplane representing the states on the threshold propagated backwards in time for an interval $t$. Combining it with eqs (3)–(5) we find that the backpropagated threshold has $t$-dependent normal and affine terms

[TABLE]

The inequality $\bm{k}^{\intercal}_{t}\bm{s}<\kappa_{t}$ determines the backpropagated half-space, which is below threshold.

4.3.2 Hypervolume boundary in parametric form

We must now find the boundary of the image of the map $E$. The latter is a closed set, being the image of a closed set under a continuous map; its boundary must therefore be the image of some points of the domain. Such points must either lie on the boundaries of the domain, $\mathbf{R}^{N-1}\times\{0\}$ and $\mathbf{R}^{N-1}\times\{h\}$, or be critical points of $E$, or both, because $E$ is differentiable. See the example of 3.

The images of the boundary are easily found from (11): one (image of $t=0$) is the threshold hyperplane, the other ($t=h$) is the hyperplane $\bm{k}^{\intercal}_{h}\bm{s}=\kappa_{h}$, with coefficients given by (12). Explicitly, in coordinates $(\bm{I},V)$,

[TABLE]

Let us find the image of the critical points of $E$. The tangent map of $E$ at a point $(\bm{x},t)$ is

[TABLE]

The tangent map is also denoted $E_{*}$ in some differential-geometry texts.

Its determinant is the inverse ratio between a volume element at that point and its image in state-space (the sign determines their relative orientation). Hence this ratio vanishes at points where volume elements are mapped onto area elements (in other words, $N$ linearly independent vectors in the domain are mapped onto $N$ linearly dependent vectors), which is a feature of the boundary. See the two-dimensional example of 3.

Let us look for points where $\det\mathrm{T}E(\bm{x},t)=0$. In eq. (15), the determinant of the exponential never vanishes, so we only have to consider the determinant of the matrix on the right. This is easily calculated by Laplace expansion along the last row, whose elements all vanish except the last. The cofactor of the last element is $\det\textbf{{{I}}}$ (modulo a sign). Hence

[TABLE]

and the coordinates $(\bm{x},t)$ of critical points satisfy

[TABLE]

This equation says that one of the coordinates $\bm{x}$ has an affine dependence on the remaining ones; let us call these $\bm{y}$. For example, if $c_{N-1}\neq 0$, the equation above has the parametric solution

[TABLE]

Denote this affine dependence by $\bm{x}(\bm{y})$. By taking the derivative of eq. (17) with respect to $\bm{y}$ we have

[TABLE]

a property we will use later.

The locus of critical points in state-space is then given parametrically by a map $\bm{\varGamma}$, found by substitution of eq. (18) in (11):

[TABLE]

This locus has four important interrelated features:

First: from the form of eqs (17) and (20), the locus of critical points is flat along $N-2$ dimensions – corresponding to a fixed value of the coordinate $t$ – and is curved normally to the direction $t$: it is an $(N-2)$-ruled surface.

Second: the points on the locus corresponding to fixed $t$ belong to a hyperplane with coefficients (12), as is easily checked by substitution. In other words, the locus of critical points is the envelope of the backpropagated threshold hyperplanes $\bm{k}^{\intercal}_{t}\bm{s}=\kappa_{t}$, eq. (12), at different times $t$, and its tangent hyperplanes have normals $\bm{k}^{\intercal}_{t}$ given by (12).

Third: considering the feature above for $t=0$, the locus of critical points is tangent to the threshold hyperplane.

Fourth: comparing the dynamics (6), the map for the threshold hyperplane (10), and the critical-point condition (17), we notice that the latter is also the condition for the trajectory (7) of a state $\bm{s}_{0}=(\bm{x},\theta)$ on the threshold hyperplane to have an extremum in the membrane potential $V$. Hence, the locus of critical points is the trajectory of the intersection between the threshold and the $V$-nullcline.

In view of the second property above, let us call the locus of critical points envelope, for brevity.

4.3.3 Hypervolume boundary in implicit form

The next step is the elimination of the parameters $(\bm{y},t)$ to express the envelope $\bm{\varGamma}$ as one or more implicit equations $\chi(\bm{s})=0$ defined on domains $D_{\chi}$ for the state-space coordinates $\bm{s}$. This step can involve a transcendental equation, therefore we cannot give a general solution for it. In this regard, a criticism might be raised, which we immediately address: it looks like we have already taken great care in avoiding to solve eq. (9), only to face another difficult transcendental equation? The current problem is more manageable for three reasons. First, we observe that the equations needed to re-express the curved surface in explicit form are generally easier to handle than (9). 5.3.3 gives a concrete example. Second, the problem of finding the explicit form of the curved surface must be solved only once, whereas the solution of the original equation (9) has to be found for each initial state $\bm{s}_{0}$. Third, the equation for the curved surface is easier to approximate numerically than the original one (9), again because the approximations do not depend on an initial state.

Suppose we have found a function $\chi$ such that $\chi(\bm{s})=0$ is the envelope:

[TABLE]

The condition above leaves $\chi$ completely undetermined (apart from smoothness requirements) outside of the envelope. Therefore, if possible, it is useful to extend the condition as follows. Since $\chi(\bm{s})=0$ is the envelope of the backpropagated threshold hyperplanes $\bm{k}^{\intercal}_{t}\bm{s}=\kappa_{t}$, eq. (12), at different times $t$, it must be tangent to each of them. The inequalities $\bm{k}^{\intercal}_{t}\bm{s}-\kappa_{t}<0$, for each $t$, correspond to the backpropagated below-threshold half-plane, and it is useful to choose $\chi$ in such a way that the inequality $\chi(\bm{s})<0$, $\bm{s}\in D_{\chi}$, determines the side corresponding to the intersection of these half-spaces. An example is given in 5.3.3. The differential $\mathrm{d}\chi$ then is a positive multiple of $\bm{k}^{\intercal}_{t}$ (note that $\mathrm{d}\chi$ and $\bm{k}^{\intercal}_{t}$ are collinear because of the tangency condition).

For the following, let us assume that $\chi$ has been chosen this way.

4.3.4 Boundary intersections and final system

If a state does not cross the threshold hyperplane at any time $t\in[0,h]$, then it must, by definition, belong for every $t$ to the image of the half-space that lies below threshold, when this half-space is backpropagated by the time $t$. Mathematically, this is simply the statement of the equivalence

[TABLE]

established in 4.2, where M and $\bm{m}$ are the coefficients of the affine evolution by time $t$.

By construction in the previous subsection, the inequality $\chi(\bm{I},V)<0$ defines the region of intersection of all such backpropagated half-spaces, bounded by the envelope $\chi(\bm{I},V)=0$. We just need to join to it the condition for the boundaries corresponding to the times $t=0$ and $t=h$, $V<\theta$, $\bm{k}^{\intercal}_{h}\bm{s}<\kappa_{h}$, discussed in 4.3.1.

The states that do not cross the threshold during the time interval $[0,h]$ belong therefore to the no-spike region defined, in coordinates, by

[TABLE]

Some inequalities in this system may turn out to be redundant, i.e. automatically satisfied if the remaining ones are, and can thus be dropped. 5.3.3 illustrates such a redundancy.

4.4 The region of missed spikes

In the previous section we mentioned that the $N$-dimensional product manifold formed by the $(N-1)$-dimensional threshold hyperplane and the time interval $[0,h]$ presents self-intersections when mapped to the state-space. Two-dimensional examples of such a region are shown in 3, lower panel, and in 5, the region called $\textit{S}_{2}$.

A state $\bm{s}$ in the self-intersection region corresponds to two or more different coordinates $(\bm{x},t)$:

[TABLE]

Note that $\bm{x}\neq\bm{x}^{\prime}$, $t=t^{\prime}$ is impossible for an affine transformation, since the first column of (15) can never vanish. The condition above simply means that $\bm{s}$ crosses the threshold at $\bm{x}$ after an interval $t$ and at $\bm{x}^{\prime}$ after an interval $t^{\prime}$. By continuity of the dynamics (6), one of the two must be a crossing from above, and one from below: this is exactly the scenario of double threshold crossing illustrated in 1.

Suppose we have a probability distribution for the initial states, for example one that is invariant under the dynamics (6). The probability of the self-intersection region $P(X)$ is the probability that the initial state will lead to a double-crossing of the threshold, and therefore be missed. This fact will be used in 5.5.2 to estimate the number of spikes missed.

4.5 Convexity, approximations, optimization

The affine dynamics (7) preserves affine combinations of solutions – and therefore convex combinations as well. If $\bm{s}_{1}$, $\bm{s}_{2}$ are two arbitrary initial states in the no-spike region, then their propagated states also satisfy $V[\bm{s}_{1}(t)]<\theta$ and $V[\bm{s}_{2}(t)]<\theta$ when $0\leqslant t\leqslant h$; and also the propagation of their convex combination $\lambda\bm{s}_{1}+(1-\lambda)\bm{s}_{2}$, $0<\lambda<1$ satisfies

[TABLE]

i.e. it lies in the no-spike region. This proves that the no-spike region is convex.

Convexity is important for approximations. If transforming the parametric equation for the envelope (20) into an implicit form (21) turns out to be analytically impossible, we can numerically find some points on the envelope and then approximate the latter by simplices constructed on these points, as finely as needed. In other words we can triangularize the envelope. Owing to convexity, the triangularization can be done completely on the side of the no-spike region (the corners of the simplices touch the envelope), or completely from that of the spike region (the barycentres of the simplices touch the envelope). This way we can formulate a test with no false positives, or one with no false negatives, or both.

The envelope is moreover a ruled surface, as shown in 4.3.2. The triangulation on the side of the no-spike region can therefore be conveniently chosen in such a way that one face of each simplex fully lies on the envelope.

5 Implementation example:

leaky integrate-and-fire neuron with exponentially decaying post-synaptic currents

5.1 The example model

In the previous section we mathematically developed the idea of propagating the threshold backward in time in order to check whether a threshold-crossing occurs in a time-stepped dynamics. The derivation, valid for an affine subthreshold dynamics, is general and therefore also quite abstract; moreover, it involves a couple of mathematical steps (III and IV in 4.1) for which no general formulae can be given.

To explain the idea in more concreteness and to give an example of how to face all its steps, we now apply the scheme to a simple but relevant model with a 2-dimensional state-space: the leaky integrate-and-fire neuron with exponentially decaying post-synaptic currents. This model has a homogeneous linear dynamics on a 3-dimensional state-space (Rotter and Diesmann, 1999), where the third coordinate is the input current. If this current is constant, the dynamics can be rewritten as a 2-dimensional affine one.

Leaky integrate-and-fire models, despite their simplicity, approximate the behavior of real neurons with high accuracy (Rauch et al., 2003). The model with exponential synaptic currents captures important properties of real neurons: The postsynaptic potential has a finite rise and decay time and the membrane potential is a continuous function of time. Continuity avoids artificial synchronization, present in simpler models. Moreover, the model is to some extent analytically tractable. For short synaptic time constants, the mean firing rate (Fourcaud and Brunel, 2002a) as well as the linear response to small inputs (Schücker et al., 2015) can be obtained analytically.

5.2 Mathematical preliminaries: terms in block form

The example model has a 2-dimensional state-space for a single neuron, defined by the post-synaptic current $I$ and the membrane potential $V$, which are also our coordinates. Its subthreshold interspike dynamics (a) in 4 is affine:

[TABLE]

where $C$ is the membrane capacitance, $\tau$ is the membrane time constant, $I$ is the synaptic input current and $I_{\text{e}}$ is the external input current. The membrane potential $V$ is subject to dissipation with time constant $\tau$ and integrates the post-synaptic current $I$. The latter decays exponentially with time constant $\tau_{\text{s}}$. Typical values of the parameters are $\tau=10\;\textrm{ms}$, $C=250\;\textrm{pF}$, $\tau_{\text{s}}=2\;\textrm{ms}$, and the threshold $\theta=20\;\textrm{mV}$.

Incoming spikes are incorporated in the equation for the current as Dirac deltas, and the external current $I_{\text{e}}$ has jump discontinuities in time. In the timestepped evolution, such discontinuous events are implemented as instantaneous changes in the initial state $\bm{s}_{0}$ at each timestep. Hence we do not need to consider them explicitly in the equations above (Rotter and Diesmann, 1999).

In terms of the block form of 4.2 we have

[TABLE]

The exponential of A is

[TABLE]

which determines the evolution of the neuron state $\bm{s}_{0}$ by an affine map as in eq. (7).

Since the state space is 2-dimensional, the “hyperplanes” and “hypersurfaces” of 4 are straight lines and curves. In particular, the threshold hyperplane is a line; when propagated by a time $t$ it maps onto a line with covector and affine term given by eq. (12), explicitly

[TABLE]

For this model, testing for threshold-crossing by checking the intersection of a propagated state and the threshold line means solving the following transcendental system in $t$, given the initial state $\bm{s}_{0}=(I,V)$:

[TABLE]

for which we cannot find a solution in analytic form; we would have to resort to bisection algorithms.

5.3 The threshold-crossing condition

5.3.1 Hypervolume in parametric form

The product manifold “threshold line${}\times{}$time interval” is in this case 2-dimensional, with coordinates $(x,t)$. Here $\bm{x}$ is the value of the current at threshold crossing and $t$ the corresponding time point. Its mapping $E$, eq. (11), to the state-space is

[TABLE]

This map is shown in 3: the $x$ isocurves are yellow and in the lower panel they represent the trajectories of states terminating on the threshold; the $t$ isolines are blue and represent the threshold line propagated at different times. Each area element $\mathrm{d}\bm{x}\land\mathrm{d}t$ in the domain – the small rectangles in the upper plane – is mapped into an area element $\mathrm{d}\bm{I}\land\mathrm{d}V$ in the image. Note how these area elements are rotated and sheared. The thicker green curve is the set of singular points where the determinant of the tangent map vanishes: $\det(\mathrm{T}E)=0$. Such points are singular because around them the images of the area elements get flattened to one dimension. In 5.4 we discuss the region of self-intersection bounded by the thick light cyan line, the thick dark blue line, and the thick green curved line.

5.3.2 Hypervolume boundary in parametric form

The boundaries of the image of the map $E$ must be, as explained in 4.3.3, a subset of the images of the boundaries of the domain, $\mathbf{R}\times\{0\}$ and $\mathbf{R}\times\{h\}$, and of the envelope.

The images of the boundaries, with general equations (13) and (14), in terms of $(I,V)$ follow

[TABLE]

The set of critical points of the map $E$ is in this case a 1-dimensional curve, given in parametric form by

[TABLE]

the coordinate $\bm{y}$ of the general form (20) do not exist in this case, because the threshold is 1-dimensional. 3 visualizes the set by the green curve.

5.3.3 Hypervolume boundary in implicit form

The next step in our procedure is to convert the parametric equation (34) of the envelope into an explicit or implicit equation for the coordinates $(I,V)$. As 4.3.3 does not provide a general algorithm, below we illustrate the process using our example model.

By equating the first component of eq. (34) to the coordinate $I$ and solving for $t$, we find

[TABLE]

required for $0\leqslant t\leqslant h$ and a real logarithm.

Substituting eq. (35a) into the $V$ coordinate of eq. (34) we find

[TABLE]

subject to the condition (35b). Let us analyse this equation, in view of its extension to an inequality of the form $\chi(I,V)<0$ as required in 4.3.3. First, we observe that

[TABLE]

The inequality can be proven by studying the derivative of the fraction with respect to $t$. The derivative is always negative in the range above and the only maximum of the fraction is the value unity assumed at $t=0$.

Inspection of equation (36) and of its parametric form (34) shows that we must consider three cases: $I_{\text{e}}\lesseqgtr\theta C/\tau$, i.e. whether the external current is smaller or larger than the rheobase current

[TABLE]

this is the current necessary to reach threshold in an infinite time starting from any state with $I\leqslant 0$.

• •

If $I_{\text{e}}<I_{\theta}$, then $I$ is restricted to

[TABLE]

In this case, using the inequality (37), the $V$ component of the envelope (34) is always smaller than the threshold:

[TABLE]

• •

If $I_{\text{e}}>I_{\theta}$, then $I$ must be negative and restricted to

[TABLE]

In this case, using the inequality (37), the $V$ component of the envelope (34) is always larger than the threshold:

[TABLE]

• •

If $I_{\text{e}}=I_{\theta}$, the envelope degenerates to a point: $\bm{\varGamma}(t)=(0,\theta)$, which is the limit point reached in infinite time from any initial point in the no-spike region; this is the geometric interpretation of the equality of external and rheobase currents.

The function representing the envelope, eq. (21) of 4.3.3, has in this case the explicit form

[TABLE]

If we calculate its differential, as discussed in 4.3.3, we find that the latter is a positive multiple of the differential of the backpropagated threshold for each $t$. This is true in each of the three cases above. Consequently, the inequality $\chi(I,V)<0$, $(I,V)\in D_{\chi}$ always includes the no-spike region.

In summary, we arrive at an analytic, implicit equation for the curved boundary (43). The expression is a transcendental function in $I$, owing to the generally irrational exponent $1-\tau_{\text{s}}/\tau$, but it is used in an inequality, hence we do not need to find its roots. This is in contrast to the original transcendental equation for the threshold-crossing condition (30), which requires a bisection algorithm to find its solution.

5.3.4 Boundary intersections and final system

We can now assemble the system of inequalities defining the no-spike region consisting of the boundaries (32) and (33), and the envelope (43).

With some rearrangements, simplifications, and the introduction of two new functions $f$, $b$, the condition reads:

[TABLE]

4 illustrates the system for the two cases $I_{\text{e}}<I_{\theta}$ and $I_{\text{e}}>I_{\theta}$. In the former case all three inequalities are necessary; in the latter, as well as for $I_{\text{e}}=I_{\theta}$, the last inequality is automatically enforced by the first because its right-hand side is larger than the threshold $\theta$ (see eq. (42)).

5.4 The region of missed spikes

In 3, lower panel, the trajectories of several states during a timestep $h$ and ending on the threshold line are represented by yellow curves. In that figure we can identify a region were such trajectories self-intersect: it is bounded by a segment of the thick light blue line, a segment of the thick dark blue line, and a portion of the thick green curved line. Trajectories with initial states in this region must therefore cross the threshold twice during the interval $]t,t+h]$. As explained in 4.4, their threshold-crossing is not detected by the sole condition $V[\bm{s}(t+h)]\geqslant\theta$. All states in this region thus generate spikes that are missed by the standard test (1).

This region is crucial for the comparison of the performances of schemes implementing the present lossless method and schemes relying on the standard test (1). This comparison is quantitatively made in 5.5.3.

5.5 Optimization, time performance, and accuracy

5.5.1 Numerical implementation and optimization

In the last section we arrived at the system of three inequalities (44) that determines whether or not the current state $(V,I)$ will cross the threshold within the timestep $h$. The state will cross the threshold if the system is not satisfied, and it will not cross the threshold if the system is satisfied. This system of inequalities constitutes the lossless method in the present model. The system requires the current timestep $h$, state $(V,I)$, and external electric current $I_{\text{e}}$ as inputs (cf. 4). The timestep $h$ is the minimum between the global timestep and the time interval up to the next input coming from other neurons, hence it can differ every time the test is called. The external electric current $I_{\text{e}}$ may also vary, stepwise, during the simulation, hence it may also be distinct at each call of the test.

If the system of inequalities is satisfied, thus predicting the absence of a spike within a timestep $h$, then the state of the neuron is evolved by applying the propagator (7) with (27)–(28), leading to a new state, and the procedure starts again. If the system is not satisfied, thus predicting the occurrence of a spike within $h$, it is then necessary to compute the time $t_{\theta}$ at which the threshold is crossed. This calculation, explained in Appendix A, is made by interpolation between the current state and time, and the state and time $t_{\text{max}}$ at which the membrane potential would reach its maximum if allowed to increase above threshold. Once $t_{\theta}$ is calculated, a spike is emitted, communicated to the postsynaptic neurons, and the membrane potential is reset to $V_{\text{reset}}$ for a refractory period. When this refractory period is over the procedure starts again.

In the evolution loop just described, the membrane potential is reset to a value below threshold as soon as it crosses the latter. Thus no initial state can have $V\geqslant\theta$. This means that the first inequality (44b) in the system is always satisfied and can be dropped. Only inequalities (44c) and (44d) have to be assessed, leading to the reduced system

[TABLE]

We first discuss the case $I_{\text{e}}<I_{\theta}\equiv\theta C/\tau$.

The geometric meaning of the inequalities above is illustrated in 5. The figure shows four regions: $\textit{NS}_{1}$, $\textit{NS}_{2}$, $\textit{S}_{1}$, $\textit{S}_{2}$. The no-spike region is the union of $\textit{NS}_{1}$ and $\textit{NS}_{2}$, the spike region the union of $\textit{S}_{1}$ and $\textit{S}_{2}$. In the figure they are separated by a thick line, partly curved and black, partly straight and blue. Subregion $\textit{S}_{2}$ is particularly important: it is the region of missed spikes discussed in 5.4, corresponding to the self-intersection region of 3, lower panel. It contains those states that lead to spikes missed by schemes that rely on the standard test (1).

Region $\textit{S}_{1}$ is separated from $\textit{S}_{2}$ by a dashed blue line, and from $\textit{NS}_{1}$ by a continuous blue line, the continuation of the dashed one. This partly dashed, partly continuous blue line corresponds to the equation $V=f_{h,I_{\text{e}}}(I)$. Hence if the inequality $V<f_{h,I_{\text{e}}}(I)$ is not satisfied then the initial state is in region $\textit{S}_{2}$ or on its blue boundary, and there will be a spike. If the inequality is satisfied the state could be in $\textit{S}_{2}$ – spike – or $\textit{NS}_{1}\cup\textit{NS}_{2}$ – no spike; an undetermined case. This inequality requires modest computational costs because it is linear in $I$ and $I_{\text{e}}$ and involves exponentials of $h$. Regions $\textit{S}_{2}$ and $\textit{NS}_{2}$ are separated by a black curve: this is the envelope, corresponding to the equation $V=b_{I_{\text{e}}}(I)$ for $I\in\left[I_{\theta}-I_{\text{e}},\mathrm{e}^{\frac{h}{\tau_{\text{s}}}}(I_{\theta}-I_{\text{e}})\right]$. Hence if the inequality $V<b_{I_{\text{e}}}(I)$ is not satisfied the initial state is in $\textit{S}_{2}$ or on its boundary, and there will be a spike. If the inequality is satisfied the initial state is either in $\textit{NS}_{2}$, or in $\textit{NS}_{1}$ with $I_{\theta}-I_{\text{e}}<I<\mathrm{e}^{\frac{h}{\tau_{\text{s}}}}(I_{\theta}-I_{\text{e}})$, and no spike will occur.

The computationally most expensive inequality is $V<b_{I_{\text{e}}}(I)$ because it involves irrational powers of $I$ and $I_{\text{e}}$. It is advisable to avoid its direct computation as often as possible by pre-testing a linear inequality. In 4.5 we discussed how such a pre-test is indeed possible thanks to the convexity of the no-spike region. There, we argued that the curved envelope can be approximated by triangular hypersurfaces, which simply reduce to one straight segment in the present two-dimensional case: this is the dot-dashed red line in 5, separating $\textit{NS}_{1}$ and $\textit{NS}_{2}$. This line has equation $V=g_{h,I_{\text{e}}}(I)$ with

[TABLE]

and the corresponding inequality

[TABLE]

has the same computational costs as $V<f_{h,I_{\text{e}}}(I)$.

If the auxiliary inequality $V<g_{h,I_{\text{e}}}(I)$ is satisfied, the initial state is in $\textit{NS}_{1}$ and $V<b_{I_{\text{e}}}(I)$ is also satisfied. It is therefore convenient to test the auxiliary inequality before the computationally costly one, which can be discarded if the test is positive. 5 suggests that this test might be positive for the majority of initial states because region $\textit{NS}_{1}$ is much wider than $\textit{NS}_{2}$. This possibility would be very advantageous, but we now argue that it should be verified by a dynamical analysis.

The system (45) can be translated into a computational algorithm in several different ways, depending on the order of evaluation of its two inequalities and of the auxiliary inequality (47). In simplified terms, such an algorithm consists in a sequence of tests – variously implemented as if, and, or constructs – for finding the initial state in space-time regions $R_{1}$, $R_{2}$, and so on. The order of these tests is important. The average time cost of the algorithm in a long simulation is given by $\mathop{\textstyle\sum}_{i}p_{i}c_{i}$, where $p_{i}$ is the frequency with which states are found in region $R_{i}$, which we call “occupation frequency”, and $c_{i}$ is the cumulative time cost of the test for region $R_{i}$. This time cost $c_{i}$ is cumulative in the sense that all tests up to the $(i-1)$th must have been performed, with false outcomes, to arrive at the test for $R_{i}$. The efficiency of an algorithm therefore depends on the mathematical form of the inequalities defining a region and on the occupation frequencies of the regions, determined by the dynamics. These two factors can be extrapolated by a theoretical analysis, or more practically measured by running long test simulations with typical network setups corresponding to the cases one is interested in.

We now try to determine the most efficient algorithm for the present case. Region $\textit{S}_{1}$ is the least costly, because bounded by one line and therefore involving one inequality linear in $I$; then region $\textit{NS}_{1}$, bounded by two lines involving two linear inequalities; and finally regions $\textit{NS}_{2}$ and $\textit{S}_{2}$, bounded by the curve that involves rational exponentiation. For this example model, we tried different orderings but show here only two possible extreme cases to illustrate that there is no significant difference in the computational cost.

Algorithm 1

bool is_spike($h$):

pre-compute $g_{h,I_{\text{e}}}(I)$

if $V\leqslant f_{h,I_{\text{e}}}(I)$ and $V<g_{h,I_{\text{e}}}(I)$ then

     return false

else if $V\geqslant g_{h,I_{\text{e}}}(I)$ then

     return true

else if $V\geqslant b_{I_{\text{e}}}(I)$ then

     return true

else

     return false

Algorithm 2

bool is_spike($h$):

if $V\geqslant g_{h,I_{\text{e}}}(I)$ or

[$V\geqslant f_{h,I_{\text{e}}}(I)$ and $V\geqslant b_{I_{\text{e}}}(I)$] then

     return true

else

     return false

Algorithm 1 is based on the assumption that the occupation frequency of a subregion is proportional to that subregion’s relative size. If we check the two largest first, in the order $\textit{NS}_{1}$, $\textit{S}_{1}$, $\textit{S}_{2}$, $\textit{NS}_{2}$, we are therefore more likely to exit the test in its first if branches. The value of $g_{h,I_{\text{e}}}$ is used in two if branches, so it is computed just once and saved before the if sequence in order to save some computations.

Algorithm 2 uses a composite or-and condition rather than several ifs. Assuming left-to-right evaluation, the algorithm corresponds to testing first $\textit{S}_{1}$ (left side of or), then $\textit{S}_{2}$ (right side of or), either case leading to a spike. If neither is true, no further tests are necessary because the state must necessarily be in the no-spike region. The test for $\textit{S}_{2}$ is made less costly on average by using the auxiliary inequality (47). This algorithm uses one test less overall than the previous one, but it may require one more test on average, if $\textit{NS}_{1}$ is the region with highest occupation frequency. The left-to-right evaluation assumption does not always hold in modern processors, which build their own statistics to optimize the test order of logical constructs.

The analysis assumed $I_{\text{e}}<I_{\theta}\equiv\theta C/\tau$. The reduced system (45) and the Algorithm 1 and 2 are, however, also valid in the case that $I_{\text{e}}\geqslant I_{\theta}$, corresponding to the lower panel of 4. In this case subregions $\textit{NS}_{2}$ and $\textit{S}_{2}$ do not exist below the threshold, and the inequalities $V<b_{I_{\text{e}}}(I)$ and $V<g_{h,I_{\text{e}}}(I)$ are always satisfied when $V<\theta$; they always evaluate to true in both algorithms. Both algorithms therefore correctly distinguish spiking from non-spiking states in this case, although they become inefficient owing to the additional superfluous evaluations of $b_{I_{\text{e}}}(I)$ and $g_{h,I_{\text{e}}}(I)$ for spiking states. We decide not to modify them in the present work because the case $I_{\text{e}}\geqslant I_{\theta}$ is unusual in real applications. More efficient algorithms for this case can be designed by interested readers following the guidelines just given in this section.

5.5.2 Occupation frequencies

We want to measure the occupation frequencies in the typical case of a neuron embedded in a recurrent network, receiving fluctuating synaptic input. The setup for this simulation is illustrated in 6A and its formulae explained in Appendix B. One neuron is coupled with strengths $J$ and $-J$ to an excitatory and an inhibitory Poisson generator and also receives a constant external current. The Poisson generators each mimic an excitatory and an inhibitory population. The neuron thus receives a fluctuating input current having average $\mu$ and variance $\sigma^{2}$. In this instance the code of the simulation includes a subroutine that informs us of the current state every time the lossless method test is called, without altering the test or the dynamics.

7A gives a visual idea of the occupation frequencies for $\mu=15\;\textrm{mV}$, $\sigma^{2}=25\;\textrm{mV}^{2}$, and $J=0.1\;\textrm{mV}$ (all expressed in volts through multiplication by a resistance of \tau/C=40\;\textrm{M\Omega}), which correspond to the case $I_{\text{e}}<I_{\theta}$. These values correspond to a composite average input of $250\,000\;\mathrm{spikes/s}$, and a total average input current of $400\;\textrm{pA}$. This presynaptic input makes the neuron fire at an average rate of $7\;\textrm{spikes/s}$. A sample of its membrane dynamics is shown in 6B. Subregion $\textit{NS}_{1}$ has the overwhelmingly largest occupation frequency. The other three subregions have actually very small areas, as clear from the axis ranges of 7B, owing to the very small value of the average timestep, $h=4\times 10^{-3}\;\textrm{ms}$, given by the inverse of the input rate in this hybrid scheme.

A more precise comparison of the occupation frequencies of the four regions $NS_{1}$, $NS_{2}$, $S_{1}$, $S_{2}$ is shown in 8 for several combinations of three network parameters, producing different dynamic regimes. The parameters are the average $\mu$, the variance $\sigma^{2}$ of the input current, and the presynaptic coupling strength $J$ (all expressed in volts through multiplication by a resistance of \tau/C=40\;\textrm{M\Omega}). The values of the parameters $(\mu,\sigma^{2},J)$ include typical realistic cases as well as some extreme cases, like unusually high coupling strengths. Each panel of 8 shows the occupation frequencies for a set of dynamical regimes with constant $J,\beta$ and several $\sigma^{2}$. The panels in the last row correspond to the case $I_{\text{e}}\equiv\mu C/\tau\geqslant I_{\theta}\equiv\theta C/\tau$, or $\mu\geqslant\theta$, in which subregions $\textit{NS}_{2}$ and $\textit{S}_{2}$ do not exist below threshold.

It is important to remember that the boundary and size of the regions of 5 vary with the timestep $h$, which is a parameter of the simulation scheme, not of the dynamics per se. In an event-driven or hybrid scheme, this step varies inversely with the event input rate, which for Poisson input generators is proportional to $\sigma^{2}/J^{2}$. As a consequence, the frequencies displayed in 8 are not determined by the Liouville distribution of the dynamics (26) alone, but also by the details of the numerical-implementation scheme. The dependence of the boundaries on $h$ is illustrated in 7B. As $h$ decreases, the line $V=f_{h,I_{\text{e}}}(I)$ and the auxiliary line $V=g_{h,I_{\text{e}}}(I)$ gets closer to the threshold, and the subregions $\textit{S}_{1}$, $\textit{S}_{2}$, $\textit{NS}_{2}$ disappear. This is plausible since in the limit of $h=0$ we are not evolving the initial state at all. As $h$ increases the point of tangency between the envelope $V=b_{I_{\text{e}}}(I)$ and the line $V=f_{h,I_{\text{e}}}(I)$ moves to increasingly lower voltages and higher currents; subregion $\textit{S}_{2}$ takes over $\textit{S}_{1}$ and subregion $\textit{NS}_{1}$ becomes wider. For typical timestep values of several milliseconds, though, subregions $\textit{S}_{1}$, $\textit{S}_{2}$, $\textit{NS}_{2}$ are still very small.

A rough estimate of the dependence of the areas of the bounded regions $\textit{S}_{2}$ and $\textit{NS}_{2}$ on the parameters $(\mu,\sigma^{2},J)$, for $\mu<\theta$, can be obtained by looking at 8 and considering that these areas together form a triangle with vertices

[TABLE]

This triangle has base $\lvert I_{f}-(I_{\theta}-I_{\text{e}})\rvert$ and height $\lvert\theta-b_{I_{\text{e}}}[\mathrm{e}^{h/\tau_{\text{s}}}\,(I_{\theta}-I_{\text{e}})]\rvert$. Expressing $h$ and $I_{\text{e}}$ in terms of $\mu$ and $\sigma^{2}$ using eqs (53), where $h$ is inversely proportional to the input rate $r_{\text{I}}+r_{\text{E}}$, we find

[TABLE]

When $\tau_{\text{s}}\lesssim\tau$ and $J^{2}\lesssim\sigma^{2}$ a Taylor expansion in $J^{2}/\sigma^{2}$ to fourth order gives a good approximation, with a relative error below $10\,\%$:

[TABLE]

This approximate formula shows that subregions $\textit{NS}_{2}$ and $\textit{S}_{2}$ grow with the square of the input mean $\mu$ and with the third or fourth power of the ratio $J^{2}/\sigma^{2}$. Recall that these regions do not exist for $\mu\geqslant\theta$. The occupation frequencies do not depend on the areas alone, however, but also on the dynamics, as explained in the previous section. We can identify several other dynamical mechanisms for their dependence on the parameters $(\mu,\sigma^{2},J)$:

• •

an increase in mean input $\mu$ leads to more frequent threshold crossings, thus frequently bringing the voltage to its reset value, underneath subregions $\textit{NS}_{2}$ and $\textit{S}_{2}$. The occupation frequencies of these subregions may therefore decrease with $\mu$ even though their areas grow with $\mu$;

• •

for low mean input $\mu$, an increase in variance $\sigma^{2}$ means a higher chance of high-$V$ regions, and thus an increase in the occupation frequencies of $\textit{NS}_{2}$ and $\textit{S}_{2}$, even though their areas shrink with $\sigma^{2}$;

• •

for mean input $\mu$ close to the threshold, an increase in the variance $\sigma^{2}$ leads to more frequent threshold crossings, and may thus increase occupation frequency of $\textit{S}_{2}$ with $\sigma^{2}$, even though its area shrinks with $\sigma^{2}$.

The occupation frequency of subregion $\textit{NS}_{1}$ (green circles) dominates all others, varying from $90\,\%$ to $100\,\%$ depending on the network parameters. Subregion $\textit{S}_{1}$ (yellow crosses) follows in order of frequency and is the most frequently visited between the two spike subregions. Subregions $\textit{NS}_{2}$ (blue squares) and $\textit{S}_{2}$ (red stars) are scarcely visited for lower synaptic amplitudes, with frequencies from [math] to $10^{-6}$; and slightly more often at higher synaptic amplitudes (frequencies from [math] to $10^{-2}$).

5.5.3 Time performance and accuracy of a hybrid scheme based on the lossless method

The average time costs of Algorithm 1 and Algorithm 2 within a hybrid scheme can be assessed by real-time simulation measurements. The average time cost of a hybrid scheme based on the standard test (1) can also be assessed in the same way for comparison. We therefore compare the two algorithms of the lossless method and the sufficiency test in this section.

The basic setup is the same as for the frequency analysis of 5.5.1, explained in Appendix B, with network parameters $(\mu,\sigma^{2},J)$. The only difference is that in the present case the code does not include the subroutine that informs us of the frequencies, which would otherwise increase and bias the real-time durations of the simulations. Three instances of the basic setup are prepared: in the first the neuron is modelled by a hybrid scheme with lossless method Algorithm 1, in the second the neuron is modelled by a hybrid scheme with lossless method Algorithm 2, and in the third the neuron is modelled by a hybrid scheme with the standard threshold-crossing test (1). The various random-number-generator seeds of the three instances are exactly the same, so that the three neurons receive exactly the same input, spike-for-spike; this is essential for a fair comparison between the three schemes. The three instances are run for a long time ($2\times 10^{6}\;\textrm{ms}$) and repeated (in parallel) for several times ($10$), enough to collect reliable statistics. The statistics are collected for the same sets of $(\mu,\sigma^{2},J)$ values as in the frequency analysis. The total real-time length of a simulation depends, in all three schemes, on how often the checkpoints and threshold-crossing tests occurs, and this in turn depends on the presynaptic input frequency, as already discussed.

Rather than showing the results for all sets of parameters $(\mu,\sigma^{2},J)$, which in this case are not very informative, we show in 9 those with the $(J,\mu)$ values that yield the slowest and fastest performances. The average computation costs of Algorithm 1 and Algorithm 2, which embody the lossless method, turn out to be very similar – within each other’s standard deviations – and basically identical in comparison with the cost of the scheme based on the standard test (1). Both are slower than the hybrid scheme with the standard test: from around $33\,\%$ slower in the case of high-activity regime with frequent incoming spikes ($\mu=18\;\textrm{mV}$, $J=0.1\;\textrm{mV}$), to around $8\,\%$ slower in the case of low activity regime and infrequent incoming spikes ($\mu=10\;\textrm{mV}$, $J=5\;\textrm{mV}$). These are the extremes shown in 9. In most other sets of network parameters the hybrid scheme with lossless method was around $20\,\%$–$25\,\%$ slower than the standard hybrid scheme with threshold-crossing test (1).

As mentioned in 5.4, subregion $\textit{S}_{2}$ contains all states for which the standard threshold-crossing test misses a spike. The occupation frequencies of this subregion are therefore a direct measure of the number of spike missed, per neuron, by the standard hybrid scheme. They are shown as red stars in 8 for the various sets of network parameters $(\mu,\sigma^{2},J)$. The frequency of missed spikes does not have a simple monotonic dependence on the three parameters, owing to the interaction of several mechanisms, discussed in 5.5.2. An increase in synaptic coupling $J$ or input current $\mu$ generally leads to more frequently missed spikes because the neuron spikes more often overall. For very high – suprathreshold – input currents, however, the frequency decreases again until no spikes are missed anymore; this change in trend happens because the checkpoints become more frequent. Regimes of low $J$ are diffusion-like processes, where frequent arrival of synaptic events does not lead to missed spikes. Regimes of high $J$ are shot-noise processes, where sudden and infrequent arrival of synaptic events leads to suprathreshold excursions and missed spikes. Separate simulations show that the hybrid scheme based on the lossless method, with either algorithm, reproduce the analytic solution of the neuron model within floating point precision.

For biologically realistic synaptic couplings, $J<1\;\textrm{mV}$, the standard hybrid scheme is $33\,\%$ faster than the scheme with the lossless method, and misses less than $1$ spike every $10^{6}$ test calls, per neuron; for very low couplings $J<1\;\textrm{mV}$ this figure even becomes less than $1$ spike every $10^{8}$ test calls. For synaptic couplings $J>1\;\textrm{mV}$ the standard hybrid scheme starts to miss more spikes, reaching even $1$ missed spike every $500$ test calls for subthreshold average input currents; and it is only $8\,\%$ faster than the hybrid scheme with the lossless method. From these figures a user can decide to use the hybrid scheme with the standard or the lossless test, depending on the desired balance of accuracy and speed.

6 Summary and discussion

We here present a general method to solve the threshold-crossing detection problem for an integrable, affine or linear neuronal dynamics. The method is based on the geometric idea of propagating the threshold plane backwards in time and to determine whether the swept volume contains the initial state, rather than propagating the initial state forward in time to check whether it crosses the threshold. These two procedures are obviously mathematically equivalent, but they distribute the computational load in different ways. The forward-propagation of the state looks for the value of the crossing time; if no such value exists, it means there is no threshold crossing. The backward-propagation of the threshold first tests whether a crossing time exists at all, without yielding its value; the latter is calculated afterwards.

The different distribution of computational load in the two procedures can be explained geometrically and algebraically. The first procedure geometrically checks for the intersection of a curve (the state trajectory) with a hypersurface (the threshold). Algebraically, this corresponds to finding the roots of a system of equations, often transcendental. The second procedure geometrically checks for the “intersection” of a point (the state) with a hypervolume (the threshold trajectory), i.e. the inclusion of the former in the latter. Algebraically, this corresponds to testing a set of inequalities. The latter procedure is more efficient for numerical computations, because it only relies on inequality tests – in which the presence of transcendental functions is much less costly than in an equation that needs to be solved for a particular variable. The determination of the exact crossing-time, which involves bisection algorithms and is the costlier part, is done only when the existence of this value is certain. No root search is unnecessarily performed.

We have calculated the system of inequalities expressing the threshold-crossing condition, for a generic affine or linear neuronal dynamics in any dimension. The result is the conjunction of inequalities (23). It consists of two affine linear inequalities in the state-space variables, voltage and currents, and a non-linear one. The numerical implementation of this system of inequalities can be further optimized, on a case-by-case basis. In order to give a concrete implementation example of the generic inequalities (23), to show their geometrical meaning, and to give an example of optimization,in the present work we apply our procedure, step-by-step, to the 2-dimensional case of a leaky integrate-and-fire neuron with exponentially decaying post-synaptic currents (Rotter and Diesmann, 1999). The generic inequalities (23) take in this case the concrete form (44).

The quantitative data in the present work are obtained by integrating these inequalities into a combined event-and-time-driven simulation framework (Morrison et al., 2007) for large-scale spiking neuronal network models as released by Bos et al. (2015). Implementation and comparison to earlier work show that:

• •

the system of inequalities, the non-linear one in particular, can be expressed analytically in terms of the state-space variables even when the original threshold-crossing condition involves a transcendental equation – and would therefore require bisection algorithms. Compare (44) with (30);

• •

the computationally expensive non-linear function in the system can be conveniently triangularized, speeding up the algorithm even further by testing a linear inequality first, ruling out the majority of initial states;

• •

the new method reproduces the analytic solution of the neuron model within floating point precision. It detects all threshold crossings, in particular those that the approximate test (1) of Hanuschkin et al. (2010), misses in some ranges of mean activity, fluctuations, and synaptic-coupling strength (8);

• •

at the default spike accuracy of $0.1\;\textrm{ms}$ of the reference simulator the new method is $8\,\%$–$33\,\%$ slower than the fastest available solver with spike loss (Hanuschkin et al., 2010). It is therefore of comparable speed as embedded event driven methods (see Hanuschkin et al. (2010) Fig 5, inset).

In practice the method of Hanuschkin et al. (2010) rarely misses spikes for biologically realistic synaptic amplitudes. At low frequencies of afferent synaptic events, however, such missed threshold-crossings can happen. Our new scheme therefore offers an alternative for users who need guaranteed spike detection and are willing to pay a price in terms of slightly longer computation time.

In state-spaces of higher dimensions it is be more difficult to derive the non-linear inequality of the system (23) in implicit form, but its set of approximating flat surfaces can still be easily calculated. In the worst-case scenario of an inequality not expressible analytically, it is still possible to construct a nested sequence of approximating flat surfaces to be tested hierarchically. Such construction only needs to be done once for any given model and detects spikes with any desired precision. Being linear, such nested inequalities likely are computationally less expensive than a bisection algorithm.

The problem of detecting some sort of threshold crossing in a system of coupled first-order linear differential equations appears in many other applications and phenomena like switching, friction, and saturation (Hiebert and Shampine, 1980). For example, in an air-conditioning unit a thermostat controls the on-off state based on a certain threshold value of the room temperature (Shampine, 1994). The dynamics of this system is similar to that described in Sections 5.1–5.2. Another example is the problem of ejecting a pilot such that collision with the aircraft stabilizer is avoided.

The present work use concepts from differential geometry, in particular extrusions (Bossavit, 2003) and critical points of maps between manifolds, and shows that these concepts have a readily understandable geometrical and visual meaning. The notion of extrusion has recently found applications in numerical and discretization techniques for partial and integral differential equations (Desbrun et al., 2005). The notion of critical points of a manifold mapping is ubiquitous in science: from the caustics of propagating seismic fronts, at which the seismic wave changes its phase (Romanowicz and Dziewonski, 2007, § 1.04), to the singularities between two coordinate charts in general relativity (Misner et al., 2003), which affect the accuracy of global navigation satellite systems (Coll et al., 2012; Sáez and Puchades, 2013).

Indeed, the neuron model analysed in 5 exhibits a similarity with the dynamics of a point mass near a black hole. If the simulation timestep $h$ is very large the curved surface separating the states that lead to a spike from those that do not acts like an event horizon in general relativity: a state evolved from the spike region can enter the no-spike region, but once there it cannot escape and will always remain a “no-spike” state. This is only true for the dynamics (6), though, with constant affine term and no resets at threshold. Inputs from other neurons lead to discontinuous changes in the affine term of the dynamics, causing a “transport” of initial states out of the event horizon, from the no-spike to the spike region. Nevertheless, maybe such similarities are more than mere coincidences. For example, the trajectory of the threshold surface of a leaky integrate-and-fire model with $\alpha$-shaped post-synaptic currents (Bernard et al., 1994) can be implicitly expressed in terms of the Lambert-$W$ function (Corless et al., 1996), as an analysis along the lines of 4.2 shows. This function also appears in the implicit expression of point-mass trajectories in $(1+1)$-dimensional general relativity (Mann and Ohta, 1997). It is surely worthwhile to bring the nascent field of neuronal dynamics closer to ideas and techniques from differential geometry and general relativity.

Appendix A Interpolation for threshold-crossing time

If we know that the neuron voltage $V$ is below threshold at times $t$ and $t+h$ but above threshold somewhere in the interval $]t,t+h]$, then by continuity it must reach a maximum above threshold at a time $t_{\text{max}}\in]t,t+h]$. This time can be obtained solving the equation $\dot{V}(t_{\text{max}})=0$, with $V(t)$ given by (26). The solution (Hanuschkin et al., 2010) is

[TABLE]

from which the potential $V(t_{\text{max}})>\theta$ can also be easily calculated.

The time $t_{\theta}\in]t,t_{\text{max}}]$ at which the first threshold crossing occurs, $V(t_{\theta})=\theta$, must lie between $t$ and $t_{\text{max}}$, with $V(t)<V(t_{\theta})<V(t_{\text{max}})$, and can thus be interpolated using a bisection algorithm (Press et al., 2007).

Appendix B Network setup

The dynamics of the membrane potential $V$ and synaptic current $I$ can be described, if the input is treated stochastically and for weak synaptic couplings, by a diffusion process with equations (Fourcaud and Brunel, 2002b)

[TABLE]

with $\xi$ a zero-mean Gaussian process. The parameters $\mu$ and $\sigma^{2}$ characterize the stochastic input and are the mean and variance of the total incoming synaptic current. They are related to the synaptic couplings $\{J_{i}\}$ and firing rates $\{r_{i}\}$ of the input neurons via $\mu=\tau\mathop{\textstyle\sum}_{i}J_{i}r_{i}$ and $\sigma^{2}=\tau\mathop{\textstyle\sum}_{i}{J_{i}}^{2}r_{i}$. If we have one excitatory and one inhibitory population of input neurons, mimicked by two Poisson generators with rates $r_{\text{E}}$ and $r_{\text{I}}$ coupled to the neuron with strengths $J_{\text{E}}$ and $J_{\text{I}}$ (where $J=\frac{\tau}{C}w$ and $w$ is the synaptic weight of the current), and by an input current $I_{\text{e}}$, then

[TABLE]

where $r$ is the output firing rate of the neuron, approximated to linear order in $\sqrt{\tau_{\text{s}}/\tau}$, $\tau_{\text{r}}$ is the refractory time, $V_{\text{r}}$ the reset voltage, and $\zeta$ the zeta function (Abramowitz and Stegun, 1972).

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Author Contributions

All authors have made significant, direct, and intellectual contribution to the work.

Funding

This work was supported by DFG Grant: GSC11; Helmholtz association: VH-NG-1028 and SMHB; EU Grant 604102 (HBP); Juelich Aachen Research Alliance (JARA).

Acknowledgments

PGLPM thanks the Forschungszentrum librarians for their always prompt and kind help, Mari & Miri for continuous encouragement and affection, Buster Keaton for filling life with awe and inspiration, and the developers and maintainers of LaTeX, Emacs, AUCTeX, MiKTeX, arXiv, biorXiv, PhilSci, Python, Inkscape, Sci-Hub for making a free and unfiltered scientific exchange possible.