# Weak synchronization and convergence in coupled genetic regulatory networks: Applications to damped oscillators and multistable circuits

Nicolas Augier, Madalena Chaves, Jean-Luc Gouzé

## ► To cite this version:

Nicolas Augier, Madalena Chaves, Jean-Luc Gouzé. Weak synchronization and convergence in coupled genetic regulatory networks: Applications to damped oscillators and multistable circuits. International Journal of Robust and Nonlinear Control, 2023, 33 (9), pp.4867-4892. 10.1002/rnc.6061. hal-03270268

HAL Id: hal-03270268

https://inria.hal.science/hal-03270268

Submitted on 24 Jun 2021

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Weak synchronization and convergence in coupled genetic
regulatory networks: Applications to damped oscillators and
multistable circuits

Nicolas Augier, Madalena Chaves, Jean-Luc Gouzé *

June 13, 2021

Abstract

The study of synchronization in coupled genetic networks is a very challenging topic that is usually analyzed on a case-by-case basis. Here we consider a general model of genetic networks and examine two forms of interconnection, either *homogeneous* or *heterogeneous coupling*, corresponding to coupling functions that are either equal or different from those governing the individual dynamics. In the case of individual subsystems having unique but different steady states, we prove that the homogeneous coupled system has a unique globally asymptotically stable steady state. Moreover, in the case of large coupling strength, we show that under suitable assumptions the network achieves *weak synchronization* in the sense that the individual steady states become arbitrarily close. In the heterogeneous case, stability conditions are more intricate and some stronger assumptions on the individual dynamics have to be made, under which we prove a similar weak synchronization result in the case of large coupling strength. We apply the results to the synchronization of damped oscillators and to the control of multistable systems.

# 1 Introduction

Developing new mathematical tools enabling to understand and control the behavior of genetic reg-
ulatory networks has been an attractive topic for engineers and biologists in the last decades. Two
characteristic types of circuits, positive or negative feedback loops, are frequently found combined in
genetic regulatory networks [6, 28], where they induce different dynamics such as multistability or os-
cillatory behavior. In parallel, many theoretical studies have been devoted to understand and classify
the conditions on the system leading to either of these behaviors, especially in the case of feedback
loops [20, 13, 10, 14]. Some systems may exhibit different dynamics in response to different kinds
of inputs, so a new direction in cellular biology is the possibility of controlling oscillatory behavior:
very efficient control strategies have been recently proposed in order to induce (respectively, suppress)
oscillations in negative loops [7]. A challenging current problem is the study of intercellular commu-
nication within a cell population and the effects of coupling similar genetic regulatory networks on
the global response at the tissue or organ level. Novel experimental observations are pointing towards
the role of intercellular communication in fundamental physiological functions: circadian clocks in
peripheral organs such as the liver may synchronize without a master pacemaker [16]; apoptotic cells
release metabolites that signal across to neighbor cells [21]. Mathematical models are thus required
to analyse the effects of cell-to-cell coupling in a network of genetic regulatory networks [2].

<footnote>*Inria, Université Côte d'Azur, INRAE, CNRS, Sorbonne Université, Biocore team, Sophia Antipolis, France</footnote>
Due to the interactions between the (possibly different) subsystems, new unexpected or undesirable behaviors may occur induced by the coupling mechanisms. For instance, a network of globally asymptotically stable subsystems connected by diffusive coupling may exhibit oscillatory behaviors [25]. In [3], we studied the diffusive coupling of bistable genetic regulatory networks and showed that the system may exhibit emergent behaviors leading to diverse synchronization patterns, similar to those described in [15]. In our recent work [3], as well as in [12, 11], coupled genetic networks are studied using piecewise affine (PWA) models, which approximate continuous sigmoidal-like activity functions by step functions and facilitate an algorithmic analysis of the dynamics. This is an interesting framework since various tools are available for analysis of PWA systems [9] and a practical method for finding Lyapunov functions for PWA systems was also recently proposed in [22].

In this paper, our goal is to improve the current understanding of the emergent behaviors that may appear through coupling, while determining which properties of the individual systems are preserved. We thus propose to focus on the behavior of a class of smoothly coupled monostable genetic regulatory networks where the smoothness hypothesis is made to enable the use of smooth functions and apply the classical Lyapunov theory. Nevertheless we conjecture that the presented results remain valid in the PWA case by singular perturbation arguments similar to those presented in [24].

## 1.1 A general genetic network model with homogeneous or heterogeneous coupling

More precisely, we focus on coupled dynamical systems whose individual dynamics are driven by the equation

$$ \dot{X}(t) = B + AH(X(t)) - \Gamma X(t), \quad (1) $$

where $X(t) \in \mathbb{R}^n$ for $n \ge 1$, where $B \in \mathbb{R}^n$, A is a $n \times n$ real matrix, and a diagonal positive matrix $\Gamma = (\gamma_j)_{j \in \{1, \dots, n\}}$, $(h_j)_{j \in \{1, \dots, n\}}$ are increasing functions from $\mathbb{R}$ to $\mathbb{R}$ which are called internal interaction functions, and for $x = (x_1, \dots, x_n) \in \mathbb{R}^n$, $H(x) = {}^t(h_1(x_1), \dots, h_n(x_n))$. A typical example is the two-dimensional $(-, +)$ damped oscillator, that can be written as Equation (2) with $n = 2$, where the functions $(h_j)_{j \in \{1, 2\}}$ are Hill functions defined, for every $j \in \{1, 2\}$, by $h_j(x) = \frac{x^{s_j}}{x^{s_j} + \theta_j^{s_j}}$ for every $x \in \mathbb{R}$, with $s_j \ge 1$ and $\theta_j > 0$ for every $j \in \{1, 2\}$,

$$ A = \begin{pmatrix} -w_1 & k_1 \\ -k_2 & -w_2 \end{pmatrix}, $$

with $w_1, w_2 \ge 0$, $k_1, k_2 > 0$, and $B = {}^t(0, k_2)$. The variables $x_1$ and $x_2$ represent two interacting proteins such that $x_1$ inhibits $x_2$ while $x_2$ promotes the production of $x_1$. Note that the diagonal terms of A which correspond to self-inhibitions of the nodes are usually taken equal to zero in the standard model (see, for instance [9]), and it is a known fact that such system converges toward a steady state both in both smooth and PWA case, which corresponds to the limit case $s_j \to +\infty$ (see, for instance [13]). Notice that a system driven by Equation (1) and associated coupled systems can be put under the general form studied in [1, Theorem 1]. However, we focus here on the case where a each variable $x_j$ acts on the dynamics of other variables through the action to a single function $h_j(x_j)$. Due to this specific structure, the stability results exposed in this paper are given under easily checkable structural properties of Equation (1) (see Remark 2.4 for more details).

In this paper, we propose to study some specific properties on the coupling of subsystems driven by Equation (1). A similar coupling mechanism has already been proposed in [12] for modelling motor circuitry in Parkinsonian tremor, and it seems to be adapted to the modelling of cell to cell communication (see [2]). We assume that the coupling involves the action of the j-th component of the q-th subsystem $X_{j,q}$ (which represents the j-th protein of the q-th cell) on the j-th component of the k-th subsystem $X_{j,k}$ which is under the form $\alpha_j L_{kq}^j (c_{j,q}(X_{j,q}) - c_{j,k}(X_{j,k}))$ for $k, q \in \{1, \dots, N\}$,
where $(c_{j,q})_{j,q}$ are functions from $\mathbb{R}$ to $\mathbb{R}$, $\alpha_j > 0$ is the *coupling strength* on the variable $j$, $L_j = (L_{kq}^j)_{k,q \in \{1, \dots, N\}}$ is the *coupling matrix* whose structure will be detailed in Section 3. Following the form of the individual dynamics given by Equation (1), the internal interaction functions are increasing real functions $(h_{j,k})_{(j,k) \in \{1, \dots, n\} \times \{1, \dots, N\}}$, and an example of graph representation for two coupled two-dimensional $(-, +)$ damped oscillator is given by Figure 1. Depending on the choice of the functions $c_{j,k}$, we will distinguish two cases, the *homogeneous coupling case* when $c_{j,k} = h_{j,k}$ for every $j, k$ (respectively, *heterogeneous coupling case* when $c_{j,k} \neq h_{j,k}$ for some $j, k$) and we will prove different dynamical results in each case. Typically, the homogeneous case (respectively, heterogeneous) corresponds to the case where cells communicate with the same (respectively, different) thresholds as (than) those which regulate their internal dynamics. In the same spirit, we will study both the case where the subsystems are *parallel*, i.e. if there exists a sequence of functions $(h_j)_{j \in \{1, \dots, n\}}$ such that $h_{j,k} = h_j$ for every $k \in \{1, \dots, N\}$, and the non-parallel case, corresponding to the coupling of subsystems having different activity thresholds then representing the case of similar cells with perturbed parameters.

<figure><img src="image_2.png" /><figcaption>Figure 1: Example of two coupled subsystems, with $n = N = 2$.</figcaption></figure>

In the homogeneous coupling case, we will tackle the coupling of subsystems having different dynamics and different steady states. It turns out that, in this case, the coupling mechanism is such that each subsystem converges globally towards a steady state which depends on the coupling. For a weak coupling (i.e. $\alpha_j$ is very small), they are close to those of the uncoupled system. A natural question is then to study the impact of the variations of the coupling on such steady states, in particular when the coupling is large. Dynamical behaviors of coupled systems in the case of large coupling strength has already been investigated for instance in [23], where the authors have shown that under suitable assumptions on the dynamics (the subsystems consider in the latter are assumed to be identical) and on the coupling (with a slightly different structure than ours) and whose strength is assumed to be large enough, the subsystems synchronize when $t \to +\infty$, in the sense that $|X_{j,q}(t) - X_{j,k}(t)| \to 0$ when $t \to +\infty$, for every $k, q \in \{1, \dots, N\}$. In our case, we will instead provide an expression for the steady states in terms of the system's parameters and prove some of their synchronization properties when the coupling strength $\alpha_j$ tends to $+\infty$, for a very general class of matrices $L_j$. In particular, we prove that for parallel subsystems, the difference between the $j$-th components of the steady states of each subsystems tends to 0 when $\alpha_j \to +\infty$. We extend
those results to the case of non-parallel subsystems, showing that under suitable assumptions on the interaction functions, the j-th components of the steady states of the individual systems converge towards a one-dimensional subspace of $\mathbb{R}^N$. The methods used to prove our convergence results in this case consist in finding suitable Lyapunov functions which take into account the interaction functions $(h_{j,k})_{(j,k)\in\{1,...,n\}\times\{1,...,N\}}$, inspiring from [1].

In the heterogeneous coupling case, the study is more intricate and some stronger assumptions on the individual dynamics have to be made in order to prove similar results. When the coupled individual dynamics are parallel and possess the same steady states, we prove that the subsystems reach synchronization by considering a well chosen Lyapunov function which depends only on the internal interaction functions. In the latter case we study weak synchronization properties, showing similar properties to those proved in the homogeneous case, under an a priori global stability assumption.

We apply our results to the synchronization and to the stabilization of two types of famous genetic circuits, namely the two-dimensional $(-, +)$ damped oscillator and the toggle switch.

## 1.2 Outline

In Section 2, we study the individual dynamics and we provide a sufficient condition on A which guarantees the convergence towards a steady state. Section 3 introduces the coupling framework and notation and then Section 4 focuses on the homogeneous coupling case. In this case, we can prove the asymptotic convergence of each subsystem towards possibly different steady states, and we propose to study the large coupling case which provides some interesting weak synchronization properties. In Section 5, we study the heterogeneous coupling case and we prove that, under some technical hypothesis, the subsystems synchronize when they possess the same steady states, whatever the coupling strength may be. Then we study weak synchronization properties in the case where the subsystems have different steady states, making an a priori global stability assumption. In Section 6, we present some applications to the coupling of different two-dimensional damped oscillators and to the control of some multistable systems. In particular we propose a control strategy ensuring the convergence coupled bistable systems towards their unstable (also known as "undifferentiated") steady states.

