Analyzing synchronized clusters in neuron networks
Matteo Lodi, Fabio Della Rossa, Francesco Sorrentino, Marco Storace
AAnalyzing synchronized clusters in neuron networks
Matteo Lodi, Fabio Della Rossa,
2, 3
Francesco Sorrentino, and Marco Storace
1, 4 DITEN, University of Genoa, Via Opera Pia 11a, I-16145, Genova, Italy Mechanical Engineering Department, University of New Mexico, Albuquerque, NM 87131, USA DEIB, Politecnico di Milano, I-20133 Milan, Italy Corresponding author: [email protected]. (Dated: July 9, 2020)The presence of synchronized clusters in neuron networks is a hallmark of information transmissionand processing. The methods commonly used to study cluster synchronization in networks of coupledoscillators ground on simplifying assumptions, which often neglect key biological features of neuronnetworks. Here we propose a general framework to study presence and stability of synchronousclusters in more realistic models of neuron networks, characterized by the presence of delays, differentkinds of neurons and synapses. Application of this framework to the directed network of the macaquecerebral cortex provides an interpretation key to explain known functional mechanisms emergingfrom the combination of anatomy and neuron dynamics. The cluster synchronization analysis iscarried out also by changing parameters and studying bifurcations. Despite some simplificationswith respect to the real network, the obtained results are in good agreement with previously reportedbiological data.
Understanding the functional mechanisms ofa given system/phenomenon and describing itthrough mathematical equations as simple as pos-sible (according to the
Occam’s razor principle)is the Holy Grail of modeling. Among the oth-ers, neuron networks are the object of many stud-ies due to their complex behaviors; understandingthe functional mechanisms of information trans-mission and processing in these networks is one ofthe most difficult and fascinating challenges facedby the scientific community, at the crossroad be-tween many disciplines. The level of abstractionused to describe neuron networks can significantlychange according to the modeling goals, complex-ity of the network to be modeled and backgroundknowledge [1]. Consequently, the basic elementsof the nervous system (neurons and synapses) aremodeled by trading off accuracy and complexity[2]. Neurons in the same network can be of differ-ent kinds and their synaptic connections, also ofdifferent kinds, can be either electrical or chemi-cal, either excitatory or inhibitory, either directedor undirected, and may transmit signals with dif-ferent delays. In this letter, we focus on determin-istic models of these networks.A commonly observed phenomenon in networksof neurons is the formation of synchronous clus-ters , i.e., groups of neurons that fulfill some syn-chrony conditions [3–5], usually expressed in termsof temporal correlation between neural signals.These clusters are strongly related to informa- tion transmission and processing [6]. Recent ef-forts have been devoted to apply nonlinear dy-namics concepts and network theory to the neu-roscience context [1, 7]. However, deterministicmodels that study the presence and the stabil-ity of synchronized clusters in networks are basedon simplifying assumptions such as identical neu-rons/synapses, weak interactions, absence of de-lays, undirected/diffusive connections. In this let-ter we propose a general method to analyze clustersynchronization (CS) in neuron networks with di-rected connections, delays, couplings that dependon both the presynaptic and the postsynaptic neu-rons, and different kinds of nodes and synapses.These networks can be described by the follow-ing set of dynamical equations, describing a multi-layer network [8], ( i = 1 , . . . , N )˙ x i = ˜ f i ( x i ( t )) + L (cid:88) k =1 σ k N (cid:88) j =1 A kij h k ( x i ( t ) , x j ( t − δ k )) , (1)where x i ∈ R n is the n -dimensional state vectorof the i -th neuron, ˜ f i : R n → R n is the vectorfield of the isolated i -th neuron, σ k ∈ R is the cou-pling strength of the k -th kind of link, A k is thepossibly weighted and directed coupling matrix (oradjacency matrix) that describes the connectivityof the network with respect to the k -th kind oflink, for which the interaction between two genericcells i and j is described by the nonlinear func-tion h k ( · ) : R n × R n → R n , and δ k is the axon a r X i v : . [ m a t h . D S ] J u l transmission delay characteristic of the k -th kindof link. For example, electrical synapses (gap junc-tions) are almost instantaneous, whereas the delayassociated with transmission of a signal through achemical synapse may be considerably longer.A neuron model is described by a state vector x i , whose first component V i typically representsthe membrane potential of the neuron. A synapsemodel can either neglect or include the neurotrans-mitter dynamics, therefore we can have instanta-neous or dynamical synapses, respectively. In bothcases, we assume that the synaptic coupling influ-ences only the dynamics of V i and not of the otherstate variables contained in x i : therefore, the firstcomponent of the vector h k ( · ) is a scalar function(called activation function ) a k ( V i ( t ) , x j ( t − δ k ))and the remaining components are null. For in-stantaneous synapses, the activation depends onthe membrane potential of the pre- and post-synaptic neurons, therefore it can be expressed as a k ( V i ( t ) , V j ( t − δ k )). By contrast, for dynamicalsynapses the activation a k is a function of a statevariable s kj (in addition to V i ), whose dynamicsusually depends on the pre-synaptic membrane po-tential V j (see Sec. 1 in [9]). For this reason, alldynamical synapses of kind k connecting the neu-ron j with other neurons share the same state s kj ,which can be added to vector x j .We further assume each individual node can beof one out of M different types (with M ≤ N ):˜ f i ( x ) = ˜ f j ( x ) if i and j are of the same type,˜ f i ( x ) (cid:54) = ˜ f j ( x ) otherwise. Often, the difference(physical or functional) between two types of neu-rons is accounted for through a different value ofone or more model parameters. Within this gen-eral framework, where all oscillators can be dif-ferent, if M << N the vector fields ˜ f i are notall different, but belong to a restricted set of M models. Assuming that all node states share thesame dimension n is not restrictive: in the caseof state vectors x i with different dimensions n i , itis sufficient to define n = max i n i and set to 0 thecomponents in excess [9].Different from most models introduced in theliterature, the set of equations (1) accounts forthe following realistic properties of neuron net-works: (i) each synapse depends (algebraically inthe case of instantaneous/fast synapses or dynam-ically in the case of slower synapses) on the stateof both the pre-synaptic and the post-synapticneuron, (ii) each synapse between two neurons is in general a direct connection that can beof different kinds (such as either chemical in-hibitory/excitatory or electrical excitatory), and(iii) the transmission of information along synapsescan be non-instantaneous, which may be due inpart to local synaptic filtering of exchanged spikes,and in part to the distribution of the axonal trans-mission delays [10]. We wish to emphasize thatcurrent methods developed to analyze CS in com-plex networks [11, 12] are unable to handle features(i), (ii) and (iii) above.Cluster synchronization of the system in Eq. (1)is defined as x i ( t ) = x j ( t ) for any t and for i, j be-longing to the same cluster of a certain partition.The set of the network nodes can be partitionedinto equitable clusters (ECs), whose presence isnecessary to achieve CS. Indeed, nodes in the sameEC receive the same amount of weighted inputs ofa certain type from the other ECs or from the ECitself. The method we propose for the analysis ofCS in networks modeled by Eq. (1) consists ofthree main steps: (S1) a coloring algorithm to findthe Q ECs C q ( q = 1 , . . . , Q ) of the network, corre-sponding to a clustering C = { C , . . . , C Q } (see theexample network in Fig. 1A, where N = 11 and Q = 4); (S2) a simplified dynamical model (called quotient network ) whose Q nodes correspond toeach one of the ECs (see Fig. 1B, which is the quo-tient network corresponding to Fig. 1A); (S3) ananalysis of the cluster stability by linearizing Eq.(1) about a state corresponding to exact synchro-nization among all the nodes within each cluster.A detailed description of steps S1, S2, and S3 isprovided in [9]. The main novelty of this methodis the analysis S3, which is tailored to Eq. (1) fol-lowing, mutatis mutandis , the guidelines defined in[11, 12]. A key step of this analysis is the construc-tion of the (unique) matrix T that transforms thecoupling matrices A k into block diagonal matri-ces, B k = T A k T T . This corresponds to a changeof perturbation coordinates that converts the nodecoordinate system to the irreducible representation [11–13] coordinate system, thus evidencing the in-terdependencies among the perturbation compo-nents. For undirected networks, the N × N matrix T can be found from the symmetry group of thenetwork, as done in [12] for the orbital case and in[14] for the equitable single-layer case. For directednetworks of two specific kinds (detailed in [9]), thematrix T can be constructed as described in Sec.2.3 in [9]. The key variational equation is reported ! ! ! " !" ! ! ! $ ! " ! ! " " $ ! " ! " ! ! $ ! " ! ! $ ! ! % % ' & & ' (! " ) ! $ " " $ " ! " ( " ) " * " + " , " "- " "" ) " ) $ ) ! ) . ! ) . " ) . * . $ %+%
22 222
22 222
22 222
AC DE
C1C2C3 C4 B FIG. 1. Example: (A) a network with N = 11 nodes, L = 2 kind of connections, and Q = 4 clusters ( C = { , , , } , C = { , } , C = { , , } , C = { , } ) and (B) the structure of the corresponding matrices T and B , illustrating their relation with the clusters. All connections are bi-directional and with weight 1, with theexception of the thick connection between nodes 5 and 7, which has weight 2. Network coloring (with a largernumber of clusters) after the breaking of the red cluster if its loss of stability is due to the MLEs correspondingto either (C) the multi-color sub-block or (D) the red sub-block of matrix B . here in compact form for ease of reference:˙ ηηη = [ ρ + ρ ( B k )] ηηη. (2)This equation describes the perturbation dynam-ics, by separating that along the synchronous man-ifold (described by the first Q components η i ) fromthat transverse to it (described by the last compo-nents η i , i ∈ [ Q + 1 , N ]). We remark that theterm ρ in Eq. (2) is a diagonal matrix, which re-lates ˙ η j only to η j . By contrast, ρ relates ˙ η j alsoto the other perturbation components through thematrix B k . Therefore, an inspection of the sub-blocks of B k allows to quickly check whether thereis coupling between the dynamics of perturbations η i and η j . To better illustrate this concept, letus consider the undirected, weighted network with N = 11 nodes, L = 2 kind of connections, and Q = 4 clusters ( C , C , C , C ) shown in Fig. 1,panel A, with nodes color coded to indicate theclusters they belong to. Note that the partitionof the network nodes is equitable and not orbital [14]. Notice also the presence of a delay δ in theconnection between nodes 5 and 6.Panel C shows the structure of the matrices T (left) and B (right) for this network. Notice thatmatrix B has the same structure as B , whosegray blocks contain only 0 entries. The upper-left Q × Q block is related to the perturbation dynam-ics along the synchronous manifold. Each whitesub-block in the lower-right ( N − Q ) × ( N − Q )sub-matrix B N − Q (with dashed black borders) de-scribes the perturbation dynamics transverse tothe synchronous manifold, thus is associated withloss of synchronization, either transient or perma-nent depending on the cluster stability. For in-stance, the 1 × C (or C or C , respectively),as pointed out in the corresponding row in matrix T , and describes the dynamics of the perturbationcomponent η (or η or η , respectively); simi-larly, the 4 × C , C , C . We remark that the structureof this sub-block implies that ˙ η , ˙ η , ˙ η , ˙ η dependon η , η , η , η but not on the other perturbations.Each transverse sub-block has an associated Max-imum Laypunov Exponent (MLE) Λ i , which canbe studied independently from each other.The stability of each cluster C q related to one ormore sub-blocks depends on the maximum MLEΛ C q among those associated to these sub-blocks: ifΛ C q is negative, the cluster C q is stable, otherwiseit is unstable. In the example, we computed theMLE associated to each sub-block: Λ (blue sub-block), Λ (multi-color sub-block), Λ (red sub-block) and Λ (yellow sub-block). The stability of C depends on the sign of Λ C = Λ = max { λ } (i.e., the maximum component of the vector λ ),whereas the stability of C depends on the sign ofΛ C = max { Λ , Λ } , the stability of C dependson the sign of Λ C = Λ and the stability of C depends on the sign of Λ C = max { Λ , Λ } .Notice that the structure of the matrix B al-lows us to state something more about the clusterstability. Indeed, the red cluster is related to twosub-blocks: the 1 × × be-comes positive, while if the MLE Λ becomes posi-tive, red, blue, and green clusters become unstabletogether (see panel D).This example clearly shows that the stability ofeach cluster in a subset of intertwined clusters [11]may depend on the stability of the other clustersthat belong to the same subset, but is decoupledfrom the clusters outside of the subset. Therefore,intertwined clusters can lose synchronization with-out causing a loss of synchronization in the clustersoutside the subset, as for the yellow cluster in panelD.We apply the proposed method to a directednetwork (shown in Fig. 2) composed of N = 29nodes, each one representing one of the 91 areasof the macaque cerebral cortex [15, 16]. The neu-ron models that represent each area are of M = 2kinds: 28 nodes are of kind i = 1 and one node(corresponding to area V1) is of kind i = 2, whichis due to this one node receiving a visual input [17].The nodes are connected through L = 2 kinds ofchemical excitatory synapses: one (for k = 1) thattransmits undelayed signals with δ = 0 (in yel-low), one (for k = 2) with delay δ > A and A provided in[9] (dataset S1). The measured connection weights[16], which range between 0 and 0.7636, have beenquantized on four levels (0, 0.1, 0.5, 1) by replacingeach original weight with the closest one accord-ing to the Euclidean distance. After that, phys-ical connections with length lower than 20 mmhave been considered instantaneous (i.e., of kind k = 1) and the corresponding quantized weightshave been stored in the matrix A , whereas thoselonger than 20mm have been considered delayed(i.e., of kind k = 2) and the corresponding quan-tized weights have been stored in the matrix A .These quantizations are justified by the fact thatexact values for the coupling strengths and the de-lays reported in the literature are inevitably sub-ject to measurement noise, and by the fact that, aswe will see, they lead to the observation of func-tional mechanisms which are in agreement withphysiological data, despite our simplifications.The network non-trivial equitable clusters (con-sisting of more than one node) are listed in Tab. I.The same information is provided in Fig. 2, wherenodes of the same color (excluding black) belongto the same cluster: green for C , red for C andblue for C . All nodes in trivial orbits are coloredblack. Obviously, the presence of a large numberof trivial clusters does not mean that the corre-sponding areas are independent: they are denselyconnected, as evidenced in Fig. 2, but they cannotbe exactly synchronized.Despite the rough quantizations applied tosynaptic weights and delays, the clusters displayed FIG. 2. Macaque cortical connectivity network: N =29 nodes, M = 2 node models, L = 2 synapse models.Trivial clusters are black. Nodes of the same (non-black) color belong to the same cluster. TABLE I. ECs of the macaque cortical networkCluster index nodes in cluster C C C 𝐵 ! 𝑄 𝑥 𝑁 𝐶 % 𝐶 & 𝐶 ’𝑄 𝑥 𝑁 𝐵 " 𝐶 ! … 𝐶 " 𝐶 ! 𝐶 " 𝐶 𝜂 "$ 𝜂 "% 𝜂 "& 𝑇 𝐶 ! 𝜂 "’ FIG. 3. Structure of the matrices T , B , B , and B for the macaque cerebral cortex. The gray blocks cor-respond to 0 entries. in Table I are consistent with some previously re-ported physiological findings. For instance, cluster C contains the nodes corresponding to visual ar-eas 8l and 9/46v in the prefrontal cortex, whichare known to be physically close and with similarconnections [16, 18]. The same holds for cluster C , which contains the nodes corresponding to theposterior and anterior portion of the inferotempo-ral cortex (TEO and TEpd, respectively).The directed connections originate from or go totrivial clusters only, therefore the cluster stabilitycould be analyzed through the proposed approach.Fig. 3 shows the structure of the matrices T (left)and B k (right) for this network.In the block upper-triangular matrices B k , thefirst Q rows are related to the perturbation dynam-ics along the synchronous manifold. Each whitesub-block in the lower-right ( N − Q ) × ( N − Q )sub-matrix B kN − Q describes the perturbation dy-namics transverse to the synchronous manifold.If we analyze the matrices B k (related to the k -th connection type), we can see that: B N − Q (re-lated to undelayed chemical excitatory synapses)has only zero entries, which implies that for thenetwork with only these synapses the dynamics ofeach perturbation component η k depends only on η k through the term ρ in Eq. (2); B N − Q (re-lated to delayed chemical excitatory synapses) hasone 1 × C , which means that for the network withonly the delayed chemical excitatory synapses thedynamics of the perturbation component η de-pends on η through both ρ and ρ in Eq. (2). In summary, if we consider the whole network, withall kinds of synapses, the three clusters C , C , C turn out to be not intertwined.The stability analysis has been carried out byvarying the delay δ between 0 and 16 ms (8 evenlyspaced values). The neurons belonging to cluster C do not receive any synaptic inputs, thereforethe cluster transverse MLE is Λ C = 0 for anyvalue of δ . Figure 4, panel A, shows the MLEsΛ C q of the other clusters C q ( q = 2 ,
3) versus thedelay δ . The green (red) regions in each plotΛ C q ( δ ) denote stability (instability) of the corre-sponding cluster C q .The vertical dotted lines mark the δ values cor-responding to the time plots shown in panel B: δ = 5 ms (a) and δ = 15 ms (b). These plotsdisplay the first state variable V i of the neurons incluster C . The panels show a window of 300 ms -202 0 300-202 (a)B (b)(a) A (b) FIG. 4. (A) MLE Λ C q of each cluster C q ( q = 2 ,
3) vs.coupling delay δ , for the macaque cortical connectivitynetwork. Horizontal dashed lines: edge of stability.Vertical dotted lines: δ values corresponding to thetime plots in panel B. (B) Time plots V i ( t ) for differentvalues of δ (5 ms (a), 15 ms (b)) for cluster C . inputV1V48m5TEO7A9/46dTEpd FIG. 5. Time responses (firing rates) to a pulse-shapedinput to area V1. after a transient of 19.5 s. The breaking of thiscluster is caused by a supercritical pitchfork bifur-cation of cycles at each transition between the redand green regions, which generates two smaller sta-ble trivial sub-clusters, each one producing one ofthe membrane voltages (black or red) in panel B,plot (b).From Fig. 4 it clearly emerges that the two neu-rons in cluster C display a phase lag for δ = 15ms. The synchronization of macaque visual cor-tex areas in response to visual stimuli has beenobserved in many experiments [17, 19]. In partic-ular, the areas 8l and 9/46v respond in a very sim-ilar way to visual inputs to area V1 [17]. We thusset δ = 5ms in order to ensure synchronization ofthese two areas.We proceeded to validate our model against thequantizations applied to the synaptic weights andaxon delays, described before. To this end, follow-ing [17] we simulated its response to a pulsed inputto the primary visual cortex (area V1). The re-sponse is propagated up the visual hierarchy, pro-gressively slowing as it proceeds, as shown in Fig.5. Early visual areas, such as V1 and V4, exhibitfast responses. By contrast, prefrontal areas, suchas 8m and 24c, exhibit slower decays to the stan-dard firing rate, with traces of the stimulus per-sisting several seconds after stimulation. This is in agreement with the results shown in [17] (comparewith Fig. 3 in [17]), which unveil a circuit mech-anism for hierarchical processing of visual stimuliin the macaque cortex. Moreover, Fig. 5 evidencesCS of the areas TEO and TEpd, corresponding tocluster C , as predicted by Fig. 4.The framework proposed in this paper is a fun-damental step towards a method that fills the gapbetween the analysis of CS and networks of neu-rons. The proposed method, based on a multi-layer network, allowed us to analyze the CS in themacaque cerebral cortex. A second example, whichillustrates symultaneous loss of stability of inter-twined clusters in a smaller network is provided in[9], showing an excellent agreement with biologicalmeasurements.The results of this paper can be extended tostudy synchronization in any network character-ized by different nodes, connections, and commu-nication delays. As a final remark, we point outthat the proposed model is completely determinis-tic and assumes that a reliable model of the net-work is available. These are quite strong model-ing assumptions, since in real neuron networks thepresence of noise is unavoidable and not alwaysneuron and synapse models can be determined ac-curately. Despite this and despite the absence ofinformation about the basins of attraction of stableclusters, our approach can provide useful informa-tion. For instance, in a real network CS will beapproximate [20], not exact, as measured by highcorrelation values between the membrane poten-tials of the neurons/nodes belonging to a given sta-ble cluster. In this perspective, the patterns foundwith the proposed method are approximations tosome more realistic solutions, which are character-ized by higher complexity. Nature is quite far fromdeterminism, therefore our analysis method is farfrom providing a general description of the dynam-ics of real neuron networks. This notwithstanding,it provides basic understanding of CS mechanisms,whose robustness can be checked by resorting toother less deterministic approaches.The authors would like to express their sincereappreciation to Maurizio Mattia, Mauro Parodiand Lou Pecora for many useful inputs and valu-able comments. [1] D. S. Bassett, P. Zurn, and J. I. Gold, NatureReviews Neuroscience , 353 (2018).[2] A. V. Herz, T. Gollisch, C. K. Machens, andD. Jaeger, Science , 80 (2006).[3] A. K. Kreiter and W. Singer, Journal of neuro-science , 2381 (1996).[4] P. E. Maldonado, S. Friedman-Hill, and C. M.Gray, Cerebral Cortex , 1117 (2000).[5] M. Glennon, M. A. Keane, M. A. Elliott, andP. Sauseng, Cerebral Cortex , 2035 (2016).[6] E. Bullmore and O. Sporns, Nature reviews neu-roscience , 186 (2009).[7] R. Guevara Erra, J. L. Perez Velazquez, andM. Rosenblum, Frontiers in computational neu-roscience , 98 (2017).[8] S. Boccaletti, G. Bianconi, R. Criado, C. I.Del Genio, J. G´omez-Gardenes, M. Romance,I. Sendina-Nadal, Z. Wang, and M. Zanin, PhysicsReports , 1 (2014).[9] See Supplemental Material for further examples,mathematical details on the proposed method,and additional datasets.[10] M. Mattia, M. Biggio, A. Galluzzi, andM. Storace, PLoS computational biology (2019).[11] L. M. Pecora, F. Sorrentino, A. M. Hagerstrom,T. E. Murphy, and R. Roy, Nature communica-tions , 4079 (2014). [12] F. Della Rossa, L. M. Pecora, K. Blaha, A. Shirin,I. Klickstein, and F. Sorrentino, Nature Commu-nications , 1 (2020).[13] M. Golubitsky, I. Stewart, and D. G. Schaef-fer, Singularities and groups in bifurcation theory ,Vol. 2 (Springer Science & Business Media, 2012).[14] A. B. Siddique, L. Pecora, J. D. Hart, and F. Sor-rentino, Physical Review E , 042217 (2018).[15] N. T. Markov, M. Ercsey-Ravasz, D. C. Van Es-sen, K. Knoblauch, Z. Toroczkai, and H. Kennedy,Science , 1238406 (2013).[16] N. T. Markov, M. Ercsey-Ravasz,A. Ribeiro Gomes, C. Lamy, L. Magrou, J. Vezoli,P. Misery, A. Falchier, R. Quilodran, M. Gariel, et al. , Cerebral cortex , 17 (2014).[17] R. Chaudhuri, K. Knoblauch, M.-A. Gariel,H. Kennedy, and X.-J. Wang, Neuron , 419(2015).[18] A. Goulas, A. Schaefer, and D. S. Margulies, BrainStructure and Function , 2939 (2015).[19] C. A. Bosman, J.-M. Schoffelen, N. Brunet,R. Oostenveld, A. M. Bastos, T. Womelsdorf,B. Rubehn, T. Stieglitz, P. De Weerd, and P. Fries,Neuron , 875 (2012).[20] F. Sorrentino and L. Pecora, Chaos: An Interdis-ciplinary Journal of Nonlinear Science26