Noise focusing in neuronal tissues : Symmetry breaking and localization in excitable networks with quenched disorder

Javier G. Orlandi1,2,* and Jaume Casademunt1,3 1Departament de Física de la Matèria Condensada, Universitat de Barcelona, E-08028 Barcelona, Spain 2Complexity Science Group, Department of Physics and Astronomy, University of Calgary, Calgary, Canada T2N 1N4 3Universitat de Barcelona Institute of Complex Systems (UBICS), Universitat de Barcelona, E-08028 Barcelona, Spain (Received 5 April 2016; revised manuscript received 16 November 2016; published 3 May 2017)


I. INTRODUCTION
The study of emergent, self-organized phenomena in neuronal networks is perceived as an important key to understand brain behavior and function [1,2].From a physical perspective, in vitro neuronal cultures have proven to be an ideal model system to search for insights that may unveil general principles of collective neuronal dynamics [3][4][5][6].One aspect of particular interest and direct relevance to actual biological tissues is the study of spontaneous activity of neuronal cultures.In particular, a puzzling phenomenon that is robustly observed at early stages of culture development is the emergence of coherence in the form of nearly periodic bursting of the whole network, out of random spontaneous firing of neurons [7][8][9][10][11].The detailed scrutiny of this phenomenon in experiments in vitro and in silico simulations has led to the concept of noise focusing [12].The idea is that the network structure amplifies and propagates the noise in such a way that spontaneous activity strongly concentrates at specific sites, which nucleate bursts that then propagate through the whole system.Remarkably, even though the occurrence of these bursts is nearly periodic, the nucleation sites alternate randomly.
It is our purpose to transcend the qualitative picture emerging from Ref. [12] by developing a theoretical framework that enables a deeper quantitative understanding of the phenomenon of noise focusing and its contextualization from the perspectives of excitable systems [13], constructive effects of noise [14], and transport in complex networks [15].The proposed model reproduces the experimental observations in cultures and at the same time provides insights on the origin of the inherent symmetry breaking emerging from the network quenched disorder.
The phenomenon of noise focusing, as introduced in Ref. [12], is characteristic of directed networks of excitable units with metric connectivity correlations, where each node is assumed to fire randomly with a mean rate, due to internal noise * javier.orlandigomez@ucalgary.ca sources.For disconnected neurons, this produces a background of incoherent spiking activity.When the neurons (nodes) are connected, each spontaneous firing transmits a simultaneous input to all nodes it is connected to that adds to their internal noise, thus facilitating their firing probability, by virtue of their integrate-and-fire dynamics.This generates cascades of activity, resulting in a distribution of avalanches with powerlaw statistics [12].Here we will refer to the spatio-temporally structured activity induced by the spontaneous firing on the network as "dressed noise." In the experiments in cultures of Ref. [12], neurons were randomly fixed in space, and they established connections to their neighbors by extending their axon in random directions.The network was thus embedded in a metric space, and it was statistically homogeneous and isotropic.The metric correlations inherent to the network generation had a finite connectivity correlation length , of the order of the size of dendritic tree of neurons.Typically, a circle with radius contained of the order of 100 neurons.yields a natural scale for the coarse-graining approach that we develop here.At this scale, the quenched relative fluctuations of homogeneously and isotropically distributed quantities are relatively small.
The key observation that we exploit here is that the amplification of noise through the network is anisotropic, as a result of the sensitivity of the dynamics to the detailed wiring, implying that at each point there is a preferential direction of avalanche growth.Directed propagation of the dressed noise may thus be expected.We postulate that, at the coarse-grained level, there exists a fixed local vector that acts as an advective carrier of the activity.This implies that, in a system that is statistically homogeneous and isotropic at the microscopic (network) level, the quenched fluctuations of the actual network realization are amplified and show up at the mesoscale in the form of an explicit symmetry breaking.

