# Propagating synchrony in feed-forward networks

^{1}Network Dynamics, Max Planck Institute for Dynamics and Self-Organization (MPIDS), Göttingen, Germany^{2}Bernstein Center for Computational Neuroscience (BCCN), Göttingen, Germany^{3}Fakultät für Physik, Georg-August-Universität Göttingen, Göttingen, Germany^{4}Department for Neuroinformatics, Donders Institute, Radboud University, Nijmegen, Netherlands

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.

## 1. 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 feed-forward 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 (van Vreeswijk and 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 connectivity—the 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.

## 2. Methods and Models

### 2.1. Neuron Model

#### 2.1.1. 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}_{k} for *t* ϵ [*t*_{0}, *t*_{0} + *t*^{ref}].

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 *m*th 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 *m*th excitatory (inhibitory) spike to neuron *k* and ϵ^{exc} > 0 (ϵ^{inh} < 0) denote the corresponding coupling strength.

#### 2.1.2. 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.

### 2.2. 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.

### 2.3. 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 (van Vreeswijk and 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}_{T} − 〈*N*^{X}_{T}〉)^{2}〉, equal ν^{X}*T* = 〈*N*^{X}_{T}〉 = 〈(*N*^{X}_{T} − 〈*N*^{X}_{T}〉)^{2}〉.

Then

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.

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

**Figure 1. Propagation of synchrony in diluted FFNs. (A,B,D,E)** Raster plots of diluted feed-forward networks (*m* = 10, ω = 200, ϵ = 0.25 mV). With increasing connection probability *p* propagation of synchrony can be enabled **(A,B)** in networks with additive (linear) and **(D,E)** in networks with non-additive (non-linear) dendritic interactions (Θ_{b} = 4 mV, κ = 11 mV). **(C,F)** Average number of synchronously active neurons in the second layer, *g*_{2}, vs. the number of synchronously active neurons in the initial layer, *g*_{1}; panel **(C)** linear, panel **(F)** non-additive dendritic interactions (average over *n* = 10,000 trials: solid line, transition probability: shading). Note that non-linear dendrites allow for sparser connectivity, **(E)** vs. **(B)** and for a sparser code, i.e., for smaller numbers of spiking neurons in an activated group, **(F)** vs. **(C)**.

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., 2010, 2011; Litvak 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 *i*th 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 ω.

## 3. 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.

### 3.1. 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.

#### 3.1.1. 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 *i*th 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 ${\sigma}{=}{{\u03f5}}^{{\text{ext}}}\sqrt{{2}{{\tau}}^{{\text{m}}}{{\nu}}^{{\text{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 *g*_{i + 1} = *g*_{i} = *g**_{i}. 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 *g*_{i} + *S*(*g**_{i} − *g*_{i}) in second order in (*g**_{i} − *g*_{i}) near *g*_{i}. The arising quadratic fixed point equation *g**_{i} = *g*_{i} + *S*(*g**_{i} − *g*_{i}) is then analytically solvable in *g**_{i}. 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 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 *n*th 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 *n*th 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 ${\omega}{\to}\frac{{\text{const}}}{{\u03f5}}$ 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}.

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

**Figure 2. Iterated map and bifurcation diagram for the average group size of a propagating synchronous pulse**. **(A)** Iterated map (Equation 19) 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.

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) which allows an direct estimate 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) with *g*_{i + 1} = *g*_{i}. At the bifurcation point (*g** = *g**_{1} = *g**_{2}) where the diagonal is touched, the function *p*_{f}(*gp*ϵ) has to be concave and monotonic increasing with respect to *g*. The definition (Equation 14) of *p*_{f}(*x*) implies that it is monotonic increasing for all *x* ≥ 0. Moreover, it is concave for all *x* ≥ Θ − μ,

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*^{(2)}_{V} (Θ − *x*_{0}); so we determine the expansion point *x**_{0} as a deflection point of *P*_{V}(·), requiring that the second derivative of *P*_{V}(*V*) vanishes for *V* = Θ − *x**_{0},

In the considered regime of low spiking rates, we find ${{x}}_{{0}}^{{\ast}}{=}{\Theta}{-}{\mu}{\pm}\frac{{\sigma}}{\sqrt{{2}}}$, cf. Equation (8). Due to Equation (35)

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

**Figure 3. Critical connectivity p*_{L} in FFNs with linear dendrites decays algebraically with coupling strength ϵ and layer size ω**. The parameters of the external inputs (emulated embedding network) are fixed (

*I*

_{0}= 5 mV, ν

^{ext}= 3 kHz, ϵ

^{ext}= 0.5 mV). Panel

**(A)**shows the critical connectivity

*p**

_{L}vs. the layer size ω for different coupling strengths (ϵ = {0.05 mV (red), 0.1 mV (cyan), 0.125 mV (green), 0.2 mV (blue), and 0.4 mV (black)}) and panel

**(B)**shows

*p**

_{L}vs. the coupling strength ϵ for different layer sizes (ω = {50 (red), 100 (cyan), 150 (green), 200 (blue), and 400 (black)}). 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) (solid lines) with

*x*

_{0}=

*x**

_{0}(cf. Equation 37) 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).

#### 3.1.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 ϵ 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 *n*th 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.

