Unique scales preserve self-similar integrate-and-fire functionality of neuronal clusters

Brains demonstrate varying spatial scales of nested hierarchical clustering. Identifying the brain’s neuronal cluster size to be presented as nodes in a network computation is critical to both neuroscience and artificial intelligence, as these define the cognitive blocks capable of building intelligent computation. Experiments support various forms and sizes of neural clustering, from handfuls of dendrites to thousands of neurons, and hint at their behavior. Here, we use computational simulations with a brain-derived fMRI network to show that not only do brain networks remain structurally self-similar across scales but also neuron-like signal integration functionality (“integrate and fire”) is preserved at particular clustering scales. As such, we propose a coarse-graining of neuronal networks to ensemble-nodes, with multiple spikes making up its ensemble-spike and time re-scaling factor defining its ensemble-time step. This fractal-like spatiotemporal property, observed in both structure and function, permits strategic choice in bridging across experimental scales for computational modeling while also suggesting regulatory constraints on developmental and evolutionary “growth spurts” in brain size, as per punctuated equilibrium theories in evolutionary biology.

www.nature.com/scientificreports/ detect a pair of experimental scales and techniques (e.g. calcium imaging combined with multi-electrode array recording) that are amenable to simultaneous modeling by virtue of sharing functional organizing principle. Such an estimate of paired similarly-behaving scales would guide the choice of imaging modalities to inform multi-scale modeling efforts aimed at maximal coverage of brain's dynamical repertoire. Some understanding of the dynamical repertoire of brain's larger than neuron structures, and particularly their repetitive and self-similar nature has been obtained: synfire chain 11,12 describes the sequentially wired sub-assemblies propagating synchronous activities; avalanche models describe the power-law statistics of sizes of events in multielectrode recordings 15 ; power-law noise statistics in brain signals have been robustly characterized [16][17][18] . Also reliably documented is the scale-free structural properties of brain network 19 . Yet the rules unifying structures and processes in brain at distinct scales are not fully described. Existing computational methods for re-scaling the neuronal network focus on enforcing known first and second order statistics of a finescale simulation to its coarse-grained version 20,21 . These approaches provide insight into statistical conditions that need to be satisfied for improved simulation runtime complexity.
The most basic functionality of brain is that of spiking neurons. Neurons, however variable, act as non-linear summing integrators, with a decay in potential and a threshold for firing; they spike when excited frequently and strongly enough. Here, we propose a scale-invariant dynamic, where a pair of similar views of the same object are separated by a specific re-scaling factor that transforms one into another. A simple analogy is provided by Sierpinski triangle 22 , where a finer stage is obtained from the previous coarser stage by halving the spatial dimension of each triangle and tripling the number of triangles, thereby characterizing each stage of re-scaling of the fractal object by two integers: (2,3). Here, we propose a functional scale-invariance or "fractality, " such that a redefinition of the spiking network's basic primitives: node, edge, and spike as a particular multiple of themselves allows one to observe the neuron-like signal-integration functionality re-emerge in the coarser structure, albeit only at unique re-scaling. The consequence of this proposition is that multi-scale models can coarse-grain for computational expediency while strategically choosing a unique set of scales specifically determined to permit translation. Thus, we expand on previous studies, which have demonstrated the feasibility of hierarchical clustering 23 and networks' scale-free features [24][25][26] with regard to brain structure to introduce self-similarity of dynamics of clusters of neurons, replicating integrate-and-fire-like functionality across scales.
To computationally test whether dynamics of signal transmission are preserved in a manner akin to "functional harmonics" of the original scale, we construct a spiking network that follows actual brain connectivity. The starting point of our study is an all-to-all resting-state human fMRI-derived functional connectivity matrix, extracted from 91,282 voxels, and providing full coverage of the human cortex 27,28 . The fMRI data describe the temporal correlation of each of these voxels with every other voxel, allowing an inference about the degree to which they are functionally connected. Such inferences have been externally validated by comparison to structural and DTI data 29 . Each voxel of this matrix is (2 mm) 3 , and thus measures a compensatory hemodynamic (blood oxygen level dependent, or BOLD) response across ~ 1000 cortical minicolumns, or about ~ 1 M neurons. Each voxel is represented as a node in our initial network, and in our computational experiment it follows the functionality of a leaky integrate-and-fire neuron. We start our computational experiment from an already non-neuronal scale of data acquisition, as we suggest that the observations made on fMRI scale may be close to one of the stages of the potential network re-scalings, along which the signal integration property is preserved.
The initial fine-grain network is reduced into clusters using a streaming hierarchical clustering algorithm 30 . We chose a streaming clustering algorithm for its computational expediency and reduced memory requirement. In deciding which nodes combine to form a cluster, we chose full-linkage clustering, requiring that only sets of nodes that have all-to-all connection weights higher than the cut-off c become a cluster of nodes (an ensemblenode). All pairs of neurons i, j belonging to the same cluster satisfy: w ij ≥ c , where w ij is the connection weight between nodes i and j . By adjusting cut-off c one can flexibly control the average number of neurons in an ensemble-node. This clustering procedure allows re-scaling to a continuous spectrum of size of networks ( Fig. 1A-C).
The set of edges running between two neighboring ensemble-nodes aggregate to give the ensemble-edge connecting these ensemble-nodes. The weight of the ensemble-edge is given by the mean weight of the edges defining it: w coarse kl = mean(w ij ) where i runs over all neurons in cluster k and j runs over all neurons in cluster l . For each ensemble-node we now define the node strength as the sum of weights of edges connecting to it: s coarse k = l w coarse kl , where superscript "coarse" indicates that the quantity pertains to the coarse-grained network. We next redefine the activity of the coarse-grained network by temporally clustering firings of an ensemblenode into an ensemble-spike. An ensemble-node consisting of N N nodes produces an ensemble-spike if its nodes produce a combined burst of above-threshold ( N S ) number of spikes in the time-step, and no ensemble-spike occurs otherwise. This characterizes the number of ensemble-spikes required to arrive in rapid succession into an ensemble-node prior to activation of its own ensemble-spike. Intuitively, only areas dense enough in spikes in fine level spike raster (Fig. 1D) become ensemble-spikes in the coarse network's ensemble-spike raster (Fig. 1E) in the corresponding ensemble-node and time step. A spiking neuron's refractory period also re-emerges in the coarse-grain network, and concomitantly guides our choice of ensemble-step: the time step re-scaling parameter. Further details are provided in the Methods section. The terms "cluster" and "ensemble-node" will be used interchangeably as well as the terms "burst and "ensemble-spike. "