## 2 Single system

Consider $B \in \mathbb{R}^n$, a $n \times n$ matrix A, and a diagonal positive matrix $\Gamma = (\gamma_j)_{j \in \{1,...,n\}}$, and let $(h_j)_{j \in \{1,...,n\}}$ be a family of strictly increasing real functions. For $x = (x_1, ..., x_n) \in \mathbb{R}^n$, define $H(x) = (h_1(x_1), ..., h_n(x_n))$, and consider the equation

$$ \dot{X}(t) = B + AH(X(t)) - \Gamma X(t), \quad (2) $$

where $X(t) \in \mathbb{R}^n$. In the following we will assume that the cone $(\mathbb{R}_+)^n$ is invariant by the dynamics of Equation (2), which is of utter importance for biological applications. This is satisfied for instance when the functions $(h_j)_{j \in \{1,...,n\}}$ are non-negative, $h_j(0) = 0$ for every $j \in \{1,...,n\}$ and B and A are such that for every $j \in \{1,...,n\}$, if $x = (x_1, ..., x_n) \ge 0$ is such that $x_j = 0$, then we have $(Ax + B)_j \ge 0$.

Denote the set of real functions which are locally Lipzchitz, bounded and increasing by $\mathcal{R}$. Concerning the existence of a steady state for Equation (2) and the definition of the solutions of Equation (2), we have the following fact.

**Lemma 2.1.** Assume that the functions $(h_j)_{j \in \{1,...,n\}}$ belong to $\mathcal{R}$. Then Equation (2) admits a steady state and the solutions of Equation (2) are uniquely globally defined.
**Proof.** We obtain easily that Equation (2) admits a steady state $\bar{X}$ if and only if the application $F: X \mapsto \Gamma^{-1}(B + AH(X))$ admits a fixed point. Under the assumption of the lemma, $X \mapsto \Gamma^{-1}(B + AH(X))$ is a continuous bounded function from $\mathbb{R}^n$ in itself. Hence $\bar{F}(\mathbb{R}^n)$, is a compact invariant set for $F$. As the functions $(h_j)_j$ are increasing, we deduce that $H(\mathbb{R}^n)$ is the direct product of real compact real intervals, and is convex. We can deduce by multiplication by the matrix $\Gamma^{-1}A$ that $F(\mathbb{R}^n)$ is convex, and this property also holds for $\bar{F}(\mathbb{R}^n)$. Hence Brouwer's Theorem (see e.g. [18]) proves the existence of a fixed point for $F$.

Concerning the definition of solutions of Equation (2), Cauchy-Lipzchitz theorem proves the existence and uniqueness of maximal solutions of Equation (2). Moreover, the boundedness of $H$ guarantees the existence of $\alpha, \beta > 0$ such that the vector field $f: X \mapsto B + AH(X) - \Gamma X$ satisfies $\|f(X)\| \le \alpha\|X\| + \beta$ for every $X \in \mathbb{R}^n$. The global definition of solutions of Equation (2) then follows by a Gronwall type argument. □

**Definition 2.2.** We say that $A$ is weakly diagonally stable if there exists a diagonal positive matrix $D = \text{diag}(d_j)_{j \in \{1, \dots, n\}}$ such that $DA + {}^tAD \le 0$, in the sense that ${}^tX(DA + {}^tAD)x \le 0$ for every $x \in \mathbb{R}^n$.

Now we prove Proposition 2.3, which is obtained under the assumption that the interaction matrix $A$ of Equation (2) is weakly diagonally stable.

**Proposition 2.3.** Assume that $A$ is weakly diagonally stable and that the functions $(h_j)_j$ belong to $\mathcal{R}$. Then Equation (2) admits a unique globally asymptotically stable steady state, which belongs to $(\mathbb{R}_+)^n$.

**Proof.** By Lemma 2.1, Equation (2) admits a steady state $\bar{X} = (\bar{X}_1, \dots, \bar{X}_n) \in \mathbb{R}^n$. Define $V_j(x) = \int_{\bar{X}_j}^x (h_j(z) - h_j(\bar{X}_j))dz$, and $V(X) = \sum_{j=1}^n d_j V_j(X_j)$. Let us show that $V$ is a Lyapunov function for Equation (2). First notice that since $h_j$ is increasing for every $j \in \{1, \dots, n\}$, each $V_j$ is non-negative and $V_j(x) = 0$ if and only if $x = \bar{X}_j$, hence the same property holds for $V$. By definition of $\bar{X}$, we have

$$ \dot{X} = A(H(X) - H(\bar{X})) - \Gamma(X - \bar{X}). $$

Hence we get

$$ \begin{aligned} \sum_{j=1}^{n} d_j \frac{d}{dt} V_j(X_j(t)) &= \sum_{j=1}^{n} d_j \dot{X}_j(t) (h_j(X_j(t)) - h_j(\bar{X}_j)) \\ &= {}^t(H(X) - H(\bar{X}))D\dot{X} \\ &= {}^t(H(X) - H(\bar{X}))DA(H(X) - H(\bar{X})) - {}^t(H(X) - H(\bar{X}))D\Gamma(X - \bar{X}). \end{aligned} $$

The diagonal stability of $A$ implies ${}^t(H(X) - H(\bar{X}))DA(H(X) - H(\bar{X})) \le 0$. The strict positivity of the diagonal matrices $D$ and $\Gamma$ and the fact that the functions $(h_j)_j$ are strictly increasing imply ${}^t(H(X) - H(\bar{X}))D\Gamma(X - \bar{X}) \ge 0$. It follows that $\frac{d}{dt}V(X(t)) \le 0$ for every $t \ge 0$. Moreover, we have ${}^t(H(X) - H(\bar{X}))D\Gamma(X - \bar{X}) = 0$ if and only if $X = \bar{X}$, because $D\Gamma$ is diagonal. Hence $V$ is a Lyapunov function for Equation (2) such that $\frac{dV}{dt}(X(t)) = 0$ implies $X(t) = \bar{X}$. The fact that the functions $(h_j)_{j \in \{1, \dots, n\}}$ belong to $\mathcal{R}$ guarantees that $V$ is radially unbounded, so that the convergence result is proved, and the assumption of forward invariance of $(\mathbb{R}_+)^n$ by the dynamics of Equation (2) ensures that $\bar{X} \in (\mathbb{R}_+)^n$. □
**Remark 2.4.** • As already mentioned in the introduction, a system driven by Equation (2) can be put under the form studied in [1, Theorem 1]. In the latter, the authors propose a very general criterion ensuring global stability, via the study of a higher dimensional dissipativity matrix associated with the interactions between the different components of the dynamics. Because of the difficulty to check such criterion in practice, and due to the specificity of the interactions occurring in Equation (2) (indeed, in our case, the action a given node j on a node k is governed by the function h_j, independently of k), we propose a more specific sufficient criterion of stability. Our criterion is easier to handle, and it turns out that it can be extended to the case of interconnected systems, as we will see in the next sections. Moreover, our framework allows to take into account self interactions, without making restrictions on the slopes of the functions (h_j)_{j \in \{1,...,n\}}, nor on the degradation rates $\gamma_j > 0$, for $j \in \{1,...,n\}$, which can be both considered as uncertain parameters of the system. Using Hill functions as in [7], defined, for every $j \in \{1,...,n\}$, by $h_j(x) = \frac{x^{s_j}}{x^{s_j}+\theta_j^{s_j}}$ for every $x \in \mathbb{R}$, with $s_j \ge 2$ where $\theta_j > 0$ for every $j \in \{1,...,n\}$, under the assumptions of Proposition 2.3, Equation (2) has a globally stable steady state (which possibly depends on $s_j \ge 2$ and on the other parameters of the system), whatever the choices of $(s_j)_{j \in \{1,...,n\}}$ and $(\gamma_j)_{j \in \{1,...,n\}}$ may be.

* Notice that the Lyapunov function used in the proof of Proposition 2.3 is the sum of a linear term w.r.t. the state x and a non-linear term involving the integrals $\int_{\bar{X}_j}^x h_j(z)dz$, and depends on $\bar{X}_j$ which depends itself strongly on $(A, B)$. Note that the computation of $\bar{X}$ can be tricky in practice, provided that fixed point or Newton type methods' precision is deteriorated for large $s_j \ge 1$.

# 3 Coupled systems

## 3.1 Individual models

In this section we consider, for $N \ge 1$, the coupling of N-subsystems whose individual dynamics is given by Equation (2), associated with different $n \times n$ matrices $(A_k)_{k \in \{1,...,N\}}$, different vectors $(B^k)_{k \in \{1,...,N\}} = (b_{jk})_{j \in \{1,...,n\}}$ of $\mathbb{R}^n$, and different functions $(h_{j,k})_{(j,k) \in \{1,...,n\} \times \{1,...,N\}}$, called internal interaction functions, belonging to $\mathcal{R}$. For the sake of notations, set $(B_j)_{j \in \{1,...,n\}} = (b_{jk})_{k \in \{1,...,N\}}$. Define, for every $k \in \{1,...,N\}$, a positive $n \times n$ diagonal matrix $\Gamma^k = \text{diag}((\gamma_{jk})_{j \in \{1,...,n\}})$. We will additionally define for $j \in \{1,...,n\}$, the $N \times N$ diagonal matrix $\Gamma_j = (\gamma_{jk})_{k \in \{1,...,N\}}$. In this setting, the individual dynamics $X^k(t) \in \mathbb{R}^n$ of the k-th subsystem reads

$$ \dot{X}^k(t) = B^k + A_k H^k(X^k(t)) - \Gamma^k X^k(t), \quad (3) $$

where for $x = (x_1, ..., x_n) \in \mathbb{R}^n$, $H^k(x) = {}^t(h_{1,k}(x_1), ..., h_{n,k}(x_n))$.

The matrices $(A_k)_{k \in \{1,...,N\}}$ are assumed to be simultaneously weakly diagonally stable, in the following sense.

**Definition 3.1.** We say that the $n \times n$ matrices $(A_k)_{k \in \{1,...,N\}}$ are simultaneously weakly diagonally stable if there exists a positive diagonal matrix $D = \text{diag}(d_j)_{j \in \{1,...,n\}}$ such that $DA_k + {}^t A_k D \le 0$ for every $k \in \{1,...,N\}$.

For the sake of notations, for every $k \in \{1,...,N\}$, denote the $(j,l)$-th coefficient of $A_k$ by $A_{jl}^k$.
## 3.2 The coupling model

### 3.2.1 Qualitative interpretation

We assume that the coupling term involves the action of the j-th variable of the q-th subsystem on the k-th subsystem which is under the form $d(X_{j,q}, X_{j,k}) = c_{j,q}(X_{j,q}) - c_{j,k}(X_{j,k})$ for $k, q \in \{1, \dots, N\}$, where $(c_{j,q})_{j,q} \text{ are strictly increasing functions from } \mathbb{R} \text{ to } \mathbb{R}, \text{ called coupling functions.}$ An example of such coupled system is given on Figure 1 in the case $n = N = 2$, on which the positive interactions are given by arrows and negative interactions by bars. Such model can be used to model the linear coupling case when the $(c_{j,k})_{(j,k)\in\{1,\dots,n\}\times\{1,\dots,N\}}$ are linear or the following case, that will be called on/off diffusion. Assuming that $c_{j,k}(x) = c_j(x) = \frac{\theta^s}{x^{s+\theta^s}}$ for every $k \in \{1, \dots, N\}$ is a Hill function with a large enough slope $s \in \mathbb{N}$ and a threshold $\theta \in \mathbb{R}$, we have