**Figure 4. Robustness of analytical estimates of the critical connectivity**. **(A–C)** We consider the critical connectivity *p**_{L} of infinitely long FFNs, that are approximated by an FFN (*m* = 20, ω = 150, ϵ = 0.2 mV) with periodic boundary conditions in direct numerical simulations (markers), for different layouts of the external network. Panels **(A,B)** show *p**_{L} vs. *I*_{0} for fixed ϵ^{ext} and panel **(C)** shows *p**_{L} vs. ϵ^{ext} for fixed *I*_{0}. The solid (colored) lines indicate the critical connectivity found by numerically determining the bifurcation point of the iterated map (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) the distribution of membrane potentials is sufficiently broad. For further explanations see text (section 3.1.3).

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 computer simulations for the entire considered range of *I*_{0} and ϵ^{ext}, (cf. Figures 4B,C; gray lines).

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 potentials *P*_{V}(*V*) is narrow [for ϵ^{ext} → 0: *P*_{V}(*V*) → δ(*V* − μ); cf. Equation (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).

#### 3.1.4. 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) the average input μ_{L} to the neurons due to a propagation of a synchronous pulse is independent of the layer size ω and coupling strength ϵ. For setups with moderate external inputs (i.e., inputs of the preceding layer influence the neurons' ground state only weakly; see also section 3.1.3) the distribution of membrane potentials *P*_{V}(·) (cf. Equation 8), the firing probability of single neurons *p*_{f}(·) (cf. Equation 14) as well as the expansion point (deflection point of *P*_{V}(·); cf. Equation 37)

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

**Figure 5. Properties of propagating synchronous pulses at the transition from the no-propagation to the propagation regime**. Panels **(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. *I*_{0}; in **(B)** the mean external input *I*_{0} is fixed and μ_{L} is plotted vs. ϵ^{ext}. The markers indicate μ_{L} for FFNs of different sizes [ω and ϵ are given by the legend in **(A)**] obtained by numerical simulations of propagating synchrony. The dashed lines shows the approximation of μ_{L} given by Equation (46) (which is independent of ω and ϵ); the solid lines indicate μ_{L} = *p**_{L}*G*_{2}ϵ; values of *p**_{L} and *G*_{2} are found semi-analytically, by numerically identifying the bifurcation point of the analytically derived iterated map (Equation 13) for the different network setups (both analytical estimates are corrected for the influence of inputs from the preceding layer up to the first order). Panel **(C)** shows the fraction *p*^{frac} of neurons in a layer that participate in the propagation of a synchronous signal vs. α [(Equation 10); main panel] and vs. the spontaneous firing rate ν (inset). Data from different network setups are plotted without distinction as black dots in the main panel and with distinction by different colors and symbols in the inset (see legend); Simulations are repeated for different layouts of the external network (*I*_{0} ∈ {1, 3, …, 11} mV; ϵ^{ext} ∈ {0.1, 0.125, …, 1.0} mV). The solid lines indicate *p*_{f} (μ_{L}) = *f*_{p} (α) as given by Equation (53). The layer size ω as well as the coupling strength ϵ influence *p*^{frac} only weakly. *p*^{frac} depends on the network setup mainly through α or, equivalently, through ν (cf. Equation 9): Measurement values from different network setups largely collapse to the graph of the function *p*_{f} (μ_{L}) = *f*_{p} (α). For further explanations see text (section 3.1.4).

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 ${\alpha}{=}\frac{{\Theta}{-}{\mu}}{{\sigma}}$ (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),

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 ${\alpha}{=}\frac{{\Theta}{-}{\mu}}{{\sigma}}$ 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.

### 3.2. 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.

#### 3.2.1. 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 supra-threshold 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 ${{\sigma}}_{{\delta}}{:}{=}\sqrt{{\delta}{(}{1}{-}{{p}}_{{\gamma}}{p}{{p}}_{{f}}{\left(}{\kappa}{\right)}{)}}$, 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 non-trivial 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)

**Figure 6. Determining the critical connectivity in FFNs with non-additive dendritic interactions. (A)** For a given setup, i.e., for a given dendritic threshold Θ_{b} and coupling strength ${\u03f5}{<}\frac{{2}{{\Theta}}_{{b}}}{{\pi}}$, the corresponding *n** (or equivalently *p*_{γ}; cf. Equation 56) is found by Equation (64). The solid line indicates *n** vs. ϵ (left vertical scale), the dashed line *p*_{γ} vs. ϵ (right vertical scale) and the markers *n**(ϵ) for ϵ = {0.075, 0.3, 2.0} mV (see legend). [Here, the dendritic threshold is Θ_{b} = 4 mV, such that the estimate (Equation 64) is valid within the range ϵ ∈ (0, 2.55] mV; Equation (69)] **(B)** Knowing *n** allows to evaluate ${\beta}{\left(}\frac{{{\Theta}}_{{b}}}{{\u03f5}}{\right)}{\in}{\left[}\frac{{1}}{{2}}{,}{1}{\right]}$ according to Equation (70). Panel **(B)** shows β (cf. Equation 70) vs. ϵ (solid line, lower horizontal axis) and β vs. *n** (dashed line, upper horizontal axis), respectively. **(C)** Finally, the critical connectivity *p**_{NL} is obtained by Equation (71) which depends on ${\beta}{\left(}\frac{{{\Theta}}_{{b}}}{{\u03f5}}{\right)}$. Panel **(C)** shows the connectivity *p*_{NL}[dashed; Equation (59)] and its approximation *p*_{NL, b} [solid; Equation (63)] vs. *n*; for ϵ ∈ (0, ϵ^{max}], *p*_{NL, b} has a local minimum which agrees with the global minimum of *p*_{NL}. The markers indicate the critical connectivity *p**_{NL} obtained by the procedure described in **(A)** and **(B)**. For further explanations see text (section 3.2.1).

However, we will now show that *p*_{NL, b} has a (local) minimum if (and only if) ${\u03f5}{\in}{\left(}{0}{,}{\text{}}\frac{{2}{{\Theta}}_{{b}}}{{\pi}}{\right]}$ which approximates the global minimum of *p*_{NL} and therefore serve as an estimate for the critical connectivity. Starting with ${\frac{{d}{{p}}_{{N}{L}{,}{\text{}}{b}}{(}{n}{)}}{{d}{n}}{|}}_{{n}{=}{{n}}^{{\ast}}}{=}{0}$ yields

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., $\sqrt{{{\Theta}}_{{b}}{/}{\u03f5}}$, 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}}^{{\ast}}{=}{{n}}^{{\ast}}{\left(}\frac{{{\Theta}}_{{b}}}{{\u03f5}}{\right)}$ as given by Equation (64). We note that ${\beta}{\left(}\frac{{{\Theta}}_{{b}}}{{\u03f5}}{\right)}$ 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 $\frac{{d}}{{d}{{n}}^{{\ast}}}{f}{\left(}{{n}}^{{\ast}}{\right)}{>}{0}$ and thus the derivative

