Synchronization of hypernetworks of coupled dynamical systems

We consider synchronization of coupled dynamical systems when different types of interactions are simultaneously present. We assume that a set of dynamical systems are coupled through the connections of two or more distinct networks (each of which corresponds to a distinct type of interaction), and we refer to such a system as a hypernetwork. Applications include neural networks formed of both electrical gap junctions and chemical synapses, the coordinated motion of shoals of fishes communicating through both vision and flow sensing, and hypernetworks of coupled chaotic oscillators. We first analyze the case of a hypernetwork formed of $m=2$ networks. We look for necessary and sufficient conditions for synchronization. We attempt at reducing the linear stability problem in a master stability function form, i.e., at decoupling the effects of the coupling functions from the structure of the networks. Unfortunately, we are unable to obtain a reduction in a master stability function form for the general case. However, we show that such a reduction is possible in three cases of interest: (i) the Laplacian matrices associated with the two networks commute; (ii) one of the two networks is unweighted and fully connected; (iii) one of the two networks is such that the coupling strength from node $i$ to node $j$ is a function of $j$ but not of $i$. Furthermore, we define a class of networks such that if either one of the two coupling networks belongs to this class, the reduction can be obtained independently of the other network. As an example of interest, we study synchronization of a neural hypernetwork for which the connections can be either chemical synapses or electrical gap junctions. We propose a generalization of our stability results to the case of hypernetworks formed of $m\geq 2$ networks.


