We study the dynamics of a general adaptive multi-layer network, which encompasses a wide variety of single-node dynamical models, different connection types, and plasticity rules. The dynamical evolution of this adaptive network composed of *N* coupled nodes can be described by the following general set of delay differential equations (*i*, *j* = 1, …, *N*):

$$\dot\boldsymbolx_i(t)={{{{{\boldsymbolF}}}}}_i(\boldsymbolx_i(t))+{\boldsymbolS}\left({{{{\boldsymbolx}}}}_i,\sum _\ell =1^L\sigma ^\ell \sum_ja_ij^\ell b_ij^\ell (t){{\boldsymbolg}}^\ell ({{{{{{{{\boldsymbolx}}}}}}}}_j(t-\tau ))\right)$$

(1a)

$$\dotb_ij^\ell (t)=-\epsilon ^\ell \left(b_ij^\ell (t)+H^\ell ({{{{{{{{\boldsymbolx}}}}}}}}_i(t),{{{{{{{{\boldsymbolx}}}}}}}}_j(t))\right),$$

(1b)

where *x*_{i} is the *m*-dimensional state vector of the node *i*, with individual dynamics described by the vector field *F*_{i} ∈ **R**^{m}, *i* = 1, …, *N*, and ** S** is an

*m*-dimensional function (e.g., a saturation function) describing the coupling of the node with the rest of the network. This network contains

*L*different types of connections, each corresponding to a different layer, indexed by

*ℓ*. For each layer (

*ℓ*= 1, 2, …,

*L*), \(A^\ell =\a_ij^\ell =a_ji^\ell \\) is the symmetric adjacency matrix that describes the connectivity (\(a_ij^\ell \in \0,1\\)),

*σ*

^{ℓ}is the overall coupling strength, and

*g*^{ℓ}is the coupling function. The adaptive nature of the coupling is taken into account by the

*N*-dimensional matrix \(B^\ell =\b_ij^\ell (t)\\) (asymmetric, in general), whose entries evolve dynamically according to the adaptive rule Eq. (1b). The parameter

*ϵ*

^{ℓ}> 0 represents the inverse of the adaptation time scale. The evolution of the adaptive coupling terms \(b_ij^\ell (t)\) is controlled by a nonlinear function

*H*

^{ℓ}, which depends on the state of the nodes

*i*and

*j*. This adaptation rule can cover a wide class of plasticity functions, including Hebbian plasticity

^{34}, Wang-Rinzel plasticity

^{35}, as well as other adaptive schemes

^{26,36,37}.

This paper deals with the emergence of CS in the general system (1). We consider a partition of the set \({\mathcalV}\) of the network nodes into clusters \(\mathcalC_1,\mathcalC_2,\ldots ,\mathcalC_Q\), \(\cup _p = 1^Q\mathcalC_p={{{{{{{\mathcalV}}}}}}}\), \(\mathcalC_p\cap \mathcalC_q=\emptyset\) for *p* ≠ *q*. In particular, we are interested in the existence and stability of CS patterns, where all the elements within each cluster have an identical time evolution, i.e., *x*_{i}(*t*) = *x*_{j}(*t*) if *i* and *j* belong to the same cluster. A first challenge is to identify possible patterns of cluster synchronization for the general system of equations (1). In the Supplementary Note 1, we show that any equitable cluster partition for the adjacency matrix *A* corresponds to a flow-invariant cluster synchronous solution for the system (1). As the considered network model can have different types of connections, a valid solution presenting CS with the minimum number *Q* of equitable clusters can be found by using a variation^{38,39} of Hasler and Belykh’s coloring algorithm^{40} (for more details see Supplementary Note 1).

A second challenge concerns the stability of CS solutions, which depends on the complete network dynamical evolution, governed by Eq. (1). Therefore, once the *Q* equitable clusters are found, the CS analysis is applied to a coarse-grained dynamical model (called quotient network) whose *Q* nodes correspond to the possible equitable clusters. In particular, we analyze the stability of the CS patterns by linearizing Eq. (1) about a state corresponding to synchronization among all the nodes within each cluster^{39}. By denoting the cluster synchronization state as *x*_{i}(*t*) = *s*_{p}(*t*), where the node *i* belongs to the cluster \(\mathcalC_p\), we can write the quotient network dynamics,

$$\dot\boldsymbols_p(t)={{{{{\boldsymbolf}}}}}_p(\boldsymbols_p(t))+{{{{{{{\boldsymbolS}}}}}}}\left(\boldsymbols_p(t),\mathop\sum _\ell =1^L\mathop\sum _q=1^Q\sigma ^\ell r_pq^\ell k_pq^\ell (t){{{{{{{{\boldsymbolg}}}}}}}}^\ell (\boldsymbols_q(t-\tau ))\right),$$

(2a)

$$\dotk_pq^\ell (t)=-\epsilon ^\ell \left(k_pq^\ell (t)+H^\ell ({{{{{\boldsymbols}}}}}_p(t),{{{{{{{{\boldsymbols}}}}}}}}_q(t))\right),$$

(2b)