Consequently β assumes its minimum

for ${\u03f5}{=}{{\u03f5}}^{{\text{max}}}{=}\frac{{2}{{\Theta}}_{{b}}}{{\pi}}$ 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 2*p*^{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).

**Figure 7. Critical connectivity in FFNs with non-linear dendrites**. The panels show **(A)** the critical connectivity *p**_{NL} vs. the layer size ω for different coupling strengths (ϵ = {0.05 mV (red), 0.1 mV (cyan), 0.125 mV (green), 0.2 mV (blue), and 0.4 mV (black)}) and **(B)** *p**_{NL} vs. the coupling strength ϵ for different layer sizes (ω = {50 (red), 100 (cyan), 150 (green), 200 (blue), and 400 (black)}). The points indicate the minimal connectivity for which a synchronous pulse propagates from the first to the last layer in an FFN with *m* = 20 layers in at least 50% of *n* = 30 trials. The critical connectivity given by Equation (71) (solid lines) is in good agreement with the computer simulations. **(C)** The critical connectivity is confined to the interval *p**_{NL} ∈ [*p*^{0}, 2*p*^{0}] [indicated by the gray area for ω = 150 (green), cf. Equation (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}}_{{N}{L}}^{{\ast}}{\propto}{{\u03f5}}^{{-}{1}}{\text{}}{\xb7}{\text{}}{1}{/}{\beta}{\left(}\frac{{{\Theta}}_{{b}}}{{\u03f5}}{\right)}$; the factor ${\beta}{\left(}\frac{{{\Theta}}_{{b}}}{{\u03f5}}{\right)}{\in}{\left[}{0.5}{,}{\text{}}{1}{\right]}$ [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.

#### 3.2.2. 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 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.

**Figure 8. Transition from linear to non-linear propagation**. The figure shows the critical connectivity *p**_{NL} vs. the parameters of the non-linear modulation function σ_{NL} (cf. Equation 4) for different network setups (color code, see **(C)**). The lines are the theoretical predictions for *p**_{NL} [solid, Equation (71)] and *p**_{L} [dashed, Equation (32)]. The markers indicate the minimal connectivity for which a synchronous pulse propagates from the first to the last layer in an FFN (*I*_{0} = 5 mV, ν^{ext} = 3 kHz, ϵ^{ext} = 0.5 mV) with *m* = 20 layers in at least 50% of *n* = 30 trials. The insets illustrate the effect of changing Θ_{b} and κ on the iterated map, cf. Equation (13), where connectivity is kept constant. **(A)** Critical connectivity vs. dendritic threshold Θ_{b} for constant enhancement level κ = μ_{L} ≈ 13.7 mV (cf. Equation 50). If the dendritic threshold Θ_{b} is sufficiently small such that *p*_{f} (Θ_{b}) « *p*_{f} (κ) (cf. 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). Between these scenarios (for moderate Θ_{b}) there is a “transition regime,” where linear and non-linear propagation mix [similarly in **(C)**]. **(B)** Critical connectivity vs. enhancement level κ for constant threshold Θ_{b} = 4 mV. For small enhancement levels κ the (maximal) spiking probability of a single neuron, *p*_{f} (κ), is small and thus the critical connectivity *p**_{NL} is large. With increasing κ, *p*_{f} (κ) increases and thus *p**_{NL} decreases; for large κ, *p*_{f} (κ) → 1 (a neuron will almost surely spike upon the receipt of a non-linearly enhanced pre-synaptic input) and the critical connectivity saturates. **(C)** Critical connectivity vs. enhancement level κ for an additive enhancement by a constant Δ = κ − Θ_{b} = 4 mV. For further explanations see text (section 3.2.2).

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

#### 3.2.3. 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 *p*_{f} (κ) → 1 it converges to Θ_{b}(ϵωβ[Θ_{b}/ϵ])^{−1}, independent of the external network.

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

**Figure 9. Dependence of the critical connectivity p*_{NL} on the layout of the external network**.

**(A,B)**The lines indicate the theoretical prediction for

*p**

_{NL}given by Equation (71) and agree well with the data from direct numerical simulations (markers; FFN with ω = 150, ϵ = 0.2 mV, Θ

_{b}= 4 mV, κ = 11 mV,

*m*= 20). Panel

**(A)**shows the critical connectivity vs. the mean external input

*I*

_{0}for fixed ϵ

^{ext}and panel

**(B)**shows the critical connectivity vs. ϵ

^{ext}for fixed mean external input

*I*

_{0}. The gray line indicates the minimal critical connectivity obtained for

*p*

_{f}(κ) = 1. With increasing mean (external) input

*I*

_{0}the distribution of membrane potentials

*P*

_{V}(

*V*) is shifted toward the somatic threshold Θ, thus the spiking probability

*p*

_{f}(κ) upon the reception of a non-linear enhanced input increases and the critical connectivity

*p**

_{NL}decreases. For

*I*

_{0}= Θ − κ,

*p*

_{f}(κ) ≈ 0.5 (cf. Equation 80) and

*p**

_{NL}is largely independent of the layout of the external network [blue solid line in

**(B)**; cf. also

**(A)**where all curves coincide]. Further explanations see text (section 3.2.3).

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 ${\sigma}{=}{{\u03f5}}^{{\text{ext}}}\sqrt{{2}{{\tau}}^{{\text{m}}}{{\nu}}^{{\text{ext}}}}$). As 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.

