Investigating structural and functional aspects of the brain’s criticality in stroke

This paper addresses the question of the brain’s critical dynamics after an injury such as a stroke. It is hypothesized that the healthy brain operates near a phase transition (critical point), which provides optimal conditions for information transmission and responses to inputs. If structural damage could cause the critical point to disappear and thus make self-organized criticality unachievable, it would offer the theoretical explanation for the post-stroke impairment of brain function. In our contribution, however, we demonstrate using network models of the brain, that the dynamics remain critical even after a stroke. In cases where the average size of the second-largest cluster of active nodes, which is one of the commonly used indicators of criticality, shows an anomalous behavior, it results from the loss of integrity of the network, quantifiable within graph theory, and not from genuine non-critical dynamics. We propose a new simple model of an artificial stroke that explains this anomaly. The proposed interpretation of the results is confirmed by an analysis of real connectomes acquired from post-stroke patients and a control group. The results presented refer to neurobiological data; however, the conclusions reached apply to a broad class of complex systems that admit a critical state.


Introduction
The concept of complexity is used to characterize natural systems consisting of large numbers of nonlinearly interacting elements, resulting in the spontaneous collective behavior of the system on the macroscopic level, called emergence.However, complex systems reveal further intriguing properties; among others, one can mention scale invariance 1 , self-organized criticality [2][3][4] , and adaptability to new conditions 5 .The variety of complex characteristics makes it impossible to describe the systems by a reductionist approach, i.e., to derive system properties as a simple consequence of a physical law.Characterizing the system structure and dynamics requires rather a holistic approach relying on describing its properties on different levels of organization.In this respect, when the exact mathematical description is unattainable, agent-based modeling 6 is especially beneficial.Simulating the system as a collection of autonomous entities allows us to explore its dynamics and helps us provide its natural description, which includes emergent phenomena.This interdisciplinary approach has been applied to study complex systems, encompassing all scientific disciplines, such as physics, chemistry, biology, and social and economic systems 7 .
A canonical example of a complex system is the human brain, whose large numbers of neuronal cells display nontrivial multiscale organization 8,9 and complex characteristics. 10,11 t has also been discovered that power-law statistics, often used to describe critical phase transitions, are present in the brain.These power laws quantify the scale-free properties of neural avalanche distributions, determine the temporal organization of the brain signals recorded from various brain imaging techniques, and characterize the dependence of the correlation length with system size.The critical brain hypothesis 3 states that neural networks evolve towards and stay most of their time 12,13 in a state around a critical phase transition (we will also use interchangeably phrases 'at a critical point', 'in a critical state' or 'at criticality '), where the competition between order and disorder emerges.Systems in that state have been argued to exhibit optimal computational properties related to information processing, such as information transmission and storage [14][15][16] , computational power 17 and maximal sensitivity to stimuli 18,19 .
However, the status of criticality and the associated properties in the brain with neurological dysfunction are still not known precisely.One possibility is that given critical state provides optimal functioning of the healthy brain, neurological dysfunctions might be associated with its loss 20,21 , which opens this area of study to clinical applications 22 .These include the study of epilepsy 23,24 , Alzheimer's and Parkinson's disease 25,26 , and analysis of cognitive processes, including human learning 27 .The ideas of critical phenomena have recently been applied to the study of changes in brain dynamics due to brain damage, both purely computationally 28,29 and with realistic connectomes of stroke patients 30 .In the latter case, with a simple computational model, the authors were able to predict critical phenomena based on the first principles.Furthermore, the presence and severity of the stroke reportedly were related to a loss of critical behavior in the brain and a possible post-stroke recovery of a patient to the recovery of the critical state.However, the opposite hypothesis that the brain remains in the critical state even in case of serious injury is also possible.
The problem of assessing whether a system is at a critical point or whether it can exhibit a critical phase transition is even more delicate due to the subtleties of measures of criticality and proper interpretation of the results, and it has been challenging in similar contexts 31 .Calculation of additional criticality-aware quantities, beyond the usual second-largest cluster size, and a comparison with an artificial system with known criticality status reveals an inconsistency.As summarized in Fig. 1, the criticality status of the stroke-affected brain becomes unclear due to an ambiguous behavior of the second-largest cluster size.Explaining this crucial observation is our main motivation.
This study is of importance from several perspectives.Brain criticality is an appealing hypothesis implying optimal neural network processing, and a comprehensive understanding of its nature is crucial to understanding brain functioning.However, estimated characteristics of the critical state should be interpreted with particular care since the analysis of the complex systems of which the brain is undoubtedly an example often exhibit nontrivial and subtle properties.Therefore, properties associated with the criticality in the brain with structural damage are still a subject of vital discussion.Moreover, potential deviation from critical neural dynamics could open the possibility for clinical application, as indicated above.In this contribution, we propose a microscopic model of brain stroke to find out the possible mechanism underlying the observed inconsistency in measures associated with neural dynamics at criticality.Such an approach allows us to reproduce the statistics observed in empirical data while controlling the system's inner organization and being able to better characterize it using graph-theoretic tools.Finally, an explanation of the ambiguous character of the criticality indicators exposes the possible difficulty of utilizing a single measure of criticality and offers a consistent data-driven argument to monitor criticality in stroke patients.
The organization of the paper is as follows; in Section 2 we describe the data set used, define the Haimovici model and provide its analysis in the context of the post-stroke brain.Section 3 presents the novel model of artificial stroke with its comprehensive analysis.The results of previous sections are jointly discussed and interpreted in Section 4. Finally, Section 5 offers conclusions.Technical details of the graph connectivity measures and Haimovici and Ising models are presented in Section 6.

Criticality in a model of brain activity 2.1 Haimovici et al. brain model
Some crucial aspects of brain dynamics are well reproduced in the critical regime of the cellular automaton-type Haimovici et al. model 32 , where simple dynamical rules are applied to the network of cells based on empirical connectome.In this section, we demonstrate the quantitative characteristics of the model and discuss their relation to critical and non-critical states.
The Haimovici et al. model is a three-state cellular automaton 33 on a connectome encoded as a network with weighted connectivity matrix W .In the context of brain activity, each network node represents a region of interest (ROI) of the brain cortex.The model dynamics are discrete.At each time step, a node is in one of three states: inactive (I), active (A), or refractory (R).The transitions between the three states for the i-th node of the network are as follows: (i) I → A always if the sum of the weights of the active neighbors of the node is greater than the activation It is readily seen that although both ρ(1) and σ A exhibit local maxima typical for systems at a critical point, the shape of S 2 is less clear-cut and varies considerably between patients (shaded areas denote standard deviations).An explanation of this crucial observation constitutes the central aim of this study.(d-e) For the sake of presentation, observables are normalized to their maximal values.For the details of the model, the numerical simulations, and the calculation of observables, see Section 2 and Methods 6.3.

