Topology-dependent coalescence controls scaling exponents in finite networks

Multiple studies of neural avalanches across different data modalities led to the prominent hypothesis that the brain operates near a critical point. The observed exponents often indicate the mean-field directed-percolation universality class, leading to the fully-connected or random network models to study the avalanche dynamics. However, the cortical networks have distinct non-random features and spatial organization that is known to affect the critical exponents. Here we show that distinct empirical exponents arise in networks with different topology and depend on the network size. In particular, we find apparent scale-free behavior with mean-field exponents appearing as quasi-critical dynamics in structured networks. This quasi-critical dynamics cannot be easily discriminated from an actual critical point in small networks. We find that the local coalescence in activity dynamics can explain the distinct exponents. Therefore, both topology and system size should be considered when assessing criticality from empirical observables.

Brain activity displays a plethora of different dynamical states, including bursts, oscillations, and irregular activity.In particular, neural activity exhibits spatiotemporal patterns compatible with the dynamics of a system close to a second-order phase transition.Operating at this regime has been linked to optimal information processing [1,2], maximal dynamic range and sensitivity to stimulus [3,4], longer timescales during selective attention [5], and better stimuli discrimination [6,7].Neural network models demonstrated that short-and long-term synaptic plasticity can self-organize brain dynamics towards a critical point [8,9].
To assess whether the brain dynamics is critical, the activity propagation between the neurons is often mapped to the branching process [1,3,[10][11][12][13][14][15][16].This mapping was motivated by in-vitro observation of outbursts of neural activity known as neuronal avalanches with sizes and durations following power-law distributions with exponents τ = 1.5 and α = 2 correspondingly [10], as expected from the branching process at the critical point.The branching process dynamics is fully characterized by a branching parameter m [17], with m = 1 being the critical value, where each neuron on average activates one other neuron creating a fluctuation-driven regime.From the statistical mechanics point of view, the critical transition at m = 1 belongs to the mean-field directed percolation (MF-DP) universality class, a non-equilibrium phase transition separating absorbing and active phases [18,19].In this universality class, the avalanche-size distribution has a power-law exponent of 1.5.
However, mapping the neural activity propagation to the branching process neglects the role of underlying network topology in shaping the dynamics.The branching process assumes non-overlapping spreading of the activity.In biological neural networks, in contrast, each neu-ron can be simultaneously excited by multiple sources.This phenomenon, known as coalescence renders the independence assumption invalid and reduces the effective branching parameter of the system since some active neurons cannot trigger spikes in already excited neighbors [16].These effects are particularly severe in structured networks, for example as observed in primate cortex [20][21][22][23][24]. Additionally, the theory of critical phenomena predicts that structured, finite-dimensional network topology affects variables such as critical exponents [18,25].At the same time, some studies of neural activity reported exponents deviating from the directed percolation universality class [26,27].Taken together, network topology should be considered when interpreting and modeling the empirical avalanche statistics.
To investigate the relationship between network topology and critical dynamics, we developed a finite-size branching network model with various connectivity structures, ranging from spatially arranged networks resembling the local connectivity structure of the cortex to random or all-to-all connectivity.Using this model, we show how the network topology affects the critical branching parameter, avalanche-size distributions, and their critical exponents.We find that the non-critical avalanche-size distributions can appear as critical in finite networks, so a quasi-critical (i.e.not truly critical but expressing some amount of scaling) finite network can be confused with a mean-field critical one.
The model consists of binary units positioned on a twodimensional (L × L) square lattice with periodic boundary conditions.The connectivity is defined by two parameters: connectivity radius k and the probability of rewiring p rw .First, each unit is connected to all units in its k-Moore neighborhood ((2k + 1) 2 − 1 neighbors).Then, with probability p rw , each connection can be se- lected for rewiring and then rewired to a uniformly chosen random location.To study the impact of different network topologies on the dynamics, we consider a network with k = 1 and p rw = 0 and then systematically increase either the radius (Fig. 1, top) or the rewiring probability (Fig. 1, bottom).For k = L/2, a network with size N = L 2 will be all-to-all connected, without any structure.For p rw = 1, we obtain a completely randomly connected network.Both limit cases correspond to the usual mean-field configuration.Each unit i in the network transitions stochastically between an active state s i = 1 and an inactive state s i = 0, depending on the connectivity and external input [28]: where p r is the probability to be excited by the active neighbor, p ext is the probability to receive external input, Ω i represents the set of neighbors of the i-th unit, and p s is the probability to maintain the active state.This model is inspired by the activity and interactions in the cortex: the units represent cortical columns, and the active and inactive states correspond to transient high and low levels of activity found in primate visual cortex [29,30].The probability p s accounts for recurrent interactions between neurons within one column and p r represents recurrence between columns (more details in [31,32]).For simulations, we assume there is no external input (p ext = 0) and take p s = 0.5.The model is a spatially structured version of a branching network (BN) [3,33].The BN with random connectivity is completely described by local branching parameter m = p s + |Ω i |p r , summing all the outgoing connection probabilities of one node.On average, the number of active units A t at time t when there was a single active unit at the previous time step is By taking p r = (m − p s )/ |Ω i | and taking p s = 0.5 we re-parameterize the model in terms of m.
We first analyze the location of the critical transition depending on network topology.The local branching parameter m is used as the control parameter.The exact location of the critical transition can be found by several methods (Appendix A).Here, we look for the critical branching parameter m c that maximizes the variance of the activity χ(ρ(m)).In the mean-field system, the absorbing-active transition happens at m c = m M F = 1 [19].However, we find that the location of the critical transition depends on the network topology, in agreement with the theory of critical phenomena [25,34] (Fig. 2).For the structured connectivity (k = 1), the phase transition happens at a larger critical branching parameter (m c = 1.109).As we move towards the mean-field connectivity (either all-to-all, or random), the critical branching parameter gradually converges to m M F = 1 (Fig. 2).In the structured networks, we refer to the dynamics at the mean-field branching parameter (m M F = 1) as quasi-critical.We will demonstrate in the following that the quasi-critical dynamics presents the classical mean-field scaling with an exponent close to 1.5 but only for the limited range of the event sizes.Thus, m M F = 1 is not an actual critical point in structured networks.
The location of m c in structured networks is modeldependent and not a universal feature.For example, in the two-dimensional contact process (CP) the critical point is located at m CP ≈ 1.6 [19].In contrast, in our structured model (k = 1) criticality appears around m c ≈ 1.1 (see Appendix D for a mapping between our model and a continuous time CP model).
Next, we compare the avalanche-size distributionsoften seen as the primary indicator for critical behavior in neural data-between the quasi-critical (with m M F ) and critical (with m c ) networks for various network structures and sizes.An avalanche is a cascade of activity propagation in the network.It starts with an external input activating a single neuron in a quiescent network and ends when the activity dies out.At criticality, avalanche sizes and durations follow a power-law distribution with an exponential cutoff whose location scales with the network size.
Quasi-critical networks exhibit apparent scale-free avalanche-size distributions with the expected MF-DP power-law exponent (τ ≈ 1.5, Fig. 3, top).In particular, small quasi-critical networks with a finite interaction radius (e.g., k = 3 and L = 8, 32) can seemingly display finite-size scaling.The apparent scaling is more visible in networks with larger connectivity radius k (Supplementary Fig. 1).However, for sufficiently large system sizes a characteristic scale becomes evident (cutoff stays independent of the system size when L > L scale ).The characteristic scale in quasi-critical networks becomes more apparent when compared to the avalanche-size distributions of critical networks (with branching parameter found in Fig. 2) with the same size and topology (Fig. 3, bottom).In critical networks, power laws extend up to much larger sizes.For nearest-neighbors connectivity the exponent shifts from the mean-field values of τ ≈ 1.5 towards τ ≈ 1.27 as expected for the two-dimensional directedpercolation (2D-DP) universality class [25] (they approach τ = 1.27 as N → ∞, Supplementary Fig. 3,  Appendix C).Changing the network topology towards random or all-to-all connectivity brings the critical point closer to m M F = 1 (Fig. 2).Hence, with larger k, the characteristic scale of avalanche sizes in quasi-critical networks increases and the exponents shifts towards 1.5 (Fig. 3, Supplementary Fig. 4).Quasi-critical avalanches follow a power law up to a cutoff that is scaling with the system size for small systems (see also Supplementary Fig. 1) but exhibits a characteristic scale for large systems (a,b).Gray lines indicate the fitted power-law distribution with the exponent τ for L = 128 (see Appendix B for fitting details).Critical branching parameters are defined by maximum variance point for each system size, as in Fig. 2.
Increasing the connectivity radius or rewiring probability in critical networks changes the critical exponents continuously from τ ≈ 1.27 (the 2D-DP universality class) to τ ≈ 1.5 (the MF-DP universality class).However, these two mechanisms affect the critical exponents in different ways.The major difference stems from their behavior in the thermodynamic limit.For any finite connectivity radius k and no rewiring, in the limit of large network sizes, connections are short-ranged and the dynamics belong to the 2D-DP universality class with the critical exponent of τ = 1.27.At the same time, finite networks with large enough k are almost fully connected, showing exponents similar to MF-DP.Thus, for fixed k, the network size affects how close the system is to a fully connected system.The combination of these two factors leads to the true scaling exponent (known for the 2D-DP and MF-DP [25]) being visible only for very large avalanches (S 1), which require very large system sizes (L → ∞).In large networks, the power-law exponent will change slowly and continuously from MF-DP for relatively small avalanches to the 2D-DP exponent for large events (Supplementary Fig. 3).The rewiring, on the other hand, conserves the network topology independently of the system size.Therefore, the critical exponents increase with increasing rewiring probability, approaching the mean-field limit.For a fixed p rw , avalanche statistics follows a well-defined power-law distribution with a single exponent, which displays the true scaling also for small event sizes and durations.The difference between both cases can be clearly understood in how the effective dimension of the network changes with increasing k or p rw (Appendix C).
Our results suggest that observing the power-law-like behavior in avalanche-size distributions from small networks cannot be a reliable signature of criticality even when the distributions scale with the system size.Close to criticality, subcritical systems can exhibit apparent scale-free behavior, where the characteristic scale is only uncovered in the limit of large systems.
We next investigate the mechanisms underlying the differences between the avalanche dynamics in mean-field and structured networks.Due to locally structured connectivity in our model, each unit can be activated by multiple sources at the same time, generating coalescence (Fig. 4a).We measure the network coalescence from the ongoing activity in the timescale-separated regime (p ext = 0), where each new avalanche is initiated after the previous one is over.Each unit i in the network can be simultaneously activated by multiple neighbors, or by the self-excitation (Fig. 4a).We define the local coalescence as the number of sources that activated unit i minus one.Let n i,t be a number of active neighbors of unit i at time t.The local coalescence of unit i at time t is a random variable where b k ∼ Bernoulli(p r ) and b 0 ∼ Bernoulli(p s ).Let A t be the number of active units in the network at time t.Then, the average normalized network coalescence C(A) is given by where average is taken over all the times with A t = A. Due to coalescence, the effective branching parameter is smaller than the local branching parameter and depends on the network topology and the number of active units.We can estimate the effective branching parameter for the given number of active units A from the simulated activity as: The estimated effective branching parameter satisfies m eff (A) < m but varies depending on the number of active units A (Fig. 4b,c, left).For small A, m eff (A) values are closer to m, but with increasing number of active units, the coalescence also increases leading to larger deviations of m eff (A) from m.Therefore, a larger local branching parameter (m = 1.109) is required to have an effective branching parameter close to 1 that creates critical dynamics.When increasing the connectivity radius or the rewiring probability, the activity can spread to a broader range of units.In addition, increasing the connectivity radius while keeping the branching parameter constant reduces the strength of individual connections, and coalescence for every unit scales as a power of the connection strength.Thus, by increasing the connectivity radius or the rewiring probability, the coalescence decreases, and the m eff (A) becomes closer to m (Supplementary Fig. 4, consistent with Fig. 2).
To capture the impact of coalescence on avalanchesize distributions, we simulated an equivalent adaptive branching process for each network model.To generate this process, we first find the m eff (A) from a long simulation of the network model using Eq. ( 5).Then we define the adaptive branching process, as a Markov process Ã(t), where each of the ancestors can generate a binomially distributed number of offsprings z i , n is a maximal number of offsprings for one ancestor that is equal to the number of neighbors (2k + 1) 2 − 1 in the original network.For the simulations and figures in the paper we took n = 8.Avalanche-size distributions of the adaptive branching processes well approximate the shape of corresponding distributions from the network dynamics (Fig. 4b,c, right).In the quasi-critical network, distributions for the branching process and the network are completely overlapping.The adaptive branching process for the critical network captures most of the distribution, except for some deviations in the tail, and it has a similar power-law exponent.The mismatch in the tail can be generated by a too scarce sampling of large avalanches to estimate the correct effective branching parameter.Overall, despite the large variability in the network coalescence, the average value of effective branching parameter is sufficient to predict the shape and exponent of the avalanche-size distribution in the structured networks.
The network coalescence can be estimated analytically for a branching network with random connectivity [16].However, analytical determination of coalescence for finite-dimensional topologies (e.g., structured networks) requires renormalization group approaches relying on the precise knowledge of the system dimension [18,19].This approach becomes especially difficult when dealing with finite system sizes.Hence, we used a simulation based approach to estimate the coalescence from the network activity.
Our results indicate that the differences between the dynamics in structured and mean-field networks arise from the coalescence created by the local network interactions.Increasing the connectivity radius or rewiring probability reduces the coalescence and the network dynamics becomes more similar to the conventional branching process.Dependence of coalescence on network topology is in line with the previous observation that the effective branching parameter reduces with increasing degree of network clustering [35].Although there have been many studies addressing the effect of long-range connections in the brain [36], and critical-like avalanches have been reported in small-world [37] and scale-free [38] net- works, there was no systematic study of the impact of varying coalescence on avalanche distributions.
We demonstrated how a range of scaling exponents can arise from the network structure and be misinterpreted in finite-size networks.This is particularly significant for interpreting observations from finite-size neural recordings, where the number of recorded neurons or electrodes are often treated as a proxy for system size, and finitesize scaling is assessed by down-sampling the recorded neurons (electrodes) [39].The precise number of neurons involved in the dynamics and their interaction radius is often unknown, and approaches which can recover the dynamical regime from subsampled networks [40] so far mainly rely on the down-sampling of available data.
We showed that while increasing the connectivity radius and rewiring affect the dimensionality of the network in different ways, both mechanisms can create exponents in between the 2D-DP universality class and the meanfield values.At the same time, critical exponents other than the MF-DP have been recently observed in neural activity [26] and linked to a phase transition at the onset of collective oscillations in a network model with excitatory and inhibitory neurons [26,41].Moreover, it was shown analytically that the absence of separation of timescales in measuring avalanches could lead to exponents between directed and standard percolation universality classes [42], or in a particular case to 1.25 [43].How mixture of different mechanisms affects observed avalanche statistics should be studied in the future to better understand the underlying mechanisms for variable exponents in neural activity.

FIG. 1 .
FIG.1.Branching network model with different types of spatial connectivity.Each unit in the network represents a cortical column (left) that can be excited by the self-excitation (probability ps, blue arrows), its neighboring units (probability pr, orange arrows), and the external input (probability pext, gray arrow).We consider networks with connectivity structures ranging from local, structured (left) to mean-field (random or all-to-all, right) generated via two pathways: by increasing the connectivity radius k (top) or rewiring local connections to random with the rewiring probability prw (bottom).

FIG. 2 .
FIG. 2. Location of the critical transition depends on the network topology.Criticality occurs at the transition to the non-zero mean activity ρ (a, b) and maximal variance χ(ρ) (c, d), represented by vertical dotted lines.With increasing connectivity radius k (a, c) or rewiring probability prw (b, d), the critical branching parameter moves to the mean-field value (vertical gray dashed line).For simulations L = 128.

FIG. 3 .
FIG.3.Avalanche-size distributions differ in quasi-critical and critical networks.System-size dependence at quasicriticality (a, b) and criticality (c, d), for k = 1 (a, c) and k = 3 (b, d).For critical systems, the cutoff location of the avalanche-size distribution shifts with system size (L 2 ) (c,d).Quasi-critical avalanches follow a power law up to a cutoff that is scaling with the system size for small systems (see also Supplementary Fig.1) but exhibits a characteristic scale for large systems (a,b).Gray lines indicate the fitted power-law distribution with the exponent τ for L = 128 (see Appendix B for fitting details).Critical branching parameters are defined by maximum variance point for each system size, as in Fig.2.

FIG. 4 .
FIG. 4. The adaptive branching process captures the shape of the avalanche-size distributions in structured networks (k = 1, L = 128).(a) Avalanches (gray frames, each row marked with an empty circle represents one network unit, filled circles are active units) are separated by quiescent moments (white frames).Two types of coalescence: simultaneous activation by multiple active neighbors (brown square), or by the self-excitation and an active neighbor (orange square).Each avalanche starts with the external input (black lightning bolts).The arrow indicates time.(b, left) The effective branching parameter (m eff (A), brown line) in the quasicritical (m = 1) network as a function of the number of active units (A, Eq. (5)) deviates from the local branching parameter (dashed line).The shading indicates ±1 s.d. of Ci,t|A.(b, right) The adaptive branching process (dashed brown line) has the same avalanche-size distribution (P (S)) as the structured branching network (BN, yellow line).The gray line shows the power-law fit with the exponent τ .(c) Same as (b) for the critical network (m = 1.109).