where *f*_{p} is the vector field of each node belonging to the cluster \(\mathcalC_p\) (i.e., *f*_{p} = *F*_{i} for any *i* ∈ *C*_{p}), \(k_pq^\ell \) represents the weight of the connection of type *ℓ* between the clusters \(\mathcalC_p\) and \(\mathcalC_q\), and the *Q* − dimensional matrix \(R^\ell =\r_pq^\ell \\) is the quotient matrix of *ℓ*-th type, such that \(r_pq^\ell =\sum _{j\in {{{{{{{{\mathcalC}}}}}}}}_q}a_ij^\ell \) (\(i\in {{{{{{{{\mathcalC}}}}}}}}_p,\quad p,q=1,2,\ldots ,Q\)). Let *n*_{p} be the number of nodes in the cluster \({{{{{{{{\mathcalC}}}}}}}}_p\), therefore \(\sum\nolimits_p = 1^Qn_p=N\). In the following, we will tacitly assume that *p* depends on *i* and *q* depends on *j*, i.e., *p* = *p*(*i*) and *q* = *q*(*j*).

To analyze the local stability of the synchronous state, we introduce infinitesimal variations about *s*_{p}(*t*): *w*_{i}(*t*) = *x*_{i}(*t*) − *s*_{p}(*t*) and \(y_ij^\ell (t)=b_ij^\ell (t)-k_pq^\ell (t)\), where \(i\in {{{{{{{{\mathcalC}}}}}}}}_p\) and \(j\in {{{{{{{{\mathcalC}}}}}}}}_q\). It is important to emphasize that the stability problem involves *N* variables *w*_{i}(*t*) and up to *N*^{2} variables \(y_ij^\ell (t)\) for each *ℓ*, which for sufficiently large *N* makes the analysis soon infeasible. This motivated us to develop a technique to reduce the dimensionality of the stability problem, which we present in the following.

First, we find the variational equations in vector form (see Supplementary Note 2), in terms of two vectors (** W**(

*t*) and \(\boldsymboly_q^\ell (t)\)) and of the

*N*×

*N*diagonal matrix

*E*

_{p}, which is the cluster indicator matrix:

*E*

_{p}has entries

*E*

_{p,ii}= 1, if node \(i\in {{{{{{{{\mathcalC}}}}}}}}_p\), 0 otherwise, i.e., this matrix identifies all the nodes that belong to cluster \({{{{{{{{\mathcalC}}}}}}}}_p\). Next, we find a coordinate transformation that separates the perturbation modes in the stability analysis as much as possible, thus reducing the stability problem into sub-problems of the lowest dimension. To this end, we compute the transformation (based on a simultaneous block diagonalization, SBD) introduced in

^{41,42,43}, which has been previously applied to analyze the stability of complete synchronization

^{44}and CS

^{45}in non-adaptive networks. The canonical transformation matrix

^{46}(see Supplementary Note 3) \(T=\left(\beginarraycT_\parallel \\ T_\perp \\ \endarray\right)\) is the orthogonal matrix that simultaneously block-diagonalizes the matrices

*A*

^{1},

*A*

^{2},…,

*A*

^{L},

*E*

_{1},

*E*

_{2},…,

*E*

_{Q}into

*M*diagonal blocks, each one of size

*d*

_{α}(

*α*= 1, …,

*M*). The

*Q*×

*N*sub-matrix

*T*

_{∥}is associated with longitudinal perturbations to the synchronization manifold, which therefore characterize the nature of the synchronized dynamics (e.g., periodic, quasi-periodic, chaotic) but not the stability of the CS. This is taken into account by transverse perturbations associated with the (

*N*−

*Q*) ×

*N*sub-matrix

*T*

_{⊥}

^{47}.

To analyse the CS stability by following the MSF approach^{30}, we now consider the transformed coordinates ** η**(

*t*) = (

*T*⊗

*I*

_{m})

**(**

*W**t*) (where ⊗ denotes the Kronecker product) and \(\boldsymbol\xi _q^\ell (t)=T{{{{{{{\boldsymboly}}}}}}}_q^\ell (t)\). In particular, we focus on the perturbations along the transverse manifold,

*η*_{⊥}(

*t*) = (

*T*

_{⊥}⊗

*I*

_{m})

**(**

*W**t*) and \(\boldsymbol\xi _q\perp ^\ell (t)=T_\perp {{{{{{{{\boldsymboly}}}}}}}}_q^\ell (t)\).

By using the matrix *T* (see details in the Supplementary Note 2), for the transverse perturbations we finally obtain,

$$\dot\boldsymbol\eta _\perp (t) \, = \, \rho _1\left(\left\{{{{{{{{{\boldsymbols}}}}}}}}_p(t)\right\}\right){{{{{\boldsymbol\eta }}}}}_\perp (t)+\rho _2\left(\left\{{{{{{{{{\boldsymbols}}}}}}}}_p(t-\tau )\right\}\right){{{{{{{{\boldsymbol\eta }}}}}}}}_\perp (t-\tau )\\ +\mathop\sum _\ell =1^L\mathop\sum _q=1^Q\rho _3^\ell \left(\left\{{{{{{{{{\boldsymbols}}}}}}}}_p(t-\tau )\right\}\right){{{{\boldsymbol\xi }}}}_q\perp ^\ell (t)$$

(3a)

