Stabilization of synchronous equilibria in regular dynamical networks with delayed coupling

We consider the synchronization problem of dynamical networks with delayed interactions. More specifically, we focus on the stabilization of synchronous equilibria in regular networks where the degrees of all nodes are equal. By studying such control near a Hopf bifurcation, we obtain necessary and sufficient conditions for stabilization. It is shown that the stabilization domains in the parameter space reappear periodically with time-delay. We find that the frequency of reappearance of the control domains is linearly proportional to the number of cycle multipartitions of the network.


Introduction
Dynamical networks with complex topology are a very powerful approach for the study of large complex systems in various application areas ranging from neuroscience, engineering, to sociology, economics or Earth sciences. From this perspective, real-world systems are modeled as networks of interacting nodes, where the nodes have their own dynamics and influence each other's behavior in complex ways [1][2][3]. Research in this area combines various application fields with theoretical approaches from dynamical system theory, statistical physics, time series analysis, or graph theory. The networks can range from a few to hundreds of billions of nodes, such as large communication networks or the brain. Connectivity in networks can vary widely, ranging from "all-to-all" to hierarchical or even sparse.
Many natural systems possess unstable states, and the stabilization of them might be of interest for applications. For example, stabilization of unstable equilibrium solution has application in semiconductor lasers [37], medicine [38], and other fields. An efficient strategy to control unstable periodic orbits, which is also applied to stabilization of equilibrium, is introducing a self-feedback delayed term. This approach is noninvasive, since the feedback term does not alter the original solutions. In particular, for stabilization of periodic orbits, the usual approach is to consider the delay equal or close to the period of the orbit, known as Pyragas control [39]. For stabilization of equilibrium, the aim can also be achieved with a suitable choice of delay [40,41].
Here we consider a network model in which the nodes are interacting with a common time-delay. Hence, the delay is not an internal property of the isolated system but it is due to network coupling. Moreover, we study diffusive coupling, i.e., for a fixed node j, the coupling function depends on the difference of the state variables x (t − τ ) − x j (t), where node is adjacent to node j and τ > 0 is the time-delay. This means that the information that arrives to node j, coming from its neighbor , arrives with the delay τ . Using ideas from [40] and the Perron-Frobenius theory [42], we study in detail the control of unstable steady states close to a Hopf bifurcation, and we show that stabilization is possible for a suitable choice of τ and network coupling. The main contributions of this paper are: -Necessary and sufficient conditions for the stability of unstable steady states in regular networks with delayed interactions. -The effects of the network structure on the control is investigated. In particular, it is shown that the reappearance of the stabilization domains with τ depends on the network topology. Namely, the reappearance of the control domains is influenced by the number of spectral roots of the graph adjacency matrix, which in particular, can be associated to the number of cycle multi-partitions of the graph. -The number of maximal stability domains is analytically provided for the case of a symmetric coupling function. -Asymmetric coupling functions imply an inclination of the stability domains.
The structure and main results of the paper is as follows: In Sect. 2 we describe the considered model, which is a diffusively and delay-coupled dynamical network. Section 3 presents the main theoretical results: the bifurcation curves representing the boundaries of the stabilization domains, and an estimate of the number of the domains. In particular, we show how the number of domains per delay interval depends on the topological properties of the network. We also refer to these domains as islands of stability, similar to some previous publications. Sections 4 and 5 give illustrative examples with networks of Stuart-Landau and FitzHugh-Nagumo oscillators. In Sects. 6 and 7 we provide the technical proofs based on linearization and the master stability formalism. In Sect. 8 we discuss the properties of the spectrum of synchronous equilibrium for the case of a large delay. Such a case allows obtaining explicit approximations for the eigenvalues, which in turn provides a different perspective on the reappearing islands of stability. Finally, Sect. 9 contains concluding remarks.

Network model
Consider the following dynamical network model (1) where f, h : R n → R n are smooth functions with h(0) = 0, κ > 0 is the overall coupling strength and N is the number of nodes in the network. The adjacency matrix A = A j N j, =1 describes the interaction structure of the network with A j = 1 if there is a link from node to node j and A j = 0 otherwise, and τ > 0 represents the time-delay for the interaction between the nodes. In the following, we consider three assumptions in order to obtain a variational equation to the overall network model and use a master stability function reduction.
The first assumption guarantees that the equilibrium is close to an oscillatory threshold.
possesses an equilibrium x * , ( f (x * ) = 0), which has a pair of two complex conjugated eigenvalues. More specifically, the associated Jacobian matrix J F = ∂ f ∂ x (x * ) has a simple pair of eigenvalues α ± iβ, and the other eigenvalues have negative real parts. Without loss of generality, we consider x * = 0.
Note that the real part α can be positive, which refers to the case of an unstable equilibrium with a twodimensional local unstable manifold and the frequency β of unstable oscillations. Assumption 2.1 holds, for example, when the vector field f undergoes a Hopf bifurcation, which is often the case in applications. A simple paradigmatic model corresponding to such a situation is the well-studied Stuart-Landau oscillator [15,21,43].
By Assumption 2.1, the equilibrium x * = 0 of the uncoupled system (κ = 0) is unstable for α > 0. For α < 0, it is stable if all remaining eigenvalues are stable. Practically interesting is the case when the eigenvalues α ± iβ are critical, i.e., the remaining spectrum having the real parts smaller than α.
The next assumption allows us to apply the master stability formalism.