3/24
threshold parameter T , i.e., ∑ j active w i j > T , and with probability r 1 otherwise; (ii) A → R with probability 1; (iii) R → I with probability r 2 .Figure 1a shows an illustrative example of these transitions.Probabilities r 1 and r 2 are small numbers, r 1 , r 2 1, chosen before the start of a simulation, and determine the timescale of the system.Brain simulations based on the model are summarized in Fig. 1 and in Methods 6.3.
The model dynamics exhibit diverse behaviors depending on the choice of the connectivity matrix W .The authors of 32 used Hagmann et al.'s empirical connectome 34 to find dynamical phase transitions in a healthy brain.The use of small-world Watts-Strogatz (WS) topology with connection weights mimicking the ones found in empirical connectomes, investigated in 35,36 , showed that depending on the network parameters, the model may find itself in different regimes, see Fig. 1d, including transience to a ground state, continuous and discontinuous dynamical phase transitions (whereby we mean dis/continuous change of the order parameter, e.g., mean neural activity, as seen in Fig. 1c), which all have distinct properties known from physical systems (e.g., in discontinuous transitions hysteresis can be observed).

Brain criticality
The activation threshold T is the parameter used to control the dynamics of the system.With Hagmann et al.'s connectome 34 as the underlying network, the model admits a dynamical phase transition, and the critical value of the threshold is T c ≈ 0.073 (cf.Methods 6.3).For very small values of T , even the weakest connections between the nodes are enough to spread the activity (supercriticality).Conversely, for large values of T , active nodes may fail to activate their neighbors, and the total activity remains very low (subcriticality; see Fig. 1c).It is in the critical regime that the brain activity simulated by the model reproduces the correlations associated with the resting-state networks (RSNs) of the brain 32 .
The state-of-the-art analysis of model dynamics has been based on the time-averaged sizes of the largest clusters (cf.Methods 6.4).The order parameter can be identified with the average largest cluster size S 1 , and the critical point T c is located near the local maximum of the size of the second largest cluster S 2 as a function of T (if the maximum exists).Such behavior is a strong indicator of a dynamical phase transition 32,35 .The phase transition can also be revealed using other quantities 28 , e.g., the variability of the total activity σ A , the variance of the largest cluster size, the first coefficient of the autocorrelation function ρ(1) 37,38 , or the eigenvalues of the correlation matrix.
Figure 1d presents three quantities for non-critical (exhibiting only transient dynamics) and critical (exhibiting a continuous phase transition) systems obtained for the Haimovici dynamics on human-connectome-based WS networks.The quantities shown are the time-averaged size of the second-largest cluster S 2 , the first autocorrelation coefficient ρ(1), and the standard deviation of the total activity σ A .Evidently, each of them has a distinct functional form that allows distinguishing between the non-critical and the critical case.In Figure 1e, the same quantities are calculated for the Haimovici dynamics run on a set of empirical brain connectomes of stroke patients studied in 30 .In this case, not all signatures of a critical phase transition are equally clear: ρ(1) and σ A exhibit a form similar to the system that has a critical point, but S 2 behaves in a way that does not correspond clearly either to the non-critical or the critical WS case, in a similar manner to how it was reported in 28 .An extended investigation of several other quantities, described in Supplementary Information 1, revealed that the second-largest cluster size was the only observable with ambiguous outcomes.

Real and artificial strokes
The empirical data on strokes were recently studied 30 with the use of a dynamical model of brain activity very similar to the one presented in Sec.2.1 (cf.Methods 6.3 for details).The main finding was the observation that in model simulations on connectomes from patients three months after a stroke the maximum in the average size of the second-largest cluster of active nodes was missing, and it reappeared in a subgroup of those patients when connectomes were acquired again twelve months after the stroke.The loss and reappearance of the peak were interpreted, respectively, as a loss and recovery of the brain's ability to reach the critical state, parallel to the behavioral post-stroke recovery of a patient.
In order to explain the origin of these findings, we propose a minimal model of artificial strokes that recreates two key features found in the empirical data: a) the signatures of criticality summarized in Fig. 1e, and in particular the anomalous behavior (the absence of a peak) in the second-largest cluster size, and b) the decrease in connectome integrity correlated with this behavior.

4/24
A minimal model of artificial strokes We introduce a model of a stroke-like modification to a healthy connectome.Given a healthy empirical connectome, an artificial stroke changes the connectivity between a particular RSN with the rest of the brain.To this end, we randomly select a fixed fraction of nodes (a proxy of stroke severity) of the RSN and completely remove connections to their neighbors not belonging to the same RSN.This way, the internal structure of the RSN remains unchanged, while effectively decreasing the connection of the RSN with the rest of the brain.The model we propose aims to reproduce global characteristics found in real connectomes affected by a stroke and does not purport to offer a biologically or neurologically plausible mechanism.In Figure 2, we compare the second-largest cluster size calculated for the Haimovici dynamics on (a) connectomes with an artificial stroke of increasing severity located in a single RSN and (b) real post-stroke connectomes.In the latter case, the connectivity matrices were normalized (see Methods 6.3) thereby shifting the critical point T c → Tc .We report a good qualitative agreement between the outputs of the proposed artificial stroke model and the real stroke dataset, as in both models the second-largest cluster sizes deviate from the standard (i.e.healthy) case with a pronounced peak.