$${\dot{{{{{{{{\boldsymbol\xi }}}}}}}}}_q\perp ^\ell (t)=-\epsilon ^\ell \left\{{{{{{{{{\boldsymbol\xi }}}}}}}}_q\perp ^\ell (t)+\rho _4\left(\left\{{{{{{{{{\boldsymbols}}}}}}}}_p(t)\right\}\right){{{{{{{{\boldsymbol\eta }}}}}}}}_\perp (t)\right\},\quad \quad q=1,\ldots ,Q$$

(3b)

where the terms *ρ*_{1}, …, *ρ*_{4} are the time-varying matrices defined in Methods. The set \(\left\{{{{{{{{{\boldsymbols}}}}}}}}_p(\cdot )\right\}\) collects all the synchronous solutions (delayed or not) corresponding to the *Q* clusters.

Through Eq. (3b), each entry of the vector \({{{{{{{{\boldsymbol\xi }}}}}}}}_q\perp ^\ell (t)\) corresponds to one and only one entry of the vector *η*_{⊥}, *q* = 1, …, *Q*. Moreover, through Eq. (3a), the entries of the vector *η*_{⊥} are organized into *M* groups, each one independent of the others, and in general with a different number of elements. *M* depends on the network topology and the analyzed cluster pattern and ranges between 1 (each component of *η*_{⊥} depends on the other components) and *N* − *Q* (all the entries of the perturbation vector *η*_{⊥,i} are independent of each other). The larger *M*, the higher the degree of decoupling and therefore the dimensional reduction. If we call *D* the sum of the distinct *d*_{α} values, the number of equations to analyze is *N*_{red} = (*m* + *L**Q*)(*Q* + *D*).

We have therefore decoupled the stability problem into independent lower-dimensional equations. We remark that the time evolution of each entry of the vector \({{{{{{{{\boldsymbol\xi }}}}}}}}_q\perp ^\ell (t)\) is decoupled from all other entries. Some of these entries are multiplied by zero coefficients on the right-hand side of Eq. (3a) and so can be safely removed from the analysis. According to the MSF approach, the stability of the cluster solution can be studied in terms of the Transverse Lyapunov Exponents (TLEs)^{48,49}.

In summary, the dimensional reduction provided by the MSF approach can be quantified as follows. The original Eqs. (1) are *N*_{orig} = *m**N* + *L**N*^{2}, which is the same number of equations of the variational system, i.e., Eqs. (2) and (3). However, the Eqs. (3) have *M* diagonal blocks, corresponding to *M* independent sets of equations, which allows the study of the stability of clusters in terms of a set of decoupled equations^{39,50,51}. All the 1-size blocks can be parameterized according to the MSF approach so that only one equation must be simulated. For example, if the network has an all-to-all topology we have the maximum reduction since all blocks are scalar, i.e., *M* = *N* − *Q*, *D* = 1, and the number of equations to analyze is \(N_red=\left.(m+LQ)(Q+1)\right)\), where we recall that *m* is the dimension of the state vector of the *i*^{th} oscillator, *Q* is the minimum number of equitable clusters, and *L* is the number of different types of connections.

### Applications of the methodology to three case studies

To validate and illustrate the potentialities of the introduced methodology, this will be applied to case studies that differ in the adaptation rules and other topological and dynamical features.

The first example refers to neural systems and to the possibility of storing stimulation patterns in memory, where the memory items are represented by synchronized clusters displaying population bursts in the *β*–*γ* range^{52,53,54,55}. The corresponding network has the following characteristics: all-to-all initial topology, uniform delay, homogeneous nodes, Hebbian adaptation rule, and a single layer.

The second example models the emergence of a disease state in a volume element of tissue represented as a two-layer network of phase oscillators, whose interactions are slower than the oscillators’ dynamics and are modeled through an adaptive process. The healthy state corresponds to synchronized dynamics, and the presence of several frequency clusters denotes a pathological regime^{56,57}. The corresponding network contains two layers with a multiplex initial topology, no delay, heterogeneous nodes, and an adaptation rule that depends on a parameter related to age.

The third case study is a small (*N* = 10) synthetic network with Erdős-Rényi initial topology, no delay, homogeneous nodes, Hebbian adaptation rule, and a single layer.

The theoretical approach previously described allows us (i) to analyze the stability of a given CS solution through system (3) and (ii) to study the behavior of the corresponding quotient network through system (2), at a reduced computational cost.

#### Population coding: a plastic network of neural masses

We consider a network of *N* identical neural mass models^{58}, arranged in a single layer and interacting via plastic gap junctions, with homogenous synaptic delay transmission *τ*. In particular, each network node corresponds to a Wilson-Cowan model, which describes the dynamics of a neural mass made up of two coupled populations of neurons, one excitatory and one inhibitory^{59}. The network plasticity obeys the Hebbian rule^{60} and involves only excitatory-excitatory connections. The model parameters (see Methods) are set to values that correspond to the emergence of a stationary regime for the dynamics of an isolated node, i.e., the excitatory and inhibitory population firing rates (represented by the state variable *E* and *I*) approach a constant value^{61}. Our goal is to study the emergence of firing rate oscillations in biologically relevant frequency ranges induced by the stimulation of a sub-group of neural populations, as in the case of visual and auditory stimuli^{62,63}. Furthermore, oscillatory behavior (population bursts) promoted by the plastic properties of the network can provide a means to store memory items, somehow in analogy to the mechanisms reported for working memory^{64}.