$$d(X_{j,q}, X_{j,k}) \approx \begin{cases} 0 & \text{if } (X_{j,q} > \theta \text{ and } X_{j,k} > \theta) \text{ or } (X_{j,q} \le \theta \text{ and } X_{j,k} \le \theta) \\ 1 & \text{if } (X_{j,q} > \theta \text{ and } X_{j,k} \le \theta) \\ -1 & \text{if } (X_{j,q} \le \theta \text{ and } X_{j,k} > \theta) \end{cases} \quad (4)$$

In other words, the diffusion mechanism activates with discrete values, and the cells decide to exchange only when they both possess a large enough concentration of a given protein. This setting is slightly different from the case $d(X_{j,q}, X_{j,k}) = c_j(X_{j,q} - X_{j,k})$, which appears more naturally when studying diffusion processes [11, 3, 26, 27]. Indeed, in our case, two coupled cells having a given protein concentration less than $\theta$ cannot cooperate in order to increase concentrations, so that they evolve independently. The dynamical study of a very similar model has been initiated in [12], in which the authors, motivated by applications in neuronal circuitry, focus on systems having a sustained oscillatory behavior.

### 3.2.2 Mathematical framework

We assume that the coupling functions $(c_{j,k})_{(j,k)\in\{1,\dots,n\}\times\{1,\dots,N\}}$ are strictly increasing functions from $\mathbb{R}$ to $\mathbb{R}$. Set a subset $C$ of $\{1, \dots, n\}$ and let $(L_j)_{j \in C \subseteq \{1, \dots, n\}}$ be a family of irreducible $N \times N$ matrices called coupling matrices, such that for every $j \in C$, $L_j$ is a non-positive off diagonal matrix such that $L_j \mathbf{1} = 0$ and $^tL_j \mathbf{1} = 0$, where $\mathbf{1}$ is the vector of $\mathbb{R}^N$ whose every components are equal to 1. For $q, k \in \{1, \dots, N\}$, denote the $(q, k)$-th coefficient of $L_j$ by $L_{qk}^j$. Notice that the coupling matrices are not necessarily symmetric, and the set $C$ corresponds to the set of variables which are subject to the coupling mechanism. By simple considerations, we show the existence of $\beta > 0$ such that $Z = \frac{1}{\beta}L_j + \text{Id}$ is a non-negative irreducible doubly stochastic matrix, and by a direct application the Perron Frobenius Theorem [5, Theorem 1.4] to $Z$, we obtain that $\text{Ker}(L_j) = \text{span}(\mathbf{1})$ for every $j \in C$.

**Coupled Dynamics:** For $X = (X_{j,k})_{(j,k)\in\{1,\dots,n\}\times\{1,\dots,N\}} \in \mathbb{R}^{nN}$, define the k-th subsystem $X^k = (X_{1,k}, \dots, X_{n,k}) \in \mathbb{R}^n$, and the vector of the j-th components $X_j = (X_{j,1}, \dots, X_{j,N}) \in \mathbb{R}^N$. Define the vectors $H_j(X_j) \in \mathbb{R}^N$, $H^k(X^k) \in \mathbb{R}^n$ and $H(X) \in \mathbb{R}^{nN}$, which are related to the interaction functions $(h_{jk})_{(j,k)\in\{1,\dots,n\}\times\{1,\dots,N\}}$ as

$$\begin{aligned} H_j(X_j) &= {}^t(h_{j,1}(X_{j,1}), \dots, h_{j,N}(X_{j,N})), \\ H^k(X^k) &= {}^t(h_{1,k}(X_{1,k}), \dots, h_{n,k}(X_{n,k})), \\ H(X) &= {}^t(H^1(X^1), \dots, H^k(X^k)). \end{aligned}$$
Similarly, define

$$
\begin{align*}
C_j(X_j) &= {}^t(c_{j,1}(X_{j,1}), \dots, c_{j,N}(X_{j,N})), \\
C^k(X^k) &= {}^t(c_{1,k}(X_{1,k}), \dots, c_{n,k}(X_{n,k})), \\
C(X) &= {}^t({}^tC^1(X^1), \dots, {}^tC^k(X^k)).
\end{align*}
$$

For $X = (X_{j,k})_{(j,k) \in \{1, \dots, n\} \times \{1, \dots, N\}}$, define

$$
Q_j(X) = {}^t(\sum_{q=1}^{n} A_{jq}^1 h_{q,1}(X_{q,1}), \dots, \sum_{q=1}^{n} A_{jq}^N h_{q,N}(X_{q,N})),
$$

and consider the equation

$$
\dot{X}_j = B_j + Q_j(X) - \Gamma_j X_j - L_j C_j(X_j), \quad (5)
$$

where $X_j(t) \in \mathbb{R}^N$ for every $j \in \{1, \dots, n\}$. Writing the system component by component, Equation (5) can be expressed, for every $(j, k) \in \{1, \dots, n\} \times \{1, \dots, N\}$, as

$$
\dot{X}_{j,k}(t) = b_{jk} - \gamma_{jk}X_{j,k} + \sum_{q=1}^{n} A_{jq}^{k}h_{q,k}(X_{q,k}) - \sum_{q=1}^{N} L_{qk}^{j}(c_{j,q}(X_{j,q}) - c_{j,k}(X_{j,k})). \quad (6)
$$

Notice that Equation (5) can be written as

$$
\dot{X} = B + \tilde{A}H(X) - \Gamma X - LC(X), \quad (7)
$$

where $X(t) \in \mathbb{R}^{nN}$, where $\tilde{A}$ and $L$ are $nN \times nN$ matrices defined as

$$
\tilde{A} = \begin{pmatrix}
A_1 & 0 & \dots & 0 \\
0 & A_2 & \ddots & \vdots \\
\vdots & \ddots & \ddots & 0 \\
0 & \dots & 0 & A_N
\end{pmatrix},
$$

whose dimension $n$ blocks are equal to $A_k$ for $k \in \{1, \dots, N\}$, and the matrix $L$ is the $nN \times nN$
block matrix

$$
L = \begin{pmatrix}
Z_{1,1} & Z_{1,2} & \cdots & Z_{1,N} \\
Z_{2,1} & Z_{2,2} & \cdots & \vdots \\
\vdots & \vdots & \ddots & Z_{N-1,N} \\
Z_{N,1} & \cdots & Z_{N,N-1} & Z_{N,N},
\end{pmatrix} \quad (8)
$$

constituted of $N^2$ blocks of dimension $n$ such that for $i,j \in \{1, \dots, N\}$, the $(i,j)$-th block is defined
as $Z_{i,j}$, where $Z_{i,j} = \text{diag}(L_{ij}^k)_{k \in \{1, \dots, n\}}$. Moreover, $\Gamma$ is a block diagonal matrix constituted of a
number $N$ of $n \times n$ positive diagonal blocks $\Gamma^k$ for $k \in \{1, \dots, N\}$, and $B \in \mathbb{R}^{nN}$ is the concatenation
of the vectors $B_j \in \mathbb{R}^N$ for $j \in \{1, \dots, n\}$. Hence one can easily prove that the assumption of
simultaneous weak diagonal stability of the matrices $(A_k)_{k \in \{1, \dots, N\}}$ imply that $A$ is weakly diagonally
stable, by using a $N$-blocks concatenation $\hat{D}$ of a diagonal positive matrix $D = \text{diag}(d_j)_{j \in \{1, \dots, n\}}$
such that $DA_k + {}^t A_k D \preceq 0$ for every $k \in \{1, \dots, N\}$. Note that the assumption that the functions
$(c_{j,k})_{(j,k) \in \{1, \dots, n\} \times \{1, \dots, N\}}$ are strictly increasing and the assumptions made of the matrices $L_j$ imply
that, if the cone $(\mathbb{R}^+)^n$ is forward invariant by the individual dynamics (3) associated with $(A_k, B^k)$ for
every $k \in \{1, \dots, N\}$, then the cone $((\mathbb{R}^+)^n)^N$ is forward invariant by the dynamics of Equation (7).

**Definition 3.2.** • System (7) is said to be homogeneous when $c_{j,k} = h_{j,k}$ for every $j, k$ (respectively, heterogeneous when $c_{j,k} \neq h_{j,k}$ for some $j, k$).

* System (7) is said to be parallel if there exists a sequence of functions $(h_j)_{j \in \{1, \dots, n\}}$ belonging to $\mathcal{R}$ such that $h_{j,k} = h_j$ for every $k \in \{1, \dots, N\}$.
### 3.2.3 Synchronization properties

In the following we will denote the Euclidean distance between $x \in \mathbb{R}^N$ and a subspace $E$ of $\mathbb{R}^N$ by $d(x, E)$.

**Definition 3.3 (Componentwise synchronization).** Set $j \in \{1, \dots, n\}$. We say that the solution $X(t)$ of Equation (7) achieves $j$-synchronization if $|X_{j,q}(t) - X_{j,k}(t)| \to 0$ when $t \to +\infty$, for every $k, q \in \{1, \dots, N\}$. In other words, this definition is equivalent to say that $d(X_j(t), \operatorname{span}(\mathbb{1})) \to 0$ when $t \to +\infty$.

**Definition 3.4 (Synchronization).** We say that the solution $X(t)$ of Equation (7) achieves synchronization if $|X_{j,q}(t) - X_{j,k}(t)| \to 0$ when $t \to +\infty$, for every $k, q \in \{1, \dots, N\}$ and every $j \in \{1, \dots, n\}$. In other words, this definition is equivalent to say that $d(X_j(t), \operatorname{span}(\mathbb{1})) \to 0$ when $t \to +\infty$, for every $j \in \{1, \dots, n\}$.

Assume that the coupling matrix $L \equiv L(\alpha)$ defined as in Section 3.2.2 depends on a parameter $\alpha > 0$ and for every $\alpha > 0$, denote the steady state of Equation (7) by $\bar{X}_\alpha \in \mathbb{R}^{nN}$.

**Definition 3.5 (Componentwise weak synchronization).** Set $j \in \{1, \dots, n\}$. We say that Equation (7) achieves $j$-weak synchronization if $\bar{X}_\alpha$ is asymptotically globally stable for Equation (7) and $d((\bar{X}^\alpha)_j, \operatorname{span}(\mathbb{1})) \to 0$ when $\alpha \to +\infty$, for every $k, q \in \{1, \dots, N\}$.

**Definition 3.6 (Weak synchronization).** We say that Equation (7) achieves weak synchronization if it achieves $j$-weak synchronization for every $j \in \{1, \dots, n\}$.

We have the following direct property concerning the dynamics of Equation (7).

**Proposition 3.7.** Assume that Equation (7) achieves $j$-weak synchronization for a given $j \in \{1, \dots, n\}$ (respectively, weak synchronization). Then, for an arbitrary $\eta > 0$, there exists $\alpha_0 > 0$ and $t_0 = t_0(\alpha_0) > 0$ such that for every $t \ge t_0$, we have $d((X_{\alpha_0})_j(t), \operatorname{Span}(\mathbb{1})) < \eta$ (respectively, $d((X_{\alpha_0})_j(t), \operatorname{Span}(\mathbb{1})) < \eta$ for every $j \in \{1, \dots, n\}$).

## 4 Homogeneous coupling results

We propose to study the homogeneous coupling case as defined in Definition 3.2, where the coupling functions coincide with the internal interaction functions, that is $c_{j,k} = h_{j,k}$ for every $(j, k) \in \{1, \dots, n\} \times \{1, \dots, N\}$. In this case, Equation (7) can be written as

