The mechanism of synchronization in feed-forward neuronal networks

Synchronization in feed-forward subnetworks of the brain has been proposed to explain the precisely timed spike patterns observed in experiments. While the attractor dynamics of these networks is now well understood, the underlying single neuron mechanisms remain unexplained. Previous attempts have captured the effects of the highly fluctuating membrane potential by relating spike intensity f(U) to the instantaneous voltage U generated by the input. This article shows that f is high during the rise and low during the decay of U(t), demonstrating that the -dependence of f, not refractoriness, is essential for synchronization. Moreover, the bifurcation scenario is quantitatively described by a simple relationship. These findings suggest as the relevant model class for the investigation of neural synchronization phenomena in a noisy environment.


Introduction
Synchronization of spike timing is an important concept in neuroscience and has been experimentally observed in many neural systems. The noisy environment of cortical neural networks counteracts synchronization. Thus, it is essential to understand the interplay between noise and synchronizing mechanisms in these systems. The phenomenology can be successfully investigated by integrate-and-fire (IAF) type models, often by simulation. However, to uncover the mechanisms underlying the observed effects it is helpful to investigate analytically tractable stochastic models. Specific systems in this context are the feed-forward subnetworks of the cortex (figure 1(a)) which have been proposed ('synfire chain', Abeles (1991)) as the generators of spatio-temporal spike patterns occurring with millisecond precision (e.g. Grün et al (2002); Prut et al (1998)). These networks are both a neurobiologically plausible substrate and structurally simple. Figure 1(b) shows an example of network activity where a packet of spikes successfully propagates, reaching a spike precision in the millisecond range. The dynamics of an individual neuron is described by the IAF model (Tuckwell 1988) with the membrane potential V governed bẏ where τ m is the time constant and C the capacity of the membrane. I (t) is the sum of all the input currents to the neuron. When V (t) reaches a threshold θ , a spike is emitted, the membrane potential is reset to V 0 , and clamped to this value for an absolute refractory period τ r . The current elicited by a synapse upon the arrival of a presynaptic spike is represented by an α-function (Rotter and Diesmann 1999)ˆ e/τ s · te −t/τ s , whereˆ is the amplitude and τ s the time constant. The neurons in the subnetwork are assumed to receive balanced Poisson spike input from a large embedding network (Brunel 2000). The resulting large fluctuations of V cause threshold crossings at a low rate consistent with the input rate at a single synapse. In a mean field approach to network dynamics, the spiking activity of a group of w neurons is replaced by its expectation value wρ(t), where ρ(t) is the spike density (figure 2(a)) of a single neuron. The approach assumes identical units receiving common input and independent realizations of the fluctuations in the system. The synchronization dynamics in feed-forward networks and its bifurcation scenarios are now well understood , based on the single neuron response, characterized by ρ(t), to synchronous activity. For this analysis, ρ(t) was obtained by simulations. In this paper, an analytical expression for ρ(t) is derived, revealing the mechanisms of synchronization. (a) IAF model neurons (gray disks; τ m = 10 ms, C = 250 pF, V 0 = −70 mV, θ = −55 mV and τ r = 2 ms) are arranged into successive groups (enclosed by dashed boundaries) and receive input (arrows; α-function current,ˆ = 46 pA, τ s = 0.3 ms and delay 2 ms) from all the w = 100 neurons in the preceding group. This structure is considered as a part of the large network of the cortex supplying each neuron with excitatory (ˆ E =ˆ and λ E = 2 Hz) and inhibitory (ˆ I = −ˆ and λ I = 12.61 Hz) independent Poisson processes at K E = 17500 and K I = 2400 synapses, respectively. (b) Dots in boxes 1 to 10 represent spike times of all neurons. In every box (labeled by group number) one vertical position is reserved for each neuron of the group. The first group (1) is stimulated by a packet of 100 spikes shown in box 0 drawn from a Gaussian distribution (center at t = −2 ms and σ = 3 ms). Activity synchronizes under an intermediate loss in the number of spikes in the packet. The neurons are spontaneously active (λ = λ E , dots occurring before the arrival of the packet) because of the voltage fluctuations induced by background.