In particular, we consider a network with a fixed, fully connected, initial topology. Since the all-to-all topology is compatible with any coloring of the network, i.e., with any CS pattern^{38,40}, the network per se admits any clustering configuration. We now aim to check if the network can encode and decode simple visual memory patterns. This, in turn, provides a unique test-bed to investigate which patterns are selected by the adaptation rule.

As a first test (see Fig. 1a) we stimulated *N*_{1} nodes of the network with a current pulse with normalized height 0.2 and 100 ms duration, mimicking an external stimulation. The other *N*_{2} nodes (with *N*_{1} + *N*_{2} = *N*) are not stimulated. We want to study the network capability of encoding/decoding information through coexisting CS solutions. Notice that the analysis depends on the ratio *N*_{1}/*N*; therefore, it is valid for any value of *N*, owing to the assumption of fully-connected initial topology. To this end, we start with the simplest case of a network with *Q* = 2 stable clusters – where *C*_{1} contains the *N*_{1} stimulated nodes and *C*_{2} contains the *N*_{2} non-stimulated nodes – to study if the nodes of the network reach a stable solution. We analyze the stability of the cluster pattern solution that emerges at regime (i.e., after the effects of the initial pulse are over), meaning that the stimulus has determined the onset of two stable clusters in the network. Once the clustering has been induced and the corresponding quotient network found, we can analyze the stability of the CS solution. The TLEs computed through Eq. (3b) are always negative (see Supplementary Note 6). Therefore, the emerging cluster solution is always stable and the nodes within each cluster are phase- and frequency-synchronized. As the nodes within each cluster are synchronized, the frequency of the collective oscillations can be computed by analyzing the behavior of the *Q* = 2 nodes of the quotient network, which is much simpler than the original network.

Figure 1b shows the color-coded asymptotic behavior of the network, for different values of the synaptic coupling strength *σ* and of the number *N*_{1} of stimulated nodes. In the dark blue region, both populations have the same low persistent activity, corresponding to a constant population firing (i.e., not corresponding to brain oscillations), as shown in Fig. 1c. In the blue region, the stimulated population produces stable brain oscillations, as shown in Fig. 1e, whereas the non-stimulated population remains in a low persistent regime. In the green region, the stimulated (non-stimulated) populations present high (low) persistent activity as shown in Fig. 1d, but again not corresponding to brain oscillations. The edges between these regions (dashed lines) have been computed semi-analytically (see Supplementary Note 5).

The oscillations in the light blue region are in the *β*–*γ* range, as illustrated in Fig. 1f, which shows the frequency of the collective oscillations at the three values of *σ* marked by colored vertical bars in Fig. 1b. Each curve in Fig. 1f has the same coloring as the corresponding bar. The curves are plotted versus the number *N*_{1} of stimulated nodes. We notice that the higher the synaptic coupling *σ*, the lower the number *N*_{1} of nodes that must be stimulated to observe brain oscillations. Moreover, the obtained range of frequencies is the same for any considered value of *σ*.

We remark that, in the absence of adaptation, this network (with the same parameters) would not be able to generate persistent brain oscillations (see Supplementary Note 6).

We have obtained proof-of-concept evidence that the considered network can code the visual stimulus into brain oscillations whose frequency depends on the global coupling *σ* and the stimulus itself. Next, we will check if the network can code the stimulus into CS solutions corresponding to *Q* = 3 clusters.

As a second test (shown in Fig. 2), we stimulated two groups of nodes for 100 ms with stimulations of different intensity: namely, *N*_{1} (*N*_{2}) nodes with a current pulse of higher (lower) normalized height. The other *N*_{3} nodes (with *N*_{1} + *N*_{2} + *N*_{3} = *N*) are not stimulated, as shown in Fig. 2a. Accordingly, to analyze the network behavior we split the nodes into *Q* = 3 clusters: *C*_{1} and *C*_{2} containing active nodes, whereas nodes in *C*_{3} are in a low persistent regime. As in the first test, the nodes within each cluster are always frequency- and phase-synchronized since the TLEs are always negative. Therefore, the activity of the whole network can be studied by analyzing the quotient network. We observe that stable oscillations in the quotient network can emerge for a suitable combination of *N*_{1} and *N*_{2} nodes. In these cases, the activity of the two stimulated clusters is phase-locked with a constant phase lag Δ (see Methods) and generates oscillations again in the range *β*–*γ*.

Figure 2 b and c show, for different values of *N*_{1} and *N*_{2}, the color-coded frequency *f* of the nodes in the clusters *C*_{1} and *C*_{2} and phase lag Δ between oscillations of nodes in *C*_{1} and in *C*_{2}. In the white regions, there are no oscillations. Moreover, the behavior along the thick diagonal corresponds to the case *N*_{3} = 0, i.e., to the case with only two populations, both stimulated. In the colored region, the studied oscillatory solutions are stable since all the TLEs are negative. Figure 2d and e show the time plots of the state variables *E*_{1} (blue) and *E*_{2} (orange) for the two sets of parameters corresponding to the red (Fig. 2d) and black (Fig. 2e) dots in Fig. 2b, c.

We observe a complete symmetry in the diagram reported in Fig. 2b by exchanging *N*_{1} with *N*_{2}, due to the homogeneity in the original network, and that *γ* oscillations emerge whenever *N*_{1} >> *N*_{2} (or equivalently for *N*_{2} >> *N*_{1}). By contrast, low frequencies in the *β*-range are observable when the composition of the two stimulated clusters is comparable.

Furthermore, we observe that the phase lag between the two clusters ranges in the interval [*π*/2, 3*π*/2]. We remark that the phase lag is *π* whenever *N*_{1} = *N*_{2} and deviates from *π* as far as ∣*N*_{1} − *N*_{2}∣ increases. Notice that, due to the symmetry of the network, the phase is symmetric with respect to the axis (*N*_{1} = *N*_{2}), i.e., Δ(*N*_{1}, *N*_{2}) = 2*π* − Δ(*N*_{2}, *N*_{1}).

