Propagating synchrony in feed-forward networks.

Coordinated patterns of precisely timed action potentials (spikes) emerge in a variety of neural circuits but their dynamical origin is still not well understood. One hypothesis states that synchronous activity propagating through feed-forward chains of groups of neurons (synfire chains) may dynamically generate such spike patterns. Additionally, synfire chains offer the possibility to enable reliable signal transmission. So far, mostly densely connected chains, often with all-to-all connectivity between groups, have been theoretically and computationally studied. Yet, such prominent feed-forward structures have not been observed experimentally. Here we analytically and numerically investigate under which conditions diluted feed-forward chains may exhibit synchrony propagation. In addition to conventional linear input summation, we study the impact of non-linear, non-additive summation accounting for the effect of fast dendritic spikes. The non-linearities promote synchronous inputs to generate precisely timed spikes. We identify how non-additive coupling relaxes the conditions on connectivity such that it enables synchrony propagation at connectivities substantially lower than required for linearly coupled chains. Although the analytical treatment is based on a simple leaky integrate-and-fire neuron model, we show how to generalize our methods to biologically more detailed neuron models and verify our results by numerical simulations with, e.g., Hodgkin Huxley type neurons.


SPIKE PATTERNS AND SIGNAL TRANSMISSION IN NEURONAL CIRCUITS
Reliable signal transmission is a core part of neuronal processing. A common hypothesis states that activity propagating along neuronal sub-populations that are connected in a feed-forward manner may support such signal transmission. Indeed, there is strong indication that activity propagation along feed-forward structures drives the generation of bird songs (Long et al., 2010) and experiments have shown propagation of synchronous and rate activity in feed-forward networks (FFNs) in vitro (Reyes, 2003;Feinerman et al., 2005;Feinerman and Moses, 2006). Sequential replay in the hippocampus and in neocortical networks also suggest underlying feed-forward mechanisms (August and Levy, 1999;Nadasdy et al., 1999;Lee and Wilson, 2002;Leibold and Kempter, 2006;Xu et al., 2011;Eagleman and Dragoi, 2012;Jahnke et al., 2012) and propagation of synchronous activity along feed-forward chains is a possible explanation for experimentally observed precise spike timing in the cortex (Riehle et al., 1997;Kilavik et al., 2009;Putrino et al., 2010). Further, the modular, hierarchical structure of many sensory and motor systems suggests propagation over sequences of areas in feedforward manner, e.g., in bottom-up signal transfer (Felleman and Van Essen, 1991;Scannell et al., 1999;Bullmore and Sporns, 2009;Kumar et al., 2010). Feed-forward structures which support the propagation of synchronous activity are termed synfire chains. The concept was introduced by Abeles (1982) as groups of neurons (layers) with dense anatomical connections between subsequent groups that are embedded in otherwise roughly randomly connected local neural circuits. Two major questions regarding the dynamical options for synfire activity include a) how synchrony may actively propagate and b) how such spatio-temporally coordinated spike timing may be robust against irregular background activity, because the synfire chains are part of a cortical network with dynamics defined by the so-called irregular balanced state Sompolinsky, 1996, 1998).
Addressing these points, theoretical studies have established conditions for stable propagation of synchrony in synfire chains (Diesmann et al., 1999;Gewaltig et al., 2001). Most synfire chain models assume functionally relevant FFNs that exhibit a very dense, often all-to-all connectivity between subsequent layers (Aviel et al., 2003;Mehring et al., 2003;Kumar et al., 2008) (see also a recent review on this topic Kumar et al., 2010). Such highly prominent feed-forward-structures, however, have not been found experimentally. Since cortical neural networks are overall sparse (e.g., Braitenberg and Schüz, 1998;Holmgren et al., 2003), we may also expect some level of dilution for embedded feed-forward chains. So far, computational model studies assumed that such chains created from existing connections in sparse recurrent networks exhibit strong synaptic efficiencies and specifically modified neuron properties to enable synchrony propagation (Vogels and Abbott, 2005).
Recently, we have shown that non-additive dendritic interactions promote propagation of synchrony (Jahnke et al., 2012). The non-additive dendritic interactions considered are mediated by fast dendritic spikes (Ariav et al., 2003;Gasparini et al., 2004;Polsky et al., 2004;Gasparini and Magee, 2006): upon stimulation within a time interval less than a few milliseconds, dendrites are capable of generating sodium spikes. These induce a strong, short and stereotypical depolarization in the soma. If this depolarization elicits a somatic spike, the spike occurs a fixed time interval after stimulation with sub-millisecond precision. This dendritic non-linearities relax the requirement of dense feed-forward anatomy and thereby allow for robust propagation of synchrony even in diluted FFNs with synapses of moderate strength within the biologically observed range.
In the present article, we analytically and numerically investigate in detail under which conditions synchronous activity may reliably propagate along the layers of an FFN where the inter-group connectivity is diluted, as may be expected when they are part of a sparse cortical network. An embedding network is mimicked by external, noisy input. We study the influence of the network setup, including the influence of the emulated embedding network, and of different types of standard linearly additive as well as non-additive dendritic interactions.
We derive analytical estimates for the critical connectivitythe minimal connectivity that allows robust propagation of synchrony. Some fundamental analytical results, in particular the ansatz for deriving a critical connectivity in the first place, have been briefly reported before (Jahnke et al., 2012). Here, we extend the approach and show how the bifurcation point, i.e., the transition point from the non-propagating to the propagating regime, can be estimated quantitatively from the neurons' ground state properties. We investigate the validity range of the analytical predictions and check them via direct numerical simulations. Furthermore, we discuss the applicability of our results to biologically more detailed neuron models and network setups. In particular, we argue that the assumptions underlying the analytical approach are met by a wide class of neuron models, including, e.g., conductance based leaky integrate-and-fire and Hodgkin-Huxley-type neurons.
The article is structured as follows: After introducing the neuron model and network setup in section 2, we study in the main part the propagation of synchrony in linearly coupled FFNs (section 3.1) and in FFNs incorporating dendritic non-linearities (section 3.2). In particular, we derive tools to study the system analytically, compare the results to computer simulations and elaborate differences of the dynamics of FFNs with and without non-additive dendritic interactions. In the final part (section 3.3), we discuss the application of our analytical results to biologically more detailed neuron models.

Linear model
Consider networks of leaky integrate-and-fire neurons that interact by sending and receiving spikes via directed connections. The state of neuron k at time t is described by its membrane potential V k (t) and its dynamics satisfy where τ m k is the membrane time constant of neuron k, I const k := I 0 k /τ m k a constant input current, I net k (t) the input current caused by spikes within the network and I ext k (t) the input current arising from spikes from external sources. When the neuron's membrane potential reaches or exceeds the threshold k its membrane potential is reset to V reset k and a spike is sent to the postsynaptic neurons n, where it changes the postsynaptic potential after a delay τ nk . After emitting a spike at t = t 0 the neuron becomes refractory for a time period t ref , i.e., V k (t) = V reset To keep the model analytically tractable, we model the fast rise of the membrane potential upon the arrival of presynaptic spikes by instantaneous jumps of the membrane potential, such that the resulting input current reads Here kl denotes the coupling strength from neuron l to neuron k, t f lm is the mth spike time of neuron l and τ kl specifies the synaptic delay. In addition to spikes from the network each neuron receives excitatory and inhibitory random inputs that emulate an embedding network. These external inputs are modeled as random Poisson spike trains with rate ν exc and ν inh , respectively. The resulting input current is given by where t ext, exc km (t ext, inh km ) is the arrival time of the mth excitatory (inhibitory) spike to neuron k and exc > 0 ( inh < 0) denote the corresponding coupling strength.

Non-linear model
In the above model all input currents are summed up linearly. To also investigate the effect of dendritic spikes we modulate the sum of synchronously arriving excitatory inputs by a non-linear dendritic modulation function σ NL (·). This can be directly read off from experimental data (Ariav et al., 2003;Gasparini et al., 2004;Polsky et al., 2004;Gasparini and Magee, 2006): If the sum of excitatory inputs is below the dendritic threshold b , the single inputs are processed linearly (σ NL (·) equals the identity). If the sum of inputs exceeds the dendritic threshold b , the depolarization is strongly non-linearly enhanced compared to that expected from linear summation. This is, in biological terms, due to a dendritic spike elicited. Larger inputs have been experimentally found to not further increase the somatic peak depolarization. The dendritic modulation function may then be modeled as The dendrites process synchronous inputs non-additively: inputs below the dendritic threshold are summed linear, inputs above this threshold are summed supra-linear and, due to the saturation, very large inputs are summed sub-linear.
If not stated otherwise, we consider only exactly simultaneous arriving spikes as sufficiently synchronous; to allow for exactly simultaneous arrivals, the synaptic delays are chosen as homogeneous τ kl ≡ τ. The input currents caused by spikes that are received from the network are then given by Here, the sum over t f denotes the sum over all times at which spike(s) are sent in the network, irrespective of which neuron(s) is (are) spiking. The sets M exc t f and M inh t f specify the set of neurons that send an excitatory or inhibitory spike at time t f , respectively. (To describe a network with linear dendrites σ NL ( ) is replaced by ).
In section 3.3.1 we consider inhomogeneous delay distributions and finite dendritic integration window t (i.e., non-linear amplification of inputs received within finite time interval t) and discuss how the results achieved for homogeneous systems can be generalized.

