Nonlinear integral coupling for synchronization in networks of nonlinear systems

This paper presents a novel approach to (controlled) synchronization of networked nonlinear systems. For classes of identical single-input–single-output nonlinear systems and networks, including oscillator networks, we propose a systematic design procedure (with generic as well as constructive conditions) for specifying nonlinear coupling functions that guarantee global asymptotic synchronization of the systems’ (oscillatory) states. The proposed coupling laws are in the form of a definite integral of a nonlinear ‘‘coupling gain’’ function. It can be fit to the system’s nonlinearities and, thus, can avoid cancelling nonlinearities by feedback or high-gain arguments commonly needed for linear (diffusive) coupling laws. As demonstrated by two examples, including a network of FitzHugh–Nagumo oscillators, this design can result in much lower synchronizing coupling gains than for the common case of linear couplings, therewith increasing energy efficiency of the coupling laws and reducing outputnoise sensitivity. The resulting coupling structure can be of a varying type, when couplings are activated/deactivated depending on the systems’ outputs without undermining overall synchronization. The approach is based on a novel notion of incremental feedback passivity with a nonlinear gain. In addition to the design contribution, these results provide a new insight into potential synchronization mechanisms in natural and artificial nonlinearly coupled systems. © 2022 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
Synchronization of oscillatory systems and, in particular, of chaotic systems is a phenomenon that received huge attention in scientific literature for the last 20 years (Strogatz, 2003). The coexistence of complex, chaotic or ''irregular'' dynamics of relatively simple systems on the one hand, and the possibility of some kind of ''order'' or synchrony in such interconnected systems, on the other hand, forms an intriguing combination for specialists in physics, mathematics, control, neuroscience and biology, thus generating a seemingly endless sequence of various results other, cf. Gambuzza and Frasca (2019), Liu and Chopra (2012), Stan and Sepulchre (2007), Stoorvogel, Saberi, and Zhang (2017), Tuna (2017) and Zhang, Trentelman, and Scherpen (2014). Although for some classes of nonlinear systems, e.g. Euler-Lagrange systems (Nijmeijer & Rodriguez-Angeles, 2003), specific structure of nonlinearities is exploited to design synchronizing linear couplings, for generic nonlinear systems, specific nonlinearities are commonly ignored, as synchronization by linear (diffusive) coupling can often only be guaranteed by suppressing them using high-gain designs (Pogromsky, 1998;Zhang et al., 2014).
Nonlinear couplings provide a greater degree of flexibility. They have been studied, in the scope of synchronization, in He and Yang (2008) and Liu and Chen (2008), in the context of mechanical oscillators (Ramirez, Olvera, Nijmeijer, & Alvarez, 2016), power networks (Dörfler, Chertkov, & Bullo, 2013), neuronal networks (Belykh & Hasler, 2011) and generic phase-oscillators (Dörfler & Bullo, 2014;Strogatz, 2000), with this list being definitely incomplete. See also Fazlyab, Dörfler, and Preciado (2017), Liu and Iwasaki (2017), Nishikawa and Motter (2006) for results on the design (or optimization) of the interaction topology for synchronization of nonlinearly coupled systems. For (firstand second-order) integrator networks, nonlinear couplings have been employed in Andreasson, Dimarogonas, Sandberg, and Johansson (2014) and Saber and Murray (2003) to achieve consensus -a form of synchronization in which all systems converge to the same constant steady state, see also (Yu, Chen, & Cao, 2011;Yu, Chen, Cao, & Kurths, 2010) for the first-and secondorder consensus problems for systems with nonlinear dynamics. Greater flexibility of nonlinear couplings makes them particularly interesting for design of couplings to achieve consensus and flocking of multi-robot systems, in particular, for collision avoidance and preserving connectivity, cf. Dimarogonas and Kyriakopoulos (2008), Ji and Egerstedt (2007) and Poonawala and Spong (2017) and the references therein.
Common ways of handling nonlinearities in synchronizing systems include changing/cancelling nonlinearities by feedback to achieve a desired system form or properties (e.g. passivity) (Arcak, 2007;Chopra & Spong, 2008), suppressing nonlinearities by high-gain couplings (Pogromsky, 1998), ensuring synchronization through absolute stability arguments (Proskurnikov, 2013) or achieving synchronization or consensus by various methods, e.g. passivity, with synchronization-error feedback couplings, linear or nonlinear, e.g. Arcak (2007) and Yu et al. (2011Yu et al. ( , 2010. Most of the works dealing with synchronization of nonlinear systems focus on attaining synchronization (or consensus). The practically important questions of transient and steady-state performance of the overall synchronizing interconnected system, energy efficiency of the coupling laws and their sensitivity to noise, remain mostly unaddressed. While for linear systems there is a well-developed machinery to address these questions, see, e.g., Liu, Saberi, Stoorvogel, and Nojavanzadeh (2020) and references therein, for nonlinear systems the problem is much more complex and very few works addressing performance of nonlinear synchronizing systems exist, see, e.g., El-Gohary (2006), Macellari, Karayiannidis, and Dimarogonas (2017) and Modares, Lewis, Kang, and Davoudi (2018). Nonequilibrium steady-state dynamics, as in the problem of synchronization of oscillatory systems, makes the problem even more challenging.
Performance optimization in nonlinear synchronizing systems is directly linked to finding the range of coupling laws that yield synchronization. This range can often be characterized by the lowest and highest (in a certain sense) coupling gains. In synchronization problems it is often only the lower bound that needs to be found, as all gains higher than that bound would yield synchronization. Couplings/controls with gains at this lower bound often lead to lower sensitivity to measurement noise and higher energy efficiency in attaining and maintaining synchronization. On the other hand, finding synchronizing coupling laws with minimal gains can increase our understanding of the system dynamics and reveal novel synchronization mechanisms.
Finding synchronizing coupling laws with minimal gains, as the first step towards performance optimization/analysis of nonlinear synchronizing systems, is the question that motivates our study. For nonlinear systems, minimal synchronizing couplings should generically be nonlinear. Nonlinear coupling laws can be considered as couplings with varying gains. They can fully utilize the intrinsic system dynamics: the gains should generically be low, or even zero, where dynamics is favorable to synchronization, and increase/be nonzero only in those parts of the state space, where system nonlinearities act against synchronization. This observation calls for research in finding minimal couplings, as it excludes most of the coupling laws studied in the literature: linear and nonlinear couplings being functions of the synchronization error e only, cf. Dörfler and Bullo (2014) and Strogatz (2000)-such couplings do not distinguish different parts of the state space, as long as the synchronization error e remains the same.
In the light of the above exposition of the existing literature, the current paper takes on the following combined open challenge: how to design nonlinear coupling laws that (1) achieve synchronization in networks of nonlinear systems (including oscillatory systems), (2) achieve this synchronization through nonlinear couplings that exploit system nonlinearities where they help synchronization and counteract nonlinearities where they act against synchronization, (3) provide performance improvements over conventional, e.g., linear, coupling laws.

Contributions
We present a systematic approach to the problem of design of nonlinear coupling functions to establish synchronization in networks of identical nonlinear single-input-single-output systems. Synchrony is here to be understood as the asymptotic match of the states of the systems, while the nonlinear systems are allowed to exhibit complex oscillatory (even chaotic) intrinsic dynamics. The presented approach allows one to find nonlinear couplings with minimal gains (understood in a nonlinear sense to be clarified later) that are tailor-made to system's nonlinearities -both helping and counteracting synchronization depending on the location of the system in its state-space.
The proposed couplings are of the form of a definite integral of some non-negative weight function (or, as we call it, a nonlinear gain -the main design parameter) with the limits being the outputs of the systems. This type of coupling has a greater degree of flexibility not only compared to linear (diffusive) couplings, but also compared to arbitrary nonlinear couplings depending only on the synchronization error. It enables variable coupling gains in different parts of the state space -the flexibility essential for finding coupling gains tailor-made to system's nonlinearities.
To formalize in which part of the state space and to what extent system nonlinearities act against synchronization, we introduce the novel notion of incremental feedback passivity with a nonlinear gain. The minimal nonlinear integral coupling is then selected to counteract system nonlinearities only in these countersynchronization parts of the state space, while being zero in the parts where system dynamics render synchronization naturally.
To extend nonlinear integral couplings from 2 to N interconnected systems with a fixed communication graph, we develop new analysis and design tools, as conventional methods based on analysis of the Laplacian matrix do not apply to couplings with gains that can be zero over large parts of the state space. Based on these tools, we present a class of networks -characterized by the existence of hierarchies of synchronization subspaces, expressed in terms of relaxed balanced coloring and sequential decoloring of the network graph -which support synchronization via nonlinear integral coupling in its full potential.
We demonstrate by two examples that the proposed nonlinear integral couplings can render synchronization by lower input energy (understood in the L 2 sense of the coupling input functions), and with less sensitivity to measurement noise compared to the minimal linear coupling laws that guarantee the same rate of convergence. For the example with FitzHugh-Nagumo oscillators, the proposed integral couplings lead to a novel type of synchronization: the couplings are zero most of the time and synchronization is achieved and maintained by spikes in the nonlinear coupling gains. This result can be interesting in neuroscience, where most of the communication happens via spikes.
The developed generic framework is supplied with constructive results leading to verifiable conditions and concrete design procedures for a class of nonlinear systems. This paper extends our initial results reported in Pavlov, Steur, and van de Wouw (2009) and Pavlov, Proskurnikov, Steur, and van de Wouw (2018): the new results are less conservative, exploit the full potential of the nonlinear coupling functions in networks, and are supplied with examples demonstrating novel phenomena.

Organization
In Section 2, we state the problem of controlled synchronization and introduce nonlinear integral couplings. We also present a simple example that demonstrates the potential of the proposed type of coupling in comparison to the conventional linear coupling functions. In Section 3 we introduce incremental feedback passivity with a nonlinear gain-the notion central in the design and analysis of systems interconnected through integral couplings. In that section we also present constructive conditions for verifying this property for a class of nonlinear systems. Section 4 covers synchronization of two systems interacting via nonlinear integral coupling. Extension of these results to synchronization of multiple networked systems is presented in Section 5, where we define the class of networks and the main result on synchronization in such networks with nonlinear integral couplings. Section 6 illustrates the developed theory with an application to synchronization of FitzHugh-Nagumo oscillators, which represent simplified models of neuron dynamics and are a benchmark in the study of synchronization of nonlinear oscillators. Conclusions are drawn in Section 7.

Controlled synchronization problem
In this paper we consider N identical nonlinear systems of the form 1 We assume that f (x, t) satisfies the following technical condition: sup t≥0 |f (0, t)| < +∞. The problem of controlled synchronization studied in this paper is to find control laws for each u i that guarantee 1 General results presented in this paper also hold for systems of the forṁ under an additional condition that |∂h/∂x(x)| ≤ C h for some C h > 0 and all x ∈ R. For simplicity of presentation we have chosen a less generic form (1).
A. boundedness of solutions of the interconnected systems. 2 B. asymptotic synchronization of the systems' states: (2) For each system i, u i is allowed to depend on the system's output y i and on the outputs of neighboring systems that can communicate to system i. In addition, it is required that for identical outputs y 1 = y 2 = · · · = y N , the controls satisfy u 1 = u 2 = · · · = u N = 0, such that in exact synchrony the systems exhibit the (oscillatory) dynamics of the unforced systems in (1). The set of systems that can communicate to system i is denoted by N i . For i = 1, . . . , N, these sets define communication graph G = (V, E), where V = {1, . . . , N} is the set of nodes, which represent the systems (1), and E ⊂ V × V is the set of edges.

Motivating example
Let us demonstrate the main idea put forward in this paper (namely that of using nonlinear integral couplings) with a simple example of synchronization of two scalar systems given bẏ with a bounded piecewise-continuous excitation term G(t). The approach developed in this paper, can be demonstrated by representing the system dynamics in an equivalent form: Then the integral coupling u 1 = −u 2 = ∫ y 2 y 1 λ(s)ds gives the following error dynamics, with the synchronization error defined as e := y 1 − y 2 : where ξ (y 1 , y 2 ) ∈ (y 1 , y 2 ) results from the mean value theorem. By designing λ(s) we can ensure that ψ(s) ≤ −ε for some ε > 0 and all s ∈ R. Thus, for V (e) = 1/2e 2 , we haveV = 2ψ (ξ )V ≤ −2εV , which guarantees exponential synchronization, provided that solutions of the closed-loop system are defined and bounded for all t ≥ 0. The latter condition needs to be proven separately. The design of λ(s) can be done in different ways: • λ(s) ≡ K := 1/2 sup s∈R (1 − s 2 + ε) = 1/2(1 + ε) gives linear coupling; • λ(s) = 1/2(1 − s 2 + ε) cancels the systems dynamics (feedback linearizing coupling) and imposes linear GES error dynamicsė = −εe; • λ(s) ≥ 1/2 max{1 − s 2 + ε, 0} cancels the system dynamics only in the area where they act against synchronization and does nothing in the areas where synchronization occurs naturally.
These three choices are illustrated in Fig. 1. 3 The last approach, while providing synchronization with the same guaranteed convergence rate as for the other two approaches, has several favorable properties. Namely, for λ(s) = 1/2 max{1 − s 2 + ε, 0}: (1) the coupling law is more energy-efficient since it avoids coupling actions in the areas where synchronization is naturally 2 Boundedness of solutions is a natural and desired property in many applications. Moreover, if overlooked, it may even undermine the validity of proof of asymptotic synchronization, as demonstrated in Proskurnikov and Cao (2017).
3 As it will be seen from technical results and examples presented further in the paper, the nonlinear integral coupling does not necessarily have to consist of a combination of a feedback-linearizing coupling and a zero coupling, as it is the case for the example presented in Fig. 2. and y 2 lie in the area where λ(s) ≡ 0: for y 1 , y 2 > (2) Since λ(s) ∈ L 1 , then the coupling actions are always bounded by |u 1,2 (t)| ≤ ∫ +∞ −∞ λ(s)ds.
(3) It has a lower coupling gain (if we define the coupling gain as g = ∫ y 2 y 1 λ(s)ds/(y 2 − y 1 ), as in Pavlov et al. (2009)) than for the other two approaches. The lower coupling gain generically gives lower sensitivity to measurement noise on the measured outputs y 1 and y 2 . Notice that these favorable properties do not hold for the other two approaches. Moreover, they cannot be achieved by couplings of the form u = φ(y 2 − y 1 ) often studied in the literature, cf. Dörfler and Bullo (2014) and Strogatz (2000), as these couplings do not distinguish between different locations of the output variables y i .
These benefits are illustrated through simulations for G(t) = 3 sin t in Fig. 2, which also includes a comparison to synchronization by the linear coupling with the minimal gain K = 1/2(1 + ε).
As can be seen from the figure, nonlinear integral coupling yields synchronization by lower coupling values, with lower sensitivity to measurement noise (both in controls and steady-state synchronization error), and with higher energy efficiency. Over the simulated period, the L 2 norms of input signals for the linear and nonlinear integral couplings are equal to ∥u Lin Note that although the nonlinear integral coupling yields the same guaranteed synchronization rate characterized by ε, its actual convergence rate for the given initial conditions can be lower than for the case of linear coupling. By designing a larger gain function λ(s) ≥ 1/2 max{1 − s 2 + ε, 0}, which is the intrinsic property of the method, one can improve the transient convergence rate, yet at the expense of increased sensitivity to the measurement noise and decreased energy efficiency. The possibility of such a trade-off between convergence rate and noise sensitivity/ energy efficiency (at lower coupling gains) is the flexibility not present in the case of linear couplings: the constant coupling gains needed for the guaranteed convergence rate ε are lower bounded by K = 1/2(1 + ε).
In the remainder of the paper we investigate the application of nonlinear integral coupling of the form u = ∫ y j y i λ(s)ds to ensure synchronization of N n-dimensional nonlinear systems coupled through a communication network with a fixed topology. Our motivation is to design synchronizing couplings with lower gains and better performance properties compared to the standard methods of linear and nonlinear couplings based on error feedback only, and nonlinear couplings based on feedback linearization. We are interested both in generic conditions under which one can find such integral nonlinear couplings (both on the level of the system dynamics, as well as on the level of communication network topology) and in constructive conditions for specific classes of nonlinear systems.

Nonlinear integral coupling
Below, we give a formal definition of nonlinear integral coupling between dynamical systems.

Fig. 2.
Simulation results for synchronization with measurement noise for nonlinear integral coupling, ε = 0.001, y 1 (0) = 2, y 2 (0) = 0 -plot 1. Comparison with synchronization by linear coupling with the gain K = 1/2(1+ε) -plots 2-4 for synchronization error e, coupling value u 1 and coupling gain g. Definition 1. Given N systems with a communication topology specified by the graph G with the corresponding sets of neighboring nodes N i , i = 1, . . . , N, we say that systems are interconnected through integral coupling with a nonlinear coupling This kind of coupling satisfies the requirement that for identical outputs the coupling must be zero, since for y i = y j ∀i, j ∈ V, we have u i = 0. Notice that for λ(s) ≡ K = Const > 0, the nonlinear integral coupling (4) represents linear diffusive

Incremental feedback passivity
In this section, we present the definition, some properties and constructive conditions for incremental feedback passivity (iFP) with a nonlinear gain, which is the key property for synchronization analysis in this paper. It is an extension of the incremental passivity property employed in various works within nonlinear systems and control, see, e.g. Pavlov and Marconi (2008) and Stan and Sepulchre (2007).

Incremental feedback passivity with nonlinear gain
Definition 2. System (1), for N = 1, is called incrementally feedback passive with nonlinear gain γ (s) ≥ 0 -denoted as iFP(−γ (s)) -if there exists a C 1 function S(x) : R n → R + and a function ρ(x) : R n → R + such that for any two inputs u a (t) and u b (t) any two solutions x a (t), x b (t) of system (1), for N = 1, corresponding to these inputs, with outputs y a and y b , satisfy the inequality .
For γ (s) ≡ 0 we obtain the standard definition of incremental passivity. For γ (s) ≥ 0, system (1) can be made incrementally passive by feedback transformation u = ∫ p(t) y λ(s)ds + v, for some integrable function λ(s) ≥ γ (s), ∀s ∈ R, arbitrary continuous function p(·) and v, which is a new input after the transformation.
In this case, for two outputs y a and y b , the corresponding u a and i.e., the system is incrementally passive with respect to the new input v. 4 For the synchronization analysis in this paper, we require an additional assumption on S(x) and ρ(x).
This assumption holds, for example, for quadratic positive (1) is given in the following lemma (see the proof in the Appendix). For notational convenience, in this lemma we omit the subscript i from the state x and control u.

Constructive conditions
All subsequent results on synchronization are formulated in terms of the iSFP(−γ (s)) property. In this section, we provide constructive results on how to verify this property and to find γ (s) in (5) for systems of the form (1). All the proofs can be found in the Appendix. (5)). If R > 0, then system (1) is iSFP(−γ (s)) with functions S(x) and ρ(x) satisfying Assumption 1.
This result can be made more constructive (see Theorem 2 below) for the following class of systems: (10) where q(z, y, t) and p(z, y, t) are continuously differentiable in (z, y) and piecewise continuous in t. Notice that this system is of the form (1) holds for all (z, y) and t ≥ 0. Letγ (y) satisfỹ for all (z, y) and some ϵ > 0 satisfying where I n−1 is the (n − 1) × (n − 1) identity matrix. Then system (10) is iSFP(−γ (s)) with γ (s) = max{0,γ (s)} and quadratic positive definite S(x) and ρ(x) satisfying Assumption 1.
A functionγ (y) satisfying (12) can be found if the right-hand side of (12) is independent of z and t or can be bounded from above by a y-dependent function.
Condition (11) guarantees that the zero dynamics of system (10) (i.e., the z-dynamics) are convergent (Pavlov, Pogromsky, van de Wouw, & Nijmeijer, 2004), which implies that for a given function y(t) all solutions of the systemż = q(z, y, t) converge to a unique bounded globally asymptotically stable steady-state solution determined only by y(t). This property of the zero dynamics can be considered as a specific minimum phase property of the overall system (10).

Synchronization of two systems
Below we apply the iSFP(−γ (s)) property to the controlled synchronization problem for two systems. (1) interconnected through the integral coupling:

Theorem 3. Consider two systems a and b of the form
then all solutions of the closed-loop system (1), (14) are bounded for all t ≥ 0 and satisfy (2) for i = a, j = b.

Synchronization: From
Since x a (t) and x b (t) belong to a compact set for all t ≥ 0, by Barbalat's lemma (Khalil, 1996, Section 8. implies (2). □ Remark 1. The case of master-slave synchronization, when one of the control inputs in (14) equals zero (e.g., u a = ∫ y b ya λ(s)ds, u b = 0) is proven in the same way. In this case, condition (15) is replaced by the condition λ(s) ≥ γ (s), ∀s ∈ R.
Remark 2. In case if S(x) = x T Px and ρ(x) = x T Rx for some P = P T > 0 and R = R T > 0, one can easily see that (16) implies exponential synchronization, i.e.: for some α > 0, β > 0. In this case, one can expect a certain degree of robustness of synchronization against variations in the right-hand sides of individual systems (1) and communication effects such as sampling and time-delays. Depending on how these variations appear in the systems' dynamics, synchronization is achieved either exactly (i.e. the asymptotic match of the systems' states) or practically (i.e. the differences in all state-trajectories evolve to a small neighborhood of zero), cf. Panteley and Loría (2017) and Steur and Nijmeijer (2011).
Remark 3. The condition that ρ(x) is positive definite in (5) can be relaxed; In particular, ρ(x a − x b ) can be substituted by ρ(y a − y b ), where ρ(y) is a positive definite function. In this case, the conditions in Theorem 3 will guarantee output synchronization: y a (t) − y b (t) → 0, as t → +∞. Full state synchronization can be achieved by the observability condition that y a (t) − y b (t) → 0 and u a (t) − u b (t) → 0 imply x a (t) − x b (t) → 0, as t → +∞.