Assumption 2.2
The coupling function h : R n → R n is smooth, bounded, and h(0) = 0. Moreover, we assume that the Jacobians evaluated at the equilibrium In particular, the assumption 2.2 holds when J H and J F commute, or in an even more special case, when J H = cI with a scalar c and identity matrix I. For instance, the scalar sine function we use in the illustrative examples below satisfies the Assumption 2.2. The simultaneous diagonalizability of J H and J F means that they share the same eigenspaces. We denote by σ p (J H ) = δ + iη the eigenvalue of J H corresponding to the eigenspace from the critical eigenvalue α + iβ of J F .
The final assumption describes the class of the considered networks.

Assumption 2.3
The graph corresponding to the adjacency matrix A is d-regular and strongly connected.
Recall that a regular (or d-regular) graph is such a graph in which all vertices have the same degree d.
In the case of a directed graph we consider the inner degree. A graph is strongly connected when there exists a path between any pair of vertices [1].
With Assumption 2.3, one can use the Perron-Frobenius theory [42] to characterize the eigenvalues of the adjacency matrix A. In particular, if σ A is the spectral radius of A, then the eigenvalues {σ 1 , · · · , σ m } on the spectral radius of A are simple and given by The number m is called the index of imprimitivity (or period) and it represents the number of partitions of the graph. The index of imprimitivity m of A plays a central role in our result. It determines the frequency of reoccurrence of the stability regions with respect to the delay. Note that an undirected non-bipartite graph possesses the index of imprimitivity m = 1, while m = 2 for a bipartite undirected graph. Recall that a graph G = (V, E) is bipartite if the set of vertices V can be partitioned into two disjoint sets such that every edge connects vertices from a different set. For a directed graph, the index of imprimitivity m can admit larger values, and it denotes the number of multi-partitions. A directed graph is called cycle m-partite, m ≥ 2, if the set of vertices V admits a partition into m classes such that every class j couples only with the class j +1 where the modulus m is applied to the classes. See Fig. 1 for an illustration. The black curves are obtained from Eqs. (3) and (4). The colored region is the numerically computed stability region from Eq. (20) via the W -Lambert function and the color scale represents the magnitude of (λ) < 0. Other parameters: α = 0.01, β = π

Proposition 1 The stability regions lie between the following two neighboring bifurcation curves
and whenever they intersect. In other words, the fixed point solution of Eq. (1) is exponentially stable inside the regions delimited by the curves τ 1 ( j, κ) and τ 2 ( j, κ).
The proof of Proposition 1 is presented in Sect. 7.1. The Proposition 2 gives sufficient conditions for the existence of the intersection points, i.e., the solutions of τ 1 ( j, κ) = τ 2 ( j, κ). Its proof is given in Sect. 7.2.1.
For η = 0, e.g., when the linear part of the coupling function is symmetric, the corresponding interval is α 2dδ , +∞ and Remark 1 The stability islands are typically considered in the parameter plane (τ, κ) (see Fig. 2). A relation between the intersection points of the curves (3) and (4) is given by where θ 1 and θ 2 are the argument angles in equations (3) and (4), respectively, and θ 1 + θ 2 = arccos δ 2 − η 2 δ 2 + η 2 is constant. If η = 0 (symmetric linearized coupling function) then τ j = π(2 j + 1) βm is the same for the same island. And, the stability islands have an inclination whenever η = 0. Section 4.2 illustrates this behaviour. The stability domains are reappearing accordingly to the imprimitivity m of the graph. In other words, the time-delay interval covering a fixed number of stable islands is shortened by a factor m. Section 4.1 illustrates this property.