This second test proves that the network can code two different stimuli into persistent brain oscillations with different frequency *f* and phase lag Δ. Therefore, we conjecture that larger numbers of stimuli can be coded by the network by generating different cluster solutions with larger numbers of clusters. Moreover, each CS solution is characterized by a specific frequency *f* and phase lag Δ, which can be used to decode the stimulus. To check the validity of this conjecture and thus verify that the plasticity rule allows storing multiple memory items in the network, we performed a third test.

The memory items we consider are black and white images (size 200 × 200) of polygons inscribed in the same circumference with different numbers of edges, as shown in Fig. 3a. Each polygon corresponds to a different item to memorize and a different network stimulation pattern. Each pixel of the image is univocally associated with a network node (i.e., we considered a network with *N* = 40000 nodes) where we stimulate the *N*_{1} (*N*_{2}) nodes corresponding to black (white) pixels with a higher (lower) current pulse for 100 ms. TLEs are always negative (see Supplementary Note 6) and we can analyze the network behavior by using the quotient network. In this case, we observe the emergence of phase-locked oscillations in the *β*–*γ* range with a constant phase-lag Δ among the two population clusters. Indeed the situation is analogous to the one reported in the previous experiment whenever *N*_{3} = 0. Figure 3b shows Δ for different values of *N*_{1}. In the interval *N*_{1}/*N* ∈ [0.4, 0.6] we observe an almost linear dependence of the phase-lag on the percentage of nodes in the cluster *C*_{1}. This dependence can be employed to encode the number of edges of each polygon. Indeed this is possible as shown in Fig. 3b, where the red circles and the number next to each of them show the phase coding of the number of edges. Figure 3c and d show the time plots of the state variables *E*_{1} (blue) and *E*_{2} (orange) for the two sets of parameters corresponding to polygons with 8 (Fig. 3c) and 4 (Fig. 3d) edges.

The results of this third test can represent an example of population coding, where the relevant quantity is the phase lag Δ among collective oscillations, somehow resembling the frequency-band coupling observed in the monkey visual cortex^{63}.

All the results reported here are obtained with the methodology previously introduced. In particular, for the third experiment, this presents the advantage that the number of equations that must be simulated can be reduced from *N*_{orig} ≈ 1.6 ⋅ 10^{9} to *N*_{red} = 12, as *m* = 2, *L* = 1, *Q* = 2, and *D* = 1. The diagonal structure of the terms *ρ*_{1}, *ρ*_{2}, *ρ*_{3} and *ρ*_{4} in Eq. (3b) (see Supplementary Note 7) implies that (i) the stability of each cluster is independent of the stability of any other clusters (i.e., the clusters are not one-way dependent and not even intertwined^{39}) and (ii) the dimension of the space we have to explore to study the stability of the synchronous solution can be strongly reduced. For this example, we notice that the term *ρ*_{2} in the variational Eq. (3a) is zero. However, the same equations are influenced by the quotient network state *s*_{p}, whose behavior is affected by *σ*, as shown in Eq. (2).

#### Propagation of tumor disease: a two-layer adaptive network model

As a second example, we analyze an adaptive two-layer network of phase oscillators with multiplex topology as a model for the emergence of pathological states induced by tumors or infections^{57}. In this simplified scheme, the healthy state corresponds to a fully frequency-synchronized regime of the network. By contrast, pathological nodes have a metabolism that is faster than healthy nodes, and the pathological state is characterized by the presence of multifrequency clusters: one with a normal metabolism, containing healthy nodes, and one or more with a faster metabolism, containing pathological nodes. The presence of multifrequency clusters leads to a decreased overall synchronization of the network. In particular, the two-layer (*L* = 2) heterogeneous network is composed of *N* nodes (of two kinds) and represents a volume element of tissue consisting of parenchyma (organ tissue), basal membrane, and stroma (immune layer). This network mimics the functional interactions between parenchyma (*N* cells in the first layer) and stroma (*N* cells in the immune layer), which lead to the propagation of tumor disease. The communication through cytokines (which mediate the interaction between the cells) is modeled by adaptive connections, as their timescale is slow compared to the cell metabolism^{65,66}.

Each node *i* of the network models a cell of both the organ tissue and the immune layer; therefore, it is described by a 2-dimensional state vector \({{{{{{{{\boldsymbolx}}}}}}}}_i=[\phi _i^(1),\phi _i^(2)]^T\) (see Methods), where \(\phi _i^(1)\) and \(\phi _i^(2)\) are phase variables related to the metabolic activity of the *i*-th cell of the organ tissue (parenchyma) and of the immune layer, respectively. The model equations are detailed in Methods.

Each parenchyma cell in isolation can be of two kinds, namely healthy or pathological, depending on the value of the natural binary frequency \(\omega _i^(1)\), which represents the metabolism velocity of an isolated cell. Accordingly, the network is non-homogeneous, because the nodes are of two different kinds. A higher (lower) natural frequency \(\omega _i^(1)\) is associated with pathological (healthy) cells, since a higher velocity corresponds to a faster metabolism in the cell, usually associated with a pathological state. However, once connected within the network, the state of a cell, either healthy or pathological, is influenced by its interactions with the other cells. Indeed, the state of the cell is related to the rate \(\dot\phi _i^(1)\) of cytokine production, which depends on both the natural frequency \(\omega _i^(1)\) and the adaptive coupling terms \(b_ij^\ell \). Once coupled, healthy and tumor cells differ by their metabolic activity, i.e., tumor cells are less energy-efficient and thus have a faster cellular metabolism, corresponding to a higher rate \(\dot\phi _i^(1)\) of cytokine production.

As the adaptive interaction mediated by cytokines can involve any pair of cells^{57}, the cells of the organ tissue are initially linked through a fully connected topology (\(a_ij^\ell =1\) in Eq. (1)), where the strength of the connections depends on the adaptation variables \(b_ij^\ell \).

The health state of the complete network depends on the rate values \(\\dot\phi _i^(1)\\)^{57} and it is monitored through the parameter \(\bar\nu ^(1)\), which provides an average measure (over different trials) of the standard deviation of the angular frequencies \(\dot\phi _i^(1)\) (see Methods). A zero value of \(\bar\nu ^(1)\) means that all the nodes are locked in frequency; in this case, the network is healthy. On the contrary, whenever \(\bar\nu ^(1)\) is non-zero, frequency clusters appear in the network, meaning that some nodes evolve with a higher frequency, corresponding to a pathological state. An equivalent parameter \(\bar\nu ^(2)\) measures the standard deviation of the angular frequencies in the second layer, but its value is not related to the onset of a pathological condition because in tumor disease mutated cells are almost always located in the parenchymal layer^{57}.

The parameter *β* governs the plasticity rule of the cytokines and it can be associated with different adaptation rules. For instance, for *β* = *π*/2, a symmetric rule is obtained where the coupling increases between any two oscillators with close-by phases^{67}; on the contrary, if *β* = 0, the coupling will be strengthened if the oscillators have a phase shift of *π*/2. *β* mimics a systemic sum parameter which may account for different influences such as physiological changes due to age in the extracellular matrix, inflammation, systemic and local inflammatory baseline, etc. For tumor disease, it can include the malignancy grade of tumor cells. For the sake of brevity, we call this parameter the age parameter.