Structural analysis
We provide a deeper analysis of the dynamics of the Haimovici model and connectome integrity using graph-theoretic methods.In Figure 3, the dynamical part is summarized on the y-axis by the area under the S 2 (T ) curve I 2 = S 2 (T )dT , which qualitatively captures the loss of a peak, explained by high values of S 2 for low T ≤ T c .The connectome integrity is probed by the normalized modularity Q where the modules are found using the Louvain algorithm (cf.Methods 6.2 for details).In Figure 3a we present results for artificial strokes with six affected RSNs and different severity levels, i.e., disconnection of between 0% and 100% of the RSN nodes from the rest of the brain, while in Fig. 3b we plot real stroke results for three patient groups.The normalization is such that mildly affected dynamics are centered near the (0, 0) point, while dynamics severely affected by a stroke tend to move away from the origin.
We provide an additional analysis of the artificial stroke in terms of the conductance, h G (cf. Methods 6.2), an alternative measure of connectivity between network subsystems 39 .This measure, unlike modularity, quantifies the connectivity between a selected RSN and the rest of the connectome.It offers an additional test of the loss of integrity between known regions rather than any changes due to Louvain algorithm reconfiguring modules after the stroke.Figure 3c presents a plot of the normalized area I 2 versus normalized conductance h G .We reveal a strong correlation between the network integrity and the anomalous behavior of the second-largest cluster size.

Clusters in a divided Ising model
In the previous section, we established that the loss of connectome integrity coincides with the anomalous behavior of the second-largest cluster size in both artificial and real strokes.In light of this relation, we continue to investigate the changes in connectome structure to understand the development of this anomalous behavior in more detail and reintroduce the question of criticality: is the brain's critical state reachable, or is it lost after the stroke?In this part, we tackle it by combining the insights from connectome integrity with the Ising model, a paradigmatic case with an existing critical phase transition.Ising spins act as the neuron nodes, and the usual two-dimensional grid takes the role of the connectome.The stroke-induced loss of integrity is taken to an edge case where the connectome is completely divided into subsystems, thus modeling a severe artificial stroke.In particular, we recreate an anomalous lack-of-peak in the second-largest cluster size and describe the underlying mechanism as a competition within the hierarchy of subsystem-wide clusters. 6/24

Model description
The Ising model is among the simplest systems that undergo a continuous phase transition as the temperature T changes [40][41][42][43] .In a similar context, it has been used to study and reproduce the behavior of neuronal populations in cultured cortical neurons, cortical slices, and visual cortex among others [44][45][46] , and in particular to study structural damage to functional networks at criticality 29 .Another study compared the fMRI correlations with the correlation matrices obtained from Ising model simulated on empirical connectomes to characterize disorders of consciousness with the model's critical temperature. 47The Ising model itself consists of a network whose sites take the values ±1 (originally representing two states of atomic spins) together with a particular definition of interactions between adjacent spins, which describes the time evolution of spin states (cf.Methods 6.5).
The usual indicator of a continuous phase transition is the correlation length, i.e., the characteristic length-scale below which the non-adjacent spins are correlated, and which diverges at the transition point.In addition, it is agreed upon that observables related to clusters, used in percolation theory, are viable probes of critical state in the two-dimensional Ising model 48 .Clusters (domains) are defined as maximal connected sets of sites with the same orientation of the spins.The size of a cluster is the number of its sites.The size distribution of the clusters depends crucially on T and is used to construct two well-established indicators of criticality: (a) an abrupt change in the time-averaged size of the largest cluster S 1 and (b) a peak in the time-averaged size of the second largest cluster S 2 ; both occur near the critical temperature T c 49 .We introduce a key modification to the model's connectivity matrix by removing certain links from the original two-dimensional square lattice to form two completely disjoint subsystems, A and B. This mimics the loss of integrity, studied previously in the case of strokes, whereas the Ising dynamics fixes the critical behavior.We consider two cases: subsystems of (a) equal N A = N B , Fig. 4, and (b) unequal N A > N B number of nodes, Fig. 5.In what follows, we inspect the average cluster sizes S 1 and S 2 as indicators of criticality in the entire system (no superscript) and in each subsystem separately (superscripts A and B).Details of the simulations are discussed in Methods 6.5.shows the largest (S 1 ) and second-largest (S 2 ) clusters of the entire system.The graphs of S 1 and S 2 have both qualitatively the same shape, saturating for small T at half of the system size.The temperature dependence of both is similar to that of the largest cluster of the unmodified system (S 0 1 in panel (a)), and in S 2 no peak is present.Panel (c) shows the largest and second-largest clusters computed for each of the subsystems separately (superscripts A/B, green and orange lines).The curves are down-scaled versions of the clusters for the unmodified system (panel (a); note the different scale of the y-axis), with the signature peak in S 2 .Corresponding curves for the subsystems fully overlap.

Inconsistent indicators of criticality
In Figure 4a, we revisit the average sizes of the two largest clusters in the standard Ising model on a 100 × 100 square lattice.The size of the largest cluster (black dashed line) saturates at the size N of the entire system; the size of the second-largest cluster (black solid line) starts at zero for small T , increases to reach a peak near the critical temperature T c , and decreases to a certain nonzero value for large T (a finite-size effect; S 2 /N → 0 for N → ∞).