$$ \dot{X} = B + (\tilde{A} - L)H(X) - \Gamma X, \quad (9) $$

where $X(t) \in \mathbb{R}^{nN}$, $\tilde{A}$, $L$, $\Gamma$, $B$ and $H$ are defined as in Section 3.2.2. The latter equation has the form of Equation (2) in dimension $nN$, with $A = \tilde{A} - L$.

This will allow us to prove Theorem 4.1 ensuring that, when the assumptions of Section 3 hold true, each subsystem $X^k$ for $k \in \{1, \dots, N\}$ converges towards a point $\bar{X}^k$ of $\mathbb{R}^n$ when $t \to +\infty$, where $\bar{X} = (\bar{X}^1, \dots, \bar{X}^N)$ is the globally attractive steady state of the corresponding Equation (9) in $\mathbb{R}^{nN}$. We prove synchronization towards a steady state in the sense of Definition 3.4 in the case where $\bar{X}^k = \bar{X}^q$ for every $k, q \in \{1, \dots, N\}$. Then we study some weak synchronization properties in the sense of Definition 3.5. when the coupling strength tends to $+\infty$.
## 4.1 Convergence result

First, one can easily prove that A is weakly diagonally stable using a N-blocks concatenation $\hat{D}$ of the matrix D as in Section 3.2.2. Indeed, we have, for every $X \in \mathbb{R}^{nN}$,

$$ 
\begin{aligned}
{}^t X \hat{D} A X &= {}^t X \hat{D} \tilde{A} X + {}^t X \hat{D} L X \\
&= \sum_{k=1}^{N} {}^t X^k D A_k X^k - \sum_{j=1}^{n} d_j {}^t X_j L_j X_j \\
&\le 0.
\end{aligned}
$$

A direct application of Lemma 2.1 and Proposition 2.3 ensures the existence of a unique globally asymptotically stable steady state $\bar{X}$ for Equation (9).

**Theorem 4.1.** Assume that $((A_k)_{k \in \{1,...,N\}}, (h_{jk})_{(j,k) \in \{1,...,n\} \times \{1,...,N\}})$ are such that the matrices $(A_k)_{k \in \{1,...,N\}}$ are simultaneously weakly diagonally stable and the functions $(h_{jk})_{(j,k) \in \{1,...,n\} \times \{1,...,N\}}$ belong to $\mathcal{R}$. Then Equation (9) has a globally asymptotically stable steady state $\bar{X} \in ((\mathbb{R}^+)^n)^N$.

In particular, if $\bar{X}^k = \bar{X}^q$ for every $k, q \in \{1, ..., N\}$, then Equation (9) achieves synchronization in the sense of Definition 3.4.

**Remark 4.2.** As in the proof of Proposition 2.3, a Lyapunov function can be built explicitely for Equation (9). For this, consider a positive diagonal matrix $D = \text{diag}(d_j)_{j \in \{1,...,n\}}$ such that $DA_k + {}^t A_k D \le 0$, for every $k \in \{1, ..., N\}$, define, for $x \in \mathbb{R}$, $V_{j,k}(x) = \int_{\bar{x}_{j,k}}^x (h_{jk}(z) - h_{jk}(\bar{X}_{j,k}))dz$, and for $X \in \mathbb{R}^{nN}$, set $V(X) = \sum_{j=1}^n d_j \sum_{k=1}^N V_{j,k}(X_{j,k})$. Then we can show that V is a globally positive definite Lyapunov function for Equation (9), which is radially unbounded and whose time derivative is globally negative definite.

## 4.2 Weak synchronization results

In the case of coupling of subsystems whose individual dynamics have different steady states, Theorem 4.1 proves that the subsystems converge asymptotically towards different points of $\mathbb{R}^n$ when $t \to +\infty$. In this section, we set $j \in \{1, ..., n\}$ such that $L_j \ne 0$, and we consider Equation (9) where $L_j$ is replaced by $\alpha_j L_j$, where $\alpha_j > 0$ corresponds to the coupling strength associated with the j-th component, and we prove a j-weak synchronization result in Theorem 4.7. In this case, we can decompose the matrix L of Equation (8) as $L(\alpha_j) = \tilde{L} + \alpha_j \tilde{L}_j$, where $\tilde{L}$ is the block matrix such that the $(q, l)$-th $n \times n$ block is $\text{diag}(L_{ql}^k(1 - \delta_{jk}))_{k \in \{1,...,n\}}$, for every $q, l \in \{1, ..., N\}$, where $\delta_{jk} = 1$ is equal to 1 if $k = j$ and $\delta_{jk} = 0$ if $k \ne j$. The matrix $\tilde{L}_j$ is the block matrix constituted of the $n \times n$ blocks $\text{diag}(L_{ql}^k \delta_{jk})_{k \in \{1,...,n\}}$.

### 4.2.1 Asymptotic behavior of the steady states

By the computations made in the proof of Lemma 2.1, we have that $\bar{X} \in \mathbb{R}^{nN}$ is a steady state of Equation (9) if and only if $\bar{X} = \Gamma^{-1}(B + AH(\bar{X})) = \Gamma^{-1}(B + (\tilde{A} - \tilde{L} - \alpha_j \tilde{L}_j)H(\bar{X}))$. For every $\alpha_j > 0$, Proposition 2.3 and Lemma 2.1 guarantee the existence of a unique steady state $\bar{X}_{\alpha_j}$ solution of

$$ \Gamma^{-1} (B + (\tilde{A} - \tilde{L} - \alpha_j \tilde{L}_j) H(\bar{X}_{\alpha_j})) - \bar{X}_{\alpha_j} = 0. $$
Setting $M = \tilde{A} - \tilde{L}$, up to a change $\epsilon = \frac{1}{\alpha_j} > 0$, we propose to study the asymptotic behavior when $\epsilon \to 0$, $\epsilon > 0$ of the solution $\bar{X}_\epsilon$ of

$$ \epsilon\Gamma^{-1}B + \epsilon\Gamma^{-1}MH(\bar{X}_\epsilon) - \epsilon\bar{X}_\epsilon - \Gamma^{-1}\tilde{L}_jH(\bar{X}_\epsilon) = 0. \quad (10) $$

Lemma 4.4 is the key technical result of this section which will allow us to prove the weak synchronization result of Proposition 4.5. Its proof requires the application of an inequality given by next Lemma that can be found in [17], and which is proved in Appendix A for completeness.

**Lemma 4.3.** Let $K = (K_{ij})_{(i,j)\in\{1,...,N\}}$ be a $N \times N$ off diagonal non-negative matrix such that $K\mathbf{1} = \mathbf{1}$, $^tK\mathbf{1} = \mathbf{1}$, and let $H$ be defined for $y = (y_1, ..., y_N) \in \mathbb{R}^N$ by $\mathcal{H}(y) = (h(y_1), ..., h(y_N))$, where $h$ is an increasing real function. Then we have, for every $y \in \mathbb{R}^N$, $^t\mathcal{H}(y)Ky \le 0$.

### 4.2.2 Weak synchronization for parallel subsystems

Assume in this section that the system is parallel in the sense of Definition 3.2, that is $h_{j,k} = h_j$ for every $(j,k) \in \{1,...,n\} \times \{1,...,N\}$, where the functions $(h_j)_{j \in \{1,...,n\}}$ belong to $\mathcal{R}$.

**Lemma 4.4.** The sequence $(\bar{X}_\epsilon)_{\epsilon>0}$ is bounded.

*Proof.* By Equation (10), we obtain

$$ \epsilon\|\bar{X}_\epsilon\|^2 = \epsilon\langle\Gamma^{-1}B, \bar{X}_\epsilon\rangle + \epsilon\langle\Gamma^{-1}MH(\bar{X}_\epsilon), \bar{X}_\epsilon\rangle - \langle\Gamma^{-1}\tilde{L}_jH(\bar{X}_\epsilon), \bar{X}_\epsilon\rangle. $$

The matrix $\Gamma^{-1}\tilde{L}_j$ is a $nN \times nN$ block matrix such that for $q, l \in \{1,...,N\}$, the $(q,l)$-th $n \times n$ block is defined as $\Gamma_q^{-1}\mathrm{diag}(L_{ql}^k\delta_{jk})_{k \in \{1,...,n\}}$. Hence we can deduce that $-\Gamma^{-1}\tilde{L}_j$ is a non-negative off diagonal matrix such that the sum of each rows and columns is equal to zero. By applying Lemma 4.3, we obtain $-\langle\Gamma^{-1}\tilde{L}_jH(\bar{X}_\epsilon), \bar{X}_\epsilon\rangle \le 0$. It follows that

$$ \|\bar{X}_\epsilon\|^2 \le \langle\Gamma^{-1}B, \bar{X}_\epsilon\rangle + \langle\Gamma^{-1}MH(\bar{X}_\epsilon), \bar{X}_\epsilon\rangle. $$

We get the result using the Cauchy-Schwarz inequality and the fact that $H$ is bounded. $\square$

As a direct consequence of Lemma 4.4, we obtain $\Gamma^{-1}\tilde{L}_jH(\bar{X}_\epsilon) \to 0$ by taking the limit $\epsilon \to 0$ in Equation (10). It follows that $\tilde{L}_jH(\bar{X}_\epsilon) \to 0$, when $\epsilon \to 0$. Using the fact that for $j \in \{1,...,n\}$ such that $L_j \ne 0$ we have $\mathrm{Ker}(L_j) = \mathrm{span}(\mathbf{1})$, we can deduce that for every $j \in \{1,...,n\}$, $d(H_j((\bar{X}_\epsilon)_j), \mathrm{span}(\mathbf{1})) \to 0$ when $\epsilon \to 0$, where $d$ is defined in Section 3.2.3.

**Proposition 4.5.** We have $d((\bar{X}_\epsilon)_j, \mathrm{span}(\mathbf{1})) \to 0$ when $\epsilon \to 0$.

*Proof.* The proof relies on the following result:

* Let $f$ be a strictly increasing continuous real function, and consider two bounded real sequences $(x_n)_{n \ge 0}$ and $(y_n)_{n \ge 0}$ such that $f(x_n) - f(y_n) \to 0$ when $n \to +\infty$. Then we have $x_n - y_n \to 0$ when $n \to +\infty$.

Up to replacing $x_n$ by $\tilde{x}_n = \max(x_n, y_n)$ and $y_n$ by $\tilde{y}_n = \min(x_n, y_n)$ for every $n \ge 0$, one can assume that $x_n \ge y_n$ for every $n \ge 0$. By contradiction, assume that there exists $\eta > 0$ such that for every $n_0 \ge 0$, there exists $n \ge n_0$ such that $|x_n - y_n| > \eta$. We can extract a subsequence $\sigma: \mathbb{N} \to \mathbb{N}$ such that $|x_\sigma(n) - y_\sigma(n)| > \eta$ for every $n \in \mathbb{N}$. By the variations of $f$, we obtain $f(x_\sigma(n)) - f(y_\sigma(n)) > f(y_\sigma(n) + \eta) - f(y_\sigma(n)) > 0$. Using the fact that $(y_\sigma(n))_n$ is a bounded real sequence, we can extract a converging subsequence $y_{\tilde{\sigma}\circ\sigma(n)}$ converging to $\bar{y} \in \mathbb{R}$ when $n \to +\infty$. The function $f$ being strictly increasing, we deduce the existence of $m > 0$ such that $f(x_{\tilde{\sigma}\circ\sigma(n)}) - f(y_{\tilde{\sigma}\circ\sigma(n)}) > m > 0$, for every $n \ge 0$. Hence we get a contradiction.