Using the MSF approach, we investigate how a small fraction *r* < 0.1 of initially pathological nodes and the age parameter *β* influence the onset of tumor disease in a network composed of *N* = 200 nodes. In particular, we analyzed the network behavior for a grid of 11 × 11 values of the parameters *r* ∈ [0, 0.1] and *β* ∈ [0.40*π*, 0.55*π*].

Due to the multiplex topology of the network, many different cluster patterns are admitted (see Supplementary Note 8). Nodes within a cluster are phase-synchronized (and therefore frequency-synchronized) and represent a group of cells in the tissue that produce cytokine with the same metabolism velocity. Depending on the evolution of the quotient network, clusters can evolve either with the same frequency (which means that all the cells in the tissue produce cytokine at the same rate and therefore the tissue is healthy) or with different frequencies (which means that some cells in the tissue produce cytokine with a faster metabolism and therefore the tissue is pathological). Therefore, to investigate if a stable cluster pattern is healthy/pathological, we evaluate if the parameter \(\bar\nu ^(1)\) in the quotient network is close to zero. In this paper, we perform two tests. The first test is to analyze the presence of multi-stability by varying *β* and *r*. In the second test, we investigate how different cluster patterns influence the onset of the pathological state.

In the first test, we focus on the simplest pattern shown in Figs. 3 and 4 of ref. ^{57}, where the nodes of the network split into two groups of phase-synchronized nodes. To analyze the stability of this pattern for different values of *β* and *r*, we initially divide the network into two phase clusters, one with *r* ⋅ 200 pathological nodes and the other one with (1 − *r*) ⋅ 200 healthy nodes. All the nodes within a cluster have the same initial phase, i.e., \({{{{{{{{\boldsymbols}}}}}}}}_i=[\phi _i^(1),\phi _i^(2)]^T={{{{{{{{\boldsymbols}}}}}}}}_j\) if *i*, *j* ∈ *C*_{p}, where *s*_{i} is the state of which we are analyzing the stability. The initial phases of each cluster are randomly selected from a uniform distribution in the range [0, 2*π*], whereas the initial coupling weights \(\b_p,q^1\\) and \(\b_p,q^2\\) between clusters are randomly selected from uniform distributions defined in the ranges [0, 2] and [ − 1, 1], respectively. In the Supplementary Note 9, we show that this 2-cluster solution synchronized in phase (and thus in frequency) exists and is stable for any considered value of *β* and *r*. Moreover, the parameter \(\bar\nu ^(1)\) is always close to zero (see Supplementary Note 9). Therefore, we found that a stable healthy state is possible in the whole considered range of the parameters. This result can hardly be obtained via the analysis performed in^{57}, which was based on extensive network simulations with different random initial conditions, and therefore led to a focus only on the pattern with the largest basin of attraction for each (*β*,*r*) pair. Our analysis suggests that there should be a coexistence among different healthy and pathological states in the whole bi-dimensional parameter plane (*β*,*r*). Therefore, a deeper analysis of the multiple stable cluster patterns admitted by the network is required.