Spike intensity in response to synchronous input
If spikes are generated by an inhomogeneous Poisson process, the neuron's instantaneous rate or intensity f (t) equals ρ(t). Arbitrarily short time intervals between spikes are allowed and occur with high probability. However, it has been indicated that refractoriness can have a strong effect on ρ(t), thereby benefiting neural precision (Berry and Meister 1998). Thus, it is a proposed synchronization mechanism in cortical feed-forward subnetworks (Gewaltig 2000;Kistler and Gerstner 2002). Here, we show that synchronization does not originate from refractoriness.
A simple way of including refractory effects is to define an inhomogeneous renewal process (Berry and Meister 1998;Cox 1962;Gerstner and Kistler 2002) by a hazard function or conditional intensity where t 1 denotes the time of the first spike, f (t) is the intensity Plesser and Gerstner 2000;Rotter 1994) when the neuron is not refractory, and r (t − t ) is a recovery function describing the reduced spike probability following the last spike at time t . In  Figure 2. Analysis of the free intensity during fast input transients. (a) Spike density ρ(t) (gray area) and free intensity f (t) (solid curve) of the IAF model in response to packets of A = 80 input spikes with temporal dispersion σ = 2 ms. ρ(t) is obtained from a histogram of spike times over 10 4 simulation trials. The dashed curve shows the membrane potential U (t) generated by the probability density of the input spikes. (b) Recovery function r (t − t ) modeling neuronal refractoriness following a spike at time t (τ a = 4 ms and σ r = 4 ms). A specific function of r (t − t ) enables the estimation of f (t) (see (a)). (c) f (t) versus U (t) (clockwise progression of time) for two different pairs (A, σ ) of input parameters. (d) Division of f byU removes ambiguity during the rise time (U > 0) of U (t). Gray curves in (c) and (d) show analytical predictions (9) for f assuming infinite refractoriness. this context, we call f (t) the free intensity, it replaces the fluctuations of the membrane potential and the threshold in the IAF model by a point process with time-dependent intensity ('escape noise').
The expected membrane potential U (t) generated by a packet of synchronous input spikes with a Gaussian probability density φ(t; σ ) is where σ is the standard deviation of φ(t; σ ), A is the number of spikes, u(t) is the voltage response to a single input spike (post-synaptic potential, PSP), and * denotes convolution (figure 2(a)). The task is to study the influence of refractoriness on the neuronal response to such synchronous activity and simultaneously find a suitable model for f . Therefore, we investigated the shape of f (t) in simulations of the IAF model. It has been shown (Bi 1989;Jones et al 1985) that the spike density ρ(t) of the process (2) is given by where is the probability of finding the neuron refractory. This implicit equation can be iteratively solved for f (t) enabling the removal of refractory effects from ρ(t). The procedure starts with f 0 (t) = ρ(t) as an initial value and then calculates f is consistent with the interval distribution of the IAF model subject to dc input (not shown, (·) is the Heaviside step function). In order not to underestimate f , refractoriness is maximized (large τ a and σ r ) under the constraint of convergence for all input packets (A, σ ). f (t) is robust with respect to variations of r (s) mainly depending on the area of 1 − r (s). Figure 2(a) illustrates a typical time course of f (t). Neuron models where f (t) only depends on the instantaneous state of the neuron provide a successful mathematical framework for the analysis of network dynamics (Gewaltig 2000;Kistler and Gerstner 2002). However, figure 2(c) demonstrates that f is not a unique function of U . The firing intensity is high while U (t) is rising and low otherwise. Note that this cannot be explained by hysteresis since the effects caused by the spike history have already been accounted for by the recovery function (2). Moreover, different input packets result in considerably different values of f at the same voltage U . Taking the IAF neuron model as a biologically plausible reference, these findings rule out f (U )-models (Gewaltig 2000;Kistler and Gerstner 2002) for the analysis of synchronization. Can we find an appropriate model of f and, if so, what are its essential ingredients? A hint is provided by figure 2(d): during the rise of U the ratio of f and the temporal derivativeU constitutes a unique function of U .