I. INTRODUCTION
Synchronization of coupled dynamical systems has been the subject of a considerable amount of research (see e.g., [1][2][3][4][5]) with applications ranging from adaptive synchronization strategies [6][7][8][9][10][11] to pinning control [12][13][14][15]. One case of interest is that of complete synchronization that occurs when the individual systems, if appropriately coupled, converge on the same time-evolution. Complete synchronization can be observed in the presence of selective coupling, i.e., the systems are coupled through the connections of a network. A common underlying assumption is that the interactions among the systems are all of the same type. For this case, it has been shown that stability of the synchronized state depends on the details of the underlying network topology.
In this framework, the master stability function (MSF) approach [2] to synchronization of networks of coupled identical dynamical systems has been widely investigated in the literature [16][17][18][19]. An outstanding problem is how to obtain a reduction of the stability problem in a MSF form when the set of coupled dynamical systems simultaneously interact through different networks, with each network being associated with a distinct coupling function.
In this paper, we will focus on complete synchronization and we will retain selective coupling but we will allow for different types of couplings between the systems. We assume that all the connections that correspond to the same type of coupling form a network and the systems are connected by more than one network. This case is relevant to any situation where the individual units are allowed to interact through different types of coupling. For example, neurons in the brain are connected through both electrical gap junctions and chemical synapses, see e.g., [20,21].
The coordinated motion of shoals of fishes depends on the sensory capabilities of each individual fish. Fishes typically use vision but also chemical/flow sensing in order to localize their mates and coordinate their individual motion with respect to the shoal [22,23] (as in other animal species, the number of neighbors that can be simultaneously sensed by each fish is typically bounded and depends on the specific kind of interaction [24]). Another example is that of interdependent networks, such as e.g., the coupled infrastructure of power stations and internet communication servers [25]. In recent years, the possibility of cascades of faults through coupled interdependent networks has been pointed out as a crucial aspect with respect to the assessment and design of critical infrastructures [26].
In this paper, we consider that a set of identical dynamical systemsẋ i = F (x i (t)), i = 1, 2, ..., N , are coupled through the connections of m different networks, and we refer to such a system as a hypernetwork, see e.g., [27][28][29], [54]. We first consider the case of m = 2 networks (a generalization to the case of m ≥ 2 networks will be presented in Sec. IV). The systems are then coupled as follows, (1) i = 1, 2, ..., N , where x i (t) = [x 1 i (t), x 2 i (t), ..., x n i (t)] T is the n-dimensional state of node i, F : R n → R n represents the dynamics of each individual unit, G : R n → R n and H : R n → R n are different coupling functions, τ g and τ h are (possibly) different interaction delays, σ A and σ B are two scalar coefficients. As can be seen from (1) An equivalent way of writing Eqs. (1) is the following, Note that if a solution belongs to I over a time interval [t 0 , t 0 + τ max ], where τ max = max τg ,τ h , then the solution will belong to I, for any time t > t 0 + τ max . In this case, the synchronized solutions x 1 (t) = x 2 (t) = ... = x N (t) = x s (t) is characterized by the same dynamics as that of an uncoupled system, The main goal of this paper is to study linear stability of the synchronous solution (3,4) for the set of equations (2). The same problem for the case that the systems are coupled through the connections of only one network, i.e., L B ij = 0 in Eq. (2) has been intensively studied in the literature, see e.g., [2,18,[30][31][32][33][34][35][36]. For this case it can be shown that linear stability of the synchronous solution can be analyzed in terms of the following low-dimensional equation where DF (DH) represents the Jacobian matrix of the function F (H). In particular, the condition for stability is that the maximum Lyapunov exponents [55] associated with Eq. (5) are negative for k = 1, ..., (N − 1). Eq. (5) for which corresponds to the linearized equation for the evolution in the synchronization manifold (3). (5) is a system of n scalar differential equations as opposed to the linearized system (2), which is described by nN scalar differential equations. Hence, system (5) is termed low-dimensional. The nice thing about this approach is that it provides necessary and sufficient conditions for synchronization. Similar conditions have been obtained for networks of groups [17], for adaptive synchronization of complex networks [37,38], for the pinning control problem applied to a complex network [39,40], and for the case that slight deviations from nominal conditions are present [19,41,42]. In this paper, we attempt at obtaining a condition in terms of a low-dimensional equation for the more complex case that the systems are coupled through the connections of two different networks (Eq. (2)). However, as we will see, our proposed problem is not easy to solve in general.
In what follows, we first consider the case that the two matrices A and B in (1) are arbitrary and we show that the stability problem does not admit a solution in a low-dimensional form. Then we focus on three examples of interest for which we show that such a reduction is possible: • The two Laplacian matrices L A and L B commute.
• One of the networks (either A or B) is unweighted and fully connected.
• One of the two networks (say e.g., A) is such that A ij = a j , i, j = 1, ..., N .
The rest of the paper is organized as follows. In Sec. II we attempt at obtaining necessary and sufficient conditions for stability of the synchronous solution for a hypernetwork (2). Yet, we show that unfortunately it is not always possible to reduce the problem in a low-dimensional form. However, we analyze three cases of interest for which such a reduction is possible. Furthermore, we define a class of networks such that if one of the two coupling networks belongs to this class, the reduction can be obtained independently of the other network.
i = 1, 2, ..., N . The set of equations (7) can be rewritten in vectorial form as follows, where δx(t) = [δx 1 (t) T , δx 2 (t) T , ..., δx N (t) T ] T and the symbol ⊗ indicates the direct product or Kronecker product. Now we proceed under the assumption that at least one of the two Laplacian matrices, say L A , is diagonalizable, where Λ A is a diagonal matrix with the elements on the main diagonal being the eigenvalues λ A 1 , λ A 2 , ..., λ A N and V is a matrix whose columns are the associated eigenvectors, v 1 , v 2 , ..., v N . Then, by introducing the change of variable, η(t) = V −1 ⊗ I n δx(t), Eq. (8) becomes, where the matrix Ξ = V −1 L B V . It would be nice if the matrix Ξ were diagonal but unfortunately there is no guarantee that this will be the case in general. Then we see from Eq. (9) that, different from the classical master stability function derivation [2], it is not possible to decouple Eq. (9) in N blocks, each one independent of the others. A. The case that the two matrices L A and L B commute A special case is when the two matrices L A and L B commute. Two matrices that commute have the property of sharing the same set of eigenvectors, i.e., assuming that they are both independently diagonalizable, it is possible to where Λ B is a diagonal matrix with the elements on the main diagonal being the eigenvalues of the matrix L B . Thus for this case the matrix Ξ coincides with the diagonal matrix Λ B as It follows that equation (9) can be decomposed in N blocks independent of each other, which corresponds to perturbations in the direction tangent to the synchronization manifold (3) and as such are not relevant in determining stability of the synchronous solution. Thus a necessary and sufficient condition for synchronization is that the Lyapunov exponents associated with Eq. (10) are negative for k = 1, 2, ..., (N − 1).
We now introduce a parametric equatioṅ FIG. 2: A hypernetwork formed of a fully connected graph (thin black arrows) and a superimposed network of 9 directed links (thick gray arrows). All the links (those associated with either one of the networks) have associated weight equal one.
where y and z are two complex parameters. We associate a master stability function with Eq. (12), which returns the maximum Lyapunov exponent of Eq. (12) as a function of the pair of complex arguments (y, z).
Then given any hypernetwork (2), stability of the synchronous solution can be evaluated by checking that M(y, z) < 0, for (y, z) = (σ A λ A k , σ B λ B k ), k = 1, 2, ..., (N − 1). Alternatively, a necessary and sufficient condition for stability of the synchronized evolution is that the pairs (σ A λ A k , σ B λ B k ), k = 1, 2, ..., (N − 1) fall in the region of the domain of the master stability function M(y, z) for which M < 0. A similar result for the case of a single network whose topology is allowed to evolve in time has been previously obtained in [43].
However, we note that the case that the two matrices L A and L B commute is quite specific and not very likely to occur in practical situations. An example of two graphs with associated commuting Laplacian matrices is shown in Fig. 1. In Sections IIB and IIC, we present two examples for which a reduction of the stability problem (7) in a low-dimensional form is possible, even if the two matrices L A and L B do not commute.
B. The case that one of the two networks is unweighted and fully connected We consider the case that one of the two networks is unweighted and fully connected. Without loss of generality we take this matrix to be A, Then where δ ij is the Kronecker delta. An example of such a hypernetwork is shown in Fig. 2.
We consider again stability of the synchronous solution (3). In what follows, we obtain a master stability function by only diagonalizing the (N − 1) dimensional subspace of transverse perturbations without worrying about the fact that these may couple into the remaining direction (which is tangent to the synchronization manifold).
The matrix L A can be diagonalized as We now look at Eq. (9). It can be shown that the matrix In fact, the matrix L B V has a column whose elements are all zero. This is due to the properties (i) that the sum of the elements in each row of the matrix L B equals zero, and (ii) that the matrix V has a column (the same column where the eigenvalue 0 of Λ A is) whose elements are all the same. It immediately follows that V −1 L B V has a column whose elements are all zero. Therefore, Eq. (9) can be re-expressed as, where the vector We note that Eq. (16) is independent from Eq. (17). Hence, we term the first as the drive system and the second as the response system. Note that η ′ corresponds to perturbations transverse to the synchronization manifold, while η N corresponds to perturbations within the synchronization manifold. Thus synchronization stability is governed by Eq.
(16), which does not involve η N . We diagonalize the matrix Ξ ′ , obtaining (N − 1) blocks of the form, are the eigenvalues of the matrix Ξ ′ . Note that the eigenvalues of the matrix Ξ ′ are the same as those of the matrix L B , except for the one eigenvalue λ B N that is equal to 0.
If the (N − 1) maximum Lyapunov exponents associated with the drive system (18) are all negative, then for large enough t, ζ k (t) → 0, k = 1, ..., (N − 1). If this happens, then Eq. (17) yields for large enough t, which corresponds to the linearized equation in the direction tangent to the synchronization manifold.
Thus we can introduce the parametric equation (12) in the pair (y, z), with the parameter z being possibly complex and a master stability function (13) which returns the maximum Lyapunov exponent of Eq. (12) as a function of the parameters y and z. For a given hypernetwork (2,14), a necessary and sufficient condition for stability of the synchronous solution (3), is that y = −σ A N and z = σ B ν k , k = 1, 2, ..., (N − 1) belong to the region of the domain of the master stability function (13), for which M(y, z) < 0.
This formulation allows to decouple the effects of the dynamical function F and the coupling functions G and H, from those of the network matrices L A and L B . In particular, for any given triplet of functions F, G, and H, the matrix B determines the parameters ν 1 , ν 2 , ..., ν N −1 , and if the master stability function M(y, z) is negative for , then the synchronization manifold is stable. An interesting thing about our derivation (18) is that we have been able to obtain a reduction of the stability problem (7) in a low-dimensional form though the two matrices L A and L B do not necessarily commute.
C. The case that Aij = aj Here we consider the case that the coupling strength from node j to node i is only a function of the source node j and not of the destination node i, that is An example of such a network is shown on the left-hand-side of Fig. 3, where the width of each link j → i represents the strength of the associated coupling A ij . The network on the right-hand-side of Fig. 3 is an outward star graph, corresponding to setting a j = 0 for j = 1, ..., (N − 1) and a N = 0 in Eq. (20). Under assumption (20), Eqs. (1) become,ẋ i = 1, 2, .., N , which can be recast in the form of Eq. (2), with the matrix L A = {L A ij } having the form, (22) has the property that it has one eigenvalue λ A N = 0 with associated eigenvector [1, 1, ..., 1] and the remaining ( It is easy to see that the matrix Ξ is in the form (15), with the entries in the N -column being all equal to zero. This allows to decouple the set of linearized equations into a drive subsystem and a response subsystem, with the response subsystem corresponding to perturbations tangent to the synchronization manifold (3) and the drive subsystem corresponding to perturbations transverse to the synchronization manifold.
Then following the analysis in Sec. IIB, it can be shown that a necessary and sufficient condition for stability of the synchronous solution for the hypernetwork (21) is that the maximum Lyapunov exponent of the low-dimensional  we note that we have been able to obtain a reduction of the stability problem (7) in a master stability function form though the two matrices L A and L B do not necessarily commute.
We wish to emphasize that the case in Sec. IIB (fully connected network) can be seen as a subcase of that in Sec.
, it coincides with the matrix L A considered in Sec. IIB up to a multiplicative factor a.

D. Necessary conditions on the matrix A
We observe here that there is a substantial difference between the conditions on the adjacency matrices A and B (the Laplacian matrices L A and L B ) discussed in Sec. IIA and those discussed in Secs. IIB and IIC. First consider the case presented in Sec. IIA, that the two Laplacian matrices L A and L B commute; then, if one of the two matrices changes, there is no guarantee that the condition would still hold. On the other hand, the conditions discussed in Sec. IIB and IIC refer essentially to one of the two matrices, allowing the other one to be freely chosen.
In sections IIB and IIC, we have found sufficient conditions on one of the two adjacency matrices, say A, that if satisfied, allow a reduction of the stability problem in a low dimensional form, irrespective of the other adjacency matrix, say B. In this section, we are interested in finding necessary condition for this to happen. We consider the set of Eqs.
(1) and we define the class C of all the networks A that satisfy the property of allowing a reduction of the stability problem in a low dimensional form, irrespective of the other network B. In what follows, we show that a network in C is such that the entries of the associated adjacency matrix A = {A ij } satisfy A ij = a j , i.e., the same condition discussed in Sec. IIC.
Hereafter, we seek to find the conditions for an adjacency matrix A (a Laplacian matrix L A ) to be in C. Based on our previous discussion in Secs. IIB and IIC, we see that the properties that the matrix L A has to satisfy are the following: The sums of the elements in the rows of the matrix L A are equal zero. This also implies that the matrix L A has one eigenvalue equal zero, with associated eigenvector [1, 1, ..., 1].
(C) The remaining (N − 1) eigenvalues are all the same.
If the three properties above are satisfied, the matrix L A can always be written as follows, where the matrix P is a diagonal matrix with all the entries on the main diagonal being equal to the same value, say p, except one entry (which, without loss of generality, we assume to be the one in the rightmost column) that is equal to zero. The matrix W is any invertible matrix with the rightmost column being equal to the vector [1, 1, ..., 1]. We note that the matrix P can be rewritten as P = p(I N − I * N ), where I N is the identity matrix and I * N is a diagonal matrix with all the entries on the main diagonal being equal to zero except the one in the rightmost column being equal to one. It follows that It is easy to see that the matrix W I * N W −1 is by construction such that the entries in each one of its columns are the same. Hence, the corresponding adjacency matrices A have to be in the form A ij = a j , discussed in Sec. IIC.
We conclude that if we are given a specific adjacency matrix B (a specific Laplacian matrix L B ), there are two possible choices of the adjacency matrix A (the Laplacian matrix L A ) for which the stability problem can be reduced in a low-dimensional form: (i) L A commutes with L B , and (ii) A belongs to C, i.e., its entries are such that A ij = a j .
Note that condition (ii) is independent of the choice of the matrix L B .

III. EXAMPLES
A. Example 1: Coordinated motion of swarms of particles.
Swarms of birds, hordes of insects, shoals of fishes, and colonies of ants have been modeled as systems of interacting self-propelled particles [23,44,45]. Here we consider a simple model of N particles moving along a fixed direction, say y, through a resistent fluid. The position (velocity) of particle i along the y direction is labeled as y i (t) (v i (t)), i = 1, ..., N . We consider the following equations of motion, i = 1, ..., N . The first term on the right hand-side of Eq. (26b) represents propulsion/friction of particle i, v r i (t) is the relative velocity along y with respect to the resistent fluid of particle i, v r ) and v f (t) is the velocity of the resistent fluid, which we model as an external input and we assume to be uniform in space. The second term on the right hand-side of Eq. (26b) represents attraction from particle j on particle i. The third term on the right hand-side of Eq. (26b) models a relative velocity adjustment between particles. m j > 0 is the mass of particle j = 1, ..., N , α, β ≥ 0, c ij (t) measures the strength of the interaction from particle j on particle i, which we set to be a function of the physical distance between particles i and j, where d ij is the distance between particles i and j in the plane orthogonal to the y direction and the exponent e determines the strength of the interaction as a function of the distance. An analogous model for particles that are allowed to move in the three-dimensional space has been considered in [46].
We note that the system of equations (26) admits a synchronous solution y 1 (t) = y 2 (t) = ... = y N (t) = y s (t), where again v f (t) is an external input. The synchronous solution corresponds to a configuration in which all the positions and the velocities of the particles along the y-direction are the same. We are interested in studying stability of this solution. In order to do that, we linearize Eq. (26) about (28), i = 1, ..., N . Equations (29) can be rewritten in matrix form, where and A ij = m j , B ij = m j (d ij ) 2e , i, j = 1, ..., N . It is easy to see that the matrix A = {A ij } belongs to class C. Hence, following Sec. IIC, the stability problem can be reduced in a low-dimensional form analogous to Eq. (23),  In what follows, we consider a hypernetwork that allows a chaotic synchronous evolution (4). We choose n = 3, is the equation of the chaotic Rössler oscillator, G(x(t)) = [x 1 (t), 0, 0] T , and H(x(t)) = [0, x 2 (t), 0] T , τ g = τ h = 0.

Stability of the synchronous solution for networks of coupled Rössler oscillators coupled via different coupling functions
has been widely investigated in the literature, see e.g., [18]. While it is known that this problem allows a lowdimensional reduction, the case of dynamical hypernetworks has not been considered. In what follows, we show that under specific conditions, a low-dimensional analysis can still be applied and based on this approach, we derive new conditions for the stability of the synchronous solution. We numerically compute the master stability function M(y, z) associated with Eq. (12) for the case that y and z are real numbers. Fig. 4 shows the results of our computations with the gray (white) area indicating a negative (positive) master stability function. We wish to emphasize that once the master stability function has been computed for a given triplet F , G, and H (as shown in Fig. 4), we are able to predict stability of the synchronous solution for any hypernetwork (2), corresponding to either one of the three cases presented in Sections IIA, IIB, and IIC. We define the average synchronization error E,

. > t indicates a time average and
x s = (x 1s , x 2s , x 3s ) T denotes the dynamics of an uncoupled system (i.e., using dynamics from Eq. (4)).
We consider the hypernetwork shown in Fig. 2. We assume that the matrix A is associated with the fully connected network (whose connections are represented as thin black arrows in the figure) and that the matrix B is associated points corresponding to all the pairs (σ A λ A k , σ B λ B k ) to Fig. 4; if all the points fall into the grey area, this ensures stability (sufficient condition for synchronization) and if only one of the points fall into the white area, this corresponds to instability (necessary condition for synchronization). However, for the network of Fig. 2 and the master stability function of Fig. 4, we note that for any fixed value of y = −σ A N , the condition for stability is that where the parameter κ is the abscissa of the intersection of the y = −σ A N line with the right profile of the gray area shown in Fig. 4. Note that κ can be either positive or negative. We define λ B max = max (λ B 1 , λ B 2 , ..., λ B (N −1) ) and 1) ). Then for this case, stability of the synchronous solution can be assessed by testing the following simple condition, In Fig. 5 we consider the following three cases: σ A = 4.5/6, σ A = 5.5/6, and σ A = 2/3 (corresponding respectively to, y = −4.5, y = −5.5, and y = −4). As can be seen from Fig. 5(b), for the first two cases κ < 0, while for the latter κ > 0. We integrate Eqs. (2) and (33)