NETWORK TOPOLOGY
We consider the propagation of synchrony in diluted Feed-Forward-Networks (FFNs, synfire-chains). They consist of a sequence of m layers, each composed of ω neurons. Neurons of one layer form excitatory projections to the neurons of the subsequent layer with probability p; the strength of an existing connection from neuron l to neuron k is denoted by kl .
For simplicity of presentation, we consider homogeneous neuronal populations, i.e., all neurons have identical properties (τ m k = τ m , k = and V reset k = V reset for all i), as well as homogeneous coupling strengths, i.e., kl = if a connection is realized, throughout this article. If not stated otherwise, we use τ m = 14 ms and = 15 mV as standard values for the membrane time constant and the neuron threshold.

GROUND STATE DYNAMICS
We consider networks, where the single neurons are placed in a "fluctuation driven regime," i.e., in the ground state the average input to each neuron is sub-threshold and spiking of neurons is caused by fluctuations of the inputs. This setup allows to emulate the dynamics of neurons which are part of a balanced network Sompolinsky, 1996, 1998). The neurons fire asynchronously and irregularly with low firing rate ν; the spike trains resemble Poissonian spike trains (Tuckwell, 1988;Brunel and Hakim, 1999;Brunel, 2000;Burkitt, 2006). Thus, the inputs to the neurons may be described by three Poissonian spike trains with rates ν exc (external, excitatory), ν inh (external, inhibitory) and ν int = νpω (inputs from the preceding layer). Since the number of inputs N X T , X ∈ {exc, inh, int}, in a time interval T is Poisson distributed, the expected number of inputs N X T and the variance N X Then μ = I 0 + τ m ν exc exc + τ m ν inh inh + τ m pων (6) is the mean of the total input to the neurons in an interval of the size of the membrane time constant, T = τ m , and is its variance. In diffusion approximation, the distribution of membrane potentials P V (V) and the mean firing rate ν can be derived analytically (Brunel and Hakim, 1999;Brunel, 2000;Helias et al., 2010). In particular, for networks with low firing rates the probability density of membrane potentials (see, e.g., Tuckwell, 1988) is Gaussian and can be expressed in terms of the input current. In this approximation the average firing rate is and depends on μ and σ only via the quotient which is the distance of the average input μ from the neurons' threshold normalized by the standard deviation σ of the input. For the analytical derivations throughout this article we focus on the regime of low spiking rates α 2; ν 1.5Hz . In the absence of synchronous activity each neuron receives a large number of inputs from the external network and only a few inputs from the previous layer of the FFN, such that the ground state dynamics of the network is mainly established by the external inputs. To keep the input balanced we choose ν exc = ν inh =: ν ext and exc = − inh =: ext throughout the article.

PROPAGATION OF SYNCHRONY
To initiate propagating synchronous activity along the considered diluted FFN, we excite in the first layer a subgroup of g 0 ≤ ω neurons to spike synchronously. This causes a synchronous input to the following layer after the synaptic delay τ and may therefore initiate synchronous spiking of a subgroup of neurons in that layer. These may again excite synchronous spiking in the next layer and so on. Depending on the ground state, i.e., the layout of the external network, on the layer size ω, and on the coupling strength , a synchronous pulse may or may not propagate along the FFN (cf. Figures 1A,B,D,E).
In addition to the triggered propagation, one might generally also expect the occurrence of spontaneous propagation of synchronous activity: Neurons of a particular layer share inputs from the previous layer and this causes correlations in their spiking  activity. Over the layers these correlations can accumulate and lead to synchronous spiking (Aviel et al., 2003;Rosenbaum et al., 2010Rosenbaum et al., , 2011Litvak et al., 2013). However, in the setups considered in this article, the effect is negligible due to two reasons: (1) each neuron receives a large number of external (uncorrelated) inputs and this background noise has a decorrelating effect, (2) we study the system near the critical point, i.e., for parameters where even synchronized spiking of all neurons of a particular layer is just sufficient to initiate a propagation of synchrony. Thus, spontaneous propagation of synchrony effectively does not occur.
We study the transition from the non-propagating to the propagating regime by means of a iterated map that yields the expectation value of the number of synchronously spiking neurons g i + 1 in layer i + 1 if g i neurons are synchronously active in layer i. There is always one trivial fixed point, G 0 , of this iterated map with 0 = G 0 = g i + 1 = g i , which corresponds to absent activity. If g i + 1 < g i for all g i > G 0 , synchronous activity will die out after a small number of layers. If g i + 1 ≥ g i for some substantial group size, g i > G 0 , a stable propagation of synchrony may be enabled (cf. Figures 1C,F). More precisely, we will show in this article that with increasing connectivity p the system undergoes a tangent bifurcation and two fixed points G 1 and G 2 ≥ G 1 appear. If existing, G 1 is always unstable (the diagonal is crossed from below; the slope of the iterated map needs to be larger than one) and G 2 is always stable [all connections within the FFN are excitatory such that the iterated map is monotonically increasing (slope larger than zero, in particular larger than −1)]; further at G 2 there is an intersection with the diagonal from above thus the slope is smaller than one and stationary propagation with group sizes around G 2 is enabled.
In computer simulations, we determine for each given network setup by the following procedure whether a propagation is possible: after some initial time t init we excite all neurons of the first layer to spike synchronously and measure the number of active neurons g i in the ith layer at the expected spiking time t exp i = t init + iτ. If g i is substantially larger than the number of active neurons arising from spontaneous activity in more than 50% of n trials (i.e., n repetition of the same simulation with different initial conditions), we denote the propagation of synchrony as successful. The critical connectivity p * , that marks the transition from a regime where propagation of synchrony is not possible to a regime where propagation of synchrony is enabled, is found by determining the lowest connection probability p for which an initial synchronous pulse propagates successfully.
As the connections within the FFN are all excitatory, it is sufficient to check whether propagation of synchrony can be initiated by inducing synchronized spiking of all ω neurons of the first layer: Stationary propagation of synchrony can be enabled if there is a non-trivial stable fixed point (G 2 ) of the iterated map for the average group size. For purely excitatory connections the basin of attraction of this fixed point is bounded from the left by an unstable fixed point (G 1 ) and from the right by the maximum group size given by the layer size ω.

RESULTS AND DISCUSSION
Under which conditions can synchronous signals propagate robustly along diluted FFNs? To answer this question in detail, we first focus on networks with linear dendrites. Afterwards we study the propagation of synchrony in networks incorporating non-additive dendritic interactions and compare with the linear case. Finally, we show that the derived results are directly applicable in biologically more detailed neuron models and network configurations.

FFNs WITH LINEAR DENDRITES
In this section, we consider linearly coupled FFNs. In the first part, we derive analytical estimates for the critical connectivity p * L that marks the transition from the non-propagating to the propagating regime; the initial steps follow the lines of Jahnke et al. (2012); Memmesheimer and Timme (2012). In the second part we investigate the influence of the external network on the propagation of synchrony and determine the parameter-region for which the analytical estimates are applicable. In particular, we show that the derived estimates are applicable in the biologically relevant parameter-region, where the spontaneous firing rate is low and the distribution of membrane potentials is sufficiently broad. Finally, we study how the properties of propagating synchronous pulses depend on different system parameters.