8/24
In the next step shown in Fig. 4b, we consider the divided Ising model with two fully disconnected subsystems of equal size and shape (rectangles 100 × 50).In this case, cluster sizes S 1 and S 2 both admit a similar functional form, as shown by the dashed and solid blue lines.Both quantities saturate similarly to the largest cluster calculated for the undivided system and have no maximum near T c .
Via an introduction of a simple division of the system into two non-interacting parts we recreate the anomalous lack of a peak in the second-largest cluster size S 2 .To gain insight into how this emerges, we plot in Fig. 4c the average sizes of the two largest clusters but restricted to each subsystem: S A 1 , S A 2 for subsystem A (green) and S B 1 , S B 2 for subsystem B (orange).Unsurprisingly, the respective curves for each subsystem overlap and their functional form is identical, up to a rescaling, as that of the cluster size curves for the undivided lattice, described above and shown in the inset of Fig. 4b.Since the two subsystems are fully disconnected, each has its own independent Ising dynamics: Both are critical around T c , and their respective cluster sizes show typical signatures of criticality.However, the lack of the characteristic peak in the size of the second-largest cluster of the entire system shown in Fig. 4b is clearly flawed by suggesting the absence of a critical phase transition in the divided system.
Cluster ordering Below we refine our understanding of the ordering of entire-system clusters in terms of subsystem clusters.As long as the subsystems are fully disconnected and thus have completely separate dynamics within each subsystem, it is the subsystem clusters that play a primary role, and the system-wide clusters only name the largest among them.The lack of a peak in the size of the second largest cluster observed at the level of the entire system comes merely from a particular ordering of the sizes of the subsystem clusters.
At each time step t, we are interested in the two largest clusters of the entire system S 1 (t), S 2 (t) and of each of the subsystems S A 1 (t), S A 2 (t), S B 1 (t), S B 2 (t).Note that in this paragraph we use momentary cluster sizes, not time averages.We can compare the sizes of the four above-mentioned subsystem clusters and write them down in a list in decreasing order.The first two entries in the list are the largest entire-system clusters: , where max i selects the i-th largest number.It is clear that regardless of the definitions of the two subsystems, S 1 (t) is either S A 1 (t) or S B 1 (t), so the functional form of the average S 1 corresponds to that of the largest subsystem cluster.Crucially, however, since the role of S 2 (t) can be assumed by any of the remaining clusters S A 1 (t), S A 2 (t), S B 1 (t), S B 2 (t), the functional form of the average S 2 depends on the relative size of the subsystem clusters and the robustness of their ordering (i.e., whether the order of the subsystem clusters on the list remains constant throughout the simulation).
If N A = N B , then {S 1 (t), S 2 (t)} = {S A 1 (t), S B 1 (t)} for most of the simulation time (where = denotes set equality).In other words, the roles of the two largest entire-system clusters are decided by the competition of the largest cluster of subsystem A and the largest cluster of subsystem B. Therefore, S 1 and S 2 both share the typical characteristics of the largest subsystem cluster, resulting in the same functional behavior, as shown by the blue lines in Fig. 4b.In the next section, we demonstrate the results of decreasing the robustness of the ordering by bringing S A 2 and S B 1 to comparable sizes.
Competition between clusters In the discussion above, we described how the momentary largest system clusters are selected from an ordered list of subsystem clusters.In this section, we apply this to the case of unequally sized subsystems N A > N B , where the ordering is less stable in the simulation time.We use a square lattice of the same dimensions as previously but with subsystem B redefined as a smaller square patch in the middle of the lattice and with the connections changed accordingly, as shown in Fig. 8.This reduces the average sizes of the subsystem-B clusters S B 1 and S B 2 .For a certain patch size, the typical size of the largest cluster in subsystem B S B 1 is comparable to the second largest cluster in subsystem A S A 2 , and these two clusters compete for the role of the second largest cluster in the entire system S 2 .8 for the geometric arrangement.(a) Temperature dependence of the largest and second-largest cluster sizes for the entire system (blue lines), subsystem A (orange lines), and subsystem B (green lines).The line styles follow Fig. 4. In subsystems, characteristic features such as saturation of S 1 at small temperatures and a peak in S 2 near the critical T c persist, but the size symmetry breaking with N A > N B results in relative growth of clusters related to subsystem A with respect to subsystem B. The vertical blue line is a near-critical region where we probe the time evolution to be shown in the right panel.(b) Interplay between momentary clusters S B 1 (t), S A 2 (t) and S 2 (t) over a number of time steps at a fixed T (the line styles are the same as for the averages in panel (a)).One can observe the competition for the place of the second largest cluster (blue solid line), with S A 2 (t) (orange solid line) dominating earlier and S B 1 (t) (green dashed line) later.

9/24
In Figure 5, this mechanism is demonstrated on a square lattice of N = 10000 sites divided into N A = 7500 and N B = 2500 sites.The loss of symmetry between the subsystem sizes results in a different temperature dependence of the cluster sizes.In the case of equally sized subsystems shown in Fig. 4b and described in previous sections, the average cluster sizes in subsystems A and B overlap: S A 1 ≈ S B 1 (dashed green and orange lines) and S A 2 ≈ S B 2 (solid green and orange lines).On the other hand, the unequally sized case discussed currently is presented in Fig. 5a and shows a clear cluster size asymmetry: S A 1 S B 1 and S A 2 S B 2 near the critical temperature indicated by the blue vertical line.In this region, moreover, S A 2 ≈ S B 1 , which signals the breaking of the cluster ordering.The competition between the clusters S A 2 and S B 1 is shown more closely in Fig. 5b.In this panel, we plot the dynamics of the clusters of interest at a fixed temperature T ≈ T c .The solid orange line and the dashed green line trace the temporal evolution of the size of the clusters that compete for the role of the entire-system cluster S 2 .In this particular time window, we observe a reversal of the roles, in which the role of S 2 is played by S A 2 in earlier time steps and by S B 1 in later time steps.The plot reveals rich intermittent dynamics that elucidate the averaged picture in Fig. 5a.The cluster order is perturbed even beyond the competition of S A 2 and S B 1 , since, at certain moments such as t = 31, the cluster S B 1 is temporarily able to take over the role of S 1 and reduce the cluster S A 1 to the role of S 2 .

Loss of peak: A mechanism common to the models of Ising and Haimovici
In Fig. 6, we juxtapose cluster analyses performed on (a) the divided Ising model of unequal sizes and (b) the Haimovici model with severe artificial stroke resulting in completely disconnected auditory RSN (different choices of base RSNs give similar results; see Supplementary Fig. 11 in Supplementary Information 7.2).For a detailed description of how the Ising model and the Haimovici brain model are different and how they allow a direct comparison, cf.Methods 6.6.In both models, the average second largest cluster size S 2 of the entire system changes its behavior in the same way: When subsystem B -or the auditory RSN -is completely disconnected, S 2 loses its characteristic critical peak.The numerical data for both models fully support our prediction: At some value of the threshold parameter T , the size of the largest cluster of the smaller subsystem S B 1 exceeds the size of the second largest cluster of the larger subsystem S A 2 .This change in the cluster size hierarchy leads to S 2 monotonically decreasing as a function of T .In both cases, the second largest cluster in subsystem A (orange solid line) exhibits a maximum around the critical point, however, just below this point, near T ≈ 0.05 in the Haimovici model (T ≈ 2.2 in the Ising case) the size of the largest cluster in subsystem B (dashed green line) becomes larger and takes over in the cluster size hierarchy.
In Figure 6, we presented the edge case of a fully disconnected subsystem, or a severe artificial stroke.In Figure 7, we extend the presentation to a comparison of gradual system modifications in (a) the Ising model with various sizes of the disconnected subsystem B and (b) the Haimovici model for varying degrees of connectivity between the auditory RSN, chosen as the smaller subsystem, and the rest of the network.Both the variation in the size of the smaller Ising subsystem and changes in connectivity between RSNs produce similar comb-like families of lines (cf.Fig. 2 for similar plots comparing connectomes affected by artificial and real strokes).These modifications are not strictly equivalent, but they provide a clear presentation of the idea.In Supplementary Information 7.2, we expand on the present analysis with complementary modifications.The small perturbation regime, where the subsystem has relatively few nodes (in the Ising model) or few interconnections are removed (in the Haimovici model), constitutes the lower part of the comb.In this region, the entire system cluster S 2 has typical critical behavior with a peak near the critical T c or T c .The large perturbation regime is, in turn, reached when the subsystem is large (the Ising model) or almost completely disconnected (the Haimovici model), and constitutes the upper part of the comb.The upper limiting cases are reached if subsystem B occupies about half of the system or if the auditory RSN is fully disconnected, as studied in Figs.4b and 6, respectively.In this regime, the cluster S 2 loses its characteristic criticality peak.