IV. GENERALIZATION TO m NETWORKS
In this section we consider synchronization of a hypernetwork formed of m ≥ 2 distinct networks. For this case, we rewrite Eq. (2) as follows,ẋ i = 1, 2, ..., N , where G k : R n → R n is the coupling function associated with the connections of network k, L k = {L k ij } is the Laplacian matrix associated with network k, σ k is a scalar measuring the overall coupling strength for network k, k = 1, ..., m. In what follows, we will generalize the main results of Sec. II to this more general case (Eq. (37)).
The delays τ k may be possibly different, i.e., τ i = τ j , i, j = 1, ..., m, i = j. The nN dimensional state space of the system described by Eqs. (37) contains the n-dimensional synchronization manifold I, defined by Eq. (3). Note that if a solution belongs to I over a time interval [t 0 , t 0 + τ max ], where τ max = max i τ i , then the solution will belong to I, for any time t > t 0 + τ max . In this case, the synchronized solutions x 1 (t) = x 2 (t) = ... = x N (t) = x s (t) is characterized by the same dynamics as that of an uncoupled system (4). In what follows, we are interested in evaluating stability of the synchronization manifold I.
As a first case, we consider that the matrices {L k }, k = 1, ..., m all commute with each other, i.e., they all share the same set of linearly independent eigenvectors. Then, similar to Sec. IIA, it can be shown that stability of the synchronous solution can be reduced in the following low-dimensional form, which returns the maximum Lyapunov exponent of the system (38) for y k = σ k λ k l . A necessary and sufficient condition for stability is that M(y 1 , y 2 , ..., y m ) < 0 for l = 1, ..., (N − 1).
We now attempt to generalize the result of Sec. IIC to a hypernetwork formed of m networks. We assume that the first (m − 1) networks, k = 1, ..., (m − 1), belong to C, while the remaining network, k = m, is arbitrary. Under these assumptions the first (m − 1) Laplacian networks are in the following form: The eigenvectors of any of these matrices can be used as a new basis, say we choose k = 1, Then it is easy to see that all the matrices V −1 Λ k V , for k = 2, ..., m, are in the form (15). It follows (similarly to Sec. IIC) that we can decouple the set of linearized equations in a drive subsystem and a response subsystem, with the response subsystem corresponding to perturbations tangent to the synchronization manifold (3) and the drive subsystem corresponding to perturbations transverse to the synchronization manifold. Moreover, it can be shown that a necessary and sufficient condition for stability of the synchronous solution for the hypernetwork (21) is that the maximum Lyapunov exponent of the low-dimensional equation, with each other, k = (m ′ + 1), ..., m. We observe that a reduction of the synchronization stability problem in a low dimensional form is possible,