Analytical derivation of critical connectivity
To access the properties of propagation of synchrony we consider average numbers of active neurons in the different layers of an FFN: for this, we derive a iterated map which yields the expected number of neurons that will spike synchronously in one layer given that in the preceding layer a certain number of neurons was synchronously active.
If in the ith layer, g i neurons spike synchronously, the number of synchronous inputs h a single neuron in layer i + 1 receives follows a binomial distribution h ∼ B g i , p . We denote the spiking probability of a single neuron due to an input of strength x by p f (x). The average or expected spiking probability p sp g i of a single neuron in layer i + 1 is then given by Here and in the following we denote the expectation value of a function f (X) of a random variable X by E[f (X)]; conditional expectations are denoted by E[f (X)|Y]. The expected number of spiking neurons in layer i + 1 is then simply If the connection probability p is low and/or the connection strengths are small, the spontaneous spiking activity in the absence of synchrony is only weakly influenced by the spiking activity within the FFN. Thus as a starting point, we assume that the ground state is exclusively governed by external inputs (effectively setting ij ≡ 0). Then, the mean input to the neurons in an interval of length τ m is μ = I 0 with standard deviation σ = ext √ 2τ m ν ext (cf. section 2.3). Using the probability density (Equation 8), we calculate the spiking probability of a single neuron, p f (x), due to an input of strength x; equals the probability of finding a neuron's membrane potential in the interval [ − x, ]. To derive a iterated map for the average number of active neurons (which maps E[g i ] → E[g i + 1 ]), we interpolate E g i + 1 g i for continuous g i and in the second step replace g i by its expectation value E g i . The fixed points, E g i + 1 E g i = E g i , qualitatively determine the propagation properties of synchronous activity. In the rest of the manuscript we are dealing with the average number of active neurons in a given layer. Therefore, for simplicity we denote the expectation value of the average number of active neurons in a given layer i by g i instead of E g i . For sufficiently small connection probabilities p the map (Equation 12) has only one (trivial) fixed point G 0 = g i + 1 = g i = 0. Any initial synchronous pulse will die out after a small number of layers (see also Figure 1). With increasing connectivity two additional fixed points G 1 (unstable) and G 2 ≥ G 1 (stable) appear via a tangent bifurcation.
For FFNs with purely excitatory couplings between the layers, the second fixed point G 2 (if it exists) is always stable: The spiking probability p f (x) is monotonically increasing with input x and thus also the iterated map (Equation 13) is monotonically increasing (i.e., the slope is larger than 0). Moreover, if G 2 exists the slope of the iterated map at this intersection point with the diagonal is smaller than 1. This implies that G 2 is stable and synchronous pulses of size g i ≥ G 1 typically initiate a propagation of synchrony with an average number of active neurons around G 2 . The critical connectivity p * L at the bifurcation point marks the minimal connectivity that allows for stable propagation of synchrony.
Although the distribution of inputs from one layer to the subsequent one and the spiking probability of a single neuron p f (·) are known, there is no analytic closed form solution to the fixed point equation In other words, we can compute the firing probability p f (x 0 ) for any x 0 , and therefore also E g i + 1 g i for any g i , but g * i = E g i + 1 g * i is transcendental. We thus derive an approximate solution. We choose some expansion point g i (see section 3.1.2 for details), and approximate the function E g i + 1 g * i by a polynomial This also allows to analytically compute the critical connectivity p * L : it is the parameter value at which the iterated map undergoes a tangent bifurcation, i.e., at which the two solutions of the fixed point equation become equal upon changing from complex-conjugate to real. Since the right hand side of Equation (13) does not offer itself for a direct series expansion in g * i , we derive g i + S(g * i − g i ) from an appropriate expansion of p f (h ) and a subsequent computation the arising expectation values.
In biologically relevant scenarios, the neurons usually receive a large number of synaptic inputs and thus the distribution of Frontiers in Computational Neuroscience www.frontiersin.org November 2013 | Volume 7 | Article 153 | 5 membrane potentials P V (V) is broad, P V (V) changes slowly with V. Then, P V (V) around some V = V 0 can be approximated by considering a series expansion with a small order and it is possible to derive an approximation for the critical connectivity p * L based on an expansion of p f (·). Expanding p f (x) into a Taylor series around some x 0 and using Equation (12) yields Here and in the following we denote the nth derivative of a function f (x) at x = x 0 by Replacing the derivatives of p f (·) by the (one order lower) derivatives of probability density of membrane potentials P V (V) according to Equation (14) yields where we defined for better readability.
We have recently shown (Jahnke et al., 2012) that it is possible to derive a scaling law for the critical connectivity using the (unknown) average input from one layer to the next during stationary synchrony propagation, as expansion point. For this choice the expectation value E (h − x 0 ) n g i in Equation (19) simplifies to where we denote by m n the nth central moment of the Binomial distribution B g i , p , specifying the distribution of inputs to the (i + 1)th layer. In the limit of large layer sizes ω and small coupling strengths keeping the maximal input ω to each layer constant (to preserve the network state), all summands for n ≥ 2 vanish, and Equation (19) simplifies to Using the implicit function theorem one can show that this implies the scaling law where λ is a constant independent of and ω (Jahnke et al., 2012). We note that for the derivation of the scaling law (Equation 24) we did not use the actual functional form of the distribution of membrane potentials P V (V). Therefore this estimate holds if P V (V) is sufficiently slow changing with V such that the Taylor expansion (cf. Equation 16) is applicable, but its validity is not restricted to the low-rate approximation. However, the dependence of the prefactor 1/λ on the layout of the external network remained unknown. Here, we present an approach that enables us to derive an approximate value for λ. We consider the expansion (Equation 19) around x 0 up to second order, The truncated series (Equation 25) is quadratic in g i such that the fixed points g * 1/2 = g i + 1 = g i can be obtained analytically, where we defined At the bifurcation point, the root in Equation (26) vanishes such that both fixed points agree (g * 1 = g * 2 ) and γ L = g * 1 = g * 2 specifies the average size of a propagating synchronous pulse. Consequently, the critical connectivity is obtained by choosing p such that which yields where we defined which is independent of the setup of the FFN and completely determined by the layout of the external network and the choice of the expansion point x 0 . As before we consider the limit of large layer sizes ω and small coupling strengths , i.e., we replace ω → const and consider the leading terms of a series expansion of Equation (29). The expansion of the square bracket in Equation (29) yields such that the critical connectivity assumes the functional form given by Equation (24), Thus λ = λ * defined by Equation (30) provides an approximation of the constant λ fully specifying the critical connectivity p * L .

Optimal expansion point
To derive Equation (30) we assumed that it is sufficient to consider the second order expansion of p f (x). It is thus necessary to choose an appropriate expansion point that results in fast convergence. In particular for the choice x 0 = x * 0 , that we will now derive, Equation (37) below, the bifurcation diagram near the bifurcation point is well approximated already for k = 2 (cf. Figure 2).
The size of a propagating group at the critical connectivity is γ L (cf. Equation 27) and thus the resulting average input is p * L γ L . Our expansion point x 0 should lie near to this value, which is, of course, unknown prior to solving the fixed point equation. We will thus compute a range in which p * L γ L has to lie and choose the expansion point appropriately within. We assume that ω is large and employ Equation (23)  ) truncated after expansion order k (color code) with x 0 = x * 0 (cf. Equation 37). (B) Fixed points of the iterated maps shown in (A); with increasing connectivity two fixed points appear by a saddle node bifurcation. We note that already a second order expansion (red), i.e., the lowest order at which a saddle node bifurcation can occur, approximates the bifurcation diagram (blue) near the bifurcation point well.
of this range as we know the functional form explicitly. Equation (23) with g i + 1 = g i is just another transcendental equation for the fixed points and it has zero, one, or two non-trivial fixed point solutions points g * 1 and g * 2 , which are then also solutions of Equation (19) where the diagonal is touched, the function p f (gp ) has to be concave and monotonic increasing with respect to g. The such that the bifurcation point satisfies The condition Equation (33) holds because P V (V) ≥ 0 is a probability density and Equation (34) is derived directly from differentiating Equation (8). To maximize the quality of the second order approximation Equation (25), we choose x 0 = x * 0 such that the contribution to the expansion (Equation 19) of the k = 3rd order term equals zero. According to Equation (19), all 3rd order terms are proportional to P In the considered regime of low spiking rates, we find For x 0 = x * 0 the bifurcation diagram near the bifurcation point is well approximated already for k = 2 (cf. Figure 2) and Equation (30) provides a good estimate of the critical connectivity p * L (cf. Figure 3).