Stuart-Landau oscillator coupled in regular rings
A simple example of a d-regular, strongly connected, and cycle multi-partite graph is the unidirectional ring of N vertices, where each node j connects only with the node j + 1 (modulo N ). In this case d = 1 and the graph is cycle N -partite. Therefore its eigenvalues are exp(2π iu/N ) with u = 1, · · · , N . See Figure 1 (right panel) for an illustration. The non-directed ring graph with 5 nodes (left panel in Fig. 1) is an unpartitionable regular graph with d = 2 and m = 1. The directed ring in Fig. 1 (right panel) is a cycle 5-partite graph with d = 1 and m = 5. Figure 2 shows the significant difference between the stability regions corresponding to the two coupling structures from Fig. 2. In particular, the frequency of the reappearance of the control islands for the directed ring is 5 times larger than those for the non-directed graph, as stated in Remark 1.
Considering the measure of area of control, one can conclude that it is easier to control the equilibrium point of the SL system in a non-directed ring network rather than the directed one illustrated in Fig. 1, this is a typical property observed in different configurations [2]. The number of stability islands is given by Eq. (5)  This example can be extended further and we see that if we let the network size N grow, then its always possible to establish control of the unstable equilibrium point in a ring non-directed network, whilst in the directed ring network, the unstable equilibrium point remains unstable for any choice of κ and τ and sufficiently large N . This illustrations are in agreement with the results in [59].

Inclination of the stability islands for non-symmetric coupling
Here we analyze the same examples as in the previous Sect. 4.1, but with a linear coupling function given by the matrix The eigenvalues of H are δ ± iη = 1 ± 0.1i. Note that H commutes with J F and, therefore, Assumption 2.2 is satisfied. As stated in Remark 1 the stability islands suffer an inclination due to the non-symmetric coupling function. Also, the stability area for each island is shortened. The reappearing property as stated in Remark 1 is not affected. The corresponding stability diagrams for the new coupling function H are shown in Fig. 3 (compare with Fig. 2).

Example: the FitzHugh-Nagumo system
Next we the network model (1) with FitzHugh-Nagumo oscillators f : R 2 → R 2 defined as follows [60]: where we use the FitzHugh's original values φ = 0.08, a = 0.7, and b = 0.8. The FitzHugh-Nagumo system  (1) with the coupling topology as in Fig. 1 and the linear coupling function (7). Remaining parameters are the same as in Hence, the equilibrium (V , W ) is the real solution of and W = V + a b . As the parameter I changes, we observe numerically a stable equilibrium (−1.1994, −0.6243) for I = 0 and a stable periodic orbit for I = 1 close to this equilibrium. Hence, we expect a Hopf bifurcation in the interval I ∈ (0, 1). To find the parameter value for the Hopf bifurcation, we calculate the Jacobian matrix of (8), evaluated at the equilibrium point where V = V (real solution of Eq. (9)). The corresponding characteristic polynomial for the eigenvalues μ reads The necessary condition for the Hopf bifurcation means that both eigenvalues have the form μ = ±iω, ω ∈ R which leads to φb+V 2 −1 = 0 and φ(1−b+bV 2 ) > 0. The second condition is satisfied since 1−b = 0.2 > 0 and φ = 0.08 > 0. The first condition leads to V = ± √ 1 − φb. With this value for V we can solve Eq. (9) with respect to I , which leads to the Hopf bifurcation condition:

Two coupled FitzHugh-Nagumo oscillators
Here we consider the simplest case for the network model (1) 20), obtained numerically using the W -Lambert function stability islands is correctly predicted by Eq. (5). Note that the value of β was strongly reduced compared to the previous example. Since we consider periodic orbits near the Hopf-bifurcation, the value of β is related to the oscillation period and its direction. The smaller β is, the slower the rotational movement is, and therefore the number of stability islands reduces as β decreases. This is described in Eq. (5) and observed in Fig. 4.

Variational equation and master stability function reduction
This section together with Sect. 7 provides main theoretical analysis and proofs of Propositions 1 and 2. The linearization of Eq. (1) at the equilibrium solution x * readṡ where d j is the degree of the node j, J H is the Jacobian matrix of f at x * , and J H is the Jacobian matrix of the coupling function h at 0. (13) can be block-diagonalized and reduced to the following set of uncoupled linear equations: (14) where j = 1, · · · , N , ξ j : R → C, σ j are the eigenvalues of the coupling matrix A, and d is the degree of the network.
Equation (15) can be block-diagonalized, if the matrices A and D are simultaneously diagonalizable. A trivial case in which it happens is when the graph is d-regular, that is, all the nodes of the network have the same degree d. This implies D = d · I N . In this paper, we restrict our analysis to regular networks (see Assumption 2.3).
By the diagonalizability of matrix A we understand that there is a decomposition A = R A R −1 such that A is the diagonal matrix with the eigenvalues σ j , j = 1, · · · , N on its diagonal, and matrix R is invertible. Applying the change of coordinates ξ = (R −1 ⊗ I n )ψ, we obtain the block-diagonal systeṁ which can be written in as Eq. (14).

