Autopoietic influence hierarchies in pancreatic β-cells
Dean Korošak, Marko Jusup, Boris Podobnik, Andraž Stožer, Jurij Dolenšek, Petter Holme, Marjan Slak Rupnik
AAutopoietic influence hierarchies in pancreatic β cells Dean Koroˇsak,
1, 2
Marko Jusup, Boris Podobnik,
4, 5, 6, 7, 8
Andraˇz Stoˇzer, Jurij Dolenˇsek,
1, 9
Petter Holme, ∗ and Marjan Slak Rupnik
1, 10, 11, † Institute of Physiology, Faculty of Medicine, University of Maribor, 2000 Maribor, Slovenia Faculty of Civil Engineering, Transportation Engineering and Architecture, University of Maribor, 2000 Maribor, Slovenia Tokyo Tech World Research Hub Initiative (WRHI), Institute of Innovative Research,Tokyo Institute of Technology, Tokyo 152-8552, Japan Faculty of Civil Engineering, University of Rijeka, 51000 Rijeka, Croatia Center for Polymer Studies, Boston University, Boston MA 02215, USA Zagreb School of Economics and Management, 10000 Zagreb, Croatia Luxembourg School of Business, 2453 Luxembourg, Luxembourg Faculty of Information Studies in Novo mesto, 8000 Novo Mesto, Slovenia Faculty of Natural Sciences and Mathematics, University of Maribor, 2000 Maribor, Slovenia Center for Physiology and Pharmacology, Medical University of Vienna, 1090 Vienna, Austria Alma Mater Europaea – European Center Maribor, 2000 Maribor, Slovenia (Dated: January 8, 2021) β cells are biologically essential for humans and other vertebrates. Because their functionalityarises from cell-cell interactions, they are also a model system for collective organization among cells.There are currently two contradictory pictures of this organization: the hub-cell idea pointing atleaders who coordinate the others, and the electrophysiological theory describing all cells as equal.We use new data and computational modeling to reconcile these pictures. We find via a networkrepresentation of interacting β cells that leaders emerge naturally (confirming the hub-cell idea),yet all cells can take the hub role following a perturbation (in line with electrophysiology). The importance of the nutrient-sensing and insulin-secreting β cells in vertebrates is hard to overstate. Thesecells reside in pancreatic islets, where they extensivelycommunicate with each other and their environment [1].The intercellular communication, in particular, serves tocoordinate and synchronize cellular operations throughwhich insulin is released in proportion to stimulation andmetabolic requirements [2, 3]. The delicate nature ofthis task is reflected in the continual need to preventoversecretion, and subsequent hypoglycemia [4], despiteintracellular stores holding sufficient insulin to exceed thelethal dose by orders of magnitude if released at once.On a microscopic scale, β cells are electrically andmetabolically coupled via gap junctions built from theprotein connexin-36 (Cx36) [5]. The key role of such cou-pling is seen in the fact that the loss of Cx36 channels hasa detrimental impact on β -cell cooperation [6, 7], lead-ing to uncoordinated plasma depolarizations, the desyn-chronization of calcium signals, and increases in basalinsulin release. Sufficient Cx36 coupling, by contrast,curbs the β -cell intrinsic heterogeneity in glucose sensi-tivity [8], and in mouse islets, limits the threshold forstimulatory glucose concentration [7] to a narrow bandaround ≈ β cells are likeindividual soldiers who fall in line when the communica-tion channels between them are open. A typical β cellis coupled to between six and eight neighbors [9, 10].According to electrophysiological theory [11], single cellslack mechanisms to become pacemakers beyond their im-mediate neighborhood. In terms of our metaphor, thereare no generals in the army.The above-mentioned lateral organization is diametri- cally opposed to the picture by functional studies of alarge number of communicating β cells (referred to as β -cell collectives). For example, intrinsic cellular hetero-geneity in glucose sensing [8], heterogeneous gap-junctioncoupling [12], and extensive paracrine signaling [13] allcontribute to an islet-wide complicated cytoplasmic Ca dynamics. This dynamics of the Ca concentration inthe cytosol [14, 15] is a key insulin-secretion driver [16]. Ithas recently become possible to record the Ca dynam-ics with a great spatial and temporal precision using thefunctional multicellular confocal imaging [17, 18]. Thehigh data content of such imaging enables mapping thefunctional organization of a β -cell collective onto a com-plex network such that a pair of cells (i.e., network nodes)is linked if the cross-correlation between the correspond-ing Ca signals is large enough [19–24].Complex-network representations of the functional or-ganization of β -cell collectives (i.e., functional networks)that arise from calcium cross-correlations in intact pan-creatic islets reveal a modular structure [25] intertwinedwith small-world [20, 26] properties. Of particular inter-est is that some studies [22, 24] point to small subsetsof highly active β cells whose connectedness and impacton synchronization across islets make these cells candi-date leaders. Because of their large degree, we call suchcandidate leaders hub cells . The proponents of the hub-cell idea list many reasons [27] why electrophysiologicalmeasurements may have missed detecting hubs. Overall,functional networks suggest that our army of β cells isled by a few generals.We will hereafter show that the hub-cell idea, ratherthan contradicting, actually complements electrophysiol- a r X i v : . [ phy s i c s . b i o - ph ] J a n ogy. Based on a remarkable agreement between empiricaland computational findings, we will argue in favor of anautopoietic influence hierarchy among β cells. Leadercells emerge through cell-cell communication, and thusneed not be genetically predisposed for leadership. In thelanguage of the army metaphor, generals do lead, but notby birthright; they get promoted by their peers. Leadercells, furthermore, remain under the radar of electrophys-iological measurements by communicating with immedi-ate neighbors just as any other cell would. A direct impli-cation is that a more-or-less arbitrary cell could emergeas a leader, which in turn makes for an ultra robust ar-chitecture of β -cell collectives. Additional wide-reachingimplications underpin some of the vital β -cell features. METHODS
We based our empirical analyses on a dataset ob-tained via Ca imaging of an acutely prepared pancre-atic tissue slice [17, 28] comprising a rodent oval-shapedislet (approx. dimensions: 370 µ m × µ m). Ca sig-nals were recorded using a functional multi-cellular imag-ing technique at 10 Hz and 256 ×
256 pixel resolution in8-bit grayscale color depth. The freely downloadabledataset [29] consisted of 65,536 Ca signals, each with23,873 data points. Our focus was on fast oscillations,which is why all signals were detrended and standard-ized before use [30].For the purpose of constructing empirical functionalnetworks, we randomly picked N = 100 signals, denoted x i ( t ), from the dataset and computed cross-correlations c ij = [ x i x j ], where [ · ] is a time-averaging operator. Thenumber of data points in each signal was L = 3 × , cor-responding to 5-minute exposures of the pancreatic tissueslice to the stimulating glucose concentration of 8 mM.We link two signals i and j in the network if c ij > c ,where c is a threshold selected to give a mean degree (cid:104) k (cid:105) = 10. We have previously shown [31] that the bulkcross-correlation spectrum of typical Ca signals followsthe predictions of random matrix theory, but the statesoutside the bulk-spectrum edges carry biological informa-tion. The delocalized states corresponding to the largesteigenvalues of the cross-correlation matrix, in particular,were found to harbor contributions from all Ca signals,thus revealing an islet-wide collective mode [32].The computational aspects of the study comprised con-structing functional networks from simulated Ca sig-nals. Simulations were based on our model [30] thathad been found to mimic empirical signals closely. Com-pared to typical electrophysiological modeling by meansof differential-algebraic systems [21], the model’s struc-ture is very simple. Each of the N = 100 nodes rep-resents a β cell arranged in a random regular networkwith the degree distribution P ( k ) = δ ( k − k ) such thattrivially the average node degree is (cid:104) k (cid:105) = k . Nodes can k P Model, h = 0.1Model, h = 0.8Empirical, randomEmpirical, per se Figure 1. Functional networks of β cells contain high-degreehubs. The plot shows the configuration-ensemble degree dis-tribution, P , as a function of the node degree, k . change their binary state, indicating the calcium activ-ity of individual β cells, from active to inactive or viceversa in two ways. Internal activation is controlled bya forcing parameter, f , that is interpretable as the in-verse of the glucose concentration to which β cells areexposed. External activation is controlled by the stateof the k nearest neighbors in conjunction with the cou-pling strength, r , that is interpretable as the intensity ofgap-junctional ion exchange. Functional networks wereconstructed from simulated Ca signals in exactly thesame way as from empirical ones. RESULTS
The degree distributions of empirical and simulatedfunctional networks are practically the same (Fig. 1). Wedistinguish between two cases. Open symbols pertain tofunctional networks extracted from empirical data perse (squares), and from simulations with strong coupling, h = 0 . P ( k ) = (cid:0) k + r − k (cid:1) p k (1 − p ) r , where p = (cid:104) k (cid:105) / ( r + (cid:104) k (cid:105) ) ≈ .
77 and r = 3 is a numerical parameter. The oppo-site case, shown using filled symbols, consists of func-tional networks extracted upon randomizing empiricaldata (squares) and running simulations with weak cou-pling, h = 0 . a ij β = β = . Model, h =0.8Empirical k i k j / N ⟨ k ⟩ Figure 2. β -cell functional networks exhibit positive assor-tative mixing. The plot compares the elements a ij of theconfiguration-ensemble adjacency matrix against the normal-ized degree product, ( k i k j / (cid:104) k (cid:105) N ) β , for two different valuesof the β exponent, β = 0 and β = 0 . tion P ( k ) = (cid:0) Nk (cid:1) p k (1 − p ) N − k with p = (cid:104) k (cid:105) /N . Thebinomial distribution in turn converges to a Poisson dis-tribution P ( k ) = (cid:104) k (cid:105) k exp( −(cid:104) k (cid:105) ) /k ! in the large networklimit N → ∞ when (cid:104) k (cid:105) is fixed.Aside from degree distributions, another informativeway of characterizing functional networks is degree cor-relations, also known as degree assortativity. Ref. [33] ex-plains that the elements a ij of a configuration-ensembleadjacency matrix should satisfy the strictly linear rela-tionship a ij = k i k j / (cid:104) k (cid:105) N in zero-assortativity situations.Otherwise, a ij = ( k i k j / (cid:104) k (cid:105) N ) β , where β < β > a ij against k i k j / (cid:104) k (cid:105) N (Fig. 2), while once again revealingan agreement between empirical and simulated functionalnetworks, primarily show that β ≈ .
1. The degree assor-tativity of the functional networks is therefore positive,meaning that network nodes are preferentially linked toother nodes of similar degree.The results so far allow us to estimate the assortativitycoefficient of functional networks in an alternative way,which we can then compare with the definition [34, 35],that is, with the degree correlation coefficient betweenpairs of linked nodes ρ = ( (cid:104) k (cid:105)(cid:104) k k nn (cid:105) − (cid:104) k (cid:105) ) / ( (cid:104) k (cid:105)(cid:104) k (cid:105) −(cid:104) k (cid:105) ), where k i nn = k − i (cid:80) j a ij k j is the average degree ofthe i th node’s nearest neighbors. Using the exponent β ≈ . (cid:28) P ( k ) = q − r k r − exp( − k/q ) / Γ( r ) for large q = (cid:104) k (cid:105) /r , theassortativity coefficient is ρ ≈ βq β ≈ .
11 (see Ref. [33]).Following the definition of ρ indeed yields a very similarvalue (Fig. 3; open symbols for N r /N = 0).We furthermore tested how two different perturba-tions [36], signal randomization and replacement, affectthe assortativity coefficient. Here, randomization andreplacement respectively mean that a fraction N r /N ofsignals used in the construction of the original configura- N r / N ρ λ – ⟨ k ⟩ ρ λ – ⟨ k ⟩ Figure 3. Assortativity and the largest adjacency-matrixeigenvalue of functional networks are robust to random per-turbations, but not to perturbations targeting hubs. Theplot shows assortativity, ρ , and the largest eigenvalue of theconfiguration-ensemble adjacency matrix, λ , (from which theaverage node degree, (cid:104) k (cid:105) , has been subtracted) both as thefunctions of the fraction of perturbed nodes, N r /N . tion ensemble were either randomly reshuffled or replacedwith other random signals from the broader dataset. Wefound that assortativity of functional networks is robustto random perturbations, but not to perturbations tar-geting hubs (Fig. 3; open symbols for N r /N > β -cells because otherwise heterogeneous net-works are expected to be disassortative [33, 35].Another measure characterizing a network’s structureis the largest eigenvalue of the adjacency matrix, λ . Sim-ilarly as with the assortativity coefficient, the value of λ is robust to random perturbations, but not to perturba-tions targeting hubs (Fig. 3). As progressively more hubsget perturbed, the maximum adjacency-matrix eigen-value falls to λ = 1 + (cid:104) k (cid:105) , which is the value character-istic of random networks with a Poisson degree distribu-tion. This last result reaffirms the importance of hubs forthe structure of functional networks and, in turn, hints atan interesting connection with recent experimental find-ings that generated much excitement.Experiments show [22, 24] that perturbing hub cellscan quickly decorrelate a β -cell collective. This impliesthat hubs carry most of the collective’s cross-correlationcontent, which could be tested using our empirical andsimulated Ca signals alike. To this end, we em-ployed the pairwise cross-correlations c ij to define c i = k − i (cid:80) j a ij c ij , that is, the i th node’s neighborhood-wideaverage cross-correlation. By pairing the c i values withthe corresponding node degree k i , we obtained a func-tion c ( k ). If the β -cell collective’s cross-correlation con-tent were equally distributed among all nodes, then we k c (a) h = 0.8 h = 0.8 h = 0.1 h = 0.1 h = 0.85678 k * k (c) k (b) c No pertubation Random perturbation Hub randomization
Figure 4. Hub cells carry most of the cross-correlation content of a β -cell collective. (a) Empirical neighborhood-wide averagecross-correlation of nodes, c , as a function of the node degree, k . Open circles, squares, and triangles respectively denoteno perturbation, random perturbation, and hub randomization. (b) Same as (a), but simulated. Open circles and trianglesrespectively denote no perturbation and random perturbation when coupling, h = 0 . h = 0 .
1, is weak. (c) The simulated average number of active neighbors, k ∗ , as a function ofthe node degree. Open and filled circles respectively denote the strong, h = 0 .
8, and the weak coupling, h = 0 . should see c ( k ) = const . , whereas if hubs carried most ofthe content, then c ( k ) should be an increasing functionof k . We find that the latter is true both for empirical(Fig. 4a) and simulated (Fig. 4b) Ca signals. As withassortativity ρ and the largest eigenvalue λ , random per-turbations have little bearing on the function c ( k ), buttargeting high-degree nodes causes decorrelation. Withthe help of the dynamical network model, we see thatperturbing about 10 % of hubs is very similar to runningsimulations with a weak coupling of h = 0 . c ( k ) couldbe due to hubs somehow communicating with more oftheir neighbors than lower-degree nodes. The number ofactive neighbors k ∗ is, as expected, independent of thedegree k when the coupling is weak, h = 0 . h = 0 . DISCUSSION
Herein, we used a combination of empirical and com-putational methods to shed new light on, among others,a fundamental tension between electrophysiological the-ory and the hub-cell idea as it pertains to pancreatic β cells. We started by constructing functional β -cell net-works from empirical Ca signals, and after that esti-mated the configuration ensemble [37, 38] of such net-works. Meanwhile, we also simulated fast β -cell activity using a dynamical network model [39, 40] whose outputsfaithfully mimic the properties of said empirical Ca sig-nals [30]. Upon repeating the construction of functionalnetworks, but now from simulated signals, we found re-markable quantitative agreement between the results.Irrespective of the signal origin, empirical or simulated,the properties of functional networks—such as assortativ-ity or the largest eigenvalue of the configuration-ensembleadjacency matrix—all critically depend on the presenceof highly connected (i.e., hub) nodes. Our results thussupport the hub-cell idea. A similar heterogeneous orga-nization of functional networks has also been discoveredin the collectives of chemosensing cells [41, 42] commu-nicating via gap junctions.The support for the hub-cell idea would be a strongblow to electrophysiological theory if nodes in the dy-namical network model had any distinctive properties,for example, different degrees or an internal structure.In the model, however, all nodes are completely iden-tical, upholding the ideas of electrophysiological theorythat no cell is predisposed for leadership. Leader cells,in fact, emerge through cell-cell communication as evi-denced by tracing the degree distributions of functionalnetworks to the coupling strength in the dynamical net-work model. In the language of the army metaphor, gen-erals do lead, but not by birthright; they get promotedby their peers. Our findings thus not only reconcile theconcept of hubs with electrophysiological theory, but alsopoint to an extremely robust architecture of collectives.Before expanding on the last remark about the robust-ness of β -cell collectives, let us examine a physiologicaladvantage provided by the heterogeneity of functionalnetworks. Cells have been shown to use cell-cell dis-tance to optimize their sensing precision [43]. In com-pact microorgans such as pancreatic islets, cells can-not easily adjust mutual distances, but could insteaduse functional-network heterogeneity to sharpen theircollective response. Let λ be the largest adjacency-matrix eigenvalue and v the corresponding eigenvec-tor [32, 44]. Because the components of the latter vectorquantify the importance of network nodes, we assumethat y i = (cid:80) j ∈ Ω i v j x j = (cid:80) j a ij v j x j represents the sig-nal that node i integrates from the nearest neighbors j ∈ Ω i , while all other signals are treated as constant-variance noise. Under this assumption, the signal-to-noise ratio of sensing cells is proportional to the vari-ance E [ y ] = (cid:80) ijk a ij v j a ik v k Cov( x j , x k ) (cid:39) c λ , where c is the cross-correlation threshold used during the con-struction of functional networks. The sensing precision isthus proportional to λ , which would equal λ = 1 + (cid:104) k (cid:105) if functional networks were random, but increases to λ = 1 + (cid:104) k (cid:105) + (cid:104) k (cid:105) /r because of heterogeneity. Accord-ingly, cell-cell communication imparts a sharp glucose-sensing acuity to β -cell collectives.The example on glucose-sensing acuity demonstrates aphysiological advantage of heterogeneous functional net-works over homogeneous ones. This advantage would,however, come at a cost of vulnerability to hub-nodefailures [36] if heterogeneity were imprinted in the un-derlying physical network of β cells. Our results insteadstrongly favor an interpretation by which influential cellsmaterialize endogenously within β -cell collectives, givingrise to influence hierarchies that are autopoietic, bothin the self-producing and self-maintaining sense of thisterm. All that is needed for hub cells to emerge amongidentical peers is cell-cell communication, and should ahub fail, there is no obstacle for another one to re-emergeas long as communication remains feasible. It would thusseem that in β -cell collectives, biology has managed tocreate an ultra-robust architecture that is physiologicallyadvantageous as well. This is an impressive feat thatcould perhaps inspire thinking about the design and en-gineering of next-generation critical infrastructure.– – – – – Acknowledgements.
DK, AS, and MSR received fi-nancial support from the Slovenian Research Agency (re-search core funding program no. P3–0396 and projectsno. N3-0048, no. N3-0133 and no. J3-9289). MSR wasfurther supported by the Austrian Science Fund / Fondszur F¨orderung der Wissenschaftlichen Forschung (bilat-eral grants I3562–B27 and I4319–B30). MJ received fi-nancial supported from the Japan Society for the Pro-motion of Science (JSPS) KAKENHI (grant no. JP20H04288) as a co-investigator. BP received financialsupport from the Slovenian Research Agency (project no.J5-8236), the University of Zagreb’s project Advancedmethods and technologies in Data Science and Coopera-tive Systems (DATACROSS; ref. KK.01.1.1.01.009), and the University of Rijeka. PH received financial supportfrom JSPS KAKENHI (grant no. JP 18H01655).
Author contributions.
All authors contributed sub-stantially to all aspects of the study.
Conflict of interest.
The authors declare no conflict ofinterest, financial or otherwise. ∗ Corresponding author: [email protected] † Corresponding author: [email protected][1] D. Pipeleers, Diabetologia , 277 (1987).[2] D. Pipeleers, E. Maes, M. Van De Winkel, et al. , Proc.Natl. Acad. Sci. USA , 7322 (1982).[3] D. Koroˇsak and M. Slak Rupnik, Front. Physiol. , 31(2018).[4] H. Kolb, K. Kempf, M. R¨ohling, and S. Martin, BMCMed. , 224 (2020).[5] P. Meda, Biochim. Biophys. Acta Biomembr. , 124(2018).[6] M. A. Ravier, M. G¨uldenagel, A. Charollais, A. Gjinovci,D. Caille, G. S¨ohl, C. B. Wollheim, K. Willecke, J.-C.Henquin, and P. Meda, Diabetes , 1798 (2005).[7] S. Speier, A. Gjinovci, A. Charollais, P. Meda, andM. Rupnik, Diabetes , 1078 (2007).[8] R. K. Benninger and D. J. Hodson, Diabetes , 537(2018).[9] Q. Zhang, J. Galvanovskis, F. Abdulkader, C. J. Par-tridge, S. O. G¨opel, L. Eliasson, and P. Rorsman, Philos.Trans. R. Soc. A , 3503 (2008).[10] M. Skelin Klemen, J. Dolenˇsek, M. Slak Rupnik, andA. Stoˇzer, Islets , 109 (2017).[11] L. S. Satin, Q. Zhang, and P. Rorsman, Diabetes , 830(2020).[12] N. L. Farnsworth, A. Hemmati, M. Pozzoli, and R. K.Benninger, J. Physiol. , 4431 (2014).[13] A. Caicedo, in Seminars in cell & developmental biology ,Vol. 24(1), edited by S. D. Roper (Elsevier, 2013) pp.11–21.[14] M. J. Berridge, P. Lipp, and M. D. Bootman, Nat. Rev.Mol. Cell Biol. , 11 (2000).[15] H. M. Colecraft, Biophys. J. , 1472 (2020).[16] O. Idevall-Hagren and A. Tengholm, in Seminars in cell& developmental biology , Vol. 103, edited by W. Han andD. Eizirik (Elsevier, 2020) pp. 20–30.[17] A. Stoˇzer, J. Dolenˇsek, and M. S. Rupnik, PLOS One ,e54638 (2013).[18] J. Dolenˇsek, M. S. Klemen, M. Gosak, L. Kriˇzanˇci´c-Bombek, V. Pohorec, M. S. Rupnik, and A. Stoˇzer(2020), bioRxiv 2020.03.11.986893.[19] D. J. Hodson, M. Schaeffer, N. Roman`o, P. Fontanaud,C. Lafont, J. Birkenstock, F. Molino, H. Christian,J. Lockey, D. Carmignac, et al. , Nat. Commun. , 605(2012).[20] A. Stoˇzer, M. Gosak, J. Dolenˇsek, M. Perc, M. Marhl,M. S. Rupnik, and D. Koroˇsak, PLOS Comput. Biol. ,e1002923 (2013).[21] C. Cherubini, S. Filippi, A. Gizzi, and A. Loppini, Phys.Rev. E , 042702 (2015).[22] N. R. Johnston, R. K. Mitchell, E. Haythorne, M. P.Pessoa, F. Semplici, J. Ferrer, L. Piemonti, P. Marchetti,M. Bugliani, D. Bosco, et al. , Cell Metab. , 389 (2016). [23] M. Gosak, R. Markoviˇc, J. Dolenˇsek, M. S. Rupnik,M. Marhl, A. Stoˇzer, and M. Perc, Phys. Life Rev. ,118 (2018).[24] V. Salem, L. D. Silva, K. Suba, E. Georgiadou, S. N. M.Gharavy, N. Akhtar, A. Martin-Alonso, D. C. Gaboriau,S. M. Rothery, T. Stylianides, et al. , Nat. Metab. , 615(2019).[25] R. Markoviˇc, A. Stoˇzer, M. Gosak, J. Dolenˇsek,M. Marhl, and M. S. Rupnik, Sci. Rep. , 7845 (2015).[26] G. A. Rutter and D. J. Hodson, Mol. Endocrinol. ,1984 (2013).[27] G. A. Rutter, N. Ninov, V. Salem, and D. J. Hodson,Diabetes , e10 (2020).[28] S. Speier and M. Rupnik, Pfl¨ugers Arch. , 553 (2003).[29] B. Podobnik, D. Koroˇsak, M. S. Klemen, A. Stoˇzer,J. Dolenˇsek, M. S. Rupnik, P. C. Ivanov, P. Holme,and M. Jusup, β -cells operate collectively to help main-tain glucose homeostasis: dataset (2020), Archived at https://doi.org/10.17605/OSF.IO/NA5H3 .[30] B. Podobnik, D. Koroˇsak, M. S. Klemen, A. Stoˇzer,J. Dolenˇsek, M. S. Rupnik, P. C. Ivanov, P. Holme, andM. Jusup, Biophys. J. , 2588 (2020).[31] M. Slak Rupnik and D. Koroˇsak, Front. Physiol. , 1194(2019).[32] V. Plerou, P. Gopikrishnan, B. Rosenow, L. A. N. Ama-ral, T. Guhr, and H. E. Stanley, Phys. Rev. E , 066126 (2002).[33] S. Johnson, J. J. Torres, J. Marro, and M. A. Munoz,Phys. Rev. Lett , 108702 (2010).[34] M. E. Newman, Phys. Rev. Lett. , 208701 (2002).[35] J. Qu, S.-J. Wang, M. Jusup, and Z. Wang, Sci. Rep. ,15450 (2015).[36] R. Albert, H. Jeong, and A.-L. Barab´asi, Nature ,378 (2000).[37] M. Molloy and B. Reed, Random Struct. Algorithms ,161 (1995).[38] R. R. Nadakuditi and M. E. Newman, Phys. Rev. E ,012803 (2013).[39] A. Majdandzic, B. Podobnik, S. V. Buldyrev, D. Y.Kenett, S. Havlin, and H. E. Stanley, Nat. Phys. , 34(2014).[40] B. Podobnik, M. Jusup, Z. Tiganj, W.-X. Wang, J. M.Buld´u, and H. E. Stanley, Proc. Natl. Acad. Sci. USA , 11826 (2017).[41] B. Sun, G. Duclos, and H. A. Stone, Phys. Rev. Lett. , 158103 (2013).[42] G. D. Potter, T. A. Byrd, A. Mugler, and B. Sun, Proc.Natl. Acad. Sci. USA , 10334 (2016).[43] S. Fancher and A. Mugler, Phys. Rev. Lett. , 078101(2017).[44] I. J. Farkas, I. Der´enyi, A.-L. Barab´asi, and T. Vicsek,Phys. Rev. E64