Synchronization of multiple systems
In this section, we extend the results from the previous section to the case of N systems (1) interconnected through integral coupling (4) with a communication network given by graph G. This task appears to be challenging because of the nonlinear nature of the coupling (4), which prevents the use of standard techniques for synchronization analysis based on analysis of the Laplacian of the network graph. The main difficulty is related to the fact that λ(s) is allowed to be zero, thereby leading to regions of output values where nodes lose connectivity. In these regions, it is not the coupling, but the original system dynamics that takes care of synchronization. Synchronization results stemming from absolute stability methods, see, e.g., Proskurnikov (2013Proskurnikov ( , 2014, cannot be applied here given the generic nonlinear nature of both the integral couplings and the system dynamics. The only applicable result in this direction was reported in Pavlov et al. (2018), stating that, for a given λ(s) and given network topology with a connected bidirectional communication graph, there exists k * > 0 such that synchronization is achieved through coupling That result is more of a qualitative nature, since there are no constructive methods for finding k * . In this section, we provide a method for analyzing synchronization of N systems interconnected through nonlinear integral couplings over directed communication graphs. We firstly explain the idea behind the method with an illustrative example, then we provide the formal method for synchronization analysis of multiple systems interconnected over a given communication graph.

Illustrative example
Consider a network of 4 systems with iSFP(−γ (s)) property satisfying Assumption 1 interconnected through nonlinear integral couplings (4) with the communication graph as in Fig. 3. To demonstrate synchronization of all nodes, let us first consider nodes 1 and 3. The nonlinear integral couplings u 1 and u 3 are given by Taking into account the fact that 4λ(s)ds. Substituting this equality in (5), gives .
for Γ = 4, then d dt S(x 1 (t) − x 3 (t)) ≤ −ρ(x 1 − x 3 ), which, by boundedness of solutions and Barbalat's lemma (Khalil, 1996, Section 8.3), implies that x 1 (t) and x 3 (t) asymptotically synchronize. Applying the same reasoning for nodes 2 and 4, we obtain that these nodes synchronize if (21) holds for Γ = 3. The difference in Γ from the case of nodes 1, 3 is due to the fact that there is only one link between nodes 2 and 4 compared to two links between 1 and 3. Thus, under condition (21) where the vanishing terms correspond to synchronizing nodes within the two clusters (provided that λ(s) is bounded on R). Applying inequality (5) to nodes 1 and 2, together with expressions (22), (23), we obtain It can be demonstrated that since ϵ 1 (t), ϵ 2 (t) → 0, as t → +∞, under Assumption 1 the last inequality implies that x 1 (t) and x 2 (t) asymptotically synchronize (for rigorous proofs see subsequent technical results). This leads to synchronization of two clusters {1, 3} and {2, 4}, i.e., synchronization of the overall network.
In our example we sequentially merged nodes and clusters into larger clusters of synchronizing nodes until all nodes in the network are merged in one synchronizing cluster. This process, together with the corresponding values of Γ that guarantee that after each merger, the new cluster will still consist of synchronizing nodes, is described below where each partition P i , i = 0, . . . , 3 consists of either individual nodes or clusters of synchronizing nodes. Note that for Γ = min k=1,2,3 Γ k , condition (21) will be satisfied also for all Γ k . Condition (21) with this Γ is our final condition on the nonlinear gain λ(s).
For general graphs, this method is applicable under certain conditions on the graph. The idea behind these conditions can be understood from the example above . The key property to ensure merging of two individual nodes a, b (a, b = 1, 3 or a, b = 2, 4 in our example) into a larger synchronizing cluster is that they have the same external neighbors, i.e. N a \{a, b} = N b \{a, b}.
Furthermore, external neighbors from the same synchronizing cluster (i.e., these have the same asymptotic behavior) can be treated as the same neighboring node (for details, see exact proofs of further technical results in this section). Each merger operation is thus performed on two nodes/clusters that have (asymptotically) the same external neighbors. Graphically, the easiest way to represent this sequential operation is by coloring the nodes and clusters, such that nodes in merged clusters get the same color. This approach is formalized in the next subsection.

Synchronization analysis by relaxed balanced coloring
Relaxed balanced coloring and related notions presented below provide both a formal framework as well as an algorithmic support for the method demonstrated in the previous subsection. The main result presented in this section provides synchronization conditions for a class of networks -so-called sequentially decolorable networks -defined below. We will arrive at this result in three steps: (1) definition of relaxed balanced coloring (it conveniently captures the necessary network partitions into nodes with the (asymptotically) the same neighboring nodes), (2) definitions of color reduction and sequentially reducible networks -the class of networks that allows for sequential merger of synchronizing clusters leading to synchronization of the whole network and (3) the main synchronization result that links nonlinear integral coupling and sequentially decolorable networks.

Relaxed balanced coloring
Definition 3. A coloring of the nodes with k ∈ {1, . . . , N} colors c 1 ,. . . , c k is called a relaxed balanced coloring if • each node is assigned a color, and • every c i -colored node receives an equal number of edges from c j -colored nodes for all j ∈ {1, 2 . . . , k}\{i}.
This notion is a relaxation of balanced coloring, where each c icolored node receives an equal number of edges from c j -colored nodes for all j ∈ {1, 2 . . . , k}, and from the nodes with its own color c i . 6 It will be shown below that this relaxation is enough for our purposes and that it covers cases not covered by the (non-relaxed) balanced coloring.
A graph can be colored according to relaxed balanced coloring in multiple ways. Each coloring defines the corresponding partitioning of the nodes denoted by P i . The two trivial colorings are given by (a) assigning each node an individual color and (b) by assigning all nodes the same color. An example of a graph and the corresponding relaxed balanced colorings are presented below.
Example 1. Consider the graph G shown in Fig. 4. All its relaxed balanced colorings are presented in Fig. 5. The corresponding partitions defined by the colors are: Suppose that a given graph G admits K additional relaxed balanced colorings besides the above two trivial ones. Of course, we require these K other relaxed balanced colorings to be distinct in the sense that the partitions corresponding to the relaxed balanced coloring are all different. We then denote the relaxed balanced coloring with N-colors by partitioning P 0 = {{1}, {2}, . . . , {N}}, the coloring with one color by partitioning P K +1 = {{1, 2, . . . , N}}, and the other relaxed balanced colorings by partitionings P j , j = 1, . . . , K . 6 (Relaxed) balanced coloring is a common terminology in the analysis of synchronization (in particular, cluster synchronization) (Belykh & Hasler, 2011;Golubitsky, Stewart, & Török, 2005;Steur, Ünal, van Leeuwen, & Michiels, 2016). For the problem considered in this paper, a relaxed balanced coloring is equivalent to an almost equitable partition (Cardoso, Delorme, & Rama, 2007), which is frequently used in the framework of consensus theory and the analysis of network controllability (Gambuzza & Frasca, 2019; Notarstefano & Parlangeli, 2013).

Color reduction
Relaxed balanced colorings can be ordered in a hierarchy by means of refinement. Here a partition P i of the set V is a refinement of a partition P j of V, denoted by P i ⊂ P j , if every element of P i is a subset of some element of P j . We say that P i is finer than P j , and P j is coarser than P i . For example, P 0 is finer than P K +1 , or, equivalently, P K +1 is coarser than P 0 . Clearly, if P i is finer (resp. coarser) than P j , then the relaxed balanced coloring associated to P i has more (resp. less) colors than the one associated to P j . Similar hierarchies have been utilized for analysis of synchrony subspaces (Aguiar & Dias, 2014).
This hierarchy defines multiple ''chains'' of refinement. For example, the network graph shown in Fig. 4 has multiple chains of refinement, e.g.
These chains of refinement can be represented in a diagram called a Hasse diagram (Birkhoff, 1948). In this diagram, there is an arrow from P i to P j (i.e. P i → P j ) if and only if P j ⊂ P i and there is no other partition P ℓ such that P j ⊂ P ℓ ⊂ P i . (The edges in a Hasse diagram are sometimes undirected. We prefer the use of directed edges.) Furthermore, if P i ⊂ P j , then P j appears lower in the diagram than P i . Consequently, P 0 is at the top of the diagram and P K +1 is at the bottom. For the network graph considered in Example 1 (Fig. 4), the corresponding Hasse diagram is presented in Fig. 6. Definition 4. A network G that admits a relaxed balanced coloring P i with k i colors from the set C i = {c 1 , . . . , c k i } is colorreducible if there exists another relaxed balanced coloring P j with k j < k i colors from the set C j such that • for all ℓ = 1, . . . , k j , every c ℓ -colored node in P i has (or can be assigned) the same color in P j .
Clearly, all possible color reductions starting from a relaxed balanced coloring P i are specified by the directed paths in the Hasse diagram with P i as the root node. The following result is immediate. Fig. 6. Hasse diagram for the network graph from Example 1 (Fig. 4) and the corresponding chains of refinement.
Lemma 2. If the Hasse diagram of relaxed balanced colorings of a graph G with the operation of refinement contains a directed path of length N −1, then there exists a sequence of color reductions from P 0 to P K +1 , where in each reduction step exactly one color is removed.

Definition 5. A network that satisfies the condition of Lemma 2 is called sequentially decolorable.
The property that the graph G is sequentially decolorable is the main condition for synchronization analysis presented further. It can be verified algorithmically. The corresponding MATLAB code is available at Steur (2019).
Note that the existence of a path of length N − 1 in the Hasse diagram is equivalent to the existence of a chain of refinements 25) of length N − 1, with i 1 , i 2 , . . . , i N ∈ {0, 1, 2, . . . , K + 1} (here the length of the chain is understood as the number of symbols ⊂ in the chain). As can be seen from the Hasse diagram presented in Fig. 6, the network graph from Example 1 (Fig. 4) is sequentially decolorable, with the corresponding chain of refinements given by (24).
The following lemma details the structure of the chain of refinements of length N −1 of a sequentially decolorable network (see the Appendix for the proof).
Lemma 3. Consider a sequentially decolorable network G with the corresponding chain of refinements (25). In each step from i k to i k+1 in the chain of refinements (25), P i k+1 is obtained from P i k by merging two elements: for some A, B, C 1 , . . . , C l . Moreover, for all a ∈ A and b ∈ B it holds Moreover, |N a ∩ B|, |N b ∩ A|, |N a ∩ C j | and |N b ∩ C j |, ∀j = 1, . . . , l, are independent of the particular choice of a ∈ A and b ∈ B.
Before formulating the main result of this section, let us give the following definition.
Definition 6. Consider a sequentially decolorable network graph G with the corresponding chain of refinement (25). For each transition P i k → P i k+1 given in (26), k = 1, . . . , N − 1, define the local coupling characteristic Γ k : for some a ∈ A, b ∈ B, with A, B defined in (26). The global coupling characteristic Γ is defined as Even though Γ k are defined for some arbitrary a ∈ A and b ∈ B, according to Lemma 3, it is independent of the specific choice of a ∈ A and b ∈ B. Thus Γ k , k = 1, . . . , N − 1 and Γ are well defined. Γ k and Γ are defined for a specific choice of the chain refinement. Thus Γ may vary depending on the selected chain of refinement (25). For the network from Example 1, one can see, by simple observation, that the global coupling characteristics equals Γ = 1.  The result of Theorem 4, in combination with previously presented constructive results, can be considered both as analysis and design tools. For design problems, as demonstrated by the example in Section 2.1, the nonlinear integral coupling with the minimal coupling ''gain'' λ(s) = γ (s)/Γ can lead to more energy efficient synchronization with lower sensitivity to measurement noise. For analysis problems, the minimal nonlinear coupling gain λ(s) can provide better understanding of the system dynamics and reveal new synchronization mechanisms, as will be demonstrated by an example in Section 6. 7 One can avoid the condition λ(·) ∈ L 1 in proving boundedness of the solutions by using the machinery of semi-passive systems (Pogromsky, 1998)

Illustrative example
Let us consider the FitzHugh-Nagumo (FHN) oscillator, which represents a simplified model of neuron dynamics (FitzHugh, 1961): where y i represents the membrane potential, z i is an internal variable related to the ionic currents, and input u i is used to establish coupling with other neurons. All other parameters are positive constants. For numerical simulations we choose the following values: a = 0.7, b = 0.8, τ = 100, and r = 0.3. We added a common periodic forcing Aext τ cos ωt to the internal z i -dynamics, which can make the neurons oscillate irregularly. (Note that for a = b = 0, we recover the dynamics of the forced Van Der Pol oscillator, which is known to behave chaotically for certain parameter values.) We take A ext = 0.03 and ω = 0.026.
We consider two cases of network topology: (a) two nodes interconnected with bi-directional coupling and (b) network topology from Fig. 4. Analysis of synchronization in a network of such oscillators with linear coupling is presented in Gorban, Jarman, Steur, van Leeuwen, and Tyukin (2015). Let us apply the theory developed in the previous sections and find a nonlinear integral coupling (4)  [A1]: System (31) is iSFP(−γ (s)) and satisfies Assumption 1: System (31) is of the form (10). Let us find γ (s) for which this system is iSFP(−γ (s)) with quadratic positive definite S(x) and ρ(x).
Notice that the latter condition guarantees that Assumption 1 is satisfied. We find this γ (s) using Theorem 2. Condition (11) where ϵ > 0 is an arbitrary parameter satisfying ϵ < b due to (13). For implementation, we select ϵ = 0.0001.
[A2]: The communication network graph is sequentially decolorable. Both networks considered in this example (two bi-directionally coupled systems and six systems coupled as in Fig. 4) are sequentially decolorable with the corresponding minimal coupling characteristics equal to Γ = 2 (for 2 systems) and Γ = 1 (for 6 systems).
[A3]: Selection of the coupling gain λ(s): As mentioned in the previous sections, we are interested in finding the smallest coupling that leads to synchronization. The structure of the smallest coupling gain λ(s) can give insights into mechanisms of synchronization present in real systems. For design problems, it can lead to lower coupling actions needed to achieve synchronization, lower energy spent on coupling actions and lower sensitivity to noise. Therefore, we select the minimal λ(s) according to the lower bound in (  Simulation results for two bi-directionally coupled FHN models: system states, controls and the variable coupling gain g 21 (t) between the 1st and 2nd system.
the generalized variable gain g 21 (t) of the nonlinear integral coupling defined as g 21 (t) = ∫ y 1 (t) y 2 (t) λ(s)ds/(y 1 (t) − y 2 (t)), Pavlov et al. (2009). We observe that the two systems rapidly synchronize. The generalized variable gain is most of the time equal to zero, indicating that no interaction is required. Only when the neuron spikes, the generalized variable gain is non-zero with a maximum of (1 + ϵ)/2. This variable gain (seen over the whole simulation period) is much lower than the best estimate of the synchronizing linear diffusive coupling gain we are aware of, which equals 0.5 (Gorban et al., 2015;Stan & Sepulchre, 2007). Lower coupling gains is a distinctive feature of the proposed method over linear diffusive couplings. Another feature of the proposed coupling is that it is bounded for all values y 1 , y 2 , which is not the case for linear diffusive coupling. Interestingly, our proposed integral coupling turns out to generate a pulse-type of interaction, as demonstrated by the spiky generalized variable gain. The gain is zero most of the time and spikes only when the y variable jumps. This is a novel type of synchronizing interaction which turns to be more energy efficient and less sensitive to measurement noise, compared to linear diffusive couplings. This will be demonstrated by the next simulation with 6 systems and measurement noise.
In Fig. 8 we depict synchronization of six FHN oscillators interconnected as in Fig. 4 through the integral coupling with the minimal gain λ(s) = γ (s)/Γ (note that in this case Γ = 1). We add noise to the outputs to simulate the effect of measurement noise. It is clearly visible in the inputs u i and the synchronization error e := max i,j∈1...6 |y i − y j | that this measurement noise only affects the dynamics when a spike in the outputs is either initiated or ended, i.e. when the generalized variable gain is nonzero. In Fig. 9 this simulation is repeated with linear coupling, λ(s) ≡ 1 + ϵ, which guarantees the same rate of convergence characterized by ϵ. Comparing the nonlinear and linear coupling results, we clearly see that the nonlinear coupling is less sensitive to measurement noise both in the control inputs u i , and the steady-state synchronization errors e i , at least, in L 1 and L 2 sense. The energy efficiency of the nonlinear integral coupling follows from comparing L 2 norms of the inputs for both the linear and nonlinear couplings: ∥u lin ∥ L 2 = 18.9, ∥u nl ∥ L 2 = 10.6. For the infinity norm, we have∥u lin ∥ L∞ = 10.0, ∥u nl ∥ L∞ = 2.7. Similar relations were observed for various initial conditions. These results demonstrate the effectiveness of the selected nonlinear integral couplings also for the case of more complex network topologies. Another example demonstrating the effectiveness of the proposed approach can be found in Pavlov et al. (2018), where the developed theory was applied to Hindmarsh-Rose oscillators,  in which not only the y-dynamics, but also the z-dynamics are nonlinear.

Conclusions
In this paper, we have considered the controlled synchronization problem for a class of identical nonlinear single-inputsingle-output systems 8 and proposed a novel nonlinear integral couplings design to achieve synchronization of systems on a class of interconnection networks. As demonstrated by the examples, the proposed nonlinear integral couplings with minimal gains (understood in a nonlinear sense) can lead to more efficient synchronization with lower input energy and lower sensitivity to measurement noise compared to conventional linear couplings. Application of the results to synchronization of FitzHugh-Nagumo oscillators revealed a novel type of synchronizing interaction, where synchronization is achieved and maintained through generically zero, but spiking coupling gains. From a design point of view, the obtained results give a constructive method to design nonlinear couplings and communication networks that render more efficient synchronization of the interconnected systems. From an analysis point of view, they can provide insight into potential nonlinear synchronization mechanisms in nature.

Acknowledgments
The authors would like to thank anonymous reviewers for their valuable comments, which helped to improve the paper, and Danila Semenov, Ph.D. student from St. Petersburg State University, Russia, for conducting initial simulations and analysis of the example from Section 6.

Appendix. Proofs
A.1. Proof of Lemma 1 Inequality (5) can be written in the form which holds for any x a , x b , u a , u b and y a = Cx a , In this derivation we have utilized inequalities (7), (8) Consider system (1) written in an equivalent form: 0γ (s)ds. With this new notation, condition (9) is equivalent to This matrix inequality implies (see, e.g. Demidovich (1967) and Pavlov et al. (2004)) that for any x 1 , x 2 it holds that Now consider the storage function S(x) = 1/2x T Px. For any two solutions x 1 (t) and x 2 (t) of system (1) with inputs u 1 , u 2 and outputs y 1 , y 2 : where in the last inequality we utilized (A.4), equality PB = C T and the additive property of definite integrals. Taking into account the definition γ (s) = max{0,γ (s)}, we obtain for all y 1 , y 2 . Substituting (A.6) into (A.5), we obtain (5) with quadratic S(x) = 1/2x T Px and ρ(x) = 1/2x T Rx. Thus, system (1) is iFP(−γ (s)). Moreover, if R is positive definite, then the system is iSFP(−γ (s)) and the functions S(x) and ρ(x) satisfy Assumption 1. □

A.3. Proof of Theorem 2
Choose the matrices P and R in (9)  ] .
Notice that this P satisfies the equality PB = C T . By combining all the terms in the first inequality in (9) in the right-hand side, one can see that for the chosen P = P T > 0 and R = R T > 0 this matrix inequality is equivalent to Due to (11) Since P i k is a refinement of P i k+1 , every element of P i k is a subset of an element of P i k+1 . Since the number of elements (colors) in P i k+1 is less than the number of elements (colors) in P i k by exactly 1, two elements from P i k will form one element of P i k+1 . Let us denote these two elements by A and B and the remaining elements of the partition P i k by C 1 , . . . , C l for some l. Then (26) holds. Let a ∈ A and b ∈ B. Since P i k+1 is a partition corresponding to relaxed balanced coloring, and A ∪ B ∈ P i k+1 , then, by the definition of relaxed balanced coloring, for any two points a, b ∈ A∪B, (27) holds. In particular, it holds for any a ∈ A and b ∈ B.
Finally, since P i k is a partition corresponding to relaxed balanced coloring, for any a 1 , a 2 ∈ A, it holds that N a 1 ∩ B = N a 2 ∩ B and N a 1 ∩ C j = N a 2 ∩ C j , for j = 1, . . . , l. Thus |N a ∩ B|, and |N a ∩ C j |, and, in the same way, |N b ∩ A| and |N b ∩ C j |, ∀j = 1, . . . , l, are independent of the particular choice of a ∈ A and b ∈ B. □

Synchronization:
To prove synchronization in Theorem 4, we only need to demonstrate that each element of each partition P i k , k = 1, . . . , N, is either an individual node, or a cluster of nodes that asymptotically synchronize. Since P i N contains only one element -the set of all nodes -this will prove asymptotic synchronization of all nodes in the network. We will show this by induction.
Induction base: Partition P i 0 , consists of individual nodes.
Induction step: Suppose partition P i k contains either individual nodes or clusters of nodes that asymptotically synchronize. Let us show that P i k+1 also consists of either individual nodes, or cluster of nodes that asymptotically synchronize. The structure of partitions P i k and P i k+1 is specified in (26), Lemma 3. Thus we only need to show that A ∪ B consists of nodes with asymptotically synchronizing states. For simplicity of presentation, we consider the case of the partition P i k consisting of 3 elements: P i k = {A, B, C}. The case of more elements C j in the partition or the case when C is empty can be proven by a minor modification. Since all nodes within A synchronize as well as all nodes within B, it is enough to demonstrate that two nodes a ∈ A and b ∈ B synchronize, i.e. |x a (t) − x b (t)| → 0 as t → +∞. Since V = A ∪ B ∪ C, by the definition of the integral coupling (4) The integral terms marked with * tend to zero, since the corresponding states (and outputs) synchronize (by the assumption that A, B and C are (individually) synchronizing clusters) and since λ(s) is bounded. Denoting the sum of these vanishing elements by ϵ a , we can rewrite u a from (A.9) and, in the same way, u b as u a = |N a ∩ B| ∫ y b ya λ(s)ds + |N a ∩ C| ∫ yc ya λ(s)ds + ϵ a , Therefore, with Γ k = (|N a ∩ B| + |N b ∩ A| + |N a ∩ C|). Notice that this Γ k coincides with the one from Definition 6, since, in our case N a ∩ C = N a ∩ (A ∪ B). In (A.10), we also utilized the fact that |N a ∩ C| = |N b ∩ C|, see condition (27) (29) implies that Γ k λ(s) ≥ Γ λ(s) ≥ γ (s), ∀s ∈ R.
(A.11) Substituting (A.10) into (5) and utilizing inequality (A.11) (note that in our case N a ∩ C = N a \(A ∪ B)), we obtain By the conditions on the function ρ(x), we have .
Therefore, for 1 2 κ(x a − x b ) ≥ |ϵ a − ϵ b | it holds that This fact, together with conditions on function S(x) implies, through an incremental ISS argument (Angeli, 2002), that vanishing ϵ a (t) − ϵ b (t) (as both ϵ a (t) and ϵ b (t) vanish) yields |x a (t) − x b (t)| → 0 as t → +∞. This completes the proof of the induction step that P i k+1 also consists of either individual nodes, or clusters of nodes that asymptotically synchronize, and thus concludes the proof of the theorem. □ seeking automatic control solutions''. His inventions within automatic optimization are implemented in several oil production fields. His current research interests include control of nonlinear systems, control and optimization of highly uncertain systems, automation and digitalization in the energy industry and various industrial applications of automatic control and optimization.  (2012)(2013)(2014)(2015). His core expertise is on the modeling and analysis of collective dynamics -in particular, cluster synchronization and pattern generation -of networked systems, with applications to neuroscience and chemistry. He has co-organized several international conferences. He currently serves as associate editor for Systems and Control Letters, and Communications in Nonlinear Science and Numerical Simulation (from January 2022).  Pavlov and H. Nijmeijer (Birkhauser, 2005) and 'Stability and Convergence of Mechanical Systems with Unilateral Constraints' with R.I. Leine (Springer-Verlag, 2008). In 2015, he received the IEEE Control Systems Technology Award "For the development and application of variable-gain control techniques for high-performance motion systems". He is an IEEE Fellow for his contributions to hybrid, data-based and networked control.