Analyzing synchronized clusters in neuron networks

The presence of synchronized clusters in neuron networks is a hallmark of information transmission and processing. The methods commonly used to study cluster synchronization in networks of coupled oscillators ground on simplifying assumptions, which often neglect key biological features of neuron networks. Here we propose a general framework to study presence and stability of synchronous clusters in more realistic models of neuron networks, characterized by the presence of delays, different kinds of neurons and synapses. Application of this framework to the directed network of the macaque cerebral cortex provides an interpretation key to explain known functional mechanisms emerging from the combination of anatomy and neuron dynamics. The cluster synchronization analysis is carried out also by changing parameters and studying bifurcations. Despite some simplifications with respect to the real network, the obtained results are in good agreement with previously reported biological data.

The presence of synchronized clusters in neuron networks is a hallmark of information transmission and processing. The methods commonly used to study cluster synchronization in networks of coupled oscillators ground on simplifying assumptions, which often neglect key biological features of neuron networks. Here we propose a general framework to study presence and stability of synchronous clusters in more realistic models of neuron networks, characterized by the presence of delays, different kinds of neurons and synapses. Application of this framework to the directed network of the macaque cerebral cortex provides an interpretation key to explain known functional mechanisms emerging from the combination of anatomy and neuron dynamics. The cluster synchronization analysis is carried out also by changing parameters and studying bifurcations. Despite some simplifications with respect to the real network, the obtained results are in good agreement with previously reported biological data.
Understanding the functional mechanisms of a given system/phenomenon and describing it through mathematical equations as simple as possible (according to the Occam's razor principle) is the Holy Grail of modeling. Among the others, neuron networks are the object of many studies due to their complex behaviors; understanding the functional mechanisms of information transmission and processing in these networks is one of the most difficult and fascinating challenges faced by the scientific community, at the crossroad between many disciplines. The level of abstraction used to describe neuron networks can significantly change according to the modeling goals, complexity of the network to be modeled and background knowledge [1]. Consequently, the basic elements of the nervous system (neurons and synapses) are modeled by trading off accuracy and complexity [2]. Neurons in the same network can be of different kinds and their synaptic connections, also of different kinds, can be either electrical or chemical, either excitatory or inhibitory, either directed or undirected, and may transmit signals with different delays. In this letter, we focus on deterministic models of these networks.
A commonly observed phenomenon in networks of neurons is the formation of synchronous clusters, i.e., groups of neurons that fulfill some synchrony conditions [3][4][5], usually expressed in terms of temporal correlation between neural signals. These clusters are strongly related to informa-tion transmission and processing [6]. Recent efforts have been devoted to apply nonlinear dynamics concepts and network theory to the neuroscience context [1,7]. However, deterministic models that study the presence and the stability of synchronized clusters in networks are based on simplifying assumptions such as identical neurons/synapses, weak interactions, absence of delays, undirected/diffusive connections. In this letter we propose a general method to analyze cluster synchronization (CS) in neuron networks with directed connections, delays, couplings that depend on both the presynaptic and the postsynaptic neurons, and different kinds of nodes and synapses. These networks can be described by the following set of dynamical equations, describing a multilayer network [8], (i = 1, . . . , N ) where x i ∈ R n is the n-dimensional state vector of the i-th neuron,f i : R n → R n is the vector field of the isolated i-th neuron, σ k ∈ R is the coupling strength of the k-th kind of link, A k is the possibly weighted and directed coupling matrix (or adjacency matrix) that describes the connectivity of the network with respect to the k-th kind of link, for which the interaction between two generic cells i and j is described by the nonlinear function h k (·) : R n × R n → R n , and δ k is the axon transmission delay characteristic of the k-th kind of link. For example, electrical synapses (gap junctions) are almost instantaneous, whereas the delay associated with transmission of a signal through a chemical synapse may be considerably longer.
A neuron model is described by a state vector x i , whose first component V i typically represents the membrane potential of the neuron. A synapse model can either neglect or include the neurotransmitter dynamics, therefore we can have instantaneous or dynamical synapses, respectively. In both cases, we assume that the synaptic coupling influences only the dynamics of V i and not of the other state variables contained in x i : therefore, the first component 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 instantaneous synapses, the activation depends on the membrane potential of the pre-and postsynaptic neurons, therefore it can be expressed as a k (V i (t), V j (t − δ k )). By contrast, for dynamical synapses the activation a k is a function of a state variable s k j (in addition to V i ), whose dynamics usually depends on the pre-synaptic membrane potential V j (see Sec. 1 in [9]). For this reason, all dynamical synapses of kind k connecting the neuron j with other neurons share the same state s k j , which can be added to vector x j .
We further assume each individual node can be of 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) =f j (x) otherwise. Often, the difference (physical or functional) between two types of neurons is accounted for through a different value of one or more model parameters. Within this general framework, where all oscillators can be different, if M << N the vector fieldsf i are not all different, but belong to a restricted set of M models. Assuming that all node states share the same dimension n is not restrictive: in the case of state vectors x i with different dimensions n i , it is sufficient to define n = max i n i and set to 0 the components in excess [9].
Different from most models introduced in the literature, the set of equations (1) accounts for the following realistic properties of neuron networks: (i) each synapse depends (algebraically in the case of instantaneous/fast synapses or dynamically in the case of slower synapses) on the state of both the pre-synaptic and the post-synaptic neuron, (ii) each synapse between two neurons is in general a direct connection that can be of different kinds (such as either chemical inhibitory/excitatory or electrical excitatory), and (iii) the transmission of information along synapses can be non-instantaneous, which may be due in part to local synaptic filtering of exchanged spikes, and in part to the distribution of the axonal transmission delays [10]. We wish to emphasize that current methods developed to analyze CS in complex 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 belonging to the same cluster of a certain partition. The set of the network nodes can be partitioned into equitable clusters (ECs), whose presence is necessary to achieve CS. Indeed, nodes in the same EC receive the same amount of weighted inputs of a certain type from the other ECs or from the EC itself. The method we propose for the analysis of CS in networks modeled by Eq. (1) consists of three main steps: (S1) a coloring algorithm to find the Q ECs C q (q = 1, . . . , Q) of the network, corresponding to a clustering C = {C 1 , . . . , C Q } (see the example network in Fig. 1A, where N = 11 and Q = 4); (S2) a simplified dynamical model (called quotient network ) whose Q nodes correspond to each one of the ECs (see Fig. 1B, which is the quotient network corresponding to Fig. 1A); (S3) an analysis of the cluster stability by linearizing Eq.
(1) about a state corresponding to exact synchronization among all the nodes within each cluster. A detailed description of steps S1, S2, and S3 is provided in [9]. The main novelty of this method is the analysis S3, which is tailored to Eq. (1) following, mutatis mutandis, the guidelines defined in [11,12]. A key step of this analysis is the construction of the (unique) matrix T that transforms the coupling matrices A k into block diagonal matrices, B k = T A k T T . This corresponds to a change of perturbation coordinates that converts the node coordinate system to the irreducible representation [11][12][13] coordinate system, thus evidencing the interdependencies among the perturbation components. For undirected networks, the N × N matrix T can be found from the symmetry group of the network, as done in [12] for the orbital case and in [14] for the equitable single-layer case. For directed networks of two specific kinds (detailed in [9]), the matrix T can be constructed as described in Sec. 2.3 in [9]. The key variational equation is reported here in compact form for ease of reference: This equation describes the perturbation dynamics, by separating that along the synchronous manifold (described by the first Q components η i ) from that transverse to it (described by the last components η i , i ∈ [Q + 1, N ]). We remark that the term ρ 1 in Eq. (2) is a diagonal matrix, which relatesη j only to η j . By contrast, ρ 2 relatesη j also to the other perturbation components through the matrix B k . Therefore, an inspection of the subblocks of B k allows to quickly check whether there is coupling between the dynamics of perturbations η i and η j . To better illustrate this concept, let us consider the undirected, weighted network with N = 11 nodes, L = 2 kind of connections, and Q = 4 clusters (C 1 , C 2 , C 3 , C 4 ) shown in Fig. 1, panel A, with nodes color coded to indicate the clusters they belong to. Note that the partition of the network nodes is equitable and not orbital [14]. Notice also the presence of a delay δ 2 in the connection between nodes 5 and 6.
Panel C shows the structure of the matrices T (left) and B 1 (right) for this network. Notice that matrix B 2 has the same structure as B 1 , whose gray blocks contain only 0 entries. The upper-left Q × Q block is related to the perturbation dynamics along the synchronous manifold. Each white sub-block in the lower-right N −Q (with dashed black borders) describes the perturbation dynamics transverse to the synchronous manifold, thus is associated with loss of synchronization, either transient or permanent depending on the cluster stability. For instance, the 1 × 1 yellow (or blue or red) sub-block, is related to cluster C 4 (or C 3 or C 1 , respectively), as pointed out in the corresponding row in matrix T , and describes the dynamics of the perturbation component η 11 (or η 5 or η 10 , respectively); similarly, the 4×4 multi-color sub-block corresponds to clusters C 1 , C 2 , C 3 . We remark that the structure of this sub-block implies thatη 6 ,η 7 ,η 8 ,η 9 depend on η 6 , η 7 , η 8 , η 9 but not on the other perturbations. Each transverse sub-block has an associated Maximum Laypunov Exponent (MLE) Λ i , which can be studied independently from each other.
The stability of each cluster C q related to one or more sub-blocks depends on the maximum MLE Λ Cq among those associated to these sub-blocks: if Λ Cq is negative, the cluster C q is stable, otherwise it is unstable. In the example, we computed the MLE associated to each sub-block: Λ 1 (blue subblock), Λ 2 (multi-color sub-block), Λ 3 (red subblock) and Λ 4 (yellow sub-block). The stability of C 4 depends on the sign of Λ C4 = Λ 4 = max{λ 11 } (i.e., the maximum component of the vector λ 11 ), whereas the stability of C 1 depends on the sign of Λ C1 = max{Λ 2 , Λ 3 }, the stability of C 2 depends on the sign of Λ C2 = Λ 2 and the stability of C 3 depends on the sign of Λ C3 = max{Λ 1 , Λ 2 }.
Notice that the structure of the matrix B 1 allows us to state something more about the cluster stability. Indeed, the red cluster is related to two sub-blocks: the 1 × 1 red sub-block and the 4 × 4 multi-color sub-block. This means what follows: it is possible for the red cluster to undergo isolated desynchronization (see panel E) if the MLE Λ 3 becomes positive, while if the MLE Λ 2 becomes positive, red, blue, and green clusters become unstable together (see panel D).
This example clearly shows that the stability of each cluster in a subset of intertwined clusters [11] may depend on the stability of the other clusters that belong to the same subset, but is decoupled from the clusters outside of the subset. Therefore, intertwined clusters can lose synchronization without causing a loss of synchronization in the clusters outside the subset, as for the yellow cluster in panel D.
We apply the proposed method to a directed network (shown in Fig. 2) composed of N = 29 nodes, each one representing one of the 91 areas of the macaque cerebral cortex [15,16]. The neuron models that represent each area are of M = 2 kinds: 28 nodes are of kind i = 1 and one node (corresponding to area V1) is of kind i = 2, which is due to this one node receiving a visual input [17]. The nodes are connected through L = 2 kinds of chemical excitatory synapses: one (for k = 1) that transmits undelayed signals with δ 1 = 0 (in yellow), one (for k = 2) with delay δ 2 > 0 (in blue). The overall network is modeled by using the neuron and synapse equations described in Sec. 3 of [9] and the coupling matrices A 1 and A 2 provided in [9] (dataset S1). The measured connection weights [16], which range between 0 and 0.7636, have been quantized on four levels (0, 0.1, 0.5, 1) by replacing each original weight with the closest one according to the Euclidean distance. After that, physical connections with length lower than 20 mm have been considered instantaneous (i.e., of kind k = 1) and the corresponding quantized weights have been stored in the matrix A 1 , whereas those longer than 20mm have been considered delayed (i.e., of kind k = 2) and the corresponding quantized weights have been stored in the matrix A 2 . These quantizations are justified by the fact that exact values for the coupling strengths and the delays reported in the literature are inevitably subject to measurement noise, and by the fact that, as we will see, they lead to the observation of functional mechanisms which are in agreement with physiological data, despite our simplifications.
The network non-trivial equitable clusters (consisting of more than one node) are listed in Tab. I. The same information is provided in Fig. 2, where nodes of the same color (excluding black) belong to the same cluster: green for C 1 , red for C 2 and blue for C 3 . All nodes in trivial orbits are colored black. Obviously, the presence of a large number of trivial clusters does not mean that the corresponding areas are independent: they are densely connected, as evidenced in Fig. 2, but they cannot be exactly synchronized.
Despite the rough quantizations applied to synaptic weights and delays, the clusters displayed   in Table I are consistent with some previously reported physiological findings. For instance, cluster C 2 contains the nodes corresponding to visual areas 8l and 9/46v in the prefrontal cortex, which are known to be physically close and with similar connections [16,18]. The same holds for cluster C 3 , which contains the nodes corresponding to the posterior and anterior portion of the inferotemporal cortex (TEO and TEpd, respectively). The directed connections originate from or go to trivial clusters only, therefore the cluster stability could 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 , the first Q rows are related to the perturbation dynamics along the synchronous manifold. Each white sub-block in the lower-right (N − Q) × (N − Q) sub-matrix B k N −Q describes the perturbation dynamics transverse to the synchronous manifold.
If we analyze the matrices B k (related to the k-th connection type), we can see that: B 1 N −Q (related to undelayed chemical excitatory synapses) has only zero entries, which implies that for the network with only these synapses the dynamics of each perturbation component η k depends only on η k through the term ρ 1 in Eq. (2); B 2 N −Q (related to delayed chemical excitatory synapses) has one 1 × 1 sub-block (with blue borders) related to cluster C 3 , which means that for the network with only the delayed chemical excitatory synapses the dynamics of the perturbation component η 29 depends on η 29 through both ρ 1 and ρ 2 in Eq. (2).
In summary, if we consider the whole network, with all kinds of synapses, the three clusters C 1 , C 2 , C 3 turn out to be not intertwined.
The stability analysis has been carried out by varying the delay δ 2 between 0 and 16 ms (8 evenly spaced values). The neurons belonging to cluster C 1 do not receive any synaptic inputs, therefore the cluster transverse MLE is Λ C1 = 0 for any value of δ 2 . Figure 4, panel A, shows the MLEs Λ Cq of the other clusters C q (q = 2, 3) versus the delay δ 2 . The green (red) regions in each plot Λ Cq (δ 2 ) denote stability (instability) of the corresponding cluster C q .
The vertical dotted lines mark the δ 2 values corresponding to the time plots shown in panel B: δ 2 = 5 ms (a) and δ 2 = 15 ms (b). These plots display the first state variable V i of the neurons in cluster C 3 . The panels show a window of 300 ms  after a transient of 19.5 s. The breaking of this cluster is caused by a supercritical pitchfork bifurcation of cycles at each transition between the red and green regions, which generates two smaller stable trivial sub-clusters, each one producing one of the membrane voltages (black or red) in panel B, plot (b). From Fig. 4 it clearly emerges that the two neurons in cluster C 3 display a phase lag for δ 2 = 15 ms. The synchronization of macaque visual cortex areas in response to visual stimuli has been observed in many experiments [17,19]. In particular, the areas 8l and 9/46v respond in a very similar way to visual inputs to area V1 [17]. We thus set δ 2 = 5ms in order to ensure synchronization of these two areas.
We proceeded to validate our model against the quantizations applied to the synaptic weights and axon delays, described before. To this end, following [17] we simulated its response to a pulsed input to the primary visual cortex (area V1). The response is propagated up the visual hierarchy, progressively slowing as it proceeds, as shown in Fig.  5. Early visual areas, such as V1 and V4, exhibit fast responses. By contrast, prefrontal areas, such as 8m and 24c, exhibit slower decays to the standard firing rate, with traces of the stimulus persisting several seconds after stimulation. This is in agreement with the results shown in [17] (compare with Fig. 3 in [17]), which unveil a circuit mechanism for hierarchical processing of visual stimuli in the macaque cortex. Moreover, Fig. 5 evidences CS of the areas TEO and TEpd, corresponding to cluster C 3 , as predicted by Fig. 4.
The framework proposed in this paper is a fundamental step towards a method that fills the gap between the analysis of CS and networks of neurons. The proposed method, based on a multilayer network, allowed us to analyze the CS in the macaque cerebral cortex. A second example, which illustrates symultaneous loss of stability of intertwined clusters in a smaller network is provided in [9], showing an excellent agreement with biological measurements.
The results of this paper can be extended to study synchronization in any network characterized by different nodes, connections, and communication delays. As a final remark, we point out that the proposed model is completely deterministic and assumes that a reliable model of the network is available. These are quite strong modeling assumptions, since in real neuron networks the presence of noise is unavoidable and not always neuron and synapse models can be determined accurately. Despite this and despite the absence of information about the basins of attraction of stable clusters, our approach can provide useful information. For instance, in a real network CS will be approximate [20], not exact, as measured by high correlation values between the membrane potentials of the neurons/nodes belonging to a given stable cluster. In this perspective, the patterns found with the proposed method are approximations to some more realistic solutions, which are characterized by higher complexity. Nature is quite far from determinism, therefore our analysis method is far from providing a general description of the dynamics of real neuron networks. This notwithstanding, it provides basic understanding of CS mechanisms, whose robustness can be checked by resorting to other less deterministic approaches.
The authors would like to express their sincere appreciation to Maurizio Mattia, Mauro Parodi and Lou Pecora for many useful inputs and valuable comments.