V. STABILITY ANALYSIS FOR A MORE GENERAL CLASS OF HYPERNETWORKS
In this section, we consider hypernetworks of coupled systems, which cannot be cast into the specific form of Eqs.
(1). We will show that under appropriate conditions, the master stability reduction studied in Sec. II can be extended to study synchronization for this more general class of hypernetworks. In particular, we focus on synchronization of neuronal networks. Global synchronization of large areas of the brain is usually associated with the onset of a pathological condition, such as Parkinson's disease or epilepsy [47].
We study a hypernetwork of neurons coupled through both chemical synapses and electrical gap junctions. Such neuronal networks of different types connecting the same set of neurons have recently been explicitly discussed in the context of the C. Elegans nervous system, which has both a gap junctional network and a chemical synaptic network [48,49]. Following [20,21,50], a neuronal hypernetwork with these characteristics can be described by the following system of differential equations,ẋ where the n-dimensional vector x i (t) = [x 1 i (t), x 2 i (t), ..., x n i (t)] is the state of neuron i, with the first variable x 1 i (t) representing its membrane potential, F : R n → R n defines the dynamics of an uncoupled neuron, the coupling matrix A = {A ij } specifies the connection topology of the network of chemical synapses j → i, while the coupling matrix B = {B ij } specifies the connection topology of the network formed of electrical gap junctions j ↔ i, k A i = j A ij , k B i = j B ij , σ A and σ B are two scalar coefficients, E j is the synaptic reverse potential of neuron j. Note that the matrix A (B) is assumed to be asymmetrical (symmetrical). The n-matrix specifies the form of the coupling, indicating that neurons are coupled through their membrane potentials, the n-vector ς = [1, 0, ..., 0] T has a similar function, that is, selecting the first state variable x 1 i of the state-vector x i ; Γ ≡ ςς T .
The dynamical variables s ij (t) represent how strongly cell j is connected to cell i and obey the following differential equation [50],ṡ i, j = 1, ..., N , where τ is the interaction delay associated with synaptic coupling (due to axonal conduction and synaptic processes), c 1 , c 2 > 0 are two scalar coefficients, S : R → R is a sigmoidal function, which we set, where v −1 sl represents the slope of the function S when its argument is small and v th is the firing threshold. As can be seen from (44), the individual neurons may simultaneously interact through two distinct networks, i.e., the network A formed of chemical synapses and the network B formed of electrical gap junctions.
The condition for the set of equations (44) to admit a synchronous solution is that E 1 = E 2 = .... = E N = E s . If this condition is satisfied, the synchronous solution x s (t) obeys, Note that differently from the case considered in Secs. I, II, and III, the synchronous solution (47) does not obey the same equation as that of an isolated system. Our goal in this section is to study stability of the synchronous solution (47) for the hypernetwork (44). In order to do that, we linearize the set of equations (44) about (47), obtaining where the matrices By multiplying (49b) by A ′ ij and summing over j, we can rewrite (49) as We can now introduce the (n + 1)-vectors δx i (t) = [δx i (t) T , δs i (t)] T , i = 1, ..., N and the (N (n + 1))-vector δx(t) = [δx 1 (t) T , δx 2 (t) T , ..., δx N (t) T ] T . Then, following Sec. II, we can rewrite the set of equations (49) in vectorial form as follows, where the Laplacian matrix L B ′ = {L B ′ ij } = {B ′ ij − δ ij } and the (n + 1)-square matrices As can be seen, the structure of the linearized equations (51) is quite similar to that of Eq. (8) in Sec. II. The main difference with Eq. (8) is that in the case above, one of the two coupling matrices, namely A ′ , is not a Laplacian matrix, as the entries along each row of the matrix A ′ sum to one and not to zero. We now wonder whether the stability problem (51) can be reduced in a low-dimensional form. As for the case of Eq. (8), the main difficulty is that in general it is impossible to decouple Eq. (51) in N independent blocks. One possibility, which we do not give further consideration in what follows, is that the two matrices A ′ and L B ′ commute. Another possibility is that the matrix L B ′ belongs to class C. If this is the case, then the matrix L B ′ can be diagonalized as in Eq. (24), i.e., from which we see that similarly to Sec. IIC, the linearized problem (51) can be decoupled into a drive subsystem and a response subsystem, with the response subsystem corresponding to perturbations tangent to the synchronization manifold (47) and the drive subsystem corresponding to perturbations transverse to the synchronization manifold.
It is known from the literature that in the visual cortex [51] and in the posterior part of the putamen [52], small groups of neurons are likely to form dense and uniform clusters of gap-junctions. Hence, assuming that the network L B ′ is of class C can be appropriate to model such agglomerates of neurons or small subsets of them. Therefore, as an example, we consider a small group of N neurons connected by a dense L B ′ network of gap junctions, with L B ′ ∈ C.
Under these assumptions, by diagonalizing the (N − 1)-dimensional subspace of transverse perturbations (see Sec. II), the high-dimensional problem (49) can be reduced into the low-dimensional form, In the more general case in which L B ′ does not belong to class C and the two matrices L B ′ and A ′ do not commute, stability of the synchronous solution results in a much more complex problem, for which (49) cannot be reduced in a low dimensional form and we expect a higher degree of complexity. The study of this case is beyond the scope of this paper.
We run numerical simulations in which each individual system is described by the FitzHugh-Nagumo model, n = 2, and we set v th = 0.3, v sl = 10 −2 , c 1 = c 2 = 10, E s = 1, σ A = 1, τ = 1. In Fig. 6(a) we plot the time evolution of the synchronous evolution, obtained by integrating Eq. (48) for this particular choice of the function F in (57) and of the parameters. We further set σ B = 0.9 and calculate the maximum Lyapunov exponent associated with the low-dimensional system (56) as a function of the parameter λ A ′ . This corresponds to a master stability function (MSF), which is plotted in Fig. 6(b) for the case that its argument is real. As can be seen from Fig. 6(  We finally run simulations of the full nonlinear hypernetwork described by Eqs. (44,45,57). We set the initial conditions for x 1 i and x 2 i , i = 1, ..., N and for s ij , i, j = 1, ..., N to be random numbers drawn from a uniform distribution in the range (0, 0.2). We consider that the network of chemical synapses is the network of N = 6 nodes and 9 directed links represented in gray in Fig. 2, i.e., the entries of the matrix A are A ij = 1 if there is a gray direct arrow from node j to node i in the figure, A ij = 0 otherwise. The spectrum of the corresponding matrix A ′ is real and λ A ′ min = − √ 2/2 > −0.74. We set the network of chemical synapses to be such that B ij = b j = j, i, j = 1, ..., 6 (note that the particular choice of the values of b j , j = 1, ..., N affects neither the spectrum of the matrix L B ′ nor the low-dimensional equation (56)). We evolve the hypernetwork (44,45,57) from t = 0 to t = 500. We monitor the quantity E(t), defined in Eq. (34). As expected, we observe that after a transient, E(t) → 0. We repeat the same experiment for the case that A ij = 1 if |i − j| = 1 and A ij = 0 otherwise. For this case, the spectrum of the corresponding matrix A ′ is real but λ A ′ min = −1 < −0.74, thus predicting that the synchronous solution is unstable. This is confirmed by our numerical experiments, showing that, when the full nonlinear system (44,45,57) is integrated from initial conditions that are close to the synchronous state (48), E(t) does not converge to 0.

VI. CONCLUSION AND DISCUSSION
In this paper we have studied synchronization of coupled dynamical systems when different types of interactions are simultaneously present. Our study applies to any situation where the individual units interact through different coupling mechanisms. For example, neurons in the brain are known to be connected through both electrical gap junctions and chemical synapses, [20,21,48,49]. Also, our study encompasses a situation where different coupling functions correspond to different interaction delays.
In our formulation, a set of identical dynamical systems are coupled through the connections of two or more distinct networks (each of which corresponds to a distinct coupling function) and we refer to such a system as a dynamical hypernetwork. We first focus on the case of a hypernetwork formed of m = 2 networks and we seek to obtain necessary and sufficient conditions for synchronization. In Sec. II we try to reduce the stability problem in a master stability function form. Though a solution in this form seems to be not available in general, we show that such a reduction is possible in three cases of interest: (i) the Laplacian matrices associated with the two networks commute; (ii) one of the two networks is unweighted and fully connected; (iii) one of the two networks is such that the coupling strength from node j to node i is a function of j but not of i, with case (ii) being a subcase of (iii). We introduce a unique master stability function that determines stability for all three cases. Also, we define the class C of networks for which the reduction is always possible, independent of the structure of the other network.
We note that in many situations, such as, e.g., in biological networks, different types of interactions are typically present, but the couplings may vary in time due to changing environmental conditions, making satisfaction of either one of conditions (i), (ii), or (iii) difficult. On the one hand, this highlights a limitation of the master stability function approach that does not seem to be applicable to situations of arbitrary complexity (see also e.g., [17]). On the other hand, it poses the fascinating challenge of defining alternative tools to addressing stability for the case of arbitrary hypernetworks. We also point out here that we cannot exclude the existence of other conditions to be satisfied simultaneously by both matrices A and B (e.g., for either the hypernetwork (1) or (44)) that allow a reduction of the stability problem in a low-dimensional form.
In Sec. IV we have proposed a generalization of our stability results to hypernetworks formed of m networks. In Sec. V we have shown the possibility of generalizing our approach to hypernetworks of coupled systems, which cannot be cast into the specific form of Eqs. (1). As an example of interest, we have studied synchronization of a neural hypernetworks for which the connections can be either chemical synapses or electrical gap junctions. The results of this paper could be also easily extended to study synchronization of dynamical hypernetworks of coupled discrete-time systems.
Appendix: The special case of hypernetworks of N = 2 nodes In this appendix we show that for hypernetworks of N = 2 nodes, the stability problem can always be reduced in a low-dimensional form. We start by considering that N is an arbitrary number and that the hypernetwork is formed of m = 2 networks (Eq. (1)). The generalization to the case of m > 2 networks is straightforward.
We look at Eq. (2). In general, a case of interest is that one of the two Laplacian matrices, say L A , can be rewritten as, where the matrix L A1 belongs to C (i.e., it is in the form of the matrix (22)) and the matrix L A2 commute with L B , that is L A2 = V Λ A V −1 and L B = V Λ B V −1 , where Λ A and Λ B are diagonal matrices. Under the condition (58), Eq.
(8) can be rewritten as, Following Sec. IIC, it can be shown that a necessary and sufficient condition for stability of the synchronous solution for the hypernetwork (59) is that the maximum Lyapunov exponent of the low-dimensional equation, is negative for k = 1, ..., (N − 1), whereā = N j=1 a j , λ A k and λ B k are respectively the (complex) eigenvalues of the matrices L A and L B that are associated with the same eigenvectors, i.e., such that Note that the one eigenvalue λ A N = λ B N = 0 is not relevant in determining stability. Now the question arises how likely it is that an arbitrary Laplacian matrix L A can be decomposed in the form (58). In general terms, an N -squared matrix is determined by its N 2 entries. At the same time, we are allowed 2N degrees of freedom in the decomposition (58), i.e., N degrees of freedom in choosing the entries a 1 , a 2 , ..., a N of the C-matrix L A1 and N degrees of freedom in choosing the eigenvalues of the matrix L A2 . It follows that only in the case that N = 2, a decomposition in the form (58) is guaranteed irrespective of the choice of the two Laplacian matrices L A and L B . This leads to the conclusion that the stability of the synchronous solution for an arbitrary N = 2-hypernetwork can always be associated with the