Influence of external network
In the previous section we derived an iterated map for the average group size (cf. Equation 13) and an approximation for the critical connectivity p * L (cf. Equations 30 and 32) that marks the transition from FFNs which do not support propagation of synchrony to FFNs that do. In this section we focus on the robustness of our results. How does the critical connectivity change with the layout of the external network? For which parameter range does the estimate of the critical connectivity (given by Equations 30 and 32) yield reasonable results?
The derivation was based on the assumption that the ground state dynamics of the neurons of the FFN is completely determined by the external inputs. This assumption holds if the spontaneous firing rate ν of the neurons and/or the coupling strengths In the main panels we use a logarithmic scale, the insets have a linear scale. The squares indicate the connectivity above which a synchronous pulse propagates from the 1st to the 20th layer of a FFN in at least 50% of n = 30 trials. The critical connectivity given by Equation (32) ) is in good agreement with computer simulations. As predicted p * L ∝ ( ω) −1 and the proportionality factor 1/λ is well approximated by the estimate 1/λ * derived in Equation (30). and/or the connectivity p are sufficiently small. We will generalize our approach and show how the impact of preceding layers on a layer's ground state can be taken into account. Thereafter we will compare the results with computer simulations, identify the regions in parameter space for which the derived approximations hold and discuss deviations between direct numerical simulations and analytics.
The first layer of an FFN receives inputs only from the external network and according to Equations (6, 7) the mean μ 1 and standard deviation σ 1 of its input is as assumed in the previous section. All following layers receive external inputs and spikes from their preceding layer(s). The mean μ n and standard deviation σ n of the input to neurons of the nth layer (with n ≥ 2) reads (cf. Equations 6 and 7) Here we denote the spontaneous firing rate (in the absence of synchrony) of neurons of the (n − 1)th layer by ν n − 1 . It is given by Equation (9) as From layer to layer, the mean input, the standard deviation as well as the firing rate increase. For setups, where the ground state of the FFN is non-pathological, i.e., the firing rates of all layers are bounded, the additional corrections X n := X n − X n − 1 for X ∈ {μ, σ, ν} decrease with n, and μ n , σ n and ν n saturate for sufficiently large n. Thus, μ ∞ and σ ∞ describe the input to the neurons of an infinitely long FFN and the single neurons of such an FFN spike with an average rate ν ∞ . Accordingly, replacing μ and σ by μ ∞ and σ ∞ in Equation (13) [where they appear as parameters of p f (·)] yields an iterated map for the average group size.
In Figure 4, we compare the critical connectivity found by numerically determining the bifurcation point of the iterated map (Equation 13) (i.e., we determined the connectivity p for which the iterated map touches the diagonal; solid lines) with computer simulations of propagating synchrony (markers). To also cover scenarios, where the input from the preceding layer is not negligible, we consider infinitely long FFNs (then, the distribution of membrane potentials is equal in all layers). In computer simulations this can be approximated by a sufficiently long FFN with periodic boundary conditions, i.e., an FFN where the last layer connects to the first layer. For moderate external inputs, i.e., moderate I 0 and ext , already the analytical results neglecting the influence of the preceding layers (using μ 1 and σ 1 ) agree well with computer simulations (cf. Figure 4A, solid lines). However, for large external inputs, i.e., large I 0 and ext , the critical connectivity is overestimated. Here, the assumption that the distribution of membrane potentials is not influenced by the connectivity of the FFN does not hold. The additional input shifts the membrane potentials to higher values and consequently a lower connectivity is required for a propagation of a synchronous pulse. The corrections given by Equations (38-42) account for these deviations to some extent (cf. Figures 4B,C; solid lines), in particular for setups where the spontaneous firing rate is low. However, for very large I 0 and ext , the critical connectivity is under-estimated. Here, the spontaneous firing rate is too high and the low-rate approximation, Equations (8-9), is not adequate to describe the system; the firing rate and thus the mean input from the previous layer are over-estimated. This becomes particularly clear in Figure 4C, where we show the critical connectivity as a function of the strength of the external inputs ext . For any given I 0 (different colors), the critical connectivity for small ext is well approximated; with increasing ext the firing rate increases [α decreases and thus ν increases; cf. Equations (9 and 10)] and when the coupling strengths ext exceed a I 0 -dependent threshold, the low-rate approximation becomes inapplicable.
Applying the methods in Brunel and Hakim (1999);Brunel (2000), the firing rate and the distribution of membrane potentials can be derived in diffusion approximation for states with higher spontaneous firing rates. Although most of the analytical considerations above are also applicable within this approximation, the determination of an optimal expansion point (cf. Equations 36 and 37) becomes more difficult and a closed form expression does not exist. However, the critical connectivity can be obtained by numerically determining the fixed points of the iterated map (Equation 13) and we find that it agrees with  (Equation 13). In panel (A) we neglect the influence of previous layers on the ground state of a considered layer in the analytical computations [i.e., we use μ 1 and σ 1 , cf. Equations (38) and (39)]. In (B,C) we employ corrections to account for their influence, cf. Equations (38-42). We show the third order correction, higher orders add only small modifications to the curves, but the numerical computations get more costly. The thick gray lines in (B,C) indicate the bifurcation point of the iterated map (Equation 13) with P V (V ) derived from the diffusion approximation of leaky integrate-and-fire neuron dynamics with Poissonian input (Brunel and Hakim, 1999;Brunel, 2000). The dashed lines are the estimates of the critical connectivity given by Equations (30 and 32). Again, in panel (A) we neglect the influence of previous groups on the ground state, in panels (B,C) we use the third order correction. The estimates agree with the data from numerical simulations within the biologically relevant parameter range, where (1) the spontaneous spiking activity is low and (2)  Analogous to the approach presented above, corrections for the influence of preceding layers can be taken into account for the analytical estimate of the critical connectivity derived in the previous section (Equations 30 and 32). Replacing the connectivity p by the approximation p * L = (λ * ω) −1 in Equations (40, 41) yields where λ * n − 1 := λ * (μ n − 1 , σ n − 1 ) is given by Equation (30) and ν n − 1 = ν (μ n − 1 , σ n − 1 ) is given by Equation (42). In Figure 4 we show the estimate of the critical connectivity p * L = λ * n ω −1 (cf. Equation 32) using λ * 1 (panel a; dashed line), i.e., neglecting the influence of the preceding layers, and using a higher correction order (panel b,c; dashed line: third order). For sufficiently large ext the critical connectivity found by numerically determining the bifurcation point agrees with the analytical estimate given by Equation (32). As discussed above, the corrections Equations (43, 44) account for the deviations from the simulated data as long as the total spontaneous firing rate is sufficiently low. However, for small ext the critical connectivity is under-estimated. Here, the standard deviation of the inputs (cf. Equation 7) is low, such that the distribution of membrane (8)], the spiking probability of one neuron, p f (·), increases steeply in a small interval [for ext → 0: p f (x) → (x − μ); cf. Equation (8)] and thus the approximation of p f (·) by the leading terms of a Taylor expansion is not sufficiently accurate. However, in the biologically plausible parameter regime, where the firing rates are small and the distribution of membrane potentials is broad, the critical connectivity is approximated well by Equation (32) together with Equation (30) (defining λ * ), Equation (37) (defining x * 0 ) and the corrections that account for the influence of the preceding layers, Equations (43, 44).

Characteristics of propagating synchronous pulses
In the previous sections, we have shown that a synchronous pulse may propagate along a diluted FFN. In this section we study the characteristics and properties of a propagating synchronous signal. We consider them at the transition to stable propagation, p * L , because there they depend only weakly on the network setup. How large is the fraction of neurons that participate in propagating synchrony? How does this fraction depend on the network setup?
To answer such questions, we consider the effect of a propagation synchronous pulse on the single layers in the network, as a measure for the effective pulse size. In other words, we consider the mean input μ L a neuron receives from the preceding layer if a synchronous pulse propagates along the FFN at the critical connectivity p * L . It is given by the product of the connection probability p * L , the connection strength and the average size of a propagating synchronous signal γ L ; using Equations (27) and (29) yields and after inserting λ * as given by Equation (30), According to Equation (46) are fully determined by the external inputs (I 0 , ν ext and ext ). Figures 5A,B illustrates the dependence of μ L on the layout of the external network and the FFN: as expected from our analytical considerations, the dependence on the layer size and coupling strength is weak when I 0 and ext are kept fixed. With increasing mean of the external input (I 0 ) the distribution of membrane potentials P V (V) is shifted toward the threshold , such that it is more likely to find the membrane potential of the neurons near the threshold and the critical connectivity decreases (cf. also Figures 4A,B). Naturally this implies a decreasing average input μ L at p * L , which is shown in Figure 5A for different external couplings ext and parameters of the FFN. Increasing the external coupling strength ext (and with it the variance of the external input) causes a broadening of the distribution of membrane potentials; the membrane potentials of some neurons are shifted toward the threshold and the membrane potentials of other neurons are shifted away from it. If the fraction of neurons that participate in the propagation of the synchronous pulse is large, this implies an increasing critical connectivity (Figure 5B; cf. also Figure 4C).
The spiking probability of a single neuron due to the mean input μ L equals the average fraction p frac of neurons of one layer that participate in a propagating synchronous pulse, Interestingly, in the considered regime of low spiking rates and sufficiently broad distribution of membrane potentials, where the approximations given in section 3.1.1 are applicable, p frac depends on the setup of the external inputs only via the quotient α = − μ σ (cf. Equation 10), or, equivalently, on the spontaneous firing rate ν of the neurons (cf. Equation 9). This can be shown by combining Equations (8, 37) and (Equation 46),  (A,B) show the mean input μ L that a layer receives due to a propagating synchronous pulse in the preceding layer. μ L measures the effective pulse size (the impact of a propagating synchronous pulse) and is mainly determined by the external inputs rather than by the setup of the FFN. In (A) the variance of the external input (measured by ext ) is fixed and μ L is plotted vs. such that In Figure 5C we compare the above predictions with direct numerical simulations: For different layer sizes ω, coupling strengths and layouts of the external networks (i.e., different values of I 0 and ext ), we detect whether propagation of a synchronous pulse is possible and if so, we numerically determine the average fraction of participating neurons as well as the spontaneous firing frequency. We find that indeed the size of the synchronous pulse is determined essentially by the quotient α = − μ σ and Equation (53) is a reasonable estimate of the average fraction of neurons spiking in each layer. With increasing α the fraction of participating neurons increases, it thus decreases with spontaneous firing rate ν see Figure 5C. For FFNs with low spontaneous spiking frequency almost all neurons of a layer participate in the propagation of a synchronous pulse.

FFNs WITH NON-LINEAR DENDRITES
In this section, we investigate propagation of synchrony mediated by dendritic non-linearities. Although the mechanism underlying the propagation is generally related to that in linear networks, the discontinuities introduced by non-additive dendritic interactions prevent a similar analytical approach. In the first part of this section, we thus derive analytical estimates for the critical connectivity p * NL in non-linearly coupled networks based on a self-consistency approach (see also Jahnke et al., 2012). In the second part, we study the transition from propagation of synchrony mediated by linear dendrites to propagation of synchrony mediated by non-additive dendritic interactions upon increasing the degree of non-linearity in the networks. In the last part, we evaluate the robustness of the analytical estimates with respect to the layout of the external network.