II. MODEL EQUATIONS AND DISCUSSION
For the coarse-grained description we propose a FitzHugh-Nagumo model [13], with a nullcline structure that enables a transition between excitable and oscillatory behavior.A simple diffusive spatial coupling with a diffusion coefficient D will be assumed to account for the spreading of the coarse-grained activity.In addition we must account for synaptic depression, namely, the fact that the strength of each network link is dynamic and decreases upon continued activity, up to complete disabling in saturation conditions.In cultures, the synaptic recovery time τ D is much larger than other time scales (see Appendix A).Accordingly, we consider a state-dependent spatial coupling by introducing an additional depression variable w(r,t) to the u(r,t)-v(r,t) FitzHugh-Nagumo variables of the form where η(r,t) is a Gaussian white noise with zero mean and η(r,t)η(r ,t ) = 0 δ(r − r )δ(t − t ), that stands for the intrinsic noise, and f (u) = au 3 + bu 2 + cu + d.With an appropriate choice of parameters a,b,c,d,κ,γ ,g, and h, we constrain the variables u and v to the interval (0,1).The field u is the coarse-grained local spiking activity, with u ∼ 0 corresponding to the activity of intrinsic noise, and u ∼ 1 to saturation of the network.The recovery field v relaxes slowly compared to u.The spatial coupling is contained in the term N [u] (see below) so it is lost at saturation, i.e., for w ∼ 1, as the synapses are depleted.The nonlinear coupling with u is such that w has an S-shaped dependence on the rate of an incoming spike train (see Appendix A).Accordingly, its stationary value takes the form of a Hill function w = u n /[(βτ D ) −1 + u n ].We take n as sufficiently large (n = 8) to ensure an appropriate delay of the growth of w in response of that of u.
The distinctive feature of our model is the term N [u] in Eq. (1), which is postulated to include the dynamics of the dressed noise, and combines both deterministic and stochastic terms, including a zero-mean multiplicative Gaussian white noise with ξ i (r,t)ξ j (r ,t ) = 2 δ ij δ(r − r )δ(t − t ).This dressed noise introduces two new fields, α(r) and V(r,t).The field α is a source term that accounts for a local amplification of the activity.The fundamental novelty is the vector that we call the avalanche field V(r,t), which enters the equation defining an activity current uV.The physical picture is that the activity that is created by the local source term αu is advected by an effective velocity field characteristic of the network.The noise term ξ (r) accounts for the stochasticity of the underlying avalanches.We choose a white noise for simplicity, implying that the avalanches are sufficiently fast in the coarse-grained time scale of observation [16].Because of the statistical isotropy and homogeneity of the network construction, the ensemble average over network realizations is V 0 (r) = 0 and α(r) = const, so the explicit forms of the 1.Phase-portrait and set of nullclines for the model for zero spatial dimensions.The trajectory of a small perturbation at t 0 is also shown (solid line).(a) Projection of the nullclines in the (u,v) coordinates.The system has a single fixed point (solid circle), and the nullcline for u = 0 is plotted at the stationary value of w = w st (dash-dot line) and also at w = 1 (dashed line).The nullcline for ẇ = 0 (not shown) corresponds to a horizontal line crossing the fixed point, since it is indendent of v. (b) Projection of the nullclines in (u,w) coordinates.The u = 0 nullcline is plotted at v = v st (dotted line).fields V 0 (r) and α(r) encode dynamical effects associated to the quenched fluctuations of each specific network.
Our model can accommodate in principle two different points of view that have been adopted in the literature to interpret the localization of burst initiation: (1) some special regions or groups of neurons (leaders) may be more active and recruit activity in their surroundings up to a critical amount that generates a propagating front, as in the one-dimensional studies [17,18], and (2) the nucleation sites may be sinks of some kind of average flow of activity originated in large basins of attraction [12].The local mechanism (1) can be made consistent with the measured nucleation probability densities (hereinafter, the nucleation maps) assuming a strong inhomogeneity of α.However, this scenario is not consistent with the detailed observation of the preburst activity in silico which, combined with independent tests of nonlocality in both simulations and experiments [12], settled the issue in favor of scenario (2).Consequently, for simplicity we assume α as constant.Upon rewriting αu The mechanism that governs the transition from excitable to pulsating behavior is already captured by the zero-dimensional version of our model.Since w is much slower than u and v, we assume that u and v define a FitzHugh-Nagumo model where the nullclines slowly drift in time following w quasistatically (see Fig. 1).The stable state defining the excitable regime is thus slowly approaching the transition to the oscillatory regime through a supercritical Hopf bifurcation [20] as the term (1 − w)αu grows.We assume that for α = 0 and V = 0 the system remains in the excitable regime.In the complete 2D model, as w relaxes, the transition may occur at the regions of sufficiently high values of α eff , i.e., of −∇ • V 0 .Which one of these peaks will end up nucleating a burst will be ultimately controlled by the multiplicative noise term ξ (r).The precise statistics of nucleation will also be affected by the presence of the advective term −V 0 • ∇u, as V 0 may include a solenoidal component.Once the excitation of u is triggered, the spatial coupling leads to a front that propagates through the whole system bringing all neurons to saturation, as u drags w.The network is thus shut down by the passing of the pulse, and a cycle of synaptic recovery resets the system.