Discussion
By considering the divided Ising model, we explained the rise of anomalous lack of peak in the size of the secondlargest cluster of active nodes, S 2 , in terms of subsystem clusters forming an ordering for the entire system.We believe this observation lies at the source of the confusion about criticality illustrated in Fig. 1.It should be emphasized that each subsystem separately preserves the typical behavior of clusters related to criticality, i.e., growth of the largest cluster, S 1 , and a maximum in the second largest cluster size, S 2 .The critical phase transition is preserved, although it becomes obfuscated by the loss of connectome integrity.Let us note, that there might be another contributory factor: a shift in T c connected to the change in the number of nodes in the network, as observed in a related model 50 .This effect, however, does not account for disparate behaviors of the criticality indicators, nor is it strong enough to explain a full loss of phase transition.
We were able to consistently translate the above explanation to the Haimovici model.There, an apparent loss of criticality was a consequence of changes in the ordering of subsystem clusters enabled by the changes in connectome integrity.In both models gradual subsystem changes led to the comb-like behavior of S 2 , Fig. 7, seen previously in the real stroke data in Fig. 2b.Due to the particularities of the Ising dynamics, instead of varying the connectivity, which would render the comb a single line, we varied the subsystem size, which reproduces the expected effect.We discuss these complementary approaches in Supplementary Information 7.2.A similar subsystem-cluster analysis was performed with each of the eight RSNs entirely disconnected, see Supplementary Information, resulting in a similar effect.
Crucially, the change in connectome integrity was the primary driver while the dynamics itself was secondary.Despite the limiting assumptions, the artificial stroke model recreated quite well the overall loss of connectome integrity in Figure 3.In both artificial and real strokes, we found that a stroke-induced loss of connectome integrity measured by an increase of the modularity Q and the decrease of conductance h G coincides with the anomalous behavior of the second-largest cluster size.
While the authors of another study simulating Ising dynamics on connectomes 47 hypothesize that structural connectivity alone cannot explain some effects observed in functional correlations of patients after severe brain 12/24 injuries, no other explanation has been provided in the literature.Our explanation of the anomalous behavior of the cluster-based indicator of criticality, on the other hand, is consistent with the distinction 51,52 between well-known structurally-driven percolation phase transitions and dynamic transitions-for which less is known. 36Reframed in these terms, our results would suggest that in brain strokes there is no change in whatever underlies dynamic transitions, but that there are structural changes that could disorganize the percolation-like transitions.

Summary and conclusions
In this study, we revisited the question of whether post-stroke brain dynamics stay at the critical point.To this end, we compared how indicators of criticality behave after real strokes and computer-simulated strokes proposed for this purpose, and we do not find evidence for the post-stroke loss of critical dynamics as previously suggested.Rather, we show that elementary indicators of criticality should be interpreted with caution.In particular, the behavior of the size of the second-largest cluster of activity in the brain affected by stroke may result solely from the loss of connectome integrity without the brain's departure from a dynamic critical transition.From this perspective, the behavior is understood in terms of subsystem clusters competing for the top rank system-wide.The results have been reproduced in classic physical models and confirmed based on graph-theoretical characteristics calculated for the empirical connectomes.Thus, when a system with an unknown subsystem structure is analyzed, stroke-induced loss of critical dynamics may be illusory, and the simpler and more plausible explanation is the described mechanism cluster competition mechanism enabled by a loss of connectome integrity alone.

Outlook
We hope that as the concept of criticality becomes increasingly relevant 22 , and the critical state of the brain is evaluated in relation to various diseases, disorders, states of consciousness, and tasks, this work offers an important consideration towards the robustness of these findings.
In addition to applications to neuroscience, our discussion is relevant in studies of artificial neural networks 53,54 where, for the specific learning dynamics of the connectivity matrix, the network self-organizes towards a critical state.In such models, typical criticality indicators are based on avalanche sizes, which are structure-agnostic like the cluster sizes considered in this work.

Connectivity matrices Connectomes of post-stroke brains
The connectivity matrices that we used in Haimovici model in Sec.2.2 were a set of 113 individual connectomes from 79 stroke patients and 47 connectomes from 28 control subjects acquired via Diffusion-Weighted Magnetic Resonance Imaging (DWI) 30,55,56 .The patients' connectomes were acquired twice, 3 months (t 1 ) and 12 months (t 2 ) after the stroke.Each connectome encodes a network of N = 324 nodes of cortical ROIs, based on Gordon's parcellation 57 .By convention, the diagonal elements of W are set to zero.

Connectomes in artificial strokes
For the minimal model of artificial strokes the Hagmann et al.'s connectome 34 served as the base connectivity matrix representing the healthy brain.

Watts-Strogatz networks
For the Haimovici model results presented in Fig. 1d, the empirical connectomes were substituted with Watts-Strogatz small-world networks 58 of comparable size (N = 2000 nodes).These networks are constructed from a ring of nodes, each symmetrically connected to its nearest neighbors with k edges whose ends are subsequently rewired to random nodes with probability π.Following 35,36 , to mimic the weight distribution of the human connectome 34 , the link weights were sampled from an exponential distribution p(w) = λ e −λ w , with λ = 12.5.The same Greenberg-Hastings 33 dynamics was used as for the Haimovici model, with t max = 10000, t init = 200, r 1 = 0.001, r 2 = 0.3, and the unnormalized symmetric weight matrices W . Node degrees and rewiring probability in WS networks were set to k = 2, π = 0.5 to obtain non-critical behavior and to k = 10, π = 0.5 to obtain the critical one. 13/24