**Figure 10. Critical connectivity and reduction factor**. Panel **(A)** shows the critical connectivity obtained from simulations of an FFN (ω = 150, ϵ = 0.2 mV, *m* = 20) incorporating non-additive dendritic interactions (Θ_{b} = 4 mV, κ = 11 mV; see also 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 83), 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).

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

#### 3.3.1. 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 *i*th layer and investigate the two-dimensional 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.

**Figure 11. Robustness against heterogeneities in the response delays**. **(A)** Critical connectivity vs. layer size for FFNs (*m* = 20, ϵ = 0.25 mV, *I*_{0} = 5 mV, ϵ^{ext} = 0.5 mV) with additive (squares) and non-additive (Θ_{b} = 4 mV, κ = 11 mV, Δ*t* = 2.5 ms; circles) dendritic interactions. Different colors indicate different widths of the delay distribution (cf. Equation 84). The solid lines indicate the critical connectivity *p**_{L} corrected for inhomogeneous delay distribution (cf. Equation 90), the dashed line *p**_{NL} for Δ*T* = 0 ms. **(B)** Critical connectivity vs. width of delay distribution Δ*T*. Different colors indicate different setups of the FFN (red: ω = 275, ϵ = 0.4 mV; green: ω = 125, ϵ = 0.25 mV; blue: ω = 200, ϵ = 0.1 mV). Solid and dashed lines are *p**_{L} and *p**_{NL} as before.

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 *l*th 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 $\frac{{\Delta}{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 interval ${\left(}{\tau}{+}\frac{{\Delta}{T}}{{2}}{\right)}$ is

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 (90) 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 *m*th 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

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.

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

**Figure 12. Coexistence of linear and non-linear propagation**. **(A)** Bifurcation diagram obtained from Equation (13) for an FFN (ω = 150, ϵ = 0.225 mV) with a non-linear modulation function σ′_{NL} with incomplete saturation [cf. Equation (94) and inset]. Panel **(B)** shows the iterated maps (Equation 13) for *p* = 0.5 with the different non-linear modulation functions considered in this article (linear coupling: green,dashed; non-linear coupling σ_{NL}: red, dashed; modified non-linear coupling σ′_{NL}: blue). Panel **(C)** depicts the development of the size of the synchronous pulse along the layers of the FFN (single trials). The blue and yellow regions are the basins of attraction of *G*_{2} and *G*_{4}, respectively, derived from the data in panel **(B)**. Panel **(D)** shows the probability *p*^{conv} of converging to the linear propagation regime (yellow area, blue line) and the non-linear propagation regime (blue area, red line) after *m* = 20 layers (*p*^{conv} is obtained from *n* = 150 runs with different networks and initial conditions).

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 non-additive 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 non-additive 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 non-linearly 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).

#### 3.3.3. 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-Huxley-type 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}.

**Figure 13. Same scaling of propagating regime for networks of biologically more detailed neuron models**. **(A,B)** Time course of the membrane potential of single neurons receiving inputs that are sufficiently strong to elicit a dendritic spike, with (non-linear model) and without (linear model) dendritic spike generation mechanism, for **(A)** a conductance based LIF-type neuron (henceforth: CB-type), and **(B)** a Hodgkin–Huxley-type neuron (HH-type). The insets show the observed peak of the induced postsynaptic potential (pEPSP) vs. the pEPSP expected from linear input summation (equivalent to the dendritic modulation function in the analytically tractable model). **(C)** Critical connectivity *p**_{L} vs. ϵω in linearly coupled networks. 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). **(F)** Critical connectivity *p**_{NL} vs. group size ω (lower horizontal scale) and coupling strength ϵ normalized by threshold Θ_{b} (upper horizontal scale), respectively. The theoretical estimate of *p**_{NL} (cf. Equation 71) is a function of ω, Θ_{b}/ϵ and *p*^{f}, therefore the predictions agree for both models and the data from direct numerical simulations are consistent with the theoretical predictions. [All simulations of FFNs in this figure are obtained for inhomogeneous delay distribution with Δ*T* = 1 ms (cf. Equation 83)].

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

## 4. 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 (van Vreeswijk and 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”) (Mehring et al., 2003). 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 non-additive 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 feed-forward 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.

Indeed, important candidate regions for the generation of propagating synchrony such as the hippocampus and other, neocortical regions exhibiting replay of activity (Nadasdy et al., 1999; Lee and Wilson, 2002; Ji and Wilson, 2007; Xu et al., 2011; Eagleman and Dragoi, 2012) are sparse and show synaptic plasticity (Debanne et al., 1998; Kobayashi and Poo, 2004). Dendritic spikes as prominently found in, e.g., the hippocampus (Ariav et al., 2003; Gasparini et al., 2004; Polsky et al., 2004; Gasparini and Magee, 2006) trigger depolarizations and calcium influx sufficient to change synaptic strengths (Golding et al., 2002; Remy and Spruston, 2007) and the dendrites itself exhibit branch “strength potentiation,” i.e., the strength of a dendritic spike on a dendritic branch exhibits experience- and activity-dependent plasticity (Losonczy et al., 2008; Makara et al., 2009; Müller et al., 2012).

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.

## Conflict of Interest Statement

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

## Acknowledgments

This work was supported by the BMBF (Grant No. 01GQ1005B) [Sven Jahnke, Marc Timme], the DFG (Grant No. TI 629/3-1) [Sven Jahnke], the Swartz Foundation [Raoul-Martin Memmesheimer], and the Max Planck Society [Marc Timme]. Simulation results of networks with biologically more complex neuron models were obtained using the simulation software NEST (Gewaltig and Diesmann, 2007). Sven Jahnke thanks Harold Gutch, Elian Moritz, and Jonna Jahnke for stimulating discussions.

## References

Abeles, M. (1982). *Local Cortical Circuits: An Electrophysiological Study*. Berlin: Springer. doi: 10.1007/978-3-642-81708-3

Ariav, G., Polsky, A., and Schiller, J. (2003). Submillisecond precision of the input-output transformation function mediated by fast sodium dendritic spikes in basal dendrites of CA1 pyramidal neurons. *J. Neurosci*. 23, 7750–7758.

August, D. A., and Levy, W. B. (1999). Temporal sequence compression by an integrate-and-fire model of hippocampal area CA3. *J. Comput. Neurosci*. 6, 71–90. doi: 10.1023/A:1008861001091

Aviel, Y., Mehring, C., Abeles, M., and Horn, D. (2003). On embedding synfire chains in a balanced network. *Neural Comp*. 15, 1321–1340. doi: 10.1162/089976603321780290

Braitenberg, V., and Schüz, A. (1998). *Cortex: Statistics and Geometry of Neuronal Connectivity*. Berlin: Springer. doi: 10.1007/978-3-662-03733-1

Brunel, N. (2000). Dynamics of sparsely connected networks of excitatory and inhibitory spiking neurons. *J. Comp. Neurosci*. 8, 183–208. doi: 10.1023/A:1008925309027

Brunel, N., and Hakim, V. (1999). Fast global oscillations in networks of integrate-and-fire neurons with low firing rates. *Neural Comp*. 11, 1621–1671. doi: 10.1162/089976699300016179

Bullmore, E., and Sporns, O. (2009). Complex brain networks: graph theoretical analysis of structural and functional systems. *Nat. Rev. Neurosci*. 10, 186–198. doi: 10.1038/nrn2618

Burkitt, A. (2006). A review of the integrate-and-fire neuron model: I. Homogeneous synaptic input. *Biol. Cybern*. 95, 1–19. doi: 10.1007/s00422-006-0082-8

Cayco-Gajic, N. A., and Shea-Brown, E. (2013). Neutral stability, rate propagation, and critical branching in feedforward networks. *Neural Comput*. 25, 1768–1806. doi: 10.1162/NECO_a_00461

Dayan, P., and Abbott, L. (2001). *Theoretical Neuroscience: Computational and Mathematical Modeling of Neural Systems*. Cambridge: MIT Press.

Debanne, D., Gähwiler, B. H., and Thompson, S. M. (1998). Long-term synaptic plasticity between pairs of individual CA3 pyramidal cells in rat hippocampal slice cultures. *J. Physiol*. 507, 237–247. doi: 10.1111/j.1469-7793.1998.237bu.x

Diesmann, M., Gewaltig, M. O., and Aertsen, A. (1999). Stable propagation of synchronous spiking in cortical neural networks. *Nature* 402, 529–533. doi: 10.1038/990101

Eagleman, S. L., and Dragoi, V. (2012). Image sequence reactivation in awake V4 networks. *Proc. Natl. Acad. Sci. U.S.A*. 109, 19450–19455. doi: 10.1073/pnas.1212059109

Feinerman, O., and Moses, E. (2006). Transport of information along unidimensional layered networks of dissociated hippocampal neurons and implications for rate coding. *J. Neurosci*. 26, 4526–4534. doi: 10.1523/JNEUROSCI.4692-05.2006

Feinerman, O., Segal, M., and Moses, E. (2005). Signal propagation along unidimensional neuronal networks. *J. Neurophysiol*. 94, 3406–3416. doi: 10.1152/jn.00264.2005

Felleman, D. J., and Van Essen, D. V. C. (1991). Distributed hierarchical processing in the primate cerebral cortex. *Cereb. Cortex* 1, 1–47. doi: 10.1093/cercor/1.1.1

Gasparini, S., and Magee, J. C. (2006). State-dependent dendritic computation in hippocampal CA1 pyramidal neurons. *J. Neurosci*. 26, 2088–2100. doi: 10.1523/JNEUROSCI.4428-05.2006

Gasparini, S., Migliore, M., and Magee, J. C. (2004). On the initiation and propagation of dendritic spikes in CA1 pyramidal neurons. *J. Neurosci*. 24, 11046–11056. doi: 10.1523/JNEUROSCI.2520-04.2004