Results
Coarse-grained structure is self-similar. Coarse-graining at various cutoffs produces networks with a range of ensemble-node sizes and number of ensemble-nodes ( Fig. 2A). Importantly, cutoff has little impact on the tail of ensemble-nodes' size distribution, which is invariant with respect to mean count and powerlaw behavior, suggesting self-similarity in agreement with a body of previous work 19,31,32 (Fig. 2B). The original functional brain network exhibits a self-similar property when it comes to its new coarsened structure: the www.nature.com/scientificreports/ connectivity network preserves the geometry of its node strength distribution after the coarse-graining procedure (Fig. 2C). We also provide a comparison to the coarse-graining of a randomly re-wired network (while preserving the strengths/degrees of nodes), and note a major difference: the distribution of cluster sizes as the randomly re-wired version of the original network is coarsened quickly loses the heavy tail of the power-law distribution within half an order of magnitude (Fig. 2E). See Methods, Re-wired Controls section for details of the randomization. The strength distribution of the network, on the other hand preserves its shape (Fig. 2F) in agreement with the original network's behavior under coarse-graining. The dependencies of the mean cluster size and the number of clusters as a function of the cut-off threshold (Fig. 2D) are included for validation and comparison with Fig. 2A. www.nature.com/scientificreports/ Coarse-grained dynamics is self-similar. Re-scaling our network to a spectrum of average ensemblenode sizes ( N N ) reveals signature of an integrate-and-fire property to a varying degree depending on the re-scaling factors ( N N and N S ) (Fig. 3A, solid lines), but with three key features. First, leading up to an outgoing ensemblespike, the cross-correlogram of ensemble-spiking activity weighted by the ensemble-edge weights (see Methods for full description) displays integration of a faster-than-exponentially rising amount of inputs from neighbors. Second, the incoming inputs peak at the time-step immediately preceding ensemble-spike output. Finally, incoming activity from a given ensemble-node's neighbors abruptly drops at the time of an outgoing ensemble-spike.
Self-similar dynamics is preserved only at specific scaling factors displaying "functional harmonics" of the original scale. This scale invariant dynamics results from the brain's organizational structure. Comparing the sharp increase in signal integration preceding the ensemble-spike observed in the clustered The intuition for integration coefficient (IC)-the degree, to which cross-correlogram of well-clustered brain's ensemble-spiking activity P(τ ) exceeds that of randomly clustered brain P random (τ ) is indicated (shaded area) and acts to quantify the neuronlike behavior of ensemble-nodes (neuronal clusters). See Methods for definition of IC and cross-correlograms. Shown for ensemble-spikes defined as N s = 2 spikes (top sub-panel) and N s = 4 spikes (bottom sub-panel). (B) (Left) Heatmap of the Integration coefficient (IC) determines the ensemble-neuron size N N and ensemblespike size N S used in re-scaling that best reproduces neuron-like behavior when clustered starting from a brain network. (Right) As a control, we show that the same signal integration behavior is lost when brain network is coarse-grained after a random shuffling (see "Methods" for procedure) that preserves strengths of all nodes (strength-sum of node's edge weights). This rules out the node strength sequence of brain network as the driver of its neuron-like behavior at the coarse scale. (C) The temporal re-scaling factor N t is obtained by calculating the auto-correlogram P auto (τ ) (see Methods for definition) and maximizing the ratio P auto (−2)/P auto (−1) to reveal refractory period in coarse-grained network. (D) Local clustering coefficient among nodes within an ensemble-node corresponds with its integrate-and-fire behavior. Positive dependence indicates that highly interconnected set of nodes can join into an ensemble-node with more neuron-like integration of input. www.nature.com/scientificreports/ network to that of a control network that was coarse grained to an identical average size of ensemble-nodes but through random clustering (ignoring the edge weights), we see a control behavior that is only exponential (Fig. 3A, dashed lines), indicating its weak input integration. Shaded areas in panels of Fig. 3A indicate the extent of the coarse network's input integration beyond that of an equivalently sized randomly clustered network. The top panel (ensemble-spike requiring ≥ 2 spikes) shows little difference from the random case, while the bottom panel (ensemble-spike requiring ≥ 4 spikes) shows that inputs of ≥ 4 spikes are integrated much more rapidly leading up to an ensemble-spike (solid line) than in a randomly clustered network (dashed line). From the starting network of 7784 nodes, by searching the space of re-scaling factors and maximizing the integration coefficient (see "Methods" section for details), we find ensemble-node size of ~ 10 nodes and ensemble-spike size of ~ 5 spikes to maximally recover integration of inputs leading up to ensemble-spikes (Fig. 3B left). Integration coefficient measures the degree to which the cross-correlogram of the coarse network exceeds what is randomly expected. We take the integral of the ratio P(τ )/P random (τ ) weighted by the inverse of absolute value of time delay τ integrated over the range of delays where the ratio exceeds 1. Here P(τ ) and P random (τ ) are the burst cross-correlogram of the coarse-grained network for well clustered (full-linkage clustered) and randomly clustered cases, respectively. In searching for the scale, to which the network coarse-grains to while preserving the most integrate-and-fire behavior, we sought to establish the triplet of re-scaling factors that maximize our measure of similarity to the integrate-and-fire functionality: (number of nodes in a cluster, number of spikes in a burst and the time rescaling factor). We establish the time-axis re-scaling factor, or ensemble-step, of ~ 4 by maximizing the refractoriness present in auto-correlogram (Fig. 3C). The refractoriness is defined as the ratio P auto (−2)/P auto (−1) of the auto-correlogram of the coarse-grained network at two time-steps prior to the burst, to that at one time-step prior to the burst. Thus, the network coarse-graining procedure returns a triplet of values extracted by maximizing the integration coefficient and the refractoriness, consisting of (i) ensemble-node = 10 nodes ( N N = 10 ), (ii) ensemble-spike = 5 spikes ( N S = 5 ), and (iii) ensemble-step = 4 time steps ( N T = 4 ) that recovers properties of spiking nodes in a network of ensemble-nodes. Searching for the culprit of the integrateand-fire behavior in the coarse-grained network, we extracted for each cluster after re-scaling i) the clustering coefficient of the fine-grain neuron-level network and ii) the integration coefficient from the neurons in the cluster. We observe that the within-cluster clustering coefficient is positively correlated with the signal integration behavior of the clusters (Fig. 3D). Thus, ensemble-nodes with greater local clustering coefficient display greater signal-integration behavior, and the dependence does not deviate from a linear trend.