Analytical derivation of critical connectivity
Neurons with non-additive dendritic interaction process excitatory input by a non-linear dendritic modulation function σ NL (see section 2.1), i.e., synchronous inputs that exceed the dendritic threshold b are amplified to an effective input of size κ (cf. Equation 4). Therefore the spiking probability of a single neuron due to a synchronous input of strength x, p f (σ NL (x)), is discontinuous and an approach based on expansion of p f (·) is inappropriate. To derive an analytical expression for the critical connectivity p * NL in FFNs incorporating dendritic non-linearities, we consider the (average) fraction of neurons of one layer, p γ , that receive an input x larger than the dendritic threshold, x ≥ b , due to the propagating synchronous pulse. If there is a stable (stationary) propagation of synchrony established, p γ is constant throughout the layers, which allows us to formulate a self-consistency equation. The basic derivations have been published recently (Jahnke et al., 2012) and will be briefly reviewed in the following for the readers convenience.
For sufficiently small dendritic thresholds b and sufficiently large κ, the spiking probability of a neuron due to a sub-threshold input is small compared to the spiking probability of a suprathreshold input. Therefore, we approximate the spiking probability of a single neuron in response to a synchronous input of strength x by i.e., we assume that somatic spikes due to the synchronous pulse are exclusively generated by dendritically enhanced inputs. We denote the fraction of neurons that receive a dendritic spike by p γ . This may be considered as constant throughout the different layers if stable propagation of synchrony is enabled. Then the probability that a neuron receives exactly k inputs from the preceding layer follows a binomial distribution k ∼ B ω, p γ p f (κ) p , where p γ p f (κ) p is the probability that (1) a neuron of the preceding layer receives a supra-threshold input (p γ ), (2) a somatic spike is elicited by that input p f (κ) and there is a connection from this spiking neuron to the considered neuron of the following layer (p). So we can formulate the self-consistency equation for p γ , To solve Equation (55) we approximate the binomial distribution by a Gaussian distribution with mean δ := ωp γ pp f (κ) and standard deviation σ δ := δ(1 − p γ pp f (κ)), which yields where we defined as the difference between the average number of inputs (δ) and the number of inputs needed to reach the dendritic threshold ( b / ) normalized by the standard deviation of the number of inputs (σ δ ). Solving definition (Equation 58) for p and replacing p γ by Equation (56) yields which is the connectivity p NL where stable propagation of synchrony with some given n (or, equivalently, some given p γ ; cf. Equation 56) is established. We note that a propagation of synchrony mediated by dendritic spikes requires (otherwise even the input caused by a synchronized spiking of all neurons of a layer in a fully connected FFN (p = 1) is not sufficient to reach the dendritic threshold b ). For parameters fulfilling the inequality (Equation 60), p NL (n) has a global minimum (see Appendix) and the critical connectivity p * NL , again defined as the smallest connectivity that allows for a stable propagation of synchrony, matches that global minimum: any connectivity p NL above the minimal connectivity p * NL has two preimages n 1 and n 2 corresponding to the both nontrivial fixed points G 1 and G 2 of the iterated map for the average group size (cf. Figure 1 and section 2.4). However, there exists smaller connectivities for which a stationary propagation can be established. At the global minima p * NL both preimages n 1 and n 2 collapse to n * = n 1 = n 2 and correspond to the fixed point G = G 1 = G 2 of the iterated map at the bifurcation point of the tangent bifurcation. Here the transition from the regime where no propagation of synchrony is possible to the regime where a propagation of synchrony is enabled takes place. For p NL smaller than p * NL there are no preimages (i.e., a stationary propagation of synchrony mediated by non-additive dendritic interactions cannot be established); this scenario correspond to the absence of the non-trivial fixed points of the iterated map for connectivities below the tangent bifurcation.
In the following we will obtain the minima of p NL (i.e., the critical connectivity p * NL ) in the limit of large layer sizes ω and small coupling strength . We first derive an approximation of Equation (59) (cf. Equation 62), determine the validity range of this approximation (cf. Equation 69) and finally obtain an estimate for the critical connectivity (cf. Equation 71). As before, we fix the maximal input ω to each neuron to preserve the network state and expand Equation (59) in a power series around → 0 and ω → ∞. Considering the leading terms yields Further a propagation mediated by dendritic spikes (as introduced above) requires that the layer size ω and the coupling strength are sufficiently large such that a sufficiently large fraction of neurons of each layer receive a total input larger than the dendritic threshold b . In particular for diluted FFNs, this requirement translates to ω b and Equation (61) simplifies further to Whereas p NL has always a global minimum for ω > b , this does not hold for the approximation p NL, b , e.g., (cf. also Figure 6C) However, we will now show that p NL, b has a (local) minimum if (and only if) ∈ 0, 2 b π which approximates the global minimum of p NL and therefore serve as an estimate for the critical connectivity. Starting with and n * specifies the extremum of p NL, b (n). The second derivative of p NL, b (n) at the extremum n * given by Equation (64) satisfies if n * > 0 such that the extremum actually is a minimum. Taken together, for a given setup, i.e., for given dendritic threshold b and coupling strength , the transcendent Equation (64) defines n * which maximizes or minimize p NL, b (n) and if additionally n * > 0 the extremum p NL, b (n * ) is a minimum. Differentiating the right hand side of Equation (64), shows that f (n * ) (as defined in Equation 64) is (1) minimal for n * = 0 and (2) monotonically increasing for n * > 0; according to Equation (64) the minimum n * = 0 corresponds to The left hand side of Equation (64), i.e., √ b / , is monotonically decreasing with from infinity to zero. Thus Equation (64) has a solution for any and p * NL := p * NL, b (n * ) is the (local) minimum of Equation (62) and provides an estimate for the critical connectivity, the (global) minimum of Equation (59).
For better readability we define the function β(·), where n * = n * b as given by Equation (64). We note that β b can also be considered as a function of n * . By combining Equations (62), (64), and (70) we obtain the critical connectivity The function β(·) itself is monotonically decreasing with in the validity range ∈ (0, max ] of the above approximation: within this interval n * > 0 and d dn * f (n * ) > 0 and thus the derivative Consequently β assumes its minimum β min = β n * = 0 = 1 2 (75) for = max = 2 b π and increases monotonically with decreasing against its asymptotic value Thus the critical connectivity is bounded by and converges to the lower bound p 0 for small and to its upper bound 2p 0 for large . In Figure 6 we visualize the determination of the critical connectivity (Equations 64, 70) and Equation (71). The critical connectivity obtained with the approach presented above agrees well with simulation data (cf. Figure 7).

Transition from linear to non-linear propagation
In the previous section we derived analytical estimates for the critical connectivity p * NL in FFNs with non-additive dendritic interactions; p * NL is determined by (1) the setup of the FFN (i.e., the layer size ω and coupling strength ; cf. Figure 7), (2) the parameters of the non-linear modulation function (i.e., the dendritic threshold b and enhancement level κ) and (3) the layout of the external network (i.e., the mean external input I 0 and its variance, which is proportional to ext ). In this section, we discuss the influence of the parameters of the non-linear modulation function and study the transition from a regime where propagation of synchrony is mediated by dendritically enhanced inputs to a regime where the majority of inputs is processed linearly.
In general, with increasing threshold b more and more inputs are needed to reach this threshold and consequently the critical connectivity p * NL increases. If b exceeds μ L , which is the average input to the neurons if a synchronous pulse propagates in linearly coupled FFNs (cf. Equation 45 and Figure 5), propagation mediated by linearly processed spikes is enabled for lower connectivities than propagation mediated by dendritic non-linearities. In this regime the linearly summed inputs (for p = p * L ) are sufficient to maintain propagation of synchrony, but are not sufficient to cross the dendritic threshold. Increasing b even further has no  (77)] and approaches its lower bound for small and its upper bound for large . Like in linearly coupled networks the critical connectivity decays inversely proportional to layer size, p * NL ∝ ω −1 , (cf. also Figure 3), but the scaling with coupling strength is more complicated, p * NL ∝ −1 · 1/β b ; the factor β b ∈ [0.5, 1] [cf. Equation (70) and Figure 6] measures the deviation from the algebraic decay (as found in linearly coupled networks). In this figure the parameters of the external network are fixed to I 0 = 5 mV, ν ext = 3 kHz, ext = 0.5 mV.
influence on the critical connectivity p * NL , here a propagation of synchrony is possible for p ≥ p * L as discussed in section 3.1. We illustrate this transition from non-linear to linear propagation in Figure 8A: We start with large b = μ L such that propagation is enabled for p ≈ p * L and also set κ = μ L . In fact, the linear critical connectivity p * L slightly under-estimates the observed critical connectivity p * NL as it does not account for the saturation of the non-linear modulation function, i.e., for the cutoff σ NL (x) = κ of inputs x ≥ κ. With decreasing b the critical connectivity is substantially reduced and well approximated by Equation (71). Propagation of synchrony is now mainly mediated by dendritically enhanced inputs as described in section 3.2.1. The inset illustrates the impact of decreasing the dendritic threshold b on the iterated map. Initially, for b = μ L = κ, the iterated map for linearly coupled and non-linearly coupled FFNs is similar; with decreasing b the jump like rise in the iterated map is shifted to lower group sizes and consequently the bifurcation point is shifted to lower connectivities.
The non-linear modulation function σ NL (·) (cf. Equation 4) saturates for strong inputs, thus the enhancement level κ defines the maximal (effective) input to a neuron and p f (κ) is an upper bound for the spiking probability of any neuron in response to incoming inputs. This implies that in contrast to linearly coupled FFNs, the average size of a propagating synchronous pulse, γ NL , given by the product of the probability of a neuron receiving sufficiently strong input to reach the dendritic threshold (p γ ; cf. Equation 56), the spiking probability due to that input [p f (κ)] and the layer size ω, is bounded from above by This bound decrease with decreasing κ as illustrated by Figure 8B (inset), where we compare the iterated maps for different values of κ. p f (κ) also influences the critical connectivity p * NL (cf. Equation  71): For small κ the spiking probability p f (κ) is low and thus p * NL is large (it may even exceed p * L ). With increasing κ also p f (κ) increases and consequently the critical connectivity p * NL decreases; for very large κ the spiking probability p f (κ) approaches 1 (cf. Equation 14) and p * NL saturates (cf. Figure 8B). In Figure 8C we show the critical connectivity for an additive enhancement by a constant , i.e., inputs exceeding the dendritic threshold b are increased by the constant value = κ − b . For small κ the critical connectivity p * NL is relatively large and may exceed p * L due to the low saturation level of the non-linear modulation function σ NL (·) (cf. also Figure 8B). As mentioned above, with increasing κ, also p f (κ) increases and the critical connectivity p * NL decreases. However, for large κ and thus large dendritic threshold b propagation of synchrony mediated by linearly processed spikes is possible for lower connectivities than propagation mediated by dendritic non-linearities. Consequently, p * NL converges toward p * L (cf. also Figure 8A).