Graph connectivity measures
First, we set the notation and some basic definitions for weighted directed graphs, such as connectomes, that we adopt throughout the paper.
For a weighted directed graph G with N nodes, we denote its binary adjacency matrix as A and a i j as its element corresponding to the directed connection from node j to i, and analogously the weighted adjacency matrix as W and w i j as the weight of the directed connection.The analog of degree in weighted directed graphs are in-degree and out-degree strength: w in i = ∑ j w i j , w out i = ∑ j w ji .The total strength of the graph is 2w = ∑ i w in i = ∑ i w out i = ∑ i ∑ j w i j , and the average in-degree strength is w in = ∑ i w in i /N.In a graph G with the set of nodes V , conductance of a node subset S ⊂ V and its complement S = V \ S is the quotient 39 h G (S, S) = |cut(S, S)| min(vol(S), vol( S)) of the weighted size |cut(S, S)| = ∑ i∈S, j∈ S w i j + w ji of the cut, i.e., the set of all edges connecting S and S, and the smaller of the total strengths vol(S) = ∑ i∈S w out i summed over all nodes belonging to these sets.In our study, the connectome was divided into the nodes of the chosen RSN and the rest of the network.
Graph modularity Q can be computed for weighted graphs both undirected (Hagmann et al.'s connectome-based networks) and directed graphs (stroke dataset connectomes) 59 using the weighted adjacency matrix w i j : where c i is the label of the module to which node i belongs, and δ c i ,c j is the Kronecker delta symbol of two such labels.In our work, the modularity was always normalized to the maximal modularity of a perfectly mixed network, The results presented for artificially modified connectomes are averages of 20 modification realizations for each value of parameters used.We used the implementation of conductance and of Louvain modularity optimization algorithm from NetworkX Python package 60 .
Normalization in Figure 3 To compare artificial and real strokes in Fig. 3, we propose a normalization of quantities defined as where A 0 denotes the value of the quantity in the unmodified system (i.e., a system not affected by stroke).Systems whose normalized parameters are close to zero are very similar to the unmodified system, whereas large values of normalized parameters indicate a large deviation from it.For real strokes, as the proxy of the unmodified system, we use the control group averages.

Haimovici model
For simulations of the Haimovici model, we use a publicly available Python code based on the Susceptible-Excited-Refractory model 61 .The simulation runs in time steps, each with the following three transitions performed on all nodes concurrently 33 : 1. Active → Refractory; All nodes activated in the previous time step become dormant.
2. Refractory → Inactive; Each dormant node may become inactive with probability r 2 . 14/24 3. Inactive → Active; Inactive nodes become active in two ways: a) due to spontaneous activation with probability r 1 and b) due to neighboring active nodes that satisfy the threshold criterion ∑ j active w i j > T .
Numerical simulations are performed for t max = 10000 time steps for 30 threshold values T , covering the critical value T c .Each simulation starts in a random state with 1% of nodes active and the rest inactive.To account for this initial randomness, the first t init = 200 time steps in all simulations are discarded from further analysis, leaving t sim = t max − t init simulation steps.The result, for each T , is a N × t sim data matrix whose entries take the values 0 (inactive nodes), 1 (active nodes) and 2 (refractory nodes).Following the literature 32,62 , after finishing the simulation, we apply a preprocessing step of conflating the inactive and refractory states by substituting 2 → 0 in the data matrix.
The states of nodes s i (t) at time t are used to calculate the total activity A(t) = ∑ N i=1 s i (t).The average activity A and the deviation of the total activity σ A are defined as follows: An approximate critical value of the threshold in the Haimovici model can be computed using the mean-field approach 62 : where r 2 is the probability of transition from the refractory state to the inactive state.In all our calculations, following previous studies 62 , we keep r 2 = r 1 0.2 , where r 1 is the probability of spontaneous activation set to r 1 = 2/N.For Hagmann et al.'s connectome, the mean-field critical threshold value is T c = 0.08.This theoretical value agrees quite well with the results of our numerical simulations.
Modifications to the Haimovici model As a way of implementing the homeostatic plasticity principle in the Haimovici model, the network excitability was balanced by normalizing the incoming node's excitatory input in the structural connectivity matrix 62 : wi j = w i j /w in i .Such equalization minimizes the variability of activity and of the position of the critical point between individual connectivity matrices.It is noteworthy that the control parameter also becomes rescaled T → T , so that the critical point moves to an empirical value Tc ≈ 0.15.The model exhibits a similar behavior of the cluster sizes as the Haimovici model, where at the critical value of the threshold the size of the second largest cluster peaks.Similarly to the Haimovici model, in numerical simulations, we use r 1 = 2/N and r 2 = r 1 0.2 and simulate the activity of the system for t max = 10000 time steps.
The normalized connectivity matrices W from the stroke data set used in 30 are available publicly 63 .

Clusters
We define clusters as maximal sets of nodes sharing the same type of activity (±1 in the Ising model and active in the Haimovici model) connected following the adjacency matrix.The size of the cluster is the total number of participating nodes.Among all cluster sizes, we focus on the two largest cluster sizes S 1 and S 2 , which are standard order parameters in percolation theory 64 .We measured the sizes averaged over the simulation time t sim : where S i (t) is the size of the momentary i-th largest cluster found at time t. 15/24

Ising model
The Ising model on a lattice with the adjacency matrix A is defined in terms of the energy function: where the spin variables s i take values ±1 and J is the coupling parameter.The probability of a particular spin configuration s = (s 1 , s 2 , . . ., s N ) depends on temperature T according to the Boltzmann distribution and the approximate value of the critical temperature is T c = 2.27, as shown in a classic paper by Kramers and Wannier 65 .
Details on numerical simulations The coupling parameter was set to unity J = 1 and a square lattice of dimensions 100 × 100 with nonperiodic boundary conditions was used.The size of the system N = 10000 was the total number of lattice sites.In the case N A = N B , the system was divided into two rectangles 100 × 50.In the case N A > N B , subsystem B was defined as a square patch in the center of the lattice and subsystem A as its complement (see Fig. 8), and modified the connections accordingly.We performed numerical simulations using Monte Carlo importance sampling 66 .Spins in the initial conditions were set to −1 with probability 0.75 and to 1 otherwise.Each time step (or sweep) comprised N spin-flips, which in turn followed the Metropolis approach 67 .First, the i-th node was flipped s i → −s i and the resulting change in energy was computed ∆E = E(spin-flip) − E 0 .The new configuration was accepted with probability min(e − 1 T ∆E , 1).Simulations were run for 30 equally spaced temperature values between T = 0.01 and T = 4.5, each simulation continuing for t max = 5000 sweeps with the t init = 200 initial sweeps discarded.The resulting data matrix for a single temperature had dimensions N × (t max − t init ).