As a second test, we extended the result obtained in^{57} by analyzing the network multi-stability, to see (i) how the fragmentation of the network in multiple clusters affects the overall health of the network, in terms of the parameter \(\bar\nu ^(1)\), and (ii) how many healthy and pathological states can coexist. To this aim, we divided initially the (1 − *r*) ⋅ 200 healthy nodes into *N*_{a} phase clusters, each one containing (1 − *r*) ⋅ 200/*N*_{a} nodes. The analyzed network structure is shown in Fig. 4a: 200 ⋅ *r* nodes (corresponding to the red node of the quotient network in Fig. 4a) are set as pathological nodes (\(\omega _i^(1)\) = 1), whereas 200 ⋅ (1 − *r*) nodes (corresponding to the green nodes of the quotient network in Fig. 4a) are set as healthy nodes (\(\omega _i^(1)\) = 0) and split into *N*_{a} phase clusters. Each cluster corresponds to a different initial phase, selected randomly, as in the first test. To study the effect of clustering on the network evolution, we considered three values of *N*_{a}, namely 6, 10, and 15, meaning that we have initially 6 (or 10 or 15) clusters of healthy nodes and we analyzed whether these clusters and the pathological nodes synchronize in frequency thus either yielding a healthy state or not.

For each parameter set, we analyze the quotient network and compute the average indicator \(\bar\nu ^(1)\) (see Methods) to discriminate between pathological (\(\bar\nu ^(1) > 0\)) and healthy (\(\bar\nu ^(1)=0\)) conditions. We take an average over \(\mathcalN\) trials with different random selections of the initial conditions for each cluster to make this indicator robust. The results are shown in Fig. 4b. For each value of *N*_{a}, a colored line marks the boundary between healthy (leftmost areas) and pathological (rightmost areas) states of the network. This implies that, for a given parameter set (*r* and *β* fixed), we can reach a healthy or pathological state depending on the initial conditions, i.e., on *N*_{a}.

The variability of tumor disease results from the initial genetic state of the tumor cells, their subsequent mutations, epithelial-mesenchymal transitions, and interactions with the innate immune system^{68}. Therefore, starting from different initial conditions allows monitoring of different possible evolutions. The model predicts that the tissue can fall into a pathological state with a high number of tumor cells when the fraction *r* of initially pathological cells or the age of the patient (*β*) increases. Moreover, the number *N*_{a} of healthy clusters influences the onset of the tumor; the higher *N*_{a}, the wider the parameter region that corresponds to the pathological state. Therefore, in the presence of smaller groups of initially healthy nodes synchronized at different phases, the network will become pathological more easily. By contrast, a single larger phase cluster of initially healthy nodes is more robust to the propagation of the tumor.

Figure 4c1 and c2 show \(\bar\nu ^(\ell )\) versus *β* for *r* = 0.07 and different values of *N*_{a}. As already shown in Fig. 4b, when *N*_{a} increases, the pathological condition (\(\bar\nu ^(1) > 0\)) is obtained in a wider range of *β*. At the same time, the value of the indicator \(\bar\nu ^(2)\) concerning the immune layer seems not to be influenced by *N*_{a} and reveals that this layer is usually not frequency-locked for *β* > 0.438*π*. Moreover, for the same set of parameters (*β*,*r*) the network exhibits different behaviors, depending on the chosen value of *N*_{a}, i.e., on the initial conditions, thus clearly showing multi-stability. Once again, this shows that knowledge of *β* and *r* is not sufficient to determine the onset of a pathological state and that the initial conditions of the network should also be taken into account.

The network dynamics is more deeply analyzed in some specific cases, by considering *r* = 0.07 and *β* = 0.55*π* for *N*_{a} = 6 (Fig. 4d) and *N*_{a} = 15 (Fig. 4e). In particular, the raster plots (see Methods) reported in the leftmost panels show the time evolution of the network. Notice that *N*_{a} = 6 corresponds to green nodes starting from 6 sets of different initial conditions (Fig. 4d1 and d2), whereas *N*_{a} = 15 corresponds to green nodes starting from 15 different sets of close initial conditions (Fig. 4e1 and e2), as better detailed in Methods. In Fig. 4d3 and d4, this leads to synchronization between green and red nodes and therefore to a healthy state. Indeed, although the pathological nodes in the parenchymal layer have different natural velocity \(\omega _i^(1)\), they synchronize their metabolic rate \(\dot\phi _i^(1)\) with the healthy nodes due to the coupling. The same information can be extracted from Fig. 4d3, which shows \(\dot\phi _i^(l)\) in the parenchymal layer. We remark that the metabolic rates in the parenchymal and immune layers are different due to the different coupling terms at regime, as pointed out in Fig. 4d5 and d6, which show the coupling matrix *B*^{1} and *B*^{2} once the asymptotic state is reached.

When the number of different sets of close initial conditions increases, multi-frequency clusters appear and the network falls into a pathological state (Fig. 4e1 and e2), where the pathological (red) nodes in the parenchyma layer evolve with a higher metabolic rate (Fig. 4e3). A corresponding cluster of lower frequencies is observed in the immune layer (Fig. 4e4). This indicates an essential change in the dynamical state of the immune layer, where all nodes have a velocity close to their velocity when isolated \(\omega _i^(2)=0\) (Fig. 4e4). Figure 4e5 and e6 confirm that the asymptotic coupling is weaker than in Fig. 4d5 and d6.

To check the robustness to heterogeneity (which is a common feature in natural networks) of the results obtained in this case study, we selected one parameter of each node randomly (similarly to what was done in refs. ^{69,70}), thus considering a fully heterogeneous network and preventing the possibility of perfect CS. In particular, the natural frequency of the isolated nodes in the first layer is selected as \(\omega _i^(1)=| \mathcalZ(0,\Sigma )|\) for healthy nodes and \(\omega _i^(1)=1+\mathcalZ(0,\Sigma )\) for pathological nodes, where \({{{{\mathcalZ}}}}(0,\Sigma )\) is a random variable with mean 0 and standard deviation Σ. Therefore, each node is associated with a different value of the parameter \(\omega _i^(1)\) and no exact CS is admitted by the network for any Σ ≠ 0. This implies that to analyze the network we have to simulate its whole dynamics and we cannot employ the quotient network. The obtained results for \(\Sigma \in \left\10^-3,10^-2,10^-1\right\\) are compared to the case Σ = 0 (solid curves in Fig. 5). For each value of Σ, the network behavior is analyzed for \({{{{{\mathcalN}}}}}=30\) different realizations of \({{{{{{{\mathcalZ}}}}}}}(0,\Sigma )\). The mean parameter \(\bar\nu ^(1)\) averaged across the \({{{{{{{\mathcalN}}}}}}}\) trials is shown in Fig. 5 for *N*_{a} = 6 (Fig. 5a), *N*_{a} = 10 (Fig. 5b), *N*_{a} = 15 (Fig. 5c).

We remark that the specific value \(\bar\nu ^(1)\) is relatively relevant, what is important is the threshold value of *β* from which \(\bar\nu ^(1)\ne 0\), corresponding to the onset of a pathological state. With this caveat in mind, for all the considered values of *N*_{a}, the results obtained with the proposed method remain valid for Σ≤10^{−2}. The results shown in Figs. 4b and c have been obtained through the analysis method introduced in this paper, which reduced the computation times by two orders of magnitude with respect to simulations of the whole network required to obtain Fig. 5. Indeed, the number of equations that must be simulated is reduced from *N*_{orig} = 80400 to *N*_{red} = 128 for *N*_{a} = 6, as *m* = 2, *L* = 2, *Q* = 7, and *D* = 1. Moreover, the reduction of the dimensionality of the state space obtained by using the quotient network enables an accurate analysis of the multi-stable solutions displayed by the model.

### Network of Lorenz systems with non-global coupling

To further investigate the potentiality of the method we analyze a synthetic network with non-global coupling. In particular, we consider a small symmetric network with one layer, no delay, *N* = 10 nodes, whose initial topology is of Erdős-Rényi kind (i.e., dense, but not global), with an edge removal probability of 0.1. This is a single-layer network, therefore we omit the index *ℓ*. Each oscillator is a Lorenz system, with coupling in the first variable (see Methods). When isolated, each node evolves toward a chaotic attractor. The adaptation law is the classical Hebb’s rule. By using the coloring method proposed in refs. ^{38,39}, we can split the network into *Q* = 3 clusters: nodes from 1 to 6 belong to cluster \({{{{{{{{\mathcalC}}}}}}}}_1\) (see Fig. 6a, dark blue nodes), nodes 7 and 8 belong to cluster \({{{{{{{{\mathcalC}}}}}}}}_2\) (light blue nodes), and the other (yellow) nodes belong to cluster \({{{{{{{{\mathcalC}}}}}}}}_3\). The second and third clusters are intertwined because their perturbations are associated with the same block of matrix \(\hatB_\perp \) (see Supplementary Note 10) and therefore share the same TLEs.

Figure 6b and c show the TLEs for a regular grid of 100 × 100 values of *σ* and *ϵ* for cluster \({{{{{{{{\mathcalC}}}}}}}}_1\) and \({{{{{{{{\mathcalC}}}}}}}}_2\) (the same result holds for cluster \({{{{{{{{\mathcalC}}}}}}}}_3\), due to the intertwining mentioned above), respectively.

TLEs lower than zero correspond to parameter sets where all the nodes within the cluster \({{{{{{{{\mathcalC}}}}}}}}_i\) are synchronized. A similar analysis can be carried out by simulating the whole network and by computing the synchronization error *E* (see Methods) between the nodes of each cluster. Figures 6d-f show the corresponding results. As expected, the two methods manifest an excellent agreement, but the proposed MSF-based method allows for reducing the computational time of one order of magnitude. Indeed, the entries of the matrices *ρ*_{1}, …, *ρ*_{4} are organized into *M* = 6 diagonal blocks (see Supplementary Note 10) and the number of equations that must be simulated is reduced from *N*_{orig} = 130 to *N*_{red} = 36, as *m* = 3, *L* = 1, *Q* = 3, and *D* = 3.

link