Using the condition $d(H_j((\bar{X}_\epsilon)_j), \mathrm{span}(\mathbf{1})) \to 0$ when $\epsilon \to 0$, a direct application of the latter result with $f = h_j$ proves the claim of the proposition. $\square$
**Remark 4.6.** Proposition 4.5 proves that $d((\bar{X}_\epsilon)_j, \operatorname{span}(\mathbb{1})) \to 0$, which guarantees that $\bar{X}_\epsilon$ has asymptotically the same j-th components when $\epsilon \to 0$. However there is no guarantee that $(\bar{X}_\epsilon)_\epsilon$ converges as $\epsilon \to 0$.

In particular, a direct consequence of this fact associated with Theorem 4.1 is the following weak synchronization result, in the sense of Definition 3.5.

**Theorem 4.7.** Equation (9) achieves j-weak synchronization when $\alpha_j \to +\infty$.

**Remark 4.8.** The same results are valid if instead of considering a unique $j \in \{1, \dots, n\}$ such that $L_j \neq 0$ and replacing $L_j$ by $\alpha_j L_j$, one would consider a subset $C$ of different $j \in \{1, \dots, n\}$ such that $L_j \neq 0$, and replacing every $L_j$ for $j \in C$ by $\alpha L_j$, where $\alpha > 0$ does not depend on $j$, then taking the limit $\alpha \to +\infty$. In this case, Equation (9) achieves weak synchronization (on all components) in the sense of Definition 3.6.

### 4.2.3 Non-parallel subsystems having shifted interaction functions

Now we propose a direct extension of Theorem 4.7 to the case of subsystems possessing shifted interaction functions $(h_{j,k})_{j,k}$ in the following sense. Consider a sequence $(h_j)_{j \in \{1, \dots, n\}}$ of functions belonging to $\mathcal{R}$, and assume that **for every** $k \in \{1, \dots, N\}$, $h_{j,k}(x) = h_j(\tau_{j,k}x)$, **for every** $x \in \mathbb{R}$, **with a given** $\tau_{j,k} > 0$. In the case where for $j \in \{1, \dots, n\}$, we have $h_j(x) = \frac{x^s}{x^s + \theta_j^s}$ for every $x \in \mathbb{R}$, where $\theta_j > 0$ for every $j \in \{1, \dots, n\}$, $s > 0$, we have $h_{j,k}(x) = \frac{x^s}{x^s + \left(\frac{\theta_j}{\tau_{j,k}}\right)^s}$, so that the latter