III. CONSTRUCTION OF THE AVALANCHE FIELD
Before showing that this phenomenological model can indeed account for the puzzling observations of Ref. [12], we will gain additional fundamental insights by pursuing a microscopic interpretation that provides a bottom-up constructive definition of V 0 .A microscopic definition of V 0 must involve not only the details of the network wiring but also of the dynamics.The dynamics are indeed quite sensitive to the detailed wiring of the network and are easily recognized in the nucleation maps.In Fig. 2 we show the different nucleation maps for a given network (same adjacency matrix) while varying the time constant τ of the excitatory currents; from 10 ms (left) to 50 ms (right).Slower decaying currents widen the nucleation sites and even allow their merger.The information contained in the adjacency matrix is thus insufficient to pinpoint the exact location and shape of the nucleation sites.It is only through its interplay with the full dynamics that the real nucleation sites are revealed.
The microscopic origin of the symmetry breaking is the anisotropic growth of large avalanches.The difference between the avalanche center of mass and its initial point defines a net displacement and, combined with the duration of the avalanche (or the inverse occurrence frequency if larger than the duration), yields a velocity that reflects the spatio-temporal structure of the avalanche.In the Supplemental Material [21], we explicitly implement this procedure to construct an average velocity field for simple, spatially homogeneous networks.For the general inhomogeneous case at hand, however, an extension of this analysis becomes impractical.Alternatively, we propose a more empiric approach within the same spirit and based on our observations in silico.We assume that the subset of ignition avalanches (IAs), that is, those that end up nucleating a burst, dominates the statistical contributions to V.These are large avalanches that peak at relatively small regions, so that they excite the local percolation fraction that excites a nucleus of the critical size [12].Smaller avalanches are more frequent but more isotropic, and large but not-igniting avalanches are more spread, so we assume that they do not contribute significantly [22].
For each IA j that is participated by a neuron i at position r i , we assign a vector pointing towards the burst ignition point r B j .Adding up for the subset of IAs I i participated by neuron i, out of a total number of IAs N I , and averaging over the N (r) neurons within a disk of radius , we have where is the probability per unit time of occurrence of an IA (for a nondepressed network), and λ is an undetermined dimensionless factor.Equation ( 6) is a hands-on implementation of where P (r ) is the probability density of r participating in a given IA, and P (r B |r ) the conditional probability density of such an avalanche to end at r B .This assumes that avalanches are sufficiently fast to neglect memory effects, and that the activity field u(r) is roughly constant within the range .Accordingly, the nucleation map reads Equation ( 8) exhibits the nonlocal character of the information contained in V 0 (r), which exhibits correlations over distances much larger than , typically of the order of the axon length.The streamlines of V 0 are plotted for an illustrative example in Fig. 3, with the corresponding nucleation map, and together with the streamlines of its irrotational I and solenoidal S components, where V 0 = I + S, with ∇ × I = 0 and ∇ • S = 0.The solenoidal component is typically smaller ( |I|/|S|) ≈ 2), but not negligible, and has to be taken into consideration.The parameter λ does not modify the streamlines but controls the degree of localization of the nucleation map, which is related to the level of network clustering [12].This can be fitted from the corresponding Lorenz curves for the area distribution of the nucleation map.
The above procedure seems to yield and optimal avalanche field to reproduce the distribution and statistics of the nucleation sites.Indeed, we can compare with other simpler alternatives to construct the avalanche fields that do not require  a detailed knowledge of all the avalanche statistics.The most naive procedure would approximate the avalanche field by some local averaging the connectivity matrix.For instance, we can associate to each connection A ij a vector r ij = r j − r i , where r i is the position of neuron i, and then define the vector field for each neuron by where the last sum goes over all the k out i connections.The vector V A i describes the average direction neuron i forms with its output connections, and its modulus is proportional to the distance and number of connections.The vector is not normalized, given that neurons with higher output connectivity must be associated to a higher activity flow.
The avalanche field generated above relies only on the connectivity, and is independent of the underlying dynamics.A different field can be generated if we take into account the information from the averaged background activity.We can generate a weighted matrix B ij from the information obtained by the analysis of the background avalanches (BAs) in the dynamical system.For instance, every time a link participates in a BA its weight is increased by a fixed amount.The full weighted matrix is then normalized so it describes the probability that a given link participates in a BA.Hence, we can define a new avalanche field as A comparison between the different definitions of the avalanche field is shown in Fig. 4 by comparing the nucleation map with the divergence of each of the possible constructions.Note how the first field we defined, the IA-based one, correlates much better with the nucleation sites.The IA-based avalanche field is constructed in such a way that it retains nonlocal information from the whole region that covers the basin of attraction of a given nucleation site, which has an extension defined by the mean axon length.This information is inaccessible to a connectivity-only construction, or from one that takes into account all avalanches.By using only the IA information, which corresponds to large avalanches, and also by using the information of the nucleation point, i.e., where the IA activity ends up concentrating, we are providing each neuron with information that is not available locally.