Discussion
We started the research looking for the cluster size that acts as a minimal cognitive building block, and our fMRI-driven interrogation of simulated spiking network data suggests that brain organization has a dense multiscale property in structure and functionality. Neurons naturally cluster across many scales, and such clusters create new networks with fewer functional units, but similar algorithmic features. This suggests a more complex answer. Our coarse-graining procedure attempts to capture only one such transition between a pair of scales by detecting similarity of dynamics, but suggests that the brain potentially preserves its dynamics across scales at many distinct levels, recycling already existing functional building blocks. It is therefore feasible that brain's least functional unit larger than a single neuron might be a cluster of neurons (an "ensemble-neuron" in our terminology) displaying a neuron-like integrate-and-fire behavior.
Testing whether neuronal clusters act as integrating neurons and, if so, at what re-scaling factor, can be established by a targeted experiment of large-scale calcium imaging and other optical recording modalities. Such an experiment would, having simultaneously recorded from large number of neurons, apply the machinery we describe, identifying a pairs of coarse and fine scales that display the most self-similar behavior. Such a constraint on the similarly behaving pair of scales would drastically narrow the space of multi-scale models of brain dynamics. Support to our work comes from previous renormalization-inspired models of neuron-level data coarsegraining. Bialek and team proposed a renormalization procedure that groups variables (neurons) into larger-scale lumped variables based on their observed pairwise correlations 33 . Their procedure reveals the self-similarity of neuron level activities in live animal recordings, but does so without identifying the re-scaling factor that separates similarly behaving scales or whether such a factor exists. They instead estimate a set of renormalizationinspired scaling parameters from the data. Shew and team carry out a renormalization of binarized membrane potential signal from both simulation and rodent data for a coarse-graining similar to ours: their procedure is also partially limited by a single renormalization, rather than a sequence, as the authors point out 34 . In addition, their procedure yields scaling laws, but not constraints on discrete ratios between similarly behaving scales. Instead of finding the scale at which coarse-graining reveals the most self-similar behavior, their procedure probes the dependence of the outcome of coarse-graining on the simulation parameters, presenting an insight complementary to ours.
In clarifying our results, we note that the slope at the tail of the distribution of cluster sizes is not fully unchanged as a function of the system size. However, as a comparison to a case where no signature of similarity of scales is present, we show that for a 1.23 decade ( log 10 (23251/1372) ≈ 1.23 ) of system sizes, the slope changes to a noticeably smaller degree than in a randomly clustered case, where the slope change is much larger during for the 0.62 decade ( log 10 (3418/812) ≈ 0.62 ) of sizes. We use the term "scale-invariance" in the sense of a pair of spatiotemporal scales that are displaying similar dynamics, rather than a stricter definition in which a system displays identical behavior across all accessed scales.
Our primary finding that there exists a particular re-scaling factor that best recovers the integrate-and-fire principle the most dynamic should not be confused with detection of the same dynamic across all accessed scales. The term scale-invariance is used in this paper to describe a pair of spatiotemporal scales that are displaying www.nature.com/scientificreports/ similar dynamics, instead of the sometimes used stricter definition where the system displays identical behavior at all scales accessible to it. We propose a future direction to strengthen the results. In the work reported here, the slope at the tail of the distribution of cluster sizes is not fully unchanged as a function of the system size. Rather as a comparison to a case where no signature of similarity of scales is present, we show that for a 1.23 decade (log 10 (23251/1372) ≈ 1.23) of system sizes, the slope changes to a noticeably smaller degree than in a randomly clustered case, where the slope change is much larger during for the 0.62 decade (log 10 (3418/812) ≈ 0.62) of sizes. An additional limitation of our work is partially a matter of terminology. Here we use the term scale-invariance, in the sense of a pair of spatiotemporal scales that are displaying similar dynamics, instead of the stricter definition where the system displays identical behavior at all scales accessible to it. The implications of our work are particularly important for neural mass models (NMMs). We add to existing experimental and analytical motivations for NMMs [35][36][37][38][39][40] and propose a mechanism for constraining scales to be modeled with NMMs. NMMs provide an attractive procedure for reducing brain dynamics to just a few quantities at the population level. Further supporting the treatment of neuronal systems as NMMs, we provide computational justification for approximating the fine-scale neural dynamics (namely integrate-and-fire rule) by a similar dynamical rule, but at a particular coarser scale out of all scales empirically available for imaging. One implication of such an approximation is that brain dynamics may be modeled using experimental data from several non-overlapping scales, as long as those scales are chosen strategically to retain self-similar dynamics. In addition, our analysis coarse-grains a network of neurons while preserving integrate-and-fire dynamics at network sizes exceeding 1000 neurons and ensemble-neurons, highlighting the potential of NMM-like models to describe interactions of a large number of neuronal populations.
From a biological perspective, the periodicity of scale invariance may offer clues into developmental and evolutionary processes. The dynamic structure of the action potential (gating, excitation, inhibition, return to baseline) reflects, at its very basis, a mechanistic solution to the biologically ubiquitous problem of allostatic regulation in response to noisy inputs. In growing the brain from unicellular to multicellular to agglomerate structure-while continuing to maintain allostatic regulation-it therefore makes sense that dynamic signatures associated with such regulation (i.e., those associated with negative feedback loops) would also be preserved. As such, the discontinuity of scales at which allostatic regulation continues to hold may impose key functional constraints on "growth spurts" in brain development [41][42][43] , as per punctuated equilibrium theories in evolutionary biology 44,45 .