Analysis of the characteristic equation
The spectrum of eigenvalues for Eq. (14) determines the stability of the solution ξ j (t) = 0 and, therefore, the local stability of the equilibrium in Eq. (1). This section reduces the problem to a set of characteristic equations, namely, the ones in which the eigenvalues σ j are on the spectral radius of the adjacency matrix.
The spectrum of eigenvalues for (14) consist of complex numbers λ corresponding to the exponential solutions ξ j (t) = e λt v j . The corresponding eigenvalue problem for λ and v j reads The simultaneous diagonalizibility of J H and J F (Assumption 2.2) implies that Eq. (17) can be diagonalized leading to the set of n scalar characteristic equations in which σ p (J F ) and σ p (J H ), p = 1, · · · , n, are the eigenvalues of J F and J H , respectively, associated to the same eigenspace. Denoting σ p (J F ) = α + iβ and σ p (J H ) = δ + iη, Eq. (18) reads −λ+α +iβ −κd(δ +iη)+κσ j e −λτ (δ +iη) = 0. (19) Equation (19) is of type −λ + a + be −λτ = 0 with a, b ∈ C. The solution λ can be obtained via the W -Lambert function, which reads λ = a + (1/τ )W (bτ exp(−aτ )). Therefore, the solutions to Eq. (19) can be formally written using the W -Lambert function: This expression determines families of eigenvalues that depend on the index j as well as on the (infinitely many) branches of the Lambert function. Equation (20) is used to obtain, numerically, the colored points inside the stability regions given in the example sections. We are interested in those eigenvalues λ that have largest real part as they determine the stability. The following Lemma shows that, under certain conditions, one can restrict the analysis to the characteristic equation with such σ j that belong to the spectral radius (2).

Lemma 2 Under the condition
the real part of the solutions λ of the characteristic equation (19) is a monotone increasing function of |σ j | in the neighborhood of (λ) = 0. Thus, for the parameter values satisfying (21), the destabilization occurs with the increasing of |σ j |, and larger |σ j | from the spectral radius of A determine the most unstable roots of Eq. (19).
is an increasing function, for sufficiently large coupling weight κ or delay τ . As we will illustrate later, the condition (21) from Lemma 2 is not very restrictive and is fulfilled in a wide range of parameters of interest.
The arguments of each part in Equation (25) must coincide. Then we find with n ∈ Z.
We can rewrite 2π j m +2π n = 2π m ( j +n·m) = 2πj m withj ∈ Z. Therefore, isolating τ from Eq. (29) we obtain The boundary of the control domain in (τ, κ)parameter plane is given by two branches corresponding to two "neighboring" bifurcations (see more geometric details in Ref. [40]). As a result, we define the two neighboring bifurcation curves given by Eqs. (3) and (4) which define the the boundary of the control islands. The arg function is restricted to [−π, π].

Existence of stability islands
A natural question is when solutions of the equation τ 1 = τ 2 exists. Here, these solutions for the general case are found only numerically. However, when the linearized coupling function is symmetric then the existence of the solutions is guaranteed.

Proof of Lemma 2
The branches τ 1 and τ 2 meet at points κ min and κ max , which are solutions of the transcendental equation Consider the function r : κ min , β dη , defined as where κ min is defined in Eq. (27). The arccos part of this function is decreasing and limited in the interval (0, π]. The square root part times τ j is growing with κ and it has a vertical asymptote at κ = β/(dη). At ψ = 0 (that is, κ given by Eq. (27)), we have and r (κ min ) > 0 for all In particular, if η = 0 then r (κ min ) > 0 for all m > 1 which is always fulfilled. If r (κ) < 0 for any interval of κ, then r (κ) changes sign twice in its domain meaning that it has two roots. For the existence of the roots, it is sufficient that r (κ) < 0 withκ being the solution of r = 0. See Fig. 5 for a typical profile of the function r .
Due to the complexity of the function r (κ) in Eq. (32), computing its derivative is not straightforward. Instead, we note that this function is continuous with respect to η and calculate r (κ) setting η = 0, which is the case of symmetrical linearized coupling function. Getting the existence conditions for η = 0, the extension for small η > 0 is straightforward due to the continuity of the function r in Eq. (32).
The derivative of r η=0 (κ) is The minimum value is attained for r η=0 (κ) = 0 that leads to Substituting this value in Eq. (32) and reorganizing some terms we have the inequality Note that the domain of τ j lies in the interval [0, 2/α], and for small enough τ j it is always possible to validate inequality (34) and, consequently, to find solutions of Eq. (31), for η = 0.
Writing y j = 1 − ατ j in Eq. (34) it becomes Setting y j = cos(θ j ) then 1 − y 2 j = sin(θ j ) and the above inequality is now θ j + sin(θ j ) < π m Considering θ j = 0 (which happens only for y j = 1) then we have sin(θ j ) θ j < π mθ j − 1 which leads to θ j < π 2m . Undoing the nested substitutions we arrive to the condition for existence of solutions given by Eq. (5)