Dependence on the voltage derivative
Consider the probability density of the membrane potential p (V, t) in the presence of a strong voltage transient U (t) as in (3). Neglecting the spike threshold and using linearity of the membrane potential dynamics, we have where p st (V ) is the stationary density. The spike density ρ(t) equals the probability current at the threshold V = θ which acts as an absorbing boundary (Tuckwell 1988). The rise time of U (t) is short compared to the membrane time constant. Thus, we assume that the shape of p(V, t) for V < θ is not changed by the interaction with the absorbing boundary. We further assume that an individual trajectory reaches the threshold at most once during the excursion described by U (t) effectively considering infinite refractoriness. Then, the spike density equals the negative rate of change of the survival probability S(t) = Pr [V (t) < θ ] = 6 and we find: The fact that a drift probability current of the form (7) contributes to the first passage time density was also found by Plesser and Gerstner (2000). Considering infinite refractoriness in (4), the spike density ρ(t) is just the instantaneous rate f (t) provided that the neuron has stayed below threshold before: ρ(t) = f (t)S(t). Solving for f and inserting (7) we obtain: Thus, f , as already observed for ρ in (7), directly depends on the speed at which the membrane potential distribution is pushed towards threshold. f does not explicitly depend on t and can be expressed as a function of U and its temporal derivativeU = dU/dt. The function factorizes into U -andU -dependent terms: f (U,U ) = U + g(U ), (10) with [x] + := (x + |x|)/2. Neglecting the threshold again, p st (V ) can be approximated by a Gaussian (Tuckwell 1988). Figure 3 compares our theory with simulations of the IAF neuron model. The free intensity f (t) resulting from (9) is in excellent agreement with the simulation results. The absorbing boundary at the threshold does not distort the Gaussian shape of the membrane potential density. Additional intensity is observed at the peak and during the descent of U (t). This can be attributed to a diffusive probability current not captured by the theory because even if the membrane potential distribution remains stationary (e.g. U (t) = 0) threshold crossings occur due to the fluctuating nature of V (t). WhenU is small or negative, fluctuations are the only cause of spiking. Note that in the example shown in figure 3 the rise time of U (t) is not much shorter than the membrane time constant. We derived an accurate model of the free intensity and related the expansive shape of the U -dependent term (cf figure 2(d)) to the shape of the stationary membrane potential density p st (V ). The drift probability current captured by a f (U,U )-model as in equation (10) is the main contribution to the observed spike densities generated by synchronous input. Is such a model sufficient to understand synchronization dynamics?