Differences between Ising and Haimovici models
Here, we describe differences between the Ising and the Haimovici models which are important for understanding the comparison between the models but are largely immaterial to our argument.The Haimovici model has starkly different dynamics than the Ising model.In the latter case, each spin, whether pointing up or down, always belongs to some cluster, unless it is surrounded by four spins of the opposite orientation.This results in large domains such that, in the low T regime, the largest and second-largest clusters can cover virtually the entire subsystems A and B, respectively (see Fig. 4b).In contrast, in the Haimovici model, it is only the clusters of active nodes that are considered.An activated node becomes refractory at the very next time step and then waits to become inactive again (with the probability of becoming inactive at a given step r 2 ≈ 0.29).This burst-like cluster formation results in the average sizes of the largest and second-largest clusters attaining much lower maximal values (S 1 ≈ 0.17N for T near 0; S 2 ≈ 0.006N at its peak) than in the case of the Ising model and rapidly decreasing to zero for increasing 16/24 values of T > T c .The unmodified networks in both models are also significantly different: the Haimovici model uses Hagmann et al.'s connectome, which belongs to the small-world class 35 , whereas the Ising model uses a regular two-dimensional square lattice.
When comparing both models, one should be mindful that the role of the threshold parameter T is inverse to that of the temperature T in the Ising model, i.e., small values of T result in supercritical states with high total activity, and large values of T result in subcritical states with low stochastically induced activity.

Investigation of various measures of criticality
We calculate an extended list of criticality measures for the empirical connectomes of stroke patients and for the critical and non-critical human-connectome-based Watts-Strogatz networks.This Supplementary Information serves as an extension to Fig. 1 where the three most relevant quantities were shown.The results of this extended study are reported in Supplementary Fig. 9.

Other modifications of the models
In this section, we apply different modifications to the Ising and Haimovici models to demonstrate the robustness of the main result presented in Fig. 7.We swap the modifications in a complementary way.We gradually disconnect subsystem B in the Ising model (in the main text, we varied the subsystem size), and we vary the subsystem size in the Haimovici model by disconnecting different single RSNs (in the main text, we gradually disconnected only the auditory RSN).Both results are shown in Supplementary Fig. 10.In this study, the Ising model runs on a 32 × 32 lattice and subsystem B is a patch of size 16 × 16.
In the Ising case, Supplementary Fig. 10a, we find a family of curves forming a degenerate comb with its upper non-critical branch consisting of a single line corresponding to all connections removed.The remaining modifications cover the lower branch of the comb and exhibit indicators of criticality.This behavior is characteristic of the model in which even a single connection between the subsystems is sufficient to connect the clusters and make the behavior qualitatively the same as in an unmodified system.
Disconnecting different RSNs serves as a proxy for varying the size of the subsystem.In this case, the resulting family of curves does form a comb-like pattern, albeit somewhat perturbed; this is not a surprise, as different brain substructures are not rescaled copies of each other.Still, the curves are size-ordered to a certain degree.The most prominent examples of partitions that preserve the critical peak (such as DorL and VisL) are the smallest RSNs, and those that appear non-critical (such as DMN) are the largest ones.This is in agreement with the Ising case with varying subsystem size, as seen in Fig. 7a.The competition between S A 2 and S B 1 and the resulting second largest cluster of the entire system when individual RSNs are fully disconnected is additionally depicted in Supplementary Fig. 11.Each panel shows the second largest cluster of the entire system, S 2 (blue line), the second largest cluster of the larger subsystem, S A 2 (orange line), and the largest cluster of the disconnected RSN (abbreviation and size give in the panel), S B 2 (dashed green line).In all the cases, at some value of the threshold parameter T , the order of sizes changes, impacting the size of the second largest entire-system cluster.Depending on the case, the peak in S 2 may persist (e.g., DorL, VisL) or be absent (e.g., DorM, Aud, and DMN).

S 2 Figure 1 .
Figure 1.Haimovici et al. brain model: description and dependence of criticality status on connectome topology.(a) Illustration of model dynamics: transitions between possible states of the network nodes (circles with numbers) with connection weights w i j .Initially, the inactive (green) central node is connected to three active (orange) nodes, and to one node in a refractory state (blue).In the next step, assuming that the sum of active neighbors' weights is larger than the threshold, i.e. w 15 + w 35 + w 45 > T , the central node is activated, while the active nodes become refractory.In the last step, the transition from the refractory to the inactive state takes place randomly at each node.(b) Example of the structural part of the model: Hagmann et al.'s connectome of healthy human subjects, represented by an adjacency matrix.The Haimovici model consists of panel (a) dynamics applied to the connectome-based network.(c) Model criticality for a healthy connectome: total activity and a time-averaged size of the largest (S 1 ) and second-largest (S 2 ) clusters of concurrently active connected nodes, for varying threshold parameter T .For small values of the threshold T T c (red), the nodes are easily activated, while for large values T T c (gray) the nodes remain mostly quiet.Near the critical threshold value T ≈ T c (purple), the activity becomes correlated, resulting in a characteristic peak in S 2 .(d-e) Beyond the Haimovici model: panels show an extension of the model resulting from a replacement of the healthy connectome by artificial networks and post-stroke connectomes.The simulated activity is used to compute various indicators of criticality: the time-averaged size of the second-largest cluster S 2 , the first autocorrelation coefficient ρ(1), and the standard deviation of the total activity σ A .(d) Known criticality status: examples of Watts-Strogatz small-world networks resulting in clearly discernible non-critical/critical dynamics (upper/lower part).(e) Stroke patients connectomes: ambiguous criticality status.It is readily seen that although both ρ(1) and σ A exhibit local maxima typical for systems at a critical point, the shape of S 2 is less clear-cut and varies considerably between patients (shaded areas denote standard deviations).An explanation of this crucial observation constitutes the central aim of this study.(d-e) For the sake of presentation, observables are normalized to their maximal values.For the details of the model, the numerical simulations, and the calculation of observables, see Section 2 and Methods 6.3.