case corresponds to the coupling of subsystems possessing different activation thresholds $\frac{\theta_j}{\tau_{j,k}}$, for $(j, k) \in \{1, \dots, n\} \times \{1, \dots, N\}$. Under those assumptions, one can easily see that $H(X) = \tilde{H}(QX)$, where for $X = (X_{j,k})_{(j,k) \in \{1, \dots, n\} \times \{1, \dots, N\}} \in \mathbb{R}^{nN}$, $\tilde{H}(X) = {}^t(H_r(X^1), \dots, {}^t(H_r(X^k))$, with $H_r(X^k) = {}^t(h_1(X_{1,k}), \dots, h_n(X_{n,k}))$ for every $k \in \{1, \dots, N\}$, and $Q$ is the block diagonal matrix composed of a number $N$ of $n \times n$ diagonal blocks $Q^k = \operatorname{diag}(\tau_{j,k})_{j \in \{1, \dots, n\}}$.

**Proposition 4.9.** The sequence of steady states $(\bar{X}_\epsilon)_\epsilon$ defined by Equation (10) satisfies

$$d((\bar{X}_\epsilon)_j, \operatorname{Span}^t(\tau_{j,1}^{-1}, \dots, \tau_{j,N}^{-1})) \to 0,$$

when $\epsilon \to 0$.

*Proof.* Following the same steps as in Section 4.2.2, Equation (10) provides

$$\epsilon \langle \bar{X}_{\epsilon}, Q \bar{X}_{\epsilon} \rangle = \epsilon \langle \Gamma^{-1} B, Q \bar{X}_{\epsilon} \rangle + \epsilon \langle \Gamma^{-1} M H(\bar{X}_{\epsilon}), Q \bar{X}_{\epsilon} \rangle - \langle \Gamma^{-1} \tilde{L}_j H(\bar{X}_{\epsilon}), Q \bar{X}_{\epsilon} \rangle.$$

By definition of $H$ and $\tilde{H}$, we have $\langle \Gamma^{-1} \tilde{L}_j H(\bar{X}_\epsilon), Q \bar{X}_\epsilon \rangle = \langle \Gamma^{-1} \tilde{L}_j \tilde{H}(Q \bar{X}_\epsilon), Q \bar{X}_\epsilon \rangle$. Hence we can immediately deduce the non-negativity of this term by the same arguments as those used in the proof of Proposition 4.4. Then it follows from the positivity of the matrix $Q$ that the sequence $(\bar{X}_\epsilon)_{\epsilon>0}$ is bounded in this case. We conclude the proof by using the same argument as in the proof of Proposition 4.5. $\square$

For every $\alpha_j > 0$, let $X_{\alpha_j}(t)$ be the solution of Equation (9) such that $X_{\alpha_j}(0) = X_0$, where $X_0 \in \mathbb{R}^{nN}$ is independent of $\alpha_j$. By the previous facts, we can state the following result, which is a consequence of the Proposition 4.9 and Theorem 4.1. Roughly speaking, it states that, in the case of similar cells with perturbed parameters, the shifts between the subsystems are approximately preserved in the large coupling regime because the coupling term are preponderant in the dynamics of Equation (9).
**Theorem 4.10.** For an arbitrary $\eta > 0$, there exists $\alpha_0 > 0$ and $t_0 = t_0(\alpha_0) > 0$ such that for every $t \ge t_0$, we have $d((X_{\alpha_0})_j(t), \text{Span}^t(\tau_{j,1}^{-1}, \dots, \tau_{j,N}^{-1})) < \eta$.

**Remark 4.11.** Assuming that there exists a subset $O \subset \{1, \dots, N\}$ such that $\tau_{j,k} = \tau_{j,q}$ for every $k, q \in O$, the latter result allows to prove that Equation (9) achieves j-weak synchronization among the subsystems belonging to $O$. This fact will be illustrated in Section 6.1.2.

# 5 Heterogeneous coupling results

In this section, we propose to study the general heterogeneous case as defined in Definition 3.2, where the coupling functions do not coincide with the internal interaction functions, that is, there exist $(j, k) \in \{1, \dots, n\} \times \{1, \dots, N\}$ such that $c_{j,k} \neq h_{j,k}$. We start by proving a synchronization result in Theorem 5.1 in the case where the individual subsystems are parallel and all possess the same steady state. Then we state a weak synchronization result in Theorem 5.3, under an a priori global stability assumption.

## 5.1 Convergence and synchronization results

We assume that for the uncoupled system, every subsystem possess the same steady state, the subsystems are parallel (i.e. they possess the same functions $(h_j)_{j \in \{1, \dots, n\}}$ belonging to $\mathcal{R}$) and (possibly unbounded) increasing functions $(c_j)_{j \in \{1, \dots, n\}}$, and that the matrices $(A_k)_{k \in \{1, \dots, N\}}$ are simultaneously weakly diagonally stable. Assume moreover that $c_{j,k} = c_j$ for every $(j, k) \in \{1, \dots, n\} \times \{1, \dots, N\}$. The main result of this section proves that, if we are in the case where the subsystems of Equation (7) have a common steady state $\bar{x} \in (\mathbb{R}_+)^n$, then the coupling model preserves the individual behavior, that is every subsystem converge towards $\bar{x}$ when $t \to +\infty$, whatever the choice of the functions $(c_j)_{j \in \{1, \dots, n\}}$ may be. This case allows for instance to consider some examples in which the coupling involves a on/off diffusion as qualitatively explained in Section 3.2.1 with thresholds which are different from the thresholds possibly occuring in the individual dynamics. The proof is based on the use of an explicit Lyapunov function as in Proposition 2.3 which requires the application of Lemma 4.3.

**Theorem 5.1.** Assume that $((A_k)_{k \in \{1, \dots, N\}}, (B^k)_{k \in \{1, \dots, N\}}, (h_j)_{j \in \{1, \dots, n\}})$ are such that Equation (2) with $A = A_k$ has a steady state $\bar{x} \in \mathbb{R}^n$ independent of $k \in \{1, \dots, N\}$, and define $\bar{X} \in \mathbb{R}^{nN}$ by $\bar{X}^k = \bar{x}$ for every $k \in \{1, \dots, N\}$. Then $\bar{X}$ is a globally asymptotically stable steady state of Equation (7). In particular, Equation (7) achieves synchronization in the sense of Definition 3.4.

*Proof.* Consider a positive diagonal matrix $D = \text{diag}(d_j)_{j \in \{1, \dots, n\}}$ such that $DA_k + {}^tA_k D \le 0$, for every $k \in \{1, \dots, N\}$. Define, for $x \in \mathbb{R}$, $V_{j,k}(x) = \int_{\bar{X}_{j,k}}^x (h_j(z) - h_j(\bar{X}_{j,k}))dz$, and for $X \in \mathbb{R}^{nN}$, $V(X) = \sum_{j=1}^n d_j \sum_{k=1}^N V_{j,k}(X_{j,k})$. We are going to show that $V$ is a globally positive definite Lyapunov function for Equation (7) which is radially unbounded, and whose time derivative is globally negative definite. First notice that since $h_j$ is increasing for every $(j, k) \in \{1, \dots, n\} \times \{1, \dots, N\}$, each $V_{j,k}$ is non-negative and $V_{j,k}(x) = 0$ if and only if $x = \bar{X}_j$, hence the same property holds for $V$. By definition of $\bar{X}$, we have

$$ \dot{X}_j = Q_j(X) - Q_j(\bar{X}) - \Gamma_j(X_j - \bar{X}_j) - L_j(C_j(X_j) - C_j(\bar{X}_j)). $$
Hence we get

$$
\begin{align*}
\sum_{j=1}^{n} d_j \sum_{k=1}^{N} \frac{d}{dt} V_{j,k}(X_{j,k}(t)) &= \sum_{j=1}^{n} d_j \sum_{k=1}^{N} \dot{X}_{j,k}(t) (h_j(X_{j,k}(t)) - h_j(\bar{X}_{j,k})) \\
&= \sum_{j=1}^{n} d_j^t (H_j(X_j) - H_j(\bar{X}_j)) \dot{X}_j \\
&= \sum_{j=1}^{n} d_j^t (H_j(X_j) - H_j(\bar{X}_j)) (Q_j(X) - Q_j(\bar{X})) - \sum_{j=1}^{n} d_j^t (H_j(X_j) - H_j(\bar{X}_j)) \Gamma_j (X_j - \bar{X}_j) \\
&\quad - \sum_{j=1}^{n} d_j^t (H_j(X_j) - H_j(\bar{X}_j)) L_j (C_j(X_j) - C_j(\bar{X}_j)) \\
&= \sum_{k=1}^{N} {}^t(H_j(X_j) - H_j(\bar{X}_j)) \Gamma_j (X_j - \bar{X}_j) \\
&\quad - \sum_{j=1}^{n} d_j^t (H_j(X_j) - H_j(\bar{X}_j)) \Gamma_j (X_j - \bar{X}_j) \\
&\quad - \sum_{j=1}^{n} d_j^t (H_j(X_j) - H_j(\bar{X}_j)) L_j (C_j(X_j) - C_j(\bar{X}_j)).
\end{align*}
$$

The simultaneous weak diagonal stability of $(A_k)_{k \in \{1, \dots, N\}}$ implies that the first term is non-positive.
Using the fact that the functions $(h_j)_{j \in \{1, \dots, n\}}$ are increasing, we get that

$$
{}^t(H_j(X_j) - H_j(\bar{X}_j)) \Gamma_j (X_j - \bar{X}_j) \geq 0
$$

for every $j \in \{1, \dots, n\}$. There remains to show that

$$
-{}^t(H_j \circ C_j^{-1}(Y_j) - H_j \circ C_j^{-1}(\bar{Y}_j)) L_j(Y_j - \bar{Y}_j) \le 0,
$$

where $Y_j = C_j(X_j)$ and $\bar{Y}_j = C_j(\bar{X}_j)$. Applying Lemma 4.3, we get that $-{}^t H_j \circ C_j^{-1}(x) L_j x \le 0$ for every $x \in \mathbb{R}^N$. By definition of $\bar{X}$, we have $\bar{X}_j \in \text{Span}(\mathbf{1})$ for every $j \in \{1, \dots, n\}$. As a consequence, we have $\bar{Y}_j \in \text{Span}(\mathbf{1})$ for every $j \in \{1, \dots, n\}$, so that we obtain $L_j \bar{Y}_j = 0$ and ${}^t L_j (H_j \circ C_j^{-1}) (\bar{Y}_j) = 0$ for every $j \in \{1, \dots, n\}$. It follows that

$$
-{}^t(H_j(X_j) - H_j(\bar{X}_j)) L_j(C_j(X_j) - C_j(\bar{X}_j)) \le 0, \quad (K)
$$

for every $j \in \{1, \dots, n\}$. It follows that $\frac{d}{dt}V(X(t)) \le 0$ for every $t \ge 0$. Moreover, we have
$${}^t(H_j(X_j) - H_j(\bar{X}_j)) \Gamma_j (X_j - \bar{X}_j) = 0$ for every $j \in \{1, \dots, n\}$ if and only if $X = \bar{X}$. The fact
that the functions $(h_j)_{j \in \{1, \dots, n\}}$ belong to $\mathcal{R}$ guarantees that $V$ is radially unbounded, so that the
result is proved.
$\square$

<figure><img src="image_3.png" /></figure>

**Remark 5.2.** • A remarkable feature is that the Lyapunov function used in the proof of Theo-
rem 5.1 only depends on the individual dynamics, and is independent of the (possibly unbounded)
coupling functions $(c_j)_{j \in \{1, \dots, n\}}$, which could be utter importance for the design of Lyapunov
functions for systems having uncertain coupling terms.

* Notice that, in the case of coupling of subsystems having different steady states, nothing guarantees that Inequality (K) is satisfied. The problem of finding suitable conditions for global stability in this case remains open.
## 5.2 Weak synchronization result

Assume in this section that the increasing functions $(c_{j,k})_{(j,k)\in\{1,...,n\}\times\{1,...,N\}}$ are bounded, so that they belong to $\mathcal{R}$. In this case, the existence of a steady state for Equation (7) is guaranteed by the same arguments as those used in the proof of Lemma 2.1. As in Section 4.2, we set $j \in \{1,...,n\}$ such that $L_j \neq 0$, and we consider Equation (9) where $L_j$ is replaced by $\alpha_j L_j$, where $\alpha_j > 0$ corresponds to the coupling strength associated with the j-th component, and define $\tilde{L}$ and $\tilde{L}_j$ similarly. By the computations made in the proof of Lemma 2.1, we have that $\bar{X} \in \mathbb{R}^{nN}$ is a steady state of Equation (9) if and only if $\bar{X} = \Gamma^{-1}(B+AH(\bar{X})-LC(X)) = \Gamma^{-1}(B + \tilde{A}H(\bar{X}) - (\tilde{L} + \alpha_j \tilde{L}_j)C(\bar{X}))$. For every $\alpha > 0$, an adaptation of Lemma 2.1 guarantees the existence of a unique steady state $\bar{X}_{\alpha_j}$ solution of

$$ \bar{X}_{\alpha} = \Gamma^{-1} (B + \tilde{A}H(\bar{X}_{\alpha_j}) - (\tilde{L} + \alpha_j \tilde{L}_j)C(\bar{X}_{\alpha_j})). $$

An application of the same arguments as in Section 4.2 proves that in the case where $c_{j,k} = c_j$ for every $(j,k) \in \{1,...,n\} \times \{1,...,N\}$, where $(c_j)_{j \in \{1,...,n\}}$ are functions belonging to $\mathcal{R}$, we have $d((\bar{X}_{\alpha_j})_j, \text{span}(\mathbb{1})) \to 0$ when $\alpha_j \to +\infty$. An immediate consequence is the following result.

**Theorem 5.3.** Assume that, for every $\alpha_j > 0$, Equation (7) admits a unique globally asymptotically stable steady state $\bar{X}_{\alpha_j} \in ((\mathbb{R}^+)^n)^N$, and that $c_{j,k} = c_j$ for every $(j,k) \in \{1,...,n\} \times \{1,...,N\}$, where $(c_j)_{j \in \{1,...,n\}}$ are functions belonging to $\mathcal{R}$. Then Equation (7) achieves j-weak synchronization when $\alpha_j \to +\infty$.

# 6 First application: coupling of two dimensional $(-,+)$ damped oscillators

We focus in this section on the example of the coupling of two dimensional oscillators, whose individual dynamics has already been given in the introduction. To this aim, we define

$$ \mathcal{O} = \{ (A, B) \in M_2(\mathbb{R}) \times \mathbb{R}^2 \mid A = \begin{pmatrix} -w_1 & k_1 \\ -k_2 & -w_2 \end{pmatrix}, B = {}^t(0, k_2), w_1, w_2, k_1, k_2 > 0 \}, $$

and

$$ \tilde{\mathcal{O}} = \{ (A, B) \in \mathcal{O} \mid w_1 w_2 \ge \frac{(k_1 - k_2)^2}{2} \}. $$

One can show that for every $j \in \{1,...,n\}$, we have $(Ax+B)_j \ge 0$ for every $x = (x_1,...,x_n) \ge 0$ such that $x_j = 0$ and $(A,B) \in \mathcal{O}$. Considering two functions $h_1, h_2$ belonging to $\mathcal{R}$ such that $h_1(0) = h_2(0) = 0$, we obtain that $(\mathbb{R}_+)^2$ is invariant by the dynamics of Equation (2) associated with $(A,B)$. Furthermore, for a sequence $(A_k, B^k)_{k \in \{1,...,N\}}$ belonging to $\tilde{\mathcal{O}}$, one can check that $A_k + {}^t A_k \le 0$, for every $k \in \{1,...,N\}$. It follows that the matrices $(A_k)_{k \in \{1,...,N\}}$ are simultaneously weakly diagonally stable.

## 6.1 Synchronization in the homogeneous coupling case

Here we consider a network of $N = 4$ coupled two dimensional oscillators given by a sequence $(A_k, B^k)_{k \in \{1,...,N\}}$ belonging to $\tilde{\mathcal{O}}$ and satisfying Equation (4), whose two components are coupled

through the symmetric matrices $L_1 = L_2 = \alpha \begin{pmatrix} -2 & 1 & 0 & 1 \\ 1 & -2 & 1 & 0 \\ 0 & 1 & -2 & 1 \\ 1 & 0 & 1 & -2 \end{pmatrix}$. Note however that our results also work with non-symmetric coupling matrices, under the assumptions stated in Section 3.2.2.
### 6.1.1 Weak synchronization in the parallel case $h_{1k} = h_1$ and $h_{2k} = h_2$ for every $k \in \{1, \dots, N\}$.

Consider non-negative functions $h_1, h_2$ belonging to $\mathcal{R}$ such that $h_1(0) = h_2(0) = 0$, and a family of functions $(h_{1k})_{k \in \{1, \dots, N\}}$ which are identically equal to $h_1$, and a family of non-negative functions $(h_{2k})_{k \in \{1, \dots, N\}}$ which are identically equal to $h_2$. By Theorem 4.1, for every $\alpha \ge 0$, each subsystem converges to a point of $(\mathbb{R}_+)^2$, and the two components of the subsystems achieve weak synchronization asymptotically when $\alpha \to +\infty$ in the sense of Theorem 4.7. We illustrate this phenomenon on Figure 2 with $N = 4$ and $h_j(x) = \frac{x^s}{x^s+\theta_j^s}$ for every $x \in \mathbb{R}$, $s \ge 2$, $j \in \{1, 2\}$, where $\theta_j > 0$, and arbitrary other parameters of the system. We plotted the subsystems that can be identified from the solution $X_\alpha(t)$ of Equation (9) such that $X_\alpha(0) = X_0 \in \mathbb{R}^8$ for $t \ge 0$. On Figure 2(a), we have plotted the trajectories of each subsystem without coupling, i.e. $\alpha = 0$ from $t=0$ to $t=T=50$. Subsystem 1 is plotted in green, subsystem 2 in red, subsystem 3 in blue, subsystem 4 in magenta, and we notice that the trajectories converge towards different points. On Figure 5(b), we have plotted the trajectories of each subsystem with $\alpha = 50$ from $t=0$ to $t=T=50$. The trajectories converge toward very close points, in accordance with the result of Theorem 4.7.

<figure><img src="image_4.png" /><figcaption>Figure 2: Trajectories of 4 coupled different two dimensional oscillators with Hill functions of order $s = 15$, $T = 50$, thresholds $\theta_1 = 1$, $\theta_2 = 1.2$.</figcaption></figure>

### 6.1.2 Partial synchronization for non-parallel shifted systems

In the case of subsystems with different interaction functions, we propose to illustrate the results of Section 4.2.3 when $N = 4$. Consider a non-negative functions $h_1, h_2$ belonging to $\mathcal{R}$ such that $h_1(0) = h_2(0) = 0$, and two families of functions $(h_{1k})_{k \in \{1, \dots, N\}}$, $(h_{2k})_{k \in \{1, \dots, N\}}$ such that $h_{11}(x) = h_1(x)$, $h_{12} = h_1(\tau_{12}x)$, $h_{13}(x) = h_1(\tau_{12}x)$, $h_{14}(x) = h_1(x)$, $h_{21}(x) = h_2(x)$, $h_{22} = h_2(\tau_{22}x)$, $h_{23}(x) = h_2(\tau_{22}x)$, $h_{24}(x) = h_2(x)$, where $\tau_{12} \ne 1$ and $\tau_{22} \ne 1$ are positive, as in Theorem 4.10. We consider the subsystems that can be identified from the solution $X_\alpha(t)$ of Equation (9) such that $X_\alpha(0) = X_0 \in \mathbb{R}^8$ for $t \ge 0$. By Theorem 4.1, for every $\alpha \ge 0$, each subsystem converges to a point of $(\mathbb{R}_+)^2$, and the two components of the subsystems 1 and 4 (respectively, 2 and 3) are close for large enough $\alpha \ge 0$ and $t$,
as claimed in Theorem 4.10. We illustrate these phenomena on Figure 3 with $N = 4$, $h_j(x) = \frac{x^s}{x^s+\theta_j^s}$ for every $x \in \mathbb{R}$, $s \ge 2$, $j \in \{1,2\}$, where $\theta_j > 0$, $\tau_{22} = 0.8$ and $\tau_{12} = 0.5$, and arbitrary other parameters of the system. Define $\tilde{\theta}_j = \theta_j \tau_{j2}$ for $j \in \{1,2\}$. We plotted the subsystems identified from the solution $X_\alpha(t)$ of Equation (9) such that $X_\alpha(0) = X_0 \in \mathbb{R}^8$ for $t \ge 0$. Subsystem 1 is plotted in green, subsystem 2 in red, subsystem 3 in blue, subsystem 4 in magenta, and we notice that subsystems 1 and 4 (respectively, 2 and 3) seem to converge towards very close points, in accordance with Theorem 4.10.

<figure><img src="image_5.png" /><figcaption>Figure 3: Trajectories of coupled different oscillators with Hill functions of order $s = 15$, $T = 50$, thresholds $\theta_1 = 1$, $\theta_2 = 1.2$, $\tau_{22} = 0.8$, $\tau_{12} = 0.5$.</figcaption></figure>

## 6.2 Synchronization in the heterogeneous coupling case

Here we consider a network of $N = 4$ coupled two dimensional oscillators given by a sequence $(A_k, B^k)_{k \in \{1,...,N\}}$ belonging to $\tilde{\mathcal{O}}$ and satisfying Equation (4), whose two components are coupled

through the symmetric matrices $L_1 = L_2 = \alpha \begin{pmatrix} -2 & 1 & 0 & 1 \\ 1 & -2 & 1 & 0 \\ 0 & 1 & -2 & 1 \\ 1 & 0 & 1 & -2 \end{pmatrix}$. We illustrate the results of

Section 5.1 on Figure 5 by taking $h_j(x) = \frac{x^s}{x^s+\theta_j^s}$, $c_j(x) = \frac{x^s}{x^s+\tilde{\theta}_j^s}$ for every $x \in \mathbb{R}$, $s \ge 2$, with $\theta_j > 0$, $\tilde{\theta}_j > 0$ for $j \in \{1,2\}$. The coupling functions $(c_j)_{j \in \{1,2\}}$ between the subsystems is assumed to depend on the thresholds $(\tilde{\theta}_j)_{j \in \{1,2\}}$, which are different from the thresholds $(\theta_j)_{j \in \{1,2\}}$ of the internal dynamics. We plotted the subsystems that can be identified from the solution $X_\alpha(t)$ of Equation (7) such that $X_\alpha(0) = X_0 \in \mathbb{R}^8$ for $t \ge 0$. Subsystem 1 is plotted in green, subsystem 2 in red, subsystem 3 in blue, subsystem 4 in magenta, and we notice that every subsystem seem to converge towards a common steady state which is very close to the point $(\theta_1, \theta_2)$, in accordance with Theorem 5.1. When there is no coupling ($\alpha = 0$), the thresholds $\tilde{\theta}_j > 0$ for $j \in \{1,2\}$ have no influence on the dynamics, while their influence can be seen when $\alpha > 0$.
<figure><img src="image_6.png" /><figcaption>Figure 4: Trajectories of coupled identical oscillators choosing Hill functions both as internal interaction and coupling functions, whose order is $s = 40$, and the thresholds are $\theta_1 = 1$, $\theta_2 = 1.2$, $\tilde{\theta}_1 = 0.7$, $\tilde{\theta}_2 = 1.4$, and the final time $T = 50$.</figcaption></figure>

# 7 Second application: stabilization of multistable systems

A challenging task for genetic regulatory networks is the stabilization of their steady states, and several control strategies have been proposed for instance in [10] in the PWA setting or in [8] in the smooth setting. For diffusively coupled bistable systems, we studied in [3] some control strategies for synchronization involving the two asymptotically stable steady states of the individual systems in the PWA case.

## 7.1 General feedback control result

Consider a sequence $(A_k, B^k)_{k \in \{1, \dots, N\}}$, and let $\bar{X}^k \in (\mathbb{R}_+)^n$ be a steady state of Equation (7) associated with $(A_k, B^k)$, set $\bar{X} = (\bar{X}^1, \dots, \bar{X}^N)$, and consider the controlled equation

$$ \dot{X} = B + \tilde{A}H(X) - \Gamma X - LC(X) + U(X), \quad (11) $$

where $X(t) \in \mathbb{R}^{nN}$, where $\tilde{A}$ and $L$ are defined as in Section 3.2.2. For every $X \in \mathbb{R}^{nN}$, the control $U(X) \in \mathbb{R}^{nN}$ is a bounded feedback control which is under the form $U(X) = P(H(X) - H(\bar{X}))$, where

$$ P = \begin{pmatrix} P_1 & 0 & \dots & 0 \\ 0 & P_2 & \dots & \vdots \\ \vdots & \ddots & \ddots & 0 \\ 0 & \dots & 0 & P_N \end{pmatrix} $$

is a concatenation of $n \times n$ real matrices $(P_k)_{k \in \{1, \dots, N\}}$.

A feedback control $U$ is said to be *admissible* when the matrices $(P_k)_{k \in \{1, \dots, N\}}$ are such that the cone $((\mathbb{R}_+)^n)^N$ is forward invariant by the dynamics of Equation (11) (i.e. the non-negativity of the variables is preserved for biological sense).
Assume that for every $k \in \{1, \dots, N\}$, the matrices $(A_k+P_k)$ are simultaneously weakly diagonally stable, and that the feedback control $U(X) = P(H(X)-H(\bar{X}))$ is admissible. Whatever the choice of the matrices $(A_k)_{k \in \{1, \dots, N\}}$ may be, the latter conditions are satisfied, for instance when for every $k \in \{1, \dots, N\}$, $P_k = -\beta \mathrm{Id}$, where $\mathrm{Id}$ is the identity $n \times n$ matrix, and $\beta > 0$ is large enough. Noticing that Equation (11) admits $\bar{X}$ as a steady state, we can apply the results of the sections 4 and 5.1 in order to deduce the following global convergence results in different frameworks:

* **Homogeneous coupling case:** if $c_{jk} = h_{jk}$ for every $j \in \{1, \dots, n\}$ and $k \in \{1, \dots, N\}$ (without supplementary assumptions on the steady state $\bar{X}$), then the solutions of Equation (11) globally converge towards $\bar{X}$.

* **Heterogeneous coupling case:** if $\bar{X}^k = \bar{X}^q$, $h_{jq} = h_{jk}$, and $c_{jq} = c_{jk}$ for every $j \in \{1, \dots, n\}$, $k, q \in \{1, \dots, N\}$, then the solutions of Equation (11) globally converge towards $\bar{X}$.

**Remark 7.1.** Such a control strategy requires the precise knowledge of the individual steady states, which is a hard task in practice. However, note that in the PWA limit, such a control strategy only depends on the regular domains to which the state $X(t)$ belongs and not of its exact value, as it would be the case for instance if $U(X)$ was a linear feedback under the form $U(X) = P(X - \bar{X})$.

## 7.2 Stabilization of coupled bistable switch systems

We focus in this section on the example of the coupling of two dimensional identical bistable switch systems, whose individual dynamics are described by

$$ 
\begin{aligned}
\dot{x}_1 &= -\gamma_1 x_1 + k_1 h_1(x_2) \\
\dot{x}_2 &= -\gamma_2 x_2 + k_2 h_2(x_1),
\end{aligned}
\qquad (B) $$

where $\gamma_j, k_j > 0$, and $h_j(x) = \frac{x^{s_j}}{x^{s_j}+\theta_j^{s_j}}$ with $s_j \ge 2$ for every $j \in \{1, 2\}$. It is known (see for instance [8]) that for $s_j \ge 2$ large enough, such system has three steady states:

* Two locally asymptotically stable steady states $\bar{x}, \tilde{x}$ belonging to $(\mathbb{R}_+)^2$.

* A saddle point $x_{\mathrm{undiff}}$ belonging to $(\mathbb{R}_+)^2$, called *undifferentiated state*.

In this section we propose to use the methodology exposed in Section 7.1 in order to synchronize coupled bistable subsystems at an arbitrary steady state. In particular, the control strategy will work to obtain synchronization of subsystems whose individual dynamics are driven by Equation (B) towards their undifferentiated steady state $x_{\mathrm{undiff}}$. To the best of the authors' knowledge, the latter problem has not been explored yet, and it is an challenging task due to the role of $x_{\mathrm{undiff}}$ in cell differentiation (see, for instance [4, 19]).

To this aim, we define

$$ S = \{ (A, B) \in M_2(\mathbb{R}) \times \mathbb{R}^2 \mid A = \begin{pmatrix} -w_1 & k_1 \\ k_2 & -w_2 \end{pmatrix}, B \ge 0, w_1, w_2, k_1, k_2 > 0 \}, $$

and

$$ \tilde{S} = \{ (A, B) \in S \mid w_1 w_2 \ge \frac{(k_1 + k_2)^2}{2} \}. $$

One can show that for every $j \in \{1, \dots, n\}$, we have $(Ax+B)_j \ge 0$ for every $x = (x_1, \dots, x_n) \ge 0$ such that $x_j = 0$ and $(A, B) \in S$. Considering two functions $h_1, h_2$ belonging to $\mathcal{R}$ such that $h_1(0) = h_2(0) = 0$, we obtain that $(\mathbb{R}_+)^2$ is invariant by the dynamics of Equation (2) associated
with $(A, B)$. Furthermore, for a sequence $(A_k, B^k)_{k \in \{1, \dots, N\}}$ belonging to $\tilde{S}$, one can check that $A_k + {}^t A_k \le 0$, for every $k \in \{1, \dots, N\}$. It follows that the matrices $(A_k)_{k \in \{1, \dots, N\}}$ are simultaneously weakly diagonally stable. Now consider a sequence $(A_k, B^k)_{k \in \{1, \dots, N\}}$, where $B^k \ge 0$ and $A_k$ has the form $\begin{pmatrix} -w_1 & k_1 \\ k_2 & -w_2 \end{pmatrix}$, with $w_1, w_2 \ge 0$, $k_1, k_2 > 0$ for every $k \in \{1, \dots, N\}$. Nevertheless, assume that there exists at least an integer $k \in \{1, \dots, N\}$ such that $(A_k, B^k) \notin \tilde{S}$. In particular, the case where the dynamics is described by Equation (B) satisfies the latter hypothesis. Define for every a $2 \times 2$ non-positive diagonal matrix $P_s = -\beta \text{Id}$ with $\beta \ge 0$ such that $(A_k + P_s, B^k) \in \tilde{S}$ for every $k \in \{1, \dots, N\}$. Select for every $k \in \{1, \dots, N\}$, a steady state $\bar{X}^k$ of Equation (3) associated with $(A_k, B^k)$, and set $\bar{X} = (\bar{X}^1, \dots, \bar{X}^N)$. Consider the corresponding controlled Equation (11) with $U(X) = P(H(X) - H(\bar{X}))$, where

$$P = \begin{pmatrix} P_s & 0 & \dots & 0 \\ 0 & P_s & \dots & \vdots \\ \vdots & \vdots & \ddots & 0 \\ 0 & \dots & 0 & P_s \end{pmatrix}$$

is a concatenation of the $n \times n$ real matrix $P_s$. Notice that the latter equation has a biological sense provided that the cone $((\mathbb{R}_+)^2)^N$ is forward invariant under its flow. In accordance with the definition of $\tilde{A}$ and $L$ given in Section 3.2.2, we consider a matrix $\tilde{A}$ depending on the sequence $(A_k)_{k \in \{1, \dots, N\}}$

and a coupling matrix $L$ depending on the matrices $L_1 = L_2 = \alpha \begin{pmatrix} -2 & 1 & 0 & 1 \\ 1 & -2 & 1 & 0 \\ 0 & 1 & -2 & 1 \\ 1 & 0 & 1 & -2 \end{pmatrix}$ with $\alpha > 0$.

On Figure 5, we propose to illustrate the corresponding global convergence results obtained in Section 7.1, by simulating Equation (11) both in the homogeneous and heterogeneous coupling cases, focusing on identical bistable switch systems that we stabilize at their common undifferentiated point $x_{\text{undiff}}$. Note that, in the homogeneous coupling case, the control strategy also works for the coupling of systems having different steady states.

# 8 Conclusion

In this work, we proposed some new structural properties that guaranteed synchronization of genetic regulatory networks, which is a crucial task for applications to real life complex systems’ modeling and control. More specifically we introduced the notion of weak synchronization in the large coupling regime, and we showed that this property holds for a large class of coupling mechanisms, then we illustrated this phenomenon for coupled damped oscillators. We also proposed a feedback control strategy which stabilizes a network of coupled toggle switch systems at their undifferentiated point. The aim of the further research in this direction is to get a suitable convergence criterion for heterogeneous networks of systems possessing different steady states, and to extend our methods to the synchronization of systems having sustained oscillatory behaviors.

# A Proof of Lemma 4.3

*Proof.* Let $K$ be defined as in the claim of the lemma. By simple considerations, we show the existence of $\alpha > 0$ such that $A = \frac{1}{\alpha}K + \text{Id}$ is a non-negative doubly stochastic matrix. We propose to prove that, for every non-negative doubly stochastic matrix $A$, we have

$$ {}^t \mathcal{H}(y) A y - {}^t \mathcal{H}(y) y \le 0, \qquad (12) $$
<figure><img src="image_7.png" /><figcaption>Figure 5: Stabilization of four identical coupled bistable switch systems towards the undifferentiated point</figcaption></figure>

for every $y \in \mathbb{R}^N$. Indeed, once this is proved, the result can be deduced using the inequality $\frac{1}{\alpha} {}^t\mathcal{H}(y)Ky \le 0$, for every $y \in \mathbb{R}^N$.

Let $A$ be a non-negative doubly stochastic matrix. According to the Birkhoff Von-Neumann Theorem, $A$ can be decomposed as a convex combination of permutation matrices according, that is, there exist $N \times N$ permutation matrices $(P_j)_{j \in \{1, \dots, k\}}$ and real numbers $(\theta_j)_{j \in \{1, \dots, k\}}$ belonging to $[0, 1]$ satisfying $\sum_{j=1}^k \theta_j = 1$, such that $A = \sum_{j=1}^k \theta_j P_j$. Hence we can reduce the problem of showing Inequality (12) to showing that, for every permutation matrix $P$, we have

$$ {}^t\mathcal{H}(y)Py - {}^t\mathcal{H}(y)y \le 0, \qquad (13) $$

for every $y \in \mathbb{R}^N$. Using the fact that permutations can be decomposed as products of disjoint support cycles, it is sufficient to show the result when $P$ is a cycle.

We propose to show the Inequality (13) by induction on the length $l$ of the cycle $\sigma$ (whose matrix is $P$ in the canonical basis of $\mathbb{R}^N$):

* For $l = 2$, consider a transposition $\sigma = (kj)$, with $k, j \in \{1, \dots, N\}$. Then we have ${}^t\mathcal{H}(y)Py - {}^t\mathcal{H}(y)y = (y_k - y_j)(h(y_j) - h(y_k))$, for every $y \in \mathbb{R}^N$, and the claim is true because $h$ is strictly increasing.

* Set $l \in \{1, \dots, N\}$, and assume that for every cycle $\sigma$ (whose matrix is $P$ in the canonical basis of $\mathbb{R}^N$) of length $l$, we have ${}^t\mathcal{H}(y)Py - {}^t\mathcal{H}(y)y \le 0$, for every $y \in \mathbb{R}^N$. Consider a cycle $\tilde{\sigma}$ of length $l+1$, associated with a $N \times N$ matrix $\tilde{P}$. Up to a reordering of the components of $y \in \mathbb{R}^N$, we can assume that $\tilde{\sigma} = (1 \cdots l+1)$. By such a choice, we have

$$ {}^t\mathcal{H}(y)\tilde{P}y - {}^t\mathcal{H}(y)y = \sum_{j \in \text{supp}(\tilde{\sigma})} h(y_j)(y_{\tilde{\sigma}(j)} - y_j). $$
For a fixed $y \in \mathbb{R}^N$, let $\bar{k} \in \{1, \dots, N\}$ be such that $y_{\bar{k}} \ge y_k$ for every $k \in \{1, \dots, N\}$, define the cycle $\sigma = (1 \cdots \bar{k}-1 \bar{k}+1 \cdots l+1)$ of length $l$, and denote the $N \times N$ matrix associated with $\sigma$ by $P$. Then we can write

$$ \begin{aligned} {}^t\mathcal{H}(y)\tilde{P}y - {}^t\mathcal{H}(y)y &= \sum_{j \in \text{supp}(\sigma)} h(y_j)(y_{\sigma(j)} - y_j) + (y_{\bar{k}+1} - y_{\bar{k}})(h(y_{\bar{k}}) - h(y_{\bar{k}-1})) \\ &= {}^t\mathcal{H}(y)Py - {}^t\mathcal{H}(y)y + (y_{\bar{k}+1} - y_{\bar{k}})(h(y_{\bar{k}}) - h(y_{\bar{k}-1})). \end{aligned} $$

By the induction hypothesis, we have ${}^t\mathcal{H}(y)Py - {}^t\mathcal{H}(y)y \le 0$. The definition of $\bar{k}$ and the fact that $h$ is strictly increasing prove that $(y_{\bar{k}+1} - y_{\bar{k}})(h(y_{\bar{k}}) - h(y_{\bar{k}-1})) \le 0$, and we obtain that ${}^t\mathcal{H}(y)\tilde{P}y - {}^t\mathcal{H}(y)y \le 0$.

We get the result by induction.

<figure><img src="image_8.png" /></figure>

## References

[1] M. Arcak and E. D. Sontag. A passivity-based stability criterion for a class of biochemical reaction networks. *Mathematical biosciences & engineering*, 5(1):1, 2008.

[2] E. Armingol, A. Officer, O. Harismendy, and N. E. Lewis. Deciphering cell-cell interactions and communication from gene expression. *Nature Reviews Genetics*, 22:71–88, 2021.

[3] N. Augier, M. Chaves, and J.-L. Gouzé. Qualitative control strategies for synchronization of bistable gene regulatory networks. https://hal.inria.fr/hal-02953502, 2020.

[4] G. Balázsi, A. van Oudenaarden, and J. J. Collins. Cellular decision making and biological noise: from microbes to mammals. *Cell*, 144(6):910–925, 2011.

[5] A. Berman and R. J. Plemmons. *Nonnegative matrices in the mathematical sciences*, volume 9 of *Classics in Applied Mathematics*. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 1994. Revised reprint of the 1979 original.

[6] O. Brandman, J. E. Ferrell Jr., R. Li, and T. Meyer. Interlinked fast and slow positive feedback loops drive reliable cell decisions. *Science*, 310:496–498, 2005.

[7] L. Chambon, I. Belgacem, and J.-L. Gouzé. Qualitative control of undesired oscillations in a genetic negative feedback loop with uncertain measurements. *Automatica*, 112:108642, 2020.

[8] L. Chambon and J.-L. Gouzé. A new qualitative control strategy for the genetic toggle switch. *IFAC-PapersOnLine*, 52(1):532–537, 2019. 12th IFAC Symposium on Dynamics and Control of Process Systems, including Biosystems DYCOPS 2019.

[9] M. Chaves and H. de Jong. *Synthetic Gene Circuits. Methods in Molecular Biology*, volume 2229, chapter Qualitative Modeling, Analysis and Control of Synthetic Regulatory Circuits, pages 1–40. Humana, New York, NY, 2021.

[10] M. Chaves and J.-L. Gouzé. Exact control of genetic networks in a qualitative framework: the bistable switch example. *Automatica J. IFAC*, 47(6):1105–1112, 2011.

[11] M. Chaves, L. Scardovi, and E. Firippi. Coupling and synchronization of piecewise linear genetic regulatory systems. In *CDC 2019 - 58th IEEE Conference on Decision and Control*, Nice, France, Dec. 2019.
[12] R. Edwards and P. Gill. On synchronization and cross-talk in parallel networks. volume 10, pages 287-300. 2003. Second International Conference on Dynamics of Continuous, Discrete and Impulsive Systems (London, ON, 2001).

[13] E. Farcot and J.-L. Gouzé. A mathematical framework for the control of piecewise-affine models of gene networks. *Automatica J. IFAC*, 44(9):2326-2332, 2008.

[14] L. Glass and J. S. Pasternack. Prediction of limit cycles in mathematical models of biological oscillations. *Bull. Math. Biology*, 40(1):27-44, 1978. Papers presented at the Society for Mathematical Biology Meeting (Univ. Pennsylvania, Philadelphia, Pa., 1976).

[15] M. Golubitsky and I. Stewart. Rigid patterns of synchrony for equilibria and periodic cycles in network dynamics. *Chaos: An Interdisciplinary Journal of Nonlinear Science*, 26:094803, 09 2016.

[16] F. S. P. Gos, V. Petrenko, C. Hagedorn, F. Kreppel, K.-F. Storch, D. Knutti, A. Liani, C. Weitz, Y. Emmenegger, P. Franken, L. Bonacina, C. Dibner, and U. Schibler. Circadian hepatocyte clocks keep synchrony in the absence of a master pacemaker in the suprachiasmatic nucleus or other extrahepatic clocks. *Genes & Development*, 35:329-334, 2021.

[17] J.-L. Gouzé. An inequality involving nonnegative matrices and increasing functions. *Inria Research Reports*, RR-1330, 1990.

[18] M. W. Hirsch. *Differential topology*. Springer-Verlag, New York-Heidelberg, 1976. Graduate Texts in Mathematics, No. 33.

[19] J.-B. Lugagne, S. S. Carrillo, M. Kirch, A. Köhler, G. Batt, and P. Hersen. Balancing a genetic toggle switch by real-time feedback control and periodic forcing. *Nature communications*, 8(1):1-8, 2017.

[20] J. Mallet-Paret and H. L. Smith. The poincaré-bendixson theorem for monotone cyclic feedback systems. *Journal of Dynamics and Differential Equations*, 2(4):367-421, 1990.

[21] C. B. Medina, P. Mehrotra, S. Arandjelovic, J. S. Perry, Y. Guo, S. Morioka, B. Barron, S. F. Walk, B. Ghesquière, A. S. Krupnick, U. Lorenz, and K. S. Ravichandran. Metabolites released from apoptotic cells act as tissue messengers. *Nature*, 580:130-135, 2020.

[22] M. Pasquini and D. Angeli. On convergence for piecewise affine models of gene regulatory networks via a lyapunov approach. *IEEE Transactions on Automatic Control*, 65(8):3333-3348, 2020.

[23] T. Pereira, J. Eldering, M. Rasmussen, and A. Veneziani. Towards a theory for diffusive coupling functions allowing persistent synchronization. *Nonlinearity*, 27(3):501-525, feb 2014.

[24] E. Plahte, T. Mestl, and S. W. Omholt. Global analysis of steady points for systems of differential equations with sigmoid interactions. *Dynamics and Stability of Systems*, 9(4):275-291, 1994.

[25] A. Pogromsky, T. Glad, and H. Nijmeijer. On diffusion driven oscillations in coupled dynamical systems. *International Journal of Bifurcation and Chaos*, 09(04):629-644, 1999.

[26] L. Scardovi, M. Arcak, and E. D. Sontag. Synchronization of interconnected systems with applications to biochemical networks: an input-output approach. *IEEE Trans. Automat. Control*, 55(6):1367-1379, 2010.
[27] L. Scardovi and R. Sepulchre. Synchronization in networks of identical linear systems. *Automatica*, 45(11):2557–2562, 2009.

[28] T. Y.-C. Tsai, Y. S. Choi, W. Ma, J. R. Pomerening, C. Tang, and J. E. Ferrell. Robust, tunable biological oscillations from interlinked positive and negative feedback loops. *Science*, 321(5885):126–129, 2008.