Synchronization dynamics
In the remainder of the paper, we analytically investigate the network dynamics induced by equation (10) using (mV) Figure 3. Free intensity f (t) and probability density of the membrane potential p (V, t) in the presence of a fast voltage transient U (t).
(a) f (t) of the IAF neuron model (solid curve, same data as in figure 2(a)). The gray curve shows an analytical f (t) based on the drift probability current generated by a Gaussian membrane potential density. The theory describes the intensity during the rise time of U (t) (see (b)). Additional intensity is observed at the peak of U (t) (dotted vertical line) and during its descent. (b) Membrane potential density p(V, t) of the IAF neuron model (color bar indicates probability density) including reset to V 0 at the spike threshold θ (white dashed line) and clamping for a refractory period τ r . U (t) (same data as in figure 2(a)) is shown as an offset to the mean µ V (solid curve) of the stationary density p st (V ) and to µ V ± 1σ V (black dashed curves), where σ V is the standard deviation of p st (V ). During the rise time of U (t) the shape of p(V, t) remains intact for V < θ .
where the membrane potential U is measured in units of the threshold relative to the mean θ − µ V , the exponent α 2 introduces a degree of nonlinearity in U , and β > 0 is a scaling constant. It is convenient to measure the number of spikes A in units of (θ − µ V )/û (amplitude of the PSPû = max u(t)). This enables us to write (3) as with the characteristic membrane potential excursion x(t; σ ) = φ(t; σ ) * u(t)/û . Then the free intensity reads wheret specifies the location of the maximum of x(t; σ ). We consider the case that refractoriness is long, allowing each neuron to spike at most once. Thus the spike density 8 ρ(t) is given by equation (8). Writing the survival probability as S(t) = exp (− t −∞ f (s) ds) (Cox 1962), the response probability a out reads withã = f (t) dt. From equation (12), we findã = β[Ax(σ )] α withx(σ ) = max t x(t; σ ). Next, we show that synchronization σ out < σ occurs under quite general conditions. Because of refractoriness (σ out σ ), it is sufficient to show thatσ < σ , whereσ is the standard deviation of f (t). For the choice x(t) =x exp (−t 2 /2ζ 2 ), a Gaussian with amplitudex =x(σ ) and standard deviation ζ = ζ (σ ),σ can be found with the help of Gaussian integrals: Therefore, synchronization only depends on the function ζ (σ ) which is well approximated by ζ (σ ) = τ 0 /2 √ 2 ln 2 + σ . The synaptic rise time τ 0 cannot destroy synchronization for input σ larger than some minimal value. Using a different argument, Abeles (1991) noticed that the spread of ρ(t) depends on the slope of U . Herrmann et al (1995) later elegantly extended the argument to an expression for the equilibrium spread. Marsalek et al (1997) derived a relation between input and output spread for the noise-free case.
We are now in a position to describe the propagation of synchronous activity in a network as specified in figure 1 by an iterative map of (a n , σ n ). With w neurons in each group, the map is constructed setting A = wa n , σ = σ n , a n+1 = a out , and σ n+1 =σ in equations (13) and (14), respectively. Figure 4 shows the state space of the iterative map for different values of w. The bifurcation scenario in the four panels agrees quantitatively with simulations of the corresponding IAF neuron model (cf Diesmann et al (1999)). At a given choice of w the trajectories of the model reproduce the synchronization dynamics exhibited by feed-forward subnetworks ( figure 1(b), see also Diesmann et al (1999); Gewaltig et al (2001)). For w > w c an attractor for synchronous activity exists. The bifurcation creating this attractor is determined by the map a n+1 = 1 − exp (−ka α n ), with k = β[wx(σ n )] α . For k k c > 0 an a-isocline with a > 0 exists and its parabolic shape is provided by α 2. The interplay of the a-isocline and the σ -isocline generates a saddle-node bifurcation. α 2 follows from (10) providing a natural explanation of bistability in the system: a nonlinear dependence of f on U is not required. The σ -isocline is determined by ζ (σ ) in (14) and does not depend on A. The synaptic rise time τ 0 limits the precision of synchronization. In simulations the σ -isocline bends towards larger σ for small A , which can be attributed to spikes generated during the descent of U (t).

Conclusions
Taken together, we have shown that in systems with strong input transients, f (U,U ) -models of spike intensity are consistent with IAF dynamics, whereas f (U )-models (Gewaltig 2000;Kistler and Gerstner 2002) are not. The mechanism of synchronization is theU -dependence of f , not refractoriness. A simple relationship guarantees synchronization under general conditions. Thus, we have uncovered the single neuron mechanism governing the attractor dynamics  of synchronous spiking activity in feed-forward subnetworks. TheU -dependence discussed here is solely induced by the noise and not by a Intersections of the σ -isocline (vertical dashed line) and the A-isocline (solid curve, panels for w 90) are fixed points of the system. The fixed point at large A is the attractor for synchronous activity (basin of attraction shaded in gray).
U -dependence of the spike generation mechanism as has been reported in experimental literature. Our findings suggest that f (U,U )-models should be considered in the investigation of neural synchronization phenomena in a noisy environment.