IV. COMPARISON WITH SIMULATIONS AND EXPERIMENTS
Our model reproduces quantitatively all the observed phenomena in experiments and network-level simulations of Ref. [12], and the additional simulations here performed.Figure 5 gives an illustrative example where both statics (the nucleation map) and dynamics [distribution of interburst intervals (IBIs)] of the network-level simulation are matched by the continuum model.We remark that, even though the preburst activity is highly inhomogeneous and anisotropic, as dictated by the avalanche field V 0 , one nucleated, bursts propagate as a circular front, restoring the isotropy and homogeneity of the medium.In Fig. 6 we show two illustrative examples of propagating fronts.
The agreement between the two nucleation maps, except in regions with poor statistics, is expected given the construction procedure proposed.The degree of focusing can be tuned by λ, whereas the mean IBI and its variance are controlled by τ D and , respectively.Figure 5 shows a series of bursts and identifies the respective initiation sites.The IBI statistics is well reproduced by the coarse-grained model, which also exhibits the random alternation between the nucleation sites.The analysis of long series of bursts shows no significant correlation between the location of successive nucleation events, provided that the IBI dispersion is large compared to the burst traveling time between sites.

V. DISCUSSION
Our results suggest that, generically, a coarse-grained description of the spontaneous activity ψ of a metric network of units with an integrate-and-fire rule and fixed firing rate, assuming no additional dynamics of neither nodes nor links, obeys a balance equation of the form ψ may designate completely different observables, be it density of neuron spikes or, for instance, density of events of rumor propagation in a social network [23].This dressed noise ψ will have to be coupled to the appropriate coarse-grained variables for each case.Our scheme defines the simplest Markovian approximation, corresponding to a Fokker-Planck truncation of a functional Kramers-Moyal expansion of the exact probability P [ψ(r),t] (see Appendix B).The explicit implementation of this procedure is highly nontrivial in the case at hand.Alternatively, here we postulate the explicit form of Eqs. ( 12) and ( 13) on the basis of experimental and simulational data from neuronal cultures.If the characteristic times of avalanches become comparable to the IBI scale, then the Markovian approximation may not be justified and time-correlated noise should be introduced.