Methods
Network structure. The starting point of our study is an all-to-all resting-state fMRI-derived (N = 68, repeated 4 times) functional connectivity matrix, extracted from 91,282 voxels, and providing full coverage of the human cortex 27,28 . Each voxel is (2 mm) 3 , and thus measures a compensatory hemodynamic (blood oxygen level dependent) response across ~ 1,000 cortical minicolumns, or about ~ 1 M neurons (in primates and most mammals, the cortex contains a 20/80 ratio of neurons to glia 46 ). Each voxel is represented as a node in a graph.
Nodes. The nodes follow a leaky integrate-and-fire neuronal model, updated via change in membrane potential: dV k /dt = −τ vVk + j w jk × I j + σ k with spike conditions: if V k (t) > V θ , then the node initiates a spike by setting: V k (t) = V r . For intuitive order of magnitude estimation, the quantities tracked and the parameters are implemented in dimensionless units, with the range from the reset membrane potential and the threshold potential equal to 1. Other scaling parameters follow the same units. Here: V k (t) is the time-dependent membrane potential (voltage) of a neuron, with k indexing the neuron integrating signal from its neighbors indexed by j.
I j indicates the input signal from pre-synaptic neurons, indexed by j . I j = 1 during the time-step neuron j fires and 0 otherwise. w jk is weight connecting the spiking neuron j to the signal-integrating neuron k . Equal to the gain ( 7.7 × 10 −6 ) times the connectivity weight.
V r is resting membrane potential or reset value, to which the neuron returns after spiking Equal to 0 due to the dimensionless units.
V θ is firing threshold for an action potential, upon reaching which the neuron initiates a spike. Equal to 1 due to the dimensionless units.
τ v is rate of leakage in membrane potential from V k (t) towards V r . Equal to 10 −4 . σ k is the normally distributed noise input to neuron k.

