Skip to main content

ORIGINAL RESEARCH article

Front. Comput. Neurosci., 15 November 2013
This article is part of the Research Topic Correlated neuronal activity and its relationship to coding, dynamics and network architecture View all 16 articles

Propagating synchrony in feed-forward networks

  • 1Network Dynamics, Max Planck Institute for Dynamics and Self-Organization (MPIDS), Göttingen, Germany
  • 2Bernstein Center for Computational Neuroscience (BCCN), Göttingen, Germany
  • 3Fakultät für Physik, Georg-August-Universität Göttingen, Göttingen, Germany
  • 4Department 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 Vk(t) and its dynamics satisfy

dVk(t)dt=Vk(t)τkm+Ikconst+Iknet(t)+Ikext(t),(1)

where τmk is the membrane time constant of neuron k, Iconstk := I0kmk a constant input current, Inetk(t) the input current caused by spikes within the network and Iextk(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 Vresetk 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 = t0 the neuron becomes refractory for a time period tref, i.e., Vk(t) = Vresetk for t ϵ [t0, t0 + tref].

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

Iknet(t)=lmϵklδ(ttlmfτkl).(2)

Here ϵkl denotes the coupling strength from neuron l to neuron k, tflm is the mth spike time of neuron l and τkl specifies the synaptic delay. In addition to spikes from the network each neuron receives excitatory and inhibitory random inputs that emulate an embedding network. These external inputs are modeled as random Poisson spike trains with rate νexc and νinh, respectively. The resulting input current is given by

Ikext(t)=mϵexcδ(ttkmext,exc)+mϵinhδ(ttkmext,inh),(3)

where text, exckm (text, inhkm) is the arrival time of the mth 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

σNL(ϵ)={ϵfor ϵ<Θbκotherwise.(4)

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

Iknet(t) =tf[σNL(lMexc(tf)ϵkl)+lMinh(tf)ϵkl ]δ(ttfτ).     (5)

Here, the sum over tf 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 Mexc(tf) and Minh(tf) specify the set of neurons that send an excitatory or inhibitory spike at time tf, 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 (τmk = τm, Θk = Θ and Vresetk = Vreset 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 NXT, X ∈ {exc, inh, int}, in a time interval T is Poisson distributed, the expected number of inputs 〈NXT〉 and the variance 〈 (NXT − 〈NXT〉)2〉, equal νXT = 〈NXT〉 = 〈(NXT − 〈NXT〉)2〉.

Then

μ=I0+τmνexcϵexc+τmνinhϵinh+τmpωνϵ(6)

is the mean of the total input to the neurons in an interval of the size of the membrane time constant, T = τm, and

σ2=τmνexc(ϵexc)2+τmνinh(ϵinh)2+τmpωνϵ2(7)

is its variance. In diffusion approximation, the distribution of membrane potentials PV(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)

PV(V)=1πσ2exp[(Vμσ)2 ](8)

is Gaussian and can be expressed in terms of the input current. In this approximation the average firing rate is

ν=1πτmΘμσexp[(Θμσ)2 ](9)

and depends on μ and σ only via the quotient

α:=Θμσ,(10)

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 g0 ≤ ω 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
www.frontiersin.org

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, g2, vs. the number of synchronously active neurons in the initial layer, g1; 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 gi + 1 in layer i + 1 if gi neurons are synchronously active in layer i. There is always one trivial fixed point, G0, of this iterated map with 0 = G0 = gi + 1 = gi, which corresponds to absent activity. If gi + 1 < gi for all gi > G0, synchronous activity will die out after a small number of layers. If gi + 1gi for some substantial group size, gi > G0, 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 G1 and G2G1 appear. If existing, G1 is always unstable (the diagonal is crossed from below; the slope of the iterated map needs to be larger than one) and G2 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 G2 there is an intersection with the diagonal from above thus the slope is smaller than one and stationary propagation with group sizes around G2 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 tinit we excite all neurons of the first layer to spike synchronously and measure the number of active neurons gi in the ith layer at the expected spiking time texpi = tinit + iτ. If gi 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 (G2) 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 (G1) 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 ith layer, gi neurons spike synchronously, the number of synchronous inputs h a single neuron in layer i + 1 receives follows a binomial distribution h ~ B(gi, p). We denote the spiking probability of a single neuron due to an input of strength x by pf(x). The average or expected spiking probability psp(gi) of a single neuron in layer i + 1 is then given by

psp(gi)=E[pf(hϵ)|gi ]=h=0gi(gih)ph(1p)gihpf(hϵ).(11)

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

E[gi+1|gi ]=ωpsp(gi)(12)
                         =ωh=0gi(gih)ph(1p)gihpf(hϵ).(13)

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 μ = I0 with standard deviation σ=ϵext2τmνext (cf. section 2.3). Using the probability density (Equation 8), we calculate the spiking probability of a single neuron, pf(x), due to an input of strength x;

pf(x)=ΘxΘPV(V)dV(14)
            =12(Erf[Θμσ ]Erf[Θμ+xσ ])(15)

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[gi] → E[gi + 1]), we interpolate E [gi + 1 | gi] for continuous gi and in the second step replace gi by its expectation value E [gi]. The fixed points, E [gi + 1 | E [gi]] = E [gi], 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 gi instead of E[gi].

For sufficiently small connection probabilities p the map (Equation 12) has only one (trivial) fixed point G0 = gi + 1 = gi = 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 G1 (unstable) and G2G1 (stable) appear via a tangent bifurcation.

For FFNs with purely excitatory couplings between the layers, the second fixed point G2 (if it exists) is always stable: The spiking probability pf(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 G2 exists the slope of the iterated map at this intersection point with the diagonal is smaller than 1. This implies that G2 is stable and synchronous pulses of size giG1 typically initiate a propagation of synchrony with an average number of active neurons around G2. 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 pf(·) are known, there is no analytic closed form solution to the fixed point equation gi + 1 = gi = g*i. In other words, we can compute the firing probability pf(x0) for any x0, and therefore also E [gi + 1 | gi] for any gi, but g*i = E [gi + 1 | g*i] is transcendental. We thus derive an approximate solution. We choose some expansion point gi (see section 3.1.2 for details), and approximate the function E [gi + 1 | g*i] by a polynomial gi + S(g*igi) in second order in (g*igi) near gi. The arising quadratic fixed point equation g*i = gi + S(g*igi) 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 gi + S(g*igi) from an appropriate expansion of pf(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 PV(V) is broad, PV(V) changes slowly with V. Then, PV(V) around some V = V0 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 pf(·). Expanding pf(x) into a Taylor series around some x0 and using Equation (12) yields

gi+1=ωE[n=0pf(n)(x0)n!(hϵx0)n|gi ](16)
          =ωn=0pf(n)(x0)n!E[(hϵx0)n|gi ].(17)

Here and in the following we denote the nth derivative of a function f(x) at x = x0 by

f(n)(x0)=ddnxf(x)|x=x0.(18)

Replacing the derivatives of pf(·) by the (one order lower) derivatives of probability density of membrane potentials PV(V) according to Equation (14) yields

gi+1=ωpf(x0)+ωn=1PV(n1)(V0)(1)n1n!E[(hϵx0)n|gi ],(19)

where we defined

V0:=Θx0(20)

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

x0=gipϵ,(21)

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ϵ − x0)n | gi] in Equation (19) simplifies to

E[(hϵx0)n|gi ]=ϵnE[(hE[ h ])n|gi ]=ϵnmn,(22)

where we denote by mn the nth central moment of the Binomial distribution B(gi, 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

gi+1=ωpf(gipϵ).(23)

Using the implicit function theorem one can show that this implies the scaling law

pL=1λϵω(24)

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 PV(V). Therefore this estimate holds if PV(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 x0 up to second order,

gi+1ωpf(x0)+ωPV(V0)·(ϵgipx0)              ωPV(1)(V0)2[(ϵgipx0)2+ϵ2gip(1p) ](25)

The truncated series (Equation 25) is quadratic in gi such that the fixed points g*1/2 = gi + 1 = gi can be obtained analytically,

g1,2=γL±γL2x0(2PV(V0)+x0PV(1)(V0))2pf(x0)p2PV(1)(V0)ϵ2,     (26)

where we defined

γL:=pϵω(2(PV(V0)+x0PV(1)(V0))+(p1)PV(1)(V0)ϵ)22p2PV(1)(V0)ϵ2ω.    (27)

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

γL2=x0(2PV(V0)+x0PV(1)(V0))2pf(x0)p2PV(1)(V0)ϵ2(28)

which yields

pL=121ϵ[λPV(1)(V0)2PV(1)(V0)ω+(ϵPV(1)(V0)2λ)24(PV(1)(V0))2 ](29)

where we defined

λ:=PV(V0)+x0PV(1)(V0)           PV(1)(V0)(x0(2PV(V0)+x0PV(1)(V0))2pf(x0))(30)

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

As before we consider the limit of large layer sizes ω and small coupling strengths ϵ, i.e., we replace ωconstϵ and consider the leading terms of a series expansion of Equation (29). The expansion of the square bracket in Equation (29) yields

λPV(1)(V0)2PV(1)(V0)ϵconst+(ϵPV(1)(V0)2λ)24(PV(1)(V0))2=[λPV(1)(V0)λPV(1)(V0) ]ϵ(1λ·const12)+O(ϵ2),(31)

such that the critical connectivity assumes the functional form given by Equation (24),

pL1λϵω.(32)

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 pf(x). It is thus necessary to choose an appropriate expansion point that results in fast convergence. In particular for the choice x0 = 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
www.frontiersin.org

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 x0 = 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 x0 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 gi + 1 = gi 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 gi + 1 = gi. At the bifurcation point (g* = g*1 = g*2) where the diagonal is touched, the function pf(gpϵ) has to be concave and monotonic increasing with respect to g. The definition (Equation 14) of pf(x) implies that it is monotonic increasing for all x ≥ 0. Moreover, it is concave for all x ≥ Θ − μ,

pf(1)(x)=PV(Θx)0for x0(33)
pf(2)(x)=PV(1)(Θx)0for xΘμ,(34)

such that the bifurcation point satisfies

x0Θμ.(35)

The condition Equation (33) holds because PV(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 x0 = 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 (Θ − x0); so we determine the expansion point x*0 as a deflection point of PV(·), requiring that the second derivative of PV(V) vanishes for V = Θ − x*0,

pf(3)(x0)=d2PV(V)dV2|V=Θx0=!0.(36)

In the considered regime of low spiking rates, we find x0=Θμ±σ2, cf. Equation (8). Due to Equation (35)

x0=Θμ+σ2.(37)

For x0 = 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
www.frontiersin.org

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 (I0 = 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 x0 = 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

μ1=I0(38)
σ1=ϵext2τmνext,(39)

as assumed in the previous section. All following layers receive external inputs and spikes from their preceding layer(s). The mean μn and standard deviation σn of the input to neurons of the nth layer (with n ≥ 2) reads (cf. Equations 6 and 7)

μn=I0+τmpωνn1ϵ(40)
σn=2νextτm(ϵext)2+pωνn1τmϵ2.(41)

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

νn1=1πτmΘμn1σn1exp[(Θμn1σn1)2 ].(42)

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 ΔXn := XnXn − 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 pf(·)] 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 I0 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 I0 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 I0 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 I0 (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 I0-dependent threshold, the low-rate approximation becomes inapplicable.

FIGURE 4
www.frontiersin.org

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. I0 for fixed ϵext and panel (C) shows p*L vs. ϵext for fixed I0. 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 PV(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 I0 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

μn=I0+τm/λn1νn1(43)
σn=2νextτm(ϵext)2+ϵνn1τm/λn1(44)

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 PV(V) is narrow [for ϵext → 0: PV(V) → δ(V − μ); cf. Equation (8)], the spiking probability of one neuron, pf(·), increases steeply in a small interval [for ϵext → 0: pf(x) → Θ (x − μ); cf. Equation (8)] and thus the approximation of pf(·) 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

μL=γLpLϵ=PV(Θx0)+PV(1)(Θx0)x0λPV(1)(Θx0)(45)

and after inserting λ* as given by Equation (30),

μL=x0(2PV(Θx0)+x0PV(1)(Θx0))2pf(x0)PV(1)(Θx0).     (46)

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 PV(·) (cf. Equation 8), the firing probability of single neurons pf(·) (cf. Equation 14) as well as the expansion point (deflection point of PV(·); cf. Equation 37)

x0=ΘI0+ϵextτmνext(47)

are fully determined by the external inputs (I0, ν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 I0 and ϵext are kept fixed. With increasing mean of the external input (I0) the distribution of membrane potentials PV(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
www.frontiersin.org

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. I0; in (B) the mean external input I0 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*LG2ϵ; values of p*L and G2 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 pfrac 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 (I0 ∈ {1, 3, …, 11} mV; ϵext ∈ {0.1, 0.125, …, 1.0} mV). The solid lines indicate pfL) = fp (α) as given by Equation (53). The layer size ω as well as the coupling strength ϵ influence pfrac only weakly. pfrac 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 pfL) = fp (α). 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 pfrac of neurons of one layer that participate in a propagating synchronous pulse,

pfrac=γLω=pf(μL).(48)

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, pfrac depends on the setup of the external inputs only via the quotient α=Θμσ (cf. Equation 10), or, equivalently, on the spontaneous firing rate ν of the neurons (cf. Equation 9). This can be shown by combining Equations (8, 37) and (Equation 46),

μL=σ(eπ2)1/4[ (2+2α)(3+2α)2eπErf(12)Erf(α) ]1/2(49)
       =:σfμ(α)(50)

such that

pf(μL)=12[Erf(Θμσ)+Erf(μLΘ+μσ) ](51)
        =12[Erf(α)+Erf(μLσα) ](52)
            =12[Erf(α)+Erf(fμ(α)α) ]=:fp(α).(53)

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 I0 and ϵext), we detect whether propagation of a synchronous pulse is possible and if so, we numerically determine the average fraction of participating neurons as well as the spontaneous firing frequency. We find that indeed the size of the synchronous pulse is determined essentially by the quotient α=Θμσ and Equation (53) is a reasonable estimate of the average fraction of neurons spiking in each layer. With increasing α the fraction of participating neurons increases, it thus decreases with spontaneous firing rate ν see Figure 5C. For FFNs with low spontaneous spiking frequency almost all neurons of a layer participate in the propagation of a synchronous pulse.

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, pfNL(x)), is discontinuous and an approach based on expansion of pf(·) 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

pf(σNL(x))={pf(κ)ifxΘb0otherwise,(54)

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γpf (κ) p), where pγpf (κ) 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 (pf (κ)) 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γ,

pγ=k=Θb/ϵ ω(ωk)(pγpf(κ)p)k(1pγpf(κ)p)ωk.(55)

To solve Equation (55) we approximate the binomial distribution by a Gaussian distribution with mean δ := ω pγppf (κ) and standard deviation σδ:=δ(1pγppf(κ)), which yields

pγ=12[1+Erf(n2) ],(56)

where we defined

n:=δΘb/ϵσδ(57)
        =ωpγppf(κ)Θb/ϵωpγppf(κ)(1pγppf(κ))(58)

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

pNL=n2ϵ+2Θb+nn2ϵ2+4Θb(ϵΘbω)pf(κ)ϵ(n2+ω)(1+Erf(n2)),(59)

which is the connectivity pNL 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

ϵω>Θb(60)

(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), pNL(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 pNL above the minimal connectivity p*NL has two preimages n1 and n2 corresponding to the both non-trivial fixed points G1 and G2 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 n1 and n2 collapse to n* = n1 = n2 and correspond to the fixed point G = G1 = G2 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 pNL 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 pNL (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

pNLpNL,a:=2Θbpf(κ)ϵω1+nϵΘb1ω1+Erf(n2).(61)

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

pNL,b:=2Θbpf(κ)ϵω1+nϵΘb1+Erf(n2).(62)

Whereas pNL has always a global minimum for ϵω > Θb, this does not hold for the approximation pNL, b, e.g., (cf. also Figure 6C)

limn(pNL,b)=.(63)
FIGURE 6
www.frontiersin.org

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 ϵ<2Θbπ, 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 β(Θbϵ)[12,1 ] 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 β(Θbϵ). Panel (C) shows the connectivity pNL[dashed; Equation (59)] and its approximation pNL, b [solid; Equation (63)] vs. n; for ϵ ∈ (0, ϵmax], pNL, b has a local minimum which agrees with the global minimum of pNL. 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 pNL, b has a (local) minimum if (and only if) ϵ(0, 2Θbπ ] which approximates the global minimum of pNL and therefore serve as an estimate for the critical connectivity. Starting with dpNL, b(n)dn|n=n=0 yields

Θbϵ=π2exp(n22)(1+Erf(n2))n=:f(n),(64)

and n* specifies the extremum of pNL, b(n). The second derivative of pNL, b(n) at the extremum n* given by Equation (64) satisfies

dpNL,b2dn2|n=n=2nΘbϵpf(κ)ω(1+Erf[n2 ])>0(65)

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 pNL, b(n) and if additionally n* > 0 the extremum pNL, b(n*) is a minimum.

Differentiating the right hand side of Equation (64),

  df(n)dn=n·en22π2(1+Erf[n2 ])(66)
d2f(n)dn2=n+(1+n2)en22π2(1+Erf[n2 ]), (67)

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

ϵmax:=Θb[f(0) ]2=2Θbπ0.64Θb.(68)

The left hand side of Equation (64), i.e., Θb/ϵ, is monotonically decreasing with ϵ from infinity to zero. Thus Equation (64) has a solution for any

ϵ(0,ϵmax ]=(0,2Θbπ ](69)

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 β(·),

β(Θbϵ):=12(1+Erf[n2 ])nen222π,(70)

where n=n(Θbϵ) as given by Equation (64). We note that β(Θbϵ) can also be considered as a function of n*. By combining Equations (62), (64), and (70) we obtain the critical connectivity

pNL=Θbpf(κ)ϵω·1β(Θbϵ).(71)

The function β(·) itself is monotonically decreasing with ϵ in the validity range ϵ ∈ (0, ϵmax] of the above approximation: within this interval n* > 0 and ddnf(n)>0 and thus the derivative

dβdϵ=dβdn·dndΘb/ϵ·dΘb/ϵdϵ(72)
        =en22n22π·(df(n)dn)1·Θb4ϵ3(73)
<0.(74)

Consequently β assumes its minimum

βmin=β(n=0)=12(75)

for ϵ=ϵmax=2Θbπ and increases monotonically with decreasing ϵ against its asymptotic value

βmax=limn[12(1+Erf[n2 ])nen222π ]=1.(76)

Thus the critical connectivity is bounded by

p0:=Θbpf(κ)ϵωpNL2·Θbpf(κ)ϵω=2·p0(77)

and converges to the lower bound p0 for small ϵ and to its upper bound 2p0 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
www.frontiersin.org

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 ∈ [p0, 2p0] [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, pNLϵ1 · 1/β(Θbϵ); the factor β(Θbϵ)[0.5, 1 ] [cf. Equation (70) and Figure 6] measures the deviation from the algebraic decay (as found in linearly coupled networks). In this figure the parameters of the external network are fixed to I0 = 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 I0 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 pp*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 pp*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
www.frontiersin.org

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 (I0 = 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 pfb) « pf (κ) (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, pf (κ), is small and thus the critical connectivity p*NL is large. With increasing κ, pf (κ) increases and thus p*NL decreases; for large κ, pf (κ) → 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 pf (κ) 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 [pf (κ)] and the layer size ω, is bounded from above by

γNL=pγpf(κ)ωωpf(κ)=:γmax.(78)

This bound decrease with decreasing κ as illustrated by Figure 8B (inset), where we compare the iterated maps for different values of κ. pf (κ) also influences the critical connectivity p*NL (cf. Equation 71): For small κ the spiking probability pf (κ) is low and thus p*NL is large (it may even exceed p*L). With increasing κ also pf (κ) increases and consequently the critical connectivity p*NL decreases; for very large κ the spiking probability pf (κ) 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 pf (κ) 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 I0 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 pf (κ) (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 pf (κ), p*NL > 1 and propagation of synchrony is not possible. With increasing pf (κ) the critical connectivity decreases and for pf (κ) → 1 it converges to Θb(ϵωβ[Θb/ϵ])−1, independent of the external network.

In the regime of low spiking rates, changing the mean external input I0 simply shifts the distribution of membrane potentials PV(V) (which is a Gaussian distribution centered at I0; cf. Equation 8). Thus, with increasing I0, pf (κ) increases and the critical connectivity p*NL decreases.

In Figure 9A we show the critical connectivity for different ϵext [which determines the width of PV(V)] vs. the mean external input I0. For I0 = Θ − κ (such that the sum of a dendritically enhanced input and the center of the distribution of membrane potentials equals the somatic threshold Θ), pf (κ) simplifies to

pf(κ)=12(Erf[ΘI0σ ]+Erf[κΘ+I0σ ])(79)
=12Erf(ΘI0σ)(80)

and thus in the regime of low spiking rates, i.e., (Θ − I0)/σ » 1, pf (κ) ≈ 0.5 independent of the width of the distribution of membrane potentials. Consequently, all curves for different ϵext coincide at this point. For I0 > Θ − κ the majority of neurons (>50%) would spike upon receipt of a dendritically enhanced input. Thus pf (κ) increases and therewith the critical connectivity decreases upon decreasing ϵext. In the limit of ϵ → 0, PV(V) converges toward a δ-distribution centered at I0 and pf becomes a step-function

pf(κ)={0κ<ΘI01κΘI0(81)

such that the critical connectivity is either constant and minimal for I0 ≥ Θ − κ or it diverges (no propagation possible) for I0 < Θ − κ (cf. Figure 9A; magenta curve).

FIGURE 9
www.frontiersin.org

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 I0 for fixed ϵext and panel (B) shows the critical connectivity vs. ϵext for fixed mean external input I0. The gray line indicates the minimal critical connectivity obtained for pf (κ) = 1. With increasing mean (external) input I0 the distribution of membrane potentials PV(V) is shifted toward the somatic threshold Θ, thus the spiking probability pf (κ) upon the reception of a non-linear enhanced input increases and the critical connectivity p*NL decreases. For I0 = Θ − κ, pf (κ) ≈ 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 I0. As discussed above for I0 = Θ − κ, pf (κ) and thus p*NL are rather independent of ϵext and for I0 > Θ − κ the critical connectivity increases with ϵext. For I0 < Θ − κ 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 pf (κ) 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., μ = I0 and σ=ϵext2τmν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

crat:=pLpNL=pf(κ)λΘbβ(Θbϵ)(82)

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 I0, 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 I0, p*NL decreases and crat increases.

FIGURE 10
www.frontiersin.org

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 crat (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

τkl[τΔT2,τ+ΔT2 ],(83)

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 (gi) also the temporal jitter (si) of the synchronous pulse in the ith layer and investigate the two-dimensional iterated map for (gi, si) (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
www.frontiersin.org

Figure 11. Robustness against heterogeneities in the response delays. (A) Critical connectivity vs. layer size for FFNs (m = 20, ϵ = 0.25 mV, I0 = 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 lth of a total number of k inputs is

tlarr=τΔT2+ΔTk1·(l1).(84)

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 (I0) with the membrane time constant τm (cf. Equations 1, 2) until the next input arrives after a time interval ΔTk1 (cf. Equation 84). Thus the total (effective) depolarization caused by the sum of these k inputs at the end of the considered time interval (τ+ΔT2) is

Δϵk=l=1kϵexp(1τmΔTk1(l1))(85)
=ϵexp(ΔTτmkk1)1exp(ΔTτm1k1)1.(86)

We consider the effective depolarization per input, ϵ′, in the limit of a large number of inputs k (k → ∞),

ϵ=limk(Δϵkk)(87)
      =τmΔT(1exp[ΔTτm ])ϵ(88)
  =:C(ΔT)ϵ.(89)

Thus the correction factor CT) ≤ 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)

pL=1C(ΔT)·1λϵω(90)

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 Δtd. 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

SkΔt(t)=lmϵχ[tΔt,t ](tlmf+τkl),(91)

where

χA(x)={1if xA0if xA(92)

is the indicator function and tflm is the mth firing time of neuron l as before. If SΔtk(t) exceeds the dendritic threshold Θb for some t = t0, neuron k is depolarized additionally (to the depolarization arising from linear spike summation) by

ϵκadd(t0)=κSkΔt(t0)(93)

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 tref,ds and does not transfer additional spikes within the interval [t0, t0 + tref,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

σNL(ϵ)={ϵforϵΘbκforΘbϵκϵforϵκ(94)

(cf. inset of Figure 12A).

FIGURE 12
www.frontiersin.org

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 G2 and G4, respectively, derived from the data in panel (B). Panel (D) shows the probability pconv 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 (pconv 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, G0 = 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 G1G2 and G3G4 appear via tangent bifurcations. The first pair of fixed points, G1 and G2, correspond to the propagation of synchrony mediated by non-additive dendritic interactions (as discussed in section 3.1), the second pair, G3 and G4, 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 G2 and G3 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 g0 < G1 typically decay to zero after a small number of layers. Pulse sizes with G1 < g0 < G3 typically initiate propagation of synchrony with an average pulse size around G2 (where the propagation is mediated by non-additive dendritic interactions) and synchronous pulses of size g0 > G3 typically initiate propagation of synchrony with average pulse sizes around G4 (linear propagation). For sufficiently large p, i.e., the fixed points G2 and G3 disappeared, a synchronized pulse of size g0G1 will initiate propagation of synchrony with pulse sizes around G4; in this parameter region the non-additive dendritic interactions essentially increase the basin of attraction of G4.

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 pf(·) 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, pf(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 pf(·) 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 pf(·) is impossible. Still even for more detailed models, pf(·) 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
www.frontiersin.org

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 pf ≈ 0.97 (solid lines; set 1) and pf ≈ 0.67 (dashed lines; set 2). pf 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 pf, 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, pfb), is significantly smaller than the spiking probability due to a suprathreshold input, pf (κ), 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), pf, due to an input exceeding the dendritic threshold level; as examples we present two different setups with pf = {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

CrossRef Full Text

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

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

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

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

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

CrossRef Full Text

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

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

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

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

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

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

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

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

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

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

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

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

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

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

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

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

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

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

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

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

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

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

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

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

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

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

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

CrossRef Full Text

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

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

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

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

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

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

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

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

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

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

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

CrossRef Full Text

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

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

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

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

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

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

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

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

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

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

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

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

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

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

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.

Pubmed Abstract | Pubmed Full Text

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

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

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

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

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

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

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

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

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

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

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

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

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

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

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.

Pubmed Abstract | Pubmed Full Text

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

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

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

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

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

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

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

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

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

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

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

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

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

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

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

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

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

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

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

CrossRef Full Text

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

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

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

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

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

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

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

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

A. Appendix

A.1. Proof of Existence of A Global Minimum of pNL(n)

We will show that pNL(n) as derived in Equation (59),

pNL(n)=n2ϵ+2Θb+nn2ϵ2+4Θb(ϵΘbω)pf(κ)ϵ(n2+ω)(1+Erf(n2))(A.1)
            =1pf(κ)2Θb+n2ϵ(1+1+αn2)(1+Erf(n2))(n2ϵ+ωϵ),(A.2)

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

α:=4Θbϵ(1Θbϵω).(A.3)

For ϵω > Θb, pNL is positive and continuous, and approaches

limn(pNL(n))=,(A.4)
   limn(pNL(n))=1pf(κ),(A.5)

in the limit of large/small n. Further, the derivative of pNL can be written as

ddnpNL(n)=(2h1(n))h2(n),(A.6)

where we defined the functions

h1(n)=1ϵω(2πen22(n2+ω)(2Θbα+n2+n+nϵ)1+Erf(n2)+4Θbnα+n2+n+αϵ(n2+ω)α+n(α+n2+n),)(A.7)
h2(n)=ω(α+n2+n)pf(κ)(Erf(n2)+1)(n2+ω)2.(A.8)

For n > 0 and ϵω > Θb,

      α>0,(A.9)
h1(n)>0,(A.10)
h2(n)>0,(A.11)

and in the limit of large n,

limnh1(n)=1ϵω(0+2Θb+αϵ2)(A.12)
           =22ΘbωϵΘb2ω2ϵ2(A.13)
limnh2(n)=0.(A.14)

For ϵω > Θb, h1(n) is smaller than two for sufficiently large n (cf. Equation A.13) and thus the derivative of pNL(n) becomes positive (cf. Equation A.6). Consequently pNL approaches 1/pf (κ) from below for large n (cf. also Equation A.5). This proves the existence of a global minimum of pNL(n), because pNL > 1/pf (κ) 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 Vl of neuron l obeys the differential equation

ClmdVl(t)dt=glL(VlrestVl(t))+glA(t)(EExVl(t))                      +glG(t)(EInVl(t))+IlDS(t)+Il0.(A.15)

Here, Cml is the membrane capacity, gLl is the resting conductance, Vrestl is the resting membrane potential, EEx and EIn are the reversal potentials, and gAl(t) and gGl(t) are the conductances of excitatory and inhibitory synaptic populations, respectively. IDSl(t) models the current pulses caused by dendritic spikes and I0l is a constant current gathering slow external and internal currents. The time course of single synaptic conductances contributing to gAl(t) and gGl(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 Vresetl and becomes refractory for a period trefl. 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 gl, Δ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,

gl,Δt(t)=jkϵljχ[tΔt,t ](tjkf+τ), (A.16)

where χ[t − Δt, t] is the characteristic function of the interval [t − Δt, t], tfjk is the kth firing time of neuron j and τ denotes the synaptic delay. We denote the peak conductance (coupling strength) for a connection from neuron j to neuron l by gmaxlj. If gl, Δt exceeds a threshold gΘ, a dendritic spike is initiated and the dendrite becomes refractory for a time window tDS,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,

IlDS(t)=c(gΔt)[AetτDS,1+BetτDS,2CetτDS,3 ],(A.17)

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 Cml = Cm = 400 pF, gLl = gL = 25 nS, Vrestl = Vrest = −65 mV, Θl = Θ = −50 mV, trefl = tref = 3 ms and Vresetl = Vreset = −65 mV. The reversal potentials are EEx = 0 mV and EIn = −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 tref, 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 (pf ≈ 0.97) we set I0l = I0 = 250 pA, νex = 2.4 kHz, νin = 0.6 kHz, ϵex = 0.6 nS and ϵin = 6.6 nS; for the second setup (pf ≈ 0.67) we set I0l = I0 = 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 Vl of neuron l obeys the differential equation

ClmdVl(t)dt=IlNa(t)+IlK(t)+IlL(t)+Il0+Ilex(t)+Ilin(t)+IlDS(t).(A.18)

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

     INa(t)=gNam(t)3h(t)[ENaV(t) ](A.19)
  IK(t)=gKn(t)4[EKV(t) ](A.20)
IL(t)=gL[ELV(t) ](A.21)

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

dm(t)dt=αm(t)[1m(t) ]βm(t)m(t)(A.22)
dh(t)dt=αh(t)[1h(t) ]βh(t)h(t)(A.23)
dn(t)dt=αn(t)[1n(t) ]βn(t)n(t),(A.24)

where the voltage dependencies are given by

  αn(t)=0.01[V˜(t)+55 ]1exp[V˜(t)+5510 ](A.25)
      βn(t)=0.125·exp[V˜(t)+6580 ](A.26)
 αm(t)=0.1[V˜(t)+40 ]1exp[V(t)+4010 ](A.27)
     βm(t)=4·exp[V˜(t)+6518 ](A.28)
      αh(t)=0.07·exp[V˜(t)+6520 ](A.29)
      βh(t)=(1+exp[V˜(t)+3510 ])1.(A.30)

In Equations (A.25–A.30) V˜(t):=V(t)1mV 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

Iex(t)=kϵkexeτexexp[tτex ]Θ[ttkex ](A.31)
Iin(t)=kϵkineτinexp[tτin ]Θ[ttkin ],(A.32)

where ϵexk > 0 (ϵink < 0) is the strength of the kth arriving excitatory (inhibitory) spike at neuron l, texk (tink) denotes the reception time of that spike and e is the Euler constant [the currents Iex(t) and Iin(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,

ϵΔt(t)=kϵkexχ[tΔt,t ](tkf+τ).(A.33)

If this sum exceeds the dendritic threshold IΘ, a dendritic spike is initiated and we model its effect is by the current pulse

IDS(t)=c(ϵΔt)[AetτDS,1+BetτDS,2CetτDS,3 ],(A.34)

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 Cm = 200 pF, EK = −77 mV, EL = −70 mV, ENa = 50 mV, gK = 3600 nS, gL = 30 nS, gNa = 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 tref,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 (pf ≈ 0.97) we set I0 = 500 pA, νex = 3 kHz, νin = 3 kHz, ϵex = 20 pA and ϵin = −20 pA; and for the second setup (pf ≈ 0.67) we set I0 = 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, Germany

Reviewed by:

Robert Rosenbaum, University of Pittsburgh, USA
Arvind 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

Disclaimer: All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.