VI. CONCLUSIONS AND OUTLOOK
In summary, we have seen that the noise focusing phenomenon can be captured in a continuum coarse-grained model of an excitable medium with unconventional properties, that reflect the presence of an underlying network.The excitable dynamics amplifies the network quenched fluctuations giving rise to an explicit symmetry breaking at the mesoscale.The emerging advective transport tends to concentrate activity in competition with the diffusive transport, which tends to homogenize it.Both mechanisms originate at the network level.Whenever advection dominates, such as in neuronal cultures, we will have noise focusing.Remarkably, the symmetry breaking affects essentially the noise-driven mechanisms, such as nucleation, but not macroscopic aspects such as front propagation.
As in Anderson localization [24], the phenomenon here elucidated results from a modification of transport properties due to disorder, even though the mechanisms are completely different.As opposed to the wave interference of multiple scattering competing with diffusive transport, in our case it is the emergence of ballistic transport that competes with diffusion.Nonetheless, in both cases the localization is a collective (nonlocal) effect of disorder in extended regions, as opposed to eventual trapping in locations with special properties of the disorder.
The phenomenon of noise focusing is generic in relatively simple excitable networks, such as in neuronal cultures, but its relevance to in vivo neuronal tissues remains to be established.From a fundamental point of view, the direct quantitative prediction of the symmetry-breaking avalanche field for a specific network, given a set of dynamical equations, remains a nontrivial open question.where t m is the spiking time, g the strength of the synapse (associated to the receptor density at the postsynaptic terminal) and τ the characteristic decay time of the postsynaptic current.D(t) accounts for STD, i.e., the temporal weakening of the synapses due to repeated stimulation [27].Its temporal evolution can be expressed as [18,28,29] where τ D is the characteristic recovery time associated to the recycling of synaptic vesicles [30,31].The release of neurotransmitters at the synapse as a consequence of firing results in a reduction of D to βD, with β < 1. Subsequent firings will induce additional currents but with a reduced amplitude if the synapse has not had enough time to recover.The total input from its neighbors a given neuron j receives is where the first summation comprises all input connections k in j on neuron j , and the second one all spikes previously generated by neuron i.Note that the subset t m is in general different from neuron to neuron.

Synaptic depression response
In the microscopic model, synaptic depression is governed by Eq. (A5), and the response of this equation to an input of varying frequency (and Poisson statistics) is shown in Fig. 7, where w accounts for the level of depression of the synapses.This kind of dynamic response is the reason behind Eq. (3).

Metric construction of the network
To construct the neuronal network we model neurons as circular cell bodies with fixed diameter φ s = 15 μm.Neurons are then randomly placed on a bidimensional area described by the coordinates (x,y) without any overlap.The total number N of neurons is given by the desired density ρ.
From each neuron on the substrate an axon grows in a random direction following a quasistraight path, as described below, with a total length that is given by a Rayleigh where σ 2 = 900 μm 2 is the variance of the distribution and its value is chosen so that the average axonal length matches the value ∼ 1.1 mm measured in experiments.To mimic axon growth we initially divide the total length into small segments = 10 μm long.The first segment is placed at the end of the neuron cell body with a random orientation.The ith segment is then placed at the end of the previous one, and oriented following a Gaussian distribution around the previous segment given by where θ i−1 is the angle between the segment i − 1 and the origin.σ θ is chosen to obtain the desired persistence length (typically σ θ ∼ 15 • ) of a few hundred microns.The growing process is then repeated until all segments are laid down.We model the dendritic tree of each neuron as a disk, centered at the cell body, with diameter φ d drawn from a Gaussian distribution with mean μ d = 300 μm and standard deviation of σ 2 d = 40 μm.Connections are formed purely by geometric considerations.A connection can be established only when the axon of a given neuron intersects the dendritic tree of any other neuron.The neurons that fulfill this geometric condition will connect with probability α.From this connectivity rule we obtain the adjacency matrix of the network A, where A ij = 1 identifies a connection.The typical range of parameters used to generate the networks are presented in Table II.