Integration coefficient and normalization by random controls.
We compare the well-clustered coarse network to a randomly clustered control to reveal effects of strong intra-cluster connectedness on coarse network's input integration. In order to quantify the effects of integration of incoming ensemble-spikes from neighbors in a coarse-grained network, we obtain the cross-correlogram of ensemble-spikes: • For a particular ensemble-node, for each observed ensemble-spike, we collect the most recent ensemblespike of each of its neighboring ensemble-nodes according to their time lag τ , by which they preceded the ensemble-spike. For each time lag τ , we then sum ensemble-edge weights for all ensemble-spikes that preceded our ensemble-spike of interest by τ time steps. Any neighboring ensemble-node contributes only one ensemble-spike to this list of sums-its most recent ensemble-spike preceding the ensemble-spike of interest, and to the sum of weights at observed time-lag τ . This tells us how much information was integrated by the ensemble-node and from how many time steps prior to producing the ensemble-spike. www.nature.com/scientificreports/ • For each time lag τ , we then sum the ensemble-edge weights collected over all ensemble-spikes of all ensemble-nodes, to obtain the (ensemble-edge-weighted) cross-correlogram P(τ ) dependent only on time-delay τ . P(τ ) tells us how much input (ensemble-spikes weighted by ensemble-edge weights) had to be received from neighboring ensemble-nodes τ time steps prior to a typical ensemble-spike in the coarse network. • The cross-correlogram P(τ ) is then normalized (divided) by the cross-correlogram P random (τ ) for the case if the original network was clustered randomly, and we integrate the ratio over interval of τ , on which the ratio exceeds 1, weighing by the inverse of the time lag τ to discount inputs received early, at large absolute value τ . We term this quantity integration coefficient ( IC): The integration coefficient then acts as a metric of the coarse network's similarity to a network of spiking nodes. The random clustering used as a baseline coarse-graining step is done by selecting a random partitioning of the nodes into the desired number of ensemble-nodes to match the average size of ensemble-nodes observed under the full-linkage clustering used in our coarse-graining. We note that this is distinct from randomly rewiring the network prior to any coarse-graining is applied in order to dissociate the network features responsible for the observed property.
Re-wired controls. As node degree sequence frequently drives emergent properties of a complex network 47 , we next examine the role of isolating the fine-grain network's node strength sequence in determining whether integrate-and-fire behavior emerges at a coarser scale. Node strength sequence is the equivalent of degree sequence for a weighted network, and is given by the sum of weights of all edges of a given node. To quantify the contribution of node strength sequence to our observations, we provide a comparison with an original network shuffled while preserving each node's strength. The shuffling follows a local rewiring rule: given 4 nodes (A, B, C and D), edges AB, CD, AC, and BD are shuffled from AB = CD = w 1 and AC = BD = w 2 to: AB = CD = w 2 and AC = BD = w 1 . This preserves the local strength of node A as: AC + AB = w 1 + w 2 before and after the rewiring, while dissolving the network's cluster-forming topology. The value of integration coefficient is lower for a shuffled network than for original brain network (Fig. 3B, right panel). It is also notable that the shuffled network, when re-scaled, is unable to produce the same average cluster size at the same connection threshold as original brain network, nor to maintain ensemble-spikes consisting of large numbers of spikes (Fig. 3B, right panel, bottom right corner where data is absent-area indicated in white).
Do clusters show neuron-like refractory periods? By means of an auto-correlogram P auto (τ ) (instead of binning ensemble-spike inputs from neighboring ensemble-nodes, we bin the ensemble-node's own most recent ensemble-spike preceding the ensemble-spike of interest), we search for the time-axis re-scaling parameters that result in the strongest degree of refractory behavior, which we measure as the ratio of ensemble-spike auto-correlogram at τ = −2 to that at τ = −1 : P auto (−2)/P auto (−1) . This helps us resolve whether a pair of consecutive ensemble-spikes of an ensemble-node had any systematic and noticeable time-gap in between, when ensemble-spiking was suppressed (Fig. 3C).

Further variations.
A particular coarse-grain version of a spiking network is realized by choosing multiple parameters: (1) c (edge-weight cut-off for full-linkage clustering), which directly influences N N -average number of nodes in an ensemble-node, (2) N S -number of spikes in an ensemble-spike and (3) N T time-axis re-scaling factor. The latter is formally an integer factor, by which the time bin used for the coarse level statistics is smaller than the mean inter-spike interval of a single node of the fine-grain network. Re-scaling the time-axis so that the number of spikes in an ensemble-node expected from the neuronal spiking frequency alone is ~ 1 (i.e. setting N T =N N ) permits singling out time bins with spiking activity elevated by a factor of N S as compared to the level expected from mean spiking frequency. These time-bins in the coarse ensemble-spike raster become the basis of the cross-correlogram and calculation of integration coefficient. Maximizing the integration coefficient with respect to N T after the combination of spatial ( N N ) and activity ( N S ) re-scaling parameters have been estimated, allows us to establish the full set of three parameters of re-scaling. This completes the analogy between spiking neurons and ensemble-spiking clusters of neurons. In search of the integrate-and-fire behavior for coarse network, we cover the space of spatial and activity re-scaling factors as follows: N N from minimum of ~ 3 neuron clusters up to ~ 12 in increments of 1; N S from minimal burst of two spikes to the maximum that can be sustained by the fine spiking network. The coarse-graining procedures described can also be iterated more than once, on already clustered network, by treating it as the fine-grain network and operating with the same machinery enabling a principled dynamics-preserving hierarchical reduction of functional models.
We note that our coarse-graining procedure only resolves a single step in what is potentially a multiple step process from neurons all the way to fMRI voxels. We consider the coarse-graining step demonstrated in the work as only one of the stages of potential renormalizations allowed in the brain with transitions between all pairs of consecutive scales preserving the dynamical principle. Therefore, we consider implanting the integrate-and-fire principle at the voxel level as a reasonable assumption, given the possibility of multiple scales with the principle emerging above the neuron level.

Data availability
The data and code that support the findings of this study are available from www.lcneu ro.org/tools . www.nature.com/scientificreports/