Figure 2 .
Figure 2. Size of the second largest cluster S 2 as a function of the threshold parameter T for artificial and real stroke.(a) Results for an artificial stroke afflicting the auditory RSN with varying stroke severity.The color lines correspond to the fraction of nodes in the RSN disconnected from the rest of the brain, ranging from 75% (the bottom-most dotted curve) to all (the topmost dashed curve).(b) Results for connectomes from the stroke dataset in both control and patient groups (cf.Fig 3b in Rocha et al. 30 ).Each line corresponds to one person.(a-b) Artificial strokes successfully recreate an anomalous loss of peak in the S 2 curve.

Figure 3 .
Figure 3. Structural analysis of stroke-affected connectomes.(a) Normalized area under the S 2 plot versus normalized modularity for the artificial strokes on Hagmann et al.'s connectomes with single RSNs gradually more disconnected from the rest of the brain, starting from 0% to 100% of nodes (from zero to high normalized modularity).The Pearson's correlation coefficient between the normalized area under S 2 and the normalized modularity is ρ = 0.88, p-value 0.001 with 95%CI[0.84,0.92].(b) Similar to (a) but for empirical connectomes from the real stroke dataset for three patient groups with correlation coefficient ρ = 0.646, p-value 0.001, 95%CI[0.55,0.73].(a-b) Both types of strokes with variable severity result in an increase in the area under the S 2 curve as well as the values of modularity.The latter is related to the loss of subsystem interconnectivity of the stroke-affected connectome.(c) Normalized area under the S 2 plot for the Haimovici model versus the normalized conductance for Hagmann et al.'s connectomes with an increasing fraction of RSN's nodes disconnected from the rest of the brain, from 0% to 100% (from low −h G (norm.) values to −h G (norm.) = 1).The correlation is ρ = 0.88 with p-value 0.001, 95%CI[0.84,0.92].

Figure 4 .
Figure 4. Divided Ising model: temperature dependence of two largest cluster sizes.A schematic representation of a square lattice and dependence of average sizes of the largest cluster S 1 (dashed lines) and of the second-largest cluster S 2 (solid lines) on temperature T .(a) In an undivided lattice (upper index 0), clusters exhibit a characteristic saturation of S 0 1 for small temperature and a maximum of S 0 2 in the vicinity of the critical temperature T c ≈ 2.27.Panels (b-c) show a system divided into two equally sized subsystems (as marked by the yellow dashed line on the lattices).Panel (b)shows the largest (S 1 ) and second-largest (S 2 ) clusters of the entire system.The graphs of S 1 and S 2 have both qualitatively the same shape, saturating for small T at half of the system size.The temperature dependence of both is similar to that of the largest cluster of the unmodified system (S 0 1 in panel (a)), and in S 2 no peak is present.Panel (c) shows the largest and second-largest clusters computed for each of the subsystems separately (superscripts A/B, green and orange lines).The curves are down-scaled versions of the clusters for the unmodified system (panel (a); note the different scale of the y-axis), with the signature peak in S 2 .Corresponding curves for the subsystems fully overlap.

Figure 5 .
Figure 5. Dependence of cluster sizes on temperature and time in the two-dimensional Ising model divided into two fully disconnected subsystems of unequal size.Subsystem A consists of N A = 7500 sites and subsystem B of N B = 2500 sites; see Fig.8for the geometric arrangement.(a) Temperature dependence of the largest and second-largest cluster sizes for the entire system (blue lines), subsystem A (orange lines), and subsystem B (green lines).The line styles follow Fig.4.In subsystems, characteristic features such as saturation of S 1 at small temperatures and a peak in S 2 near the critical T c persist, but the size symmetry breaking with N A > N B results in relative growth of clusters related to subsystem A with respect to subsystem B. The vertical blue line is a near-critical region where we probe the time evolution to be shown in the right panel.(b) Interplay between momentary clusters S B 1 (t), S A 2 (t) and S 2 (t) over a number of time steps at a fixed T (the line styles are the same as for the averages in panel (a)).One can observe the competition for the place of the second largest cluster (blue solid line), with S A 2 (t) (orange solid line) dominating earlier and S B 1 (t) (green dashed line) later.

Figure 6 .
Figure 6.Cluster sizes in fully disconnected subsystems in the Ising and Haimovici models.(a) The Ising model with subsystems of unequal size N A = 7500, N B = 2500.(This example was shown in Fig. 5.) The legend is valid for both plots.(b) The Haimovici model with subsystems of sizes N A = 879 and N B = 119, where B denotes the auditory RSN, and A is the rest of the network.(a-b) In the unmodified systems (insets), S 2 shows a characteristic maximum near the critical point.With subsystem B disconnected, the second largest cluster saturates (panel (a) Ising) or grows monotonically (panel (b) Haimovici) (blue solid line) as we lower T .In both cases, the second largest cluster in subsystem A (orange solid line) exhibits a maximum around the critical point, however, just below this point, near T ≈ 0.05 in the Haimovici model (T ≈ 2.2 in the Ising case) the size of the largest cluster in subsystem B (dashed green line) becomes larger and takes over in the cluster size hierarchy.

Figure 7 .
Figure 7.The size of the second largest cluster in gradually modified systems.(a) The Ising model with varying subsystem size.The family of color lines represents S 2 for different sizes of subsystem B ranging from N B = 400 (the bottom dotted dark blue line) to N B = 3600 (the topmost dashed yellow line).Above a certain size of subsystem B, the critical maximum is lost.The red arrow marks the critical temperature.(b) The Haimovici model with changing RSN connectivity.The color lines represent S 2 for connectomes with a varying number of connections between the auditory RSN and the rest of the brain removed, ranging from 75% of the connections removed (the bottom-most dotted light blue curve), to all connections between the RSN and the rest of the brain removed (the topmost dashed yellow curve).The red arrow indicates the critical value of the threshold.(a-b) In both models, above a certain proportion of removed connections, the maximum in S 2 is lost.

Figure 8 .
Figure 8. Unequal Ising subsystems.Division of the Ising model lattice into unequal sizes N A > N B .

Figure 9 .Figure 10 . 1 Figure 11 .
Figure 9. Criticality measures in the Haimovici model.Left and middle: noncritical and critical connectome-based Watts-Strogatz networks.Right: stroke dataset connectomes.The measures include: S 2 , the size of the second largest cluster of activity; ρ(1), the first coefficient of the autocorrelation function; σ A , the standard deviation of the total activity; λ 2 , the second largest eigenvalue of the node activity cross-correlation matrix; CC , the average of cross-correlation matrix elements.