Spectrum for large delay
This section is devoted to the analysis of the spectrum, i.e., the zeros of Eq. (18) in the limit of large delay. In such a limit, one can obtain more explicit approximations for the eigenvalues, in contrast to the Lambert function representation by Eq. (20). Moreover, we will see how the index of imprimitivity of the network affects the eigenvalues.

Lemma 3 The solutions of the characteristic equation Eq. (18) can be estimated as
and Proof Using the scaling λ = γ /τ + iω in Eq. (19) we have The asymptotic continuous spectrum is computed for large τ by neglecting the term −γ /τ in Eq. (37) (for a rigorous justification, we refer to [62]). Taking absolute values, we solve Eq. (37) with respect to γ = γ a obtaining Eq. (36). Also, the arguments in Eq. (37) must coincide. Then Isolating ω we end up with Eq. (35). (37) is independent on m and j which means that asymptotic stability properties of different j modes are the same. The maximal value of γ a is attained for ω = β − κdη obtaining

Remark 3
The positions of the eigenvalues on the spectral lines depend on m. The vertical distance between two consecutive eigenvalues is In other words, the asymptotic spectral lines for different j = 1, ..., m coincide; on each curve, the eigenvalues are shifted by 2π/τ (which is standard for large delay [72]), but theyare shifted by 2π mτ with respect to each other for different j curves. Figure 6 illustrate the results for the example of Stuart-Landau oscillators. Figure 6 shows the following information:  (37)). Figure 6 clarifies the relationship between the frequency of recurrence of the stability islands and the graph structure. The fact that the distance between neighboring eigenvalues on the spectral curve is inversely proportional to m explains the proportionality of the frequency of recurrence of the bifurcation lines to m and, therefore, the control domains in the bifurcation diagrams from Figs. 2 and 3. With the increasing of τ , the eigenvalues move along the spectral curve and therefore, cross the imaginary axis more frequently.

Conclusions
In the context of regular networks with delayed interactions, we have proved that a stabilization of a synchronous equilibrium point is possible when the system is near the Hopf-bifurcation. Moreover, we have derived explicit expressions for the control domains (islands), showed that the control occurs periodically with time-delay, and that the frequency of this periodicity is characterized by the number of multi-partitions  37) and (36) for two coupled Stuart-Landau oscillators with identity coupling function (δ = 1 and η = 0) and with parameters α = 0.01, β = π , κ = 1 and delay τ = 7.5. The blue and red dots are the eigenvalues obtained for different j = 1, 2. The black curve is the asymptotic continuous spectrum given by the Eq. (36) of the graph. We have also approximated the maximal delay for which stabilization is observable and the maximum number of stable islands, again closely related to the graph structure.
This paper contributes to the major question of how the network structure effects its dynamical (functional) properties. In d-regular and strongly connected networks where all nodes have the same degree d, the index of imprimitivity m determines the number components that are circularly coupled to each other in a feed-forward manner. The number of the stability islands within the same delay interval in such networks is shown to be linearly proportional to m.
The intersection points of the curves that delineate the stability boundaries could not be determined analytically here, and perhaps a deeper analysis can be used in a future extension of this paper to overcome this point. We have only determined them numerically, and even this approach was challenging because these points are near the bottom of the curve where it becomes nearly vertical. We have solved this problem by setting a smaller step size in the code, but it must be adjusted accordingly when the parameters are changed.
Interestingly, we show the role of the asymmetry of the coupling function: for a symmetric function, the stability islands are "straight", while an asymmetry leads to an inclination of the islands with respect to the delay parameter. How much the asymmetric coupling functions affect the stability islands is not investigated here, and this point may be a subject for a future research.
The obtained results not only shed light on the properties of the synchronization islands in delay-coupled systems, but also open new perspectives. An important practical question is how a deviation of the graph topology from the d-regular structure affects the stability islands. For example, one could test whether a weak perturbation of the d-regular structure by a few link modifications would lead to a small perturbation of the stability islands.