APPENDIX B: FOKKER-PLANCK APPROXIMATION OF THE MESOSCALE MODEL
In the derivation of stochastic equations at the mesoscale description of systems with many degrees of freedom, in it customary to search for a Markovian models with the idea the time scales associated to internal degrees of freedom that show up in the form of noise terms are negligible against the slow variables that are kept in the effective description.In our system, the spatio-temporal structure of the microscopic activity is complex and includes many time scales and power laws.However, for simplicity, we may also assume that there is a sufficient separation of time scales so that a Markovian approximation is valid.This is what is implicitly assumed in the equation ψ = α(r)ψ − ∇ • J + η(r,t), (B1) The only consistent truncation of this expansion is to second order, which yields the Fokker-Planck approximation, and ensures continuity of the fields.This procedure defines formally the terms that we should put in the Langevin equation (B1) and (B2).What we do in our model is to postulate an explicit form for the functionals F 1 and F 2 .Using the Ito interpretation of the multiplicative noise, this means that we postulate In the case of the Stratonovich interpretation of the multiplicative noise, the Fokker-Planck equation must be rewritten to account for the so-called noise-induced drift.In our case, we do assume this interpretation, because this corresponds to the case where the internal noise has a finite correlation time that is taken in the limit of being vanishingly small.In our case, we assume that finite but small correlations of the noise are present due to the avalanche structure, that contains indeed a large range of time scales.Consistently, if we take the Stratonovich interpretation, our the definition of the drift terms in Eqs.(B1) and (B2) is the same as in the right-hand side of Eq. (B4), which corresponds to F 1 once the noise-induced drift is removed, since this is effectively incorporated in the noise terms.
FIG.1.Phase-portrait and set of nullclines for the model for zero spatial dimensions.The trajectory of a small perturbation at t 0 is also shown (solid line).(a) Projection of the nullclines in the (u,v) coordinates.The system has a single fixed point (solid circle), and the nullcline for u = 0 is plotted at the stationary value of w = w st (dash-dot line) and also at w = 1 (dashed line).The nullcline for ẇ = 0 (not shown) corresponds to a horizontal line crossing the fixed point, since it is indendent of v. (b) Projection of the nullclines in (u,w) coordinates.The u = 0 nullcline is plotted at v = v st (dotted line).

FIG. 2 .
FIG. 2. Sensitivity of the nucleation maps to the dynamics.Left to right (a-e): nucleation maps for a given network for different values of the time constant of the excitatory currents τ from 10 to 50 ms.

FIG. 3 .
FIG. 3. Nucleation probability distribution function for a given network (nucleation map) and its associated avalanche field.(a) Nucleation map.(b) Streamlines of the associated avalanche field.(c) Irrotational and (d) solenoidal parts of the field.Note that the irrotational part of the field strongly correlates with the nucleation map, identifying the nucleation sites as sinks of activity.The spacial average of the modulus of the irrotational part is typically a factor 2 larger than the solenoidal part.

FIG. 4 .
FIG. 4. Divergence maps for each of the avalanche field definitions.(a) Nucleation PDF of the network for reference.(b) Divergence map from the IA-based avalanche field (−∇ • V I ).(c) Divergence map from the BA-based avalanche field (−∇ • V B ).(d) Divergence map from the connectivity-based avalanche field (−∇ • V A ).The correlation between the IA-based avalanche field divergence map and the nucleation sites is clear, where the zones of maximum (negative) divergence indicate the nucleation sites.The other maps show many regions of maximum divergence that do not correlate with the nucleation sites.

FIG. 6 .
FIG. 6.Two examples of traveling pulses for a system with periodic boundary conditions, S = 5 × 5 mm 2 and density ρ = 400 neurons/mm 2 .The onset time is defined as the first time a point crosses an arbitrary large threshold u th = 0.8.(a) Pulse that originates in the middle right area.(b) Pulse that originates in the middle left area.The waves present an average velocity of V = 120 ± 20 mm/s.

FIG. 7 .
FIG. 7. Depression response of the microscopic dynamics to a Poisson spike train.Steady-state value of the depression of a synapse w = 1 − D, where D(t) is governed by Eq. (A5) in (a) linear and (b) semilogarithmic scale.

TABLE II .
Parameters used to generate the neuronal networks.
r) + ξ (r,t)], (B2) where ψ designate the relevant field and the noise terms η and ξ are Gaussian white noises.Formally, the most general equation for the probability functional P[ψ,t] of a Markovian field variable ψ(r) takes the form of a functional Kramers-Moyal expansion