Influence of external network
In section 3.2.1 we derived an estimate of the critical connectivity p * NL for FFNs with non-additive dendritic interactions. So far we discussed the influence of the setup of the FFN (layer size ω and coupling strength ) as well as the parameters of the non-linear modulation function σ NL (dendritic threshold b and enhancement level κ). In the current section, we focus on the remaining determining factor, the layout of the external network. How does the critical connectivity change with the mean external input I 0 and external coupling strength ext and how well are these changes covered by our analytics?
For the derivation of p * NL we assumed that somatic spikes are elicited exclusively by dendritically enhanced inputs (cf. Equation 54) and thus the critical connectivity depends on the layout of the external network only via p f (κ) (cf. also Equation 71), i.e., on the average spiking probability of a neuron receiving an input larger than the dendritic threshold x ≥ b . For sufficiently small p f (κ), p * NL > 1 and propagation of synchrony is not possible. With increasing p f (κ) the critical connectivity decreases and for In the regime of low spiking rates, changing the mean external input I 0 simply shifts the distribution of membrane potentials P V (V) (which is a Gaussian distribution centered at I 0 ; cf. Equation 8). Thus, with increasing I 0 , p f (κ) increases and the critical connectivity p * NL decreases. In Figure 9A we show the critical connectivity for different ext [which determines the width of P V (V)] vs. the mean external input I 0 . For I 0 = − κ (such that the sum of a dendritically enhanced input and the center of the distribution of membrane lin. see ( Equation 54), the propagation of synchrony is mainly mediated by non-linear enhanced inputs and the critical connectivity can be estimated by Equation (71). For large b the probability that an input from the preceding layer exceeds the dendritic threshold is very low, propagation of synchrony is mainly mediated by linearly processed inputs and the critical connectivity is given by Equation (32) potentials equals the somatic threshold ), p f (κ) simplifies to and thus in the regime of low spiking rates, i.e., ( − I 0 ) /σ 1, p f (κ) ≈ 0.5 independent of the width of the distribution of membrane potentials. Consequently, all curves for different ext coincide at this point. For I 0 > − κ the majority of neurons (>50%) would spike upon receipt of a dendritically enhanced input. Thus p f (κ) increases and therewith the critical connectivity decreases upon decreasing ext . In the limit of → 0, P V (V) converges toward a δ-distribution centered at I 0 and p f becomes a step-function such that the critical connectivity is either constant and minimal for I 0 ≥ − κ or it diverges (no propagation possible) for I 0 < − κ (cf. Figure 9A; magenta curve).
In Figure 9B we illustrate the effect of changing ext on the critical connectivity for constant I 0 . As discussed above for I 0 = − κ, p f (κ) and thus p * NL are rather independent of ext and for I 0 > − κ the critical connectivity increases with ext . For I 0 < − κ an increase of the width of the distribution of membrane potentials shifts the membrane potential of more and more neurons toward the relevant interval [ − κ, ] and thus p f (κ) increases and the critical connectivity p * NL decreases. For the derivation of p * NL we have assumed that the ground state dynamics is essentially not influenced by the spontaneous activity of the FFN itself (i.e., μ = I 0 and σ = ext √ 2τ m ν ext ). As Frontiers in Computational Neuroscience www.frontiersin.org November 2013 | Volume 7 | Article 153 | 15 discussed in section 3.1.3, we can correct the results for such influences. However, since in non-linearly coupled FFNs the impact of (non-linearly enhanced) synchronous activity is much stronger than the impact of spontaneous activity (which is irregular and not amplified by non-additive dendritic interactions), we find that the deviations between the corrected and uncorrected version of p * NL is negligible. Finally, we compare the critical connectivity for networks with and without non-additive dendritic interactions: The factor measures how much the connectivity within the FFN can be reduced by introducing non-additive dendritic interactions. It is independent of the layer size ω and becomes maximal in the limit of small coupling strengths as β ( b / ) → β max = 1 for → 0 (cf. Equation 76). It increases with decreasing b and increasing κ (see discussion in section 3.2.2). In Figure 10 we show the influence of the external network. As discussed above, for small I 0 , propagation of synchrony is not possible (the non-linear enhanced input is insufficient to elicit sufficiently many spikes in the layers of the FFN; white areas in Figure 10). With increasing I 0 , p * NL decreases and c rat increases.

GENERALIZATIONS
In the final section we discuss generalizations of the methods and results we derived. Compared to biological neurons, our models have simplifications which enable the analytical treatment, but might be suspected to be influential on the final result. These simplifications are the homogeneous delay distribution, the simplified initiation and impact of dendritic spikes, the limit of short synaptic currents and the sub-threshold leaky integrate-and-fire  Figure 9). Within the white area, propagation of synchrony is impossible because even for a fully coupled chain the input to the next layer (limited by the saturation of the non-linear modulation function and the layer size) is insufficient . Panel (B) shows the reduction factor c rat (cf. Equation 82), the quotient between the critical connectivity in FFNs without and with non-additive dendritic interactions. The lines enclose the area for which the spontaneous firing is between ν ∈ [0.5, 1.5] Hz obtained from simulations (solid) and low-rate approximation (cf. Equation 9; dashed).
dynamics. Here, we verify that our results generalize to biologically more detailed neurons without these simplifications. In particular, we show that the estimates for the critical connectivity hold. Further, we consider a qualitatively different dendritic interaction function which assumes that the saturation is incomplete, i.e., beyond a region of saturation the impact of larger inputs increases. We show that the tools developed in the article are still applicable and reveal a new phenomenon, the coexistence of linear and non-linear propagation of synchrony.
In the first part (section 3.3.1), we discuss the influence of inhomogeneous delay distribution and finite dendritic integration windows. In the second part (section 3.3.2), we consider the non-linear modulation function with incomplete saturation. Finally, we consider biologically more detailed neuron models (section 3.3.3).