Gewaltig, M. O., and Diesmann, M. (2007). (NEural Simulation Tool). *Scholarpedia* 2:1430. doi: 10.4249/scholarpedia.1430

Gewaltig, M. O., Diesmann, M., and Aertsen, A. (2001). Propagation of cortical synfire activity: survival probability in single trials and stability in the mean. *Neural Netw*. 14, 657–673. doi: 10.1016/S0893-6080(01)00070-3

Goedeke, S., and Diesmann, M. (2008). The mechanism of synchronization in feed-forward neuronal networks. *New J. Phys*. 10:015007. doi: 10.1088/1367-2630/10/1/015007

Golding, N. L., Staff, N. P., and Spruston, N. (2002). Dendritic spikes as a mechanism for cooperative long-term potentiation. *Nature* 418, 326–331. doi: 10.1038/nature00854

Helias, M., Deger, M., Rotter, S., and Diesmann, M. (2010). Instantaneous non-linear processing by pulse-coupled threshold units. *PLoS Comput. Biol*. 6:e1000929. doi: 10.1371/journal.pcbi.1000929

Holmgren, C., Harkany, T., Svennenfors, B., and Zilberter, Y. (2003). Pyramidal cell communication within local networks in layer 2/3 of rat neocortex. *J. Physiol*. 551, 139–153. doi: 10.1113/jphysiol.2003.044784

Jahnke, S., Timme, M., and Memmesheimer, R.-M. (2012). Guiding synchrony through random networks. *Phys. Rev. X* 2:041016. doi: 10.1103/PhysRevX.2.041016

Ji, D., and Wilson, M. A. (2007). Coordinated memory replay in the visual cortex and hippocampus during sleep. *Nat. Neurosci*. 10, 100–107. doi: 10.1038/nn1825

Kilavik, B. E., Roux, S., Ponce-Alvarez, A., Confais, J., Grün, S., and Riehle, A. (2009). Long-term modifications in motor cortical dynamics induced by intensive practice. *J. Neurosci*. 29, 12653–126563. doi: 10.1523/JNEUROSCI.1554-09.2009

Kobayashi, K., and Poo, M.-M. (2004). Spike train timing-dependent associative modification of Hippocampal CA3 recurrent synapses by mossy fibers. *Neuron* 41, 445–454. doi: 10.1016/S0896-6273(03)00873-0

Kumar, A., Rotter, S., and Aertsen, A. (2008). Conditions for propagating synchronous spiking and asynchronous firing rates in a cortical network model. *J. Neurosci*. 28, 5268–5280. doi: 10.1523/JNEUROSCI.2542-07.2008

Kumar, A., Rotter, S., and Aertsen, A. (2010). Spiking activity propagation in neuronal networks: reconciling differentperspectives on neural coding. *Nat. Rev. Neurosci*. 11, 615–627. doi: 10.1038/nrn2886

Lee, A. K., and Wilson, M. A. (2002). Memory of sequential experience in the Hippocampus during slow wave sleep. *Neuron* 36, 1183–1194. doi: 10.1016/S0896-6273(02)01096-6

Leibold, C., and Kempter, R. (2006). Memory capacity for sequences in a recurrent network with biological constraints. *Neural Compt*. 18, 904–941. doi: 10.1162/neco.2006.18.4.904

Litvak, V., Sompolinsky, H., Segev, I., and Abeles, M. (2013). On the transmission of rate code in long feedforward networks with excitatory-inhibitory balance. *J. Neurosci*. 23, 3006–3015.

Long, M. A., Jin, D. Z., and Fee, M. S. (2010). Support for a synaptic chain model of neuronal sequence generation. *Nature* 468, 394–399. doi: 10.1038/nature09514

Losonczy, A., Makara, J. K., and Magee, J. C. (2008). Compartmentalized dendritic plasticity and input feature storage in neurons. *Nature* 452, 436–441. doi: 10.1038/nature06725

Makara, J. K., Losonczy, A., Wen, Q., and Magee, J. C. (2009). Experience-dependent compartmentalized dendritic plasticity in rat hippocampal CA1 pyramidal neurons. *Nat. Neurosci*. 12, 1485–1487. doi: 10.1038/nn.2428

Mehring, C., Hehl, U., Kubo, M., Diesmann, M., and Aertsen, A. (2003). Activity dynamics and propagation of synchronous spiking in locally connected random networks. *Biol. Cybern*. 88, 395–408. doi: 10.1007/s00422-002-0384-4

Memmesheimer, R.-M. (2010). Quantitative prediction of intermittent high-frequency oscillations in neural networks with supralinear dendritic interactions. *Proc. Natl. Acad. Sci. U.S.A*. 107, 11092–11097. doi: 10.1073/pnas.0909615107

Memmesheimer, R.-M., and Timme, M. (2012). Non-additive coupling enables propagation of synchronous spiking activity in purely random networks. *PLoS Comput. Biol*. 8:e1002384. doi: 10.1371/journal.pcbi.1002384

Müller, C., Beck, H., Coulter, D., and Remy, S. (2012). Inhibitory control of linear and supralinear dendritic excitation in CA1 pyramidal neurons. *Neuron* 75, 851–864. doi: 10.1016/j.neuron.2012.06.025

Nadasdy, Z., Hirase, H., Czurko, A., Csicsvari, J., and Buzsaki, G. (1999). Replay and time compression of recurring spike sequences in the Hippocampus. *J. Neurosci*. 19, 9497–9507.

Nowotny, T., and Huerta, R. (2003). Explaining synchrony in feed-forward networks: are McCulloch-Pitts neurons good enough? *Biol. Cybern*. 89, 237–241. doi: 10.1007/s00422-003-0431-9

Polsky, A., Mel, B. W., and Schiller, J. (2004). Computational subunits in thin dendrites of pyramidal cells. *Nat. Neurosci*. 7, 621–627. doi: 10.1038/nn1253

Putrino, D., Brown, E. N., Mastaglia, F. L., and Ghosh, S. (2010). Differential involvement of excitatory and inhibitory neurons of cat motor cortex in coincident spike activity related to behavioral context. *J. Neurosci*. 30, 8048–8056. doi: 10.1523/JNEUROSCI.0770-10.2010

Remy, S., and Spruston, N. (2007). Dendritic spikes induce single-burst long-term potentiation. *Proc. Natl. Acad. Sci. U.S.A*. 104, 17192–17197. doi: 10.1073/pnas.0707919104

Reyes, A. D. (2003). Synchrony-dependent propagation of firing rate in iteratively constructed networks *in vitro*. *Nat. Neurosci*. 6, 593–599. doi: 10.1038/nn1056

Riehle, A., Grün, S., Diesmann, M., and Aertsen, A. (1997). Spike synchronization and rate modulation differentially involved in motor cortical function. *Science* 278, 1950–1953. doi: 10.1126/science.278.5345.1950

Rosenbaum, R. J., Trousdale, J., and Josic, K. (2010). Pooling and correlated neural activity. *Front. Comput. Neurosci*. 4:9. doi: 10.3389/fncom.2010.00009

Rosenbaum, R., Trousdale, J., and Josic, K. (2011). The effect of pooling on spike train correlations. *Front. Neurosci*. 5:58. doi: 10.3389/fnins.2011.00058

Scannell, J. W., Burns, G. A., Hilgetag, C. C., O'Neil, M. A., and Young, M. P. (1999). The connectional organization of the cortico-thalamic system of the cat. *Cereb. Cortex* 9, 277–299. doi: 10.1093/cercor/9.3.277

Tuckwell, H. (1988). *Introduction to Theoretical Neurobiology*. Cambridge: Cambridge University Press. doi: 10.1017/CBO9780511623202

van Vreeswijk, C., and Sompolinsky, H. (1996). Chaos in neuronal networks with balanced excitatory and inhibitory activity. *Science* 274, 1724–1726. doi: 10.1126/science.274.5293.1724

van Vreeswijk, C., and Sompolinsky, H. (1998). Chaotic balanced state in a model of cortical circuits. *Neural Comp*. 10, 1321–1371. doi: 10.1162/089976698300017214

Vogels, T. P., and Abbott, L. F. (2005). Signal propagation and logic gating in networks of integrate-and-fire neurons. *J. Neurosci*. 25, 10786–10795. doi: 10.1523/JNEUROSCI.3508-05.2005

Xu, S., Jiang, W., Poo, M.-M., and Dan, Y. (2012). Activity recall in a visual cortical ensemble. *Nat. Neurosci*. 15, 449–455. doi: 10.1038/nn.3036

## A. Appendix

### A.1. Proof of Existence of A Global Minimum of *p*_{NL}(*n*)

We will show that *p*_{NL}(*n*) as derived in Equation (59),

has a global minimum for ϵω > Θ_{b}. In Equation (A.2) we defined

For ϵω > Θ_{b}, *p*_{NL} is positive and continuous, and approaches

in the limit of large/small *n*. Further, the derivative of *p*_{NL} can be written as

where we defined the functions

For *n* > 0 and ϵω > Θ_{b},

and in the limit of large *n*,

For ϵω > Θ_{b}, *h*_{1}(*n*) is smaller than two for sufficiently large *n* (cf. Equation A.13) and thus the derivative of *p*_{NL}(*n*) becomes positive (cf. Equation A.6). Consequently *p*_{NL} approaches 1/*p*_{f} (κ) from below for large *n* (cf. also Equation A.5). This proves the existence of a global minimum of *p*_{NL}(*n*), because *p*_{NL} > 1/*p*_{f} (κ) for sufficiently small *n* (cf. Equation A.4).

### A.2. Biological More Detailed Neuron Models

In section 3.3.3 we consider biologically more detailed neuron models. In this appendix we present descriptions of these models including the parameters used for the numerical simulations in Figure 13. These simulations were done using NEST (Gewaltig and Diesmann, 2007), a simulator for spiking neural network models (available at http://www.nest-initiative.org). We implemented new model classes within the NEST framework to handle conductance-based leaky integrate-and-fire neurons with double exponential input conductances as well as non-linear dendritic interactions (source code available from Sven Jahnke).

#### A.2.1. CB-type model

The CB-type model is a leaky integrate-and-fire neuron with conductance based synapses, augmented with a mechanism for the generation of current pulses mimicking the effect of a dendritic spike (see also Memmesheimer, 2010; Jahnke et al., 2012). The subthreshold dynamics of the membrane potential *V*_{l} of neuron *l* obeys the differential equation

Here, *C*^{m}_{l} is the membrane capacity, *g*^{L}_{l} is the resting conductance, *V*^{rest}_{l} is the resting membrane potential, *E*^{Ex} and *E*^{In} are the reversal potentials, and *g*^{A}_{l}(*t*) and *g*^{G}_{l}(*t*) are the conductances of excitatory and inhibitory synaptic populations, respectively. *I*^{DS}_{l}(*t*) models the current pulses caused by dendritic spikes and *I*^{0}_{l} is a constant current gathering slow external and internal currents. The time course of single synaptic conductances contributing to *g*^{A}_{l}(*t*) and *g*^{G}_{l}(*t*) is given by the difference between two exponential functions (e.g., Dayan and Abbott, 2001) with time constants τ^{A, 1} and τ^{A, 2} for the excitatory and τ^{G, 1} and τ^{G, 2} for the inhibitory conductances. Whenever the membrane potential reaches the spike threshold Θ_{l}, the neuron sends a spike to its postsynaptic neurons, is reset to *V*^{reset}_{l} and becomes refractory for a period *t*^{ref}_{l}. Additionally to inputs from the preceding layer each neuron receives excitatory and inhibitory Poissonian input spike trains with rates ν^{ex} and ν^{in}; single inputs have coupling strength ϵ^{ex} and ϵ^{in}, respectively.

To account for dendritic spike generation, we consider the sum *g*_{l, Δt} of excitatory input strengths (characterized by the coupling strengths), arriving at an excitatory neuron *l* within the time window Δ*t* for non-linear dendritic interactions,

where χ_{[t − Δt, t]} is the characteristic function of the interval [*t* − Δ*t, t*], *t*^{f}_{jk} is the *k*th 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.

#### A.2.2 HH-type model

We employ a standard model provided by NEST (“hh_psc_alpha”; Hodgkin–Huxley type neuron with alpha-function shaped postsynaptic currents) and incorporated a dendritic spike current as in the CB-Model. The membrane potential *V*_{l} of neuron *l* obeys the differential equation

For clarity we drop the index *l* in the following; all quantities refer to some neuron *l*. In Equation (A.18),

specify the Na^{+} current, the K^{+} current and leak current. The dynamics of the gating variables *m, n* and *h* are governed by

where the voltage dependencies are given by

In Equations (A.25–A.30) $\tilde{{V}}{(}{t}{)}{:}{=}\frac{{V}{(}{t}{)}}{{1}{\text{mV}}}$ is the value of membrane potential normalized by 1 mV. Spikes are detected by a combined threshold-and-local-maximum search, if there is a local maximum above a certain threshold of the membrane potential, *U*_{Θ} = 0 mV, it is considered a spike (for more details see the NEST manual and the model implementation available at http://www.nest-initiative.org). After a synaptic delay time τ a spike initiates an alpha-function shaped current pulse at the post-synaptic neurons. The total excitatory and inhibitory input to neuron *l* is given by

where ϵ^{ex}_{k} > 0 (ϵ^{in}_{k} < 0) is the strength of the *k*th arriving excitatory (inhibitory) spike at neuron *l, t*^{ex}_{k} (*t*^{in}_{k}) denotes the reception time of that spike and *e* is the Euler constant [the currents *I*^{ex}(*t*) and *I*^{in}(*t*) are normalized such that an input of strength ϵ = 1 pA 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*,

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 suprathreshold inputs until the effects of linearly summed input exceed the effects of the dendritic spike (cf. inset of Figure 13B).

#### A.2.3 Parameters for Figure 13

As before, we consider homogeneous neuronal properties. The single neuron parameters for the numerical simulations are *C*^{m} = 200 pF, *E*^{K} = −77 mV, *E*^{L} = −70 mV, *E*^{Na} = 50 mV, *g*^{K} = 3600 nS, *g*^{L} = 30 nS, *g*^{Na} = 12000 nS, τ^{ex} = 2 ms and τ^{in} = 2 ms. The parameters of the dendritic spike current are Δ*t* = 3.5 ms, *I*^{Θ} = 270 pA, τ^{DS} = 2.7 ms, *A* = 27.5 nA, *B* = 32 nA, *C* = 4.5 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*(ϵ) = max {1.54 − ϵ · 0.002 pA^{−1}, 0}. For the first setup (*p*^{f} ≈ 0.97) we set *I*^{0} = 500 pA, ν^{ex} = 3 kHz, ν^{in} = 3 kHz, ϵ^{ex} = 20 pA and ϵ^{in} = −20 pA; and for the second setup (*p*^{f} ≈ 0.67) we set *I*^{0} = 250 pA, ν^{ex} = 10 kHz, ν^{in} = 10 kHz, ϵ^{ex} = 20 pA and ϵ^{in} = −20 pA.

Keywords: synchrony, networks, synfire chains, spike pattern, mathematical neuroscience, non-additive coupling, non-linear dendrites

Citation: Jahnke S, Memmesheimer R-M and Timme M (2013) Propagating synchrony in feed-forward networks. *Front. Comput. Neurosci*. **7**:153. doi: 10.3389/fncom.2013.00153

Received: 23 June 2013; Accepted: 11 October 2013;

Published online: 15 November 2013.

Edited by:

Tatjana Tchumatchenko, Max Planck Institute for Brain Research, GermanyReviewed by:

Robert Rosenbaum, University of Pittsburgh, USAArvind Kumar, University of Freiburg, Germany

Raul C. Muresan, Romanian Institute of Science and Tehnology, Romania

Copyright © 2013 Jahnke, Memmesheimer and Timme. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Sven Jahnke, Network Dynamics, Max Planck Institute for Dynamics and Self-Organization (MPIDS), Am Faßberg 17, 37077 Göttingen, Germany e-mail: sjahnke@nld.ds.mpg.de