Heterogeneous delays
So far we considered FFNs with homogeneous delay distribution and dendritic modulation functions with integration window of zero length, i.e., only exactly synchronized inputs were possibly non-linearly amplified. Are these assumptions crucial for the obtained results? How does the critical connectivity change in the presence of heterogeneous delay distributions?
To answer this question, we consider synaptic delays τ kl (specifying the synaptic delay between neuron l and k) uniformly drawn from where τ is the mean delay. A direct consequence of heterogeneous delay distribution is that the spikes of the propagating synchronous signal are not simultaneous (i.e., exactly synchronized) anymore. To describe the system accurately one has to consider additionally to the size (g i ) also the temporal jitter (s i ) of the synchronous pulse in the ith layer and investigate the twodimensional iterated map for (g i , s i ) (e.g., Diesmann et al., 1999;Gewaltig et al., 2001;Goedeke and Diesmann, 2008). However, even if the synchronous pulse is blurred out to a pulse packet with finite width, for sufficiently large connectivity stable propagation still can be obtained (see e.g., Gewaltig et al., 2001). For linearly coupled FFNs, with increasing width of the delay distribution, T, the propagating pulse becomes broader and thus the critical connectivity p * L increases (cf. Figures 11A,B; squares). However, the scaling with layer size (cf. Figure 11A) and coupling strength (data not shown) is the same.
Under the assumption that the width of the pulse packet stays bounded, one can derive a lower bound for the critical connectivity. We assume that a pulse in layer i is perfectly synchronized and calculate the effective peak of the depolarization in the (i + 1)th layer. Replacing the coupling strength by the effective depolarization (derived below, cf. Equation 89) in the estimate of the critical connectivity (cf. Equation 32) one gains an estimate of the critical connectivity for systems with heterogeneous delays [Equation (90); shown in Figure 11]. Consider a perfectly synchronized pulse in layer i. Due to inhomogeneities in the delay, the inputs arriving at the (i + 1)th layer are distributed uniformly in an interval of size T (Equation 83). We assume that all inputs arriving at a neuron of layer i + 1 are equidistantly distributed over [− T/2, T/2], i.e., the arrival time of the lth of a total number of k inputs is We consider the subthreshold dynamics only. Each single input depolarizes the neuron by an amount and afterwards the membrane potential V(t) decays exponentially toward its asymptotic value (I 0 ) with the membrane time constant τ m (cf. Equations 1, 2) until the next input arrives after a time interval T k − 1 (cf. Equation 84). Thus the total (effective) depolarization caused by the sum of these k inputs at the end of the considered time We consider the effective depolarization per input, , in the limit of a large number of inputs k (k → ∞), Thus the correction factor C ( T) ≤ 1 defined in Equation (89) relates the coupling strength to the effective coupling strength in the presence of inhomogeneous delays. The critical connectivity is then given by (cf. Equation 32) and this estimate agrees well with direct numerical simulations (cf. Figure 11). For FFNs with dendritic non-linearities and inhomogeneous delays τ kl , one has to consider a finite dendritic integration window t d . Instead of amplifying only simultaneously received spikes (cf. Equation 5), the sum of spikes within the time interval t is considered. We denote the sum of inputs to a neuron within the time interval [t − t, t] by where is the indicator function and t f lm is the mth firing time of neuron l as before. If S t k (t) exceeds the dendritic threshold b for some t = t 0 , neuron k is depolarized additionally (to the depolarization arising from linear spike summation) by add such that the total (effective) depolarization caused by an input x ≥ b equals κ, modeling the effect of a dendritic spike; cf. also section 3.3.3. After such an additional depolarization the dendrite becomes refractory for a time t ref,ds and does not transfer additional spikes within the interval t 0 , t 0 + t ref,ds . For t = 0 we recover the non-linear modulation function σ NL (·) given by Equation (4). Due to the finite dendritic interaction window, a delay distribution with T ≤ t affects the critical connectivity only weakly (cf. Figure 11B). For T > t, some of the inputs received from the preceding layer upon a propagation of synchrony fall out of the dendritic interaction window T and thus the critical connectivity increases. However, the scaling with layer size ω (cf. Figure 11B) and coupling strength (data not shown) is practically identical with the scenario T = 0. Before we discuss propagation of synchrony in biologically more plausible neuron models in section 3.3.3, we consider generalization of the non-linear modulation function in the following section.

Coexistence of linear and non-linear propagation
In this article, we employed a non-linear modulation function σ NL ( ) that is linear for dendritic stimulation smaller than the dendritic threshold, < b , and constant (i.e., saturates) for supra-threshold stimulation, ≥ b (cf. Equation 4).
Biologically, if the linear inputs are transmitted despite the dendritic sodium spike and are not shadowed by, e.g., an NMDA spike, they may lead to a second, later peak depolarization after the one generated by the sodium spike. Since our models replace depolarizations by jumps to the peak depolarization, we have to account for the later peak as soon as it exceeds the earlier one. In this part, we thus assume that if the synchronous input is so large that the depolarization it generates upon linear summation exceeds the depolarization κ generated by the dendritic spike, this former is considered as the effect of the input. In other words, we assume that the dendritic modulation function continues linearly beyond κ, i.e., we define (cf. inset of Figure 12A). The iterated map, mapping the number of active neurons in layer i to the average number of active neurons in layer i + 1 may now have (depending on the system parameters) between one and five fixed points (cf. Figure 12). As before, G 0 = 0 is a trivial fixed point corresponding to the level of absent activity and the only fixed point of the iterated map for small connectivity p. With increasing connectivity p, two additional pairs of fixed points G 1 ≤ G 2 and G 3 ≤ G 4 appear via tangent bifurcations. The first pair of fixed points, G 1 and G 2 , correspond to the propagation of synchrony mediated by non-additive dendritic interactions (as discussed in section 3.1), the second pair, G 3 and G 4 , correspond to propagation of synchrony mediated by linearly processed inputs (as discussed in section 3.2). By further increasing the connectivity p, the fixed points G 2 and G 3 disappear via a tangent bifurcation (cf. Figure 12A). Within the region, where five fixed points exists, both types of propagation of synchrony coexists (illustrated in Figures 12B-D): Synchronized pulses of size g 0 < G 1 typically decay to zero after a small number of layers. Pulse sizes with G 1 < g 0 < G 3 typically initiate propagation of synchrony with an average pulse size around G 2 (where the propagation is mediated by nonadditive dendritic interactions) and synchronous pulses of size g 0 > G 3 typically initiate propagation of synchrony with average pulse sizes around G 4 (linear propagation). For sufficiently large p, i.e., the fixed points G 2 and G 3 disappeared, a synchronized pulse of size g 0 ≥ G 1 will initiate propagation of synchrony with pulse sizes around G 4 ; in this parameter region the nonadditive dendritic interactions essentially increase the basin of attraction of G 4 .
Within the framework of our analytical tractable model, we neglect, e.g., the initiation time of a dendritic spike (in our model non-linear amplifications are instantaneous) or the different shapes of potential deflections caused by linearly and nonlinearly processed inputs. Therefore, propagating synchronous signals mediated either by linear or non-linear dendrites differ only in their size. In biological more detailed models (briefly discussed in section 3.3.3 below) both propagation types will be more distinct, e.g., the propagation frequency (speed) and the quality of synchrony of the propagating pulses are different (see also Jahnke et al., 2012).

Biological more detailed models
The model we mainly consider in this article has the advantage of being analytically tractable. Here we ask whether it over-simplifies the considered systems. More precisely, we study whether the results derived above, in particular the analytical estimates for the critical connectivity, generalize to biologically more detailed models. The main assumption underlying our analysis of linearly coupled networks is a very general one, namely that synchronous single inputs sum up linearly: we assumed that the spiking probability p f (·) of a neuron due to the reception of x synchronous inputs of size equals the spiking probability due to the reception of one single input of size y = x . Therefore, the results will hold also for more complex neuron models, as long as the effect of a synchronous input pulse is approximately the sum of the effects of single inputs. In particular, if the spiking probability due to an input of strength x, p f (x), is sufficiently slowly changing with x, according to Equation (24) the critical connectivity scales like p * L ∝ ( ω) −1 for sufficiently large layer sizes and small coupling strengths. To fully compute the critical connectivity, the actual form of p f (·) has to be known. Our leaky integrate-and-fire neuron with infinitesimally short current pulses approximates the behavior of a wide class of neuron models for which an analytical derivation of p f (·) is impossible. Still even for more detailed models, p f (·) is accessible for measurements in single neuron (computer) experiments.
In Figure 13 we verify our predictions exemplary for two types of neuron models: We employ a model of conductance based leaky integrate-and-fire-type neurons with exponential input conductances (CB-type; see Appendix) and a Hodgkin-Huxleytype neuron model with alpha-function shaped input currents (HH-type; see Appendix). The post-synaptic potential induced by single excitatory inputs is shown in panels (a) and (b) and the scaling of the critical connectivity p * L with ω in panel (c): the scaling of p * L is well described by p * L ∝ ( ω) −1 . The main assumptions underlying our analysis of non-linearly coupled networks are (1) that the maximal spiking probability due to inputs which are subthreshold relative to the dendritic threshold, p f ( b ), is significantly smaller than the spiking probability due to a suprathreshold input, p f (κ), and (2) that the temporal jitter of somatic spikes evoked by suprathreshold inputs is small such that synchronized inputs stay synchronized. Both conditions have been found to be satisfied in biological neurons (e.g., Ariav et al., 2003). Therefore, Equation (71) specifying the critical connectivity p * NL also holds for more detailed neuron models if these models incorporate biologically plausible features of fast dendritic spikes. To obtain a quantitative prediction of p * NL , it is sufficient to estimate (a) the number of inputs needed to elicit a dendritic spike, b / , (b) the layer size ω, and (c) the spiking probability due to the reception of a total input that is sufficiently strong to elicit a dendritic spike.
To investigate the scaling of the critical connectivity p * NL in direct numerical simulations, we account for the effects of dendritic spikes in the CB-type and HH-type: When the total excitatory input within the dendritic integration window exceeds the dendritic threshold level, a current pulse modeling the effect of a dendritic spike is initiated and causes an additional   For each value ω, we evaluated the critical connectivity for four different group sizes ω = 100, 300, 500, 700 and four different coupling strengths = 0.3, 0.6, 0.9, 1.2 nS (CB-type; squares; lower horizontal axis) and = 9, 18, 27, 36 pA (HH-type; crosses; upper horizontal axis), respectively. The lines are fitted functions of the form (λ ω) −1 . The analytical estimate given by Equation (24) holds in the limit of large layer sizes ω and small couplings , therefore we exclude data points from the fitting where a single input yields an EPSP larger than 0.6 mV (CB-type: ≥ 1.4 nS; HH-type: ≥ 46 pA; these points are marked in gray). (D,E) Probability distribution of somatic spike times after stimulation of the neuron by an input which is sufficiently strong to generate a dendritic spike (D: CB-type, E: HH-type). We show exemplary two different configurations for the external inputs, which result in a total somatic spiking probability after dendritic spike generation of p f ≈ 0.97 (solid lines; set 1) and p f ≈ 0.67 (dashed lines; set 2). p f equals the saturation level of the corresponding cumulative distribution function (shown in the insets). depolarization of the soma of the post-synaptic neuron (see Appendix for details; cf. also section 3.3.1). In Figure 13 we compare the results of direct numerical simulations with the estimate given by Equation (71). The post-synaptic potential induced by single excitatory inputs is shown in panels (A) and (B). Panel (D) and (E) shows the spiking probability of a single neuron (in the ground state of the FFN), p f , due to an input exceeding the dendritic threshold level; as examples we present two different setups with p f = {0.67, 0.97}. Panel (F) shows the scaling of p * NL with layer size and coupling strength and the good agreement of the analytical estimate with direct numerical simulations.

SUMMARY AND CONCLUSIONS
Propagation of synchrony in feed-forward sub-structures that are embedded in randomly connected recurrent networks has been a research topic for more than two decades now [see, e.g., review on this topic (Kumar et al., 2010)] and it is hypothesized that such propagation possibly explain the emergence of spatio-temporal spike patterns and information transmission.
In this article, we have analyzed diluted FFNs and investigated their capability to propagate synchrony. In addition to conventional additive (linear) input processing at single neurons, we considered non-additive dendritic interactions modeling the impact of fast dendritic spikes (Ariav et al., 2003;Gasparini et al., 2004;Polsky et al., 2004;Gasparini and Magee, 2006). We emulated the influence of the embedding recurrent network which establishes the irregular ground state in the FFN, by random Poissonian inputs Sompolinsky, 1996, 1998;Brunel, 2000). This approach does not account for back-reactions of activity within the FFN on the embedding network. It is justified as long as the connectivity and connection strength between the neurons of the FFN and the embedding network is low and weak compared to the feed-forward connectivity and connection strength. The back-reaction then influences the activity of the embedding network only weakly and a robust propagation of synchrony can be achieved (Vogels and Abbott, 2005;Kumar et al., 2008;Jahnke et al., 2012). Yet, if the condition is not met, synchronous activity within the FFN may spread out over the embedding network and potentially cause pathological activity ("synfire-explosions") . For specifically structured networks also more complex interactions are possible, such as an enhancement of propagating synchrony (manuscript in preparation).
In the main part of the article, we studied the propagation of synchrony employing leaky integrate-and-fire neurons in the limit of temporally short synaptic inputs and homogeneous synaptic delays. Synchronous pulses consist of exactly synchronized (simultaneous) spikes. This allows to investigate propagation of synchrony by considering the size of a synchronized pulse only, so that the analysis becomes analytically tractable. Nevertheless, in the second part of our article we also consider systems with heterogeneous coupling delays and temporally extended interactions. In agreement with the literature (e.g., Diesmann et al., 1999;Gewaltig et al., 2001;Goedeke and Diesmann, 2008), we observe that pulse packets tend to synchronize along the layers of the FFN so that the results of our simplified description are directly applicable.
We derived scaling laws as well as quantitative estimates for the critical connectivity marking the bifurcation point between the regime where robust propagation of synchrony is possible and where it is not. In particular, based on a suitable series expansion we have shown that for linearly coupled FFNs the critical connectivity decays inversely proportional to layer size and coupling strength. Moreover, the proportionality factor can be estimated from the ground state properties of the single neurons. The estimate agrees with direct numerical simulations within the biologically relevant parameter regime where (a) the spontaneous firing rate of the neurons is low and (b) the distribution of membrane potentials is broad (each neuron receives a huge number of almost random presynaptic inputs). If a synchronous pulse propagates along the layers of a linearly coupled FFN, most of the neurons of each layer participate in the propagation of synchrony, independent of the actual layer size, coupling strength or layout of the external network.
For neurons incorporating non-additive dendritic interactions, the spiking probability as a function of the dendritic stimulation becomes discontinuous. Therefore, the analytical estimation of the critical connectivity in non-linearly coupled FFNs required a different approach than the treatment of linearly coupled FFNs. We have shown that the critical connectivity decays inversely proportional to the layer size (as in linearly coupled FFNs), and we have derived the dependence on the coupling strength which is more complicated. The critical connectivity is completely determined by layer size, spiking probability of the single neuron upon the reception of a non-linearly enhanced presynaptic input and the number of inputs required to reach the dendritic threshold. Our results indicate that in presence of non-linear dendrites, neurons process synchronous inputs similar to threshold units. Such units have been previously used as simplified rate neuron models to study activity propagation in discrete time, e.g., in Nowotny and Huerta (2003); Leibold and Kempter (2006); Cayco-Gajic and Shea-Brown (2013). Because the non-linear modulation function saturates, FFNs with nonadditive dendritic interactions allow for a sparser coding, i.e., only a sub-fraction of each layer (the actual size depends on the non-linear enhancement level) participates in the propagation of synchrony. Whereas stable propagation of synchrony is possible in systems with and without dendritic non-linearities, it occurs in non-linearly coupled FFNs with substantially reduced feed-forward anatomy (reduced connectivity or reduced coupling strength) compared to linearly coupled FFNs.
The analytic derivation of the critical connectivity is based on rather general assumptions: (a) the effect of a synchronous input pulse is approximately the sum of the effects of single inputs and (b) for networks with non-additive dendritic interactions the spiking probability due to non-linearly enhanced input is substantially larger than due to a non-enhanced input. Therefore the predictions and estimates are directly applicable to networks of biologically more detailed neuron models.
In our article we have shown that even highly diluted feedforward structures are suitable to reliably support the directed and constrained propagation of synchronous activity. Such structures occur naturally in sparse, random recurrent networks which are typical for the cortex. These structures might be enhanced by simple synaptic plasticity to enable synchrony propagation. Fast dendritic spikes promote this propagation, as they selectively amplify synchronous inputs and are only weakly influenced by irregular background activity.
Our work indicates that fast dendritic spikes reduce the required synaptic strength and connection density for replay of spike patterns. Moreover, their saturation and the resulting sparse coding might explain the observed variability during replay. Thus, in particular, our understanding of propagation along diluted feed-forward chains may now be combined with knowledge on synaptic plasticity and generation of activity accompanying replay (e.g., sharp wave/ripples) to gain an integrated mechanistic understanding for encoding, replay and memory transfer.
where χ [t − t, t] is the characteristic function of the interval [t − t, t], t f jk is the kth firing time of neuron j and τ denotes the synaptic delay. We denote the peak conductance (coupling strength) for a connection from neuron j to neuron l by g max lj . If g l, t exceeds a threshold g , a dendritic spike is initiated and the dendrite becomes refractory for a time window t DS,ref . The effect of the dendritic spike is incorporated into the model by the current pulse that reaches the soma a time τ DS thereafter. This current pulse is modeled as the sum of three exponential functions, with prefactors A > 0, B > 0, C > 0, decay time constants τ DS,1 , τ DS,2 , τ DS,3 and a dimensionless correction factor c g t , where g t is the summed excitatory input at the initiation time of the dendritic spike as given by Equation (A.16). The factor c g t modulates the pulse strength, ensuring that the peak of the excitatory postsynaptic potential (pEPSP) reaches the experimentally observed region of saturation. At very high excitatory inputs, the conventionally generated depolarization exceeds the level of saturation, c g t is zero and the pEPSP increases (cf. inset of Figure 13A).

Parameters for Figure 13
The single neuron parameters for the numerical simulations are C m l = C m = 400 pF, g L l = g L = 25 nS, V rest l = V rest = −65 mV, l = = −50 mV, t ref l = t ref = 3 ms and V reset l = V reset = −65 mV. The reversal potentials are E Ex = 0 mV and E In = −75 mV and the time constants for the excitatory and inhibitory conductances are τ A,1 = τ G,1 = 2.5 ms and τ A,2 = τ G,2 = 0.5 ms. The parameters of the dendritic spike current are t = 2 ms, g = 8.65 nS, τ DS = 2.7 ms, A = 55 nA, B = 64 nA, C = 9 nA, τ DS,1 = 0.2 ms, τ DS,2 = 0.3 ms, τ DS,3 = 0.7 ms and t ref, DS = 5.2 ms and the dimensionless correction factor is given by c(g) = max 1.5 − g · 0.053nS −1 , 0 . For the first setup (p f ≈ 0.97) we set I 0 l = I 0 = 250 pA, ν ex = 2.4 kHz, ν in = 0.6 kHz, ex = 0.6 nS and in = 6.6 nS; for the second setup (p f ≈ 0.67) we set I 0 l = I 0 = 0 pA, ν ex = 20 kHz, ν in = 5 kHz, ex = 0.6 nS and in = −6.6 nS. causes a peak current of 1 pA]. The time constants τ ex and τ in are the synaptic time constants. As before, we account for dendritic spike generation by considering the sum of excitatory input strengths received by neuron l within the time window t, (A.33) If this sum exceeds the dendritic threshold I , a dendritic spike is initiated and we model its effect is by the current pulse starting after a delay time τ DS after the initiation time of the dendritic spike. The correction factor c ( t ) modulates the pulse strength such that the depolarization saturates for suprathresh-old inputs until the effects of linearly summed input exceed the effects of the dendritic spike (cf. inset of Figure 13B).