M A Aziz-Alaoui

Guillaume Cantin
email: guillaume.cantin@univ-nantes.fr

Alexandre Thorel

Synchronization of Turing patterns in complex networks of reaction-diffusion systems set in distinct domains

Keywords: Reaction-diffusion, complex network, synchronization, Turing pattern, Robin boundary condition

We present an innovative complex network of reaction-diffusion systems set in distinct domains, with boundary couplings. The complex network models the evolution of interacting populations living in a heterogeneous and fragmented habitat, whose biological individuals migrate from one patch to another. In our model, the displacements of individuals are described by mixed boundary couplings, involving both the Neumann and Robin boundary conditions, which improve the modeling of migrations by point-wise couplings. We investigate the cases of diffusion in isotropic and anisotropic habitats and establish original sufficient conditions of synchronization in this complex network model, for complete graphs, cyclic graphs and rings of nearest neighbors. In parallel, we apply our theoretical framework to a nonlinear predator-prey model with Leslie-Gower-type functional response and explore numerically the emergence of synchronization on heterogeneous Turing patterns.

Introduction

In this paper, we propose to study the dynamics of complex networks of reaction-diffusion systems of the form

                     ∂U i ∂t = D i ∆ Ωi U i + f (U i ), (t, x i ) ∈ (0, ∞) × Ω i , ( 1
)
∂U i ∂ν i = 0, (t, x i ) ∈ (0, ∞) × Γ N i , ( 2
)
∂U i ∂ν i = -µ i (x i ) j∈Ni (U i -U j ) , (t, x i ) ∈ (0, ∞) × Γ R i , ( 3
)
U i (0, x i ) = U i,0 (x i ), x i ∈ Ω i , (4)
with i ∈ {1, . . . , n}. Here, n is a positive integer; the domains (Ω i ) 1≤i≤n are open, bounded and connected sets included in R M with M ≥ 1. For each i ∈ {1, . . . , n}, we assume that the boundary ∂Ω i of Ω i is regular and can be split into two disjoint boundaries:

∂Ω i = Γ N i ∪ Γ R i , Γ N i ∩ Γ R i = ∅.
The outward unit normal vector at point x i of ∂Ω i is denoted by ν i (x i ) or simply ν i . The boundaries Γ N i and Γ R i are associated with the boundary conditions [START_REF] Aguirre | Three limit cycles in a Leslie-Gower predator-prey model with additive Allee effect[END_REF] of Neumann type and (3) of Robin type, respectively. The unknown functions (U i ) 1≤i≤n are defined in (0, +∞) × Ω 1 , . . . , (0, +∞) × Ω n , respectively, with values in R m , where m is a positive integer; the functions (U i ) 1≤i≤n model population densities of interacting biological individuals; each domain Ω i , 1 ≤ i ≤ n, models the habitat of the interacting species U i = (U i,1 , . . . , U i,m ) . These interacting species are subject to spatial mobility within the domains (Ω i ) 1≤i≤n and between those domains, as will be soon explained. The matrices (D i ) 1≤i≤n are diagonal and we assume that their diagonal coefficients D i,k , 1 ≤ k ≤ m, are positive. The operators (∆ Ωi ) 1≤i≤n are multiple instances of the Laplace operator, defined as

∆ Ωi = ∂ 2 ∂x 2 i,1 + • • • + ∂ 2 ∂x 2 i,M , 1 ≤ i ≤ n, ( 5
)
where the space variable in Ω i is denoted by x i = (x i,1 , . . . , x i,M ) . In equation ( 1), the diffusion term D i ∆ Ωi U i models the spatial diffusion of U i in the interior of the domain Ω i (1 ≤ i ≤ n). In addition, the domains (Ω i ) 1≤i≤n are coupled in such a way that each domain Ω i , 1 ≤ i ≤ n, admits a finite number of neighbors; we denote by N i the set of indices j ∈ {1, . . . , n} such that Ω i is coupled with Ω j . The connections between the domains (Ω i ) 1≤i≤n are symmetric, that is, j ∈ N i if and only if i ∈ N j (1 ≤ i, j ≤ n); those connections determine the edges of a graph G , whose vertices are the domains (Ω i ) 1≤i≤n (see Figure 1 below). If a domain Ω i is coupled with another domain Ω j (1 ≤ i, j ≤ n), then we introduce a boundary coupling function µ i , defined on the boundary Γ R i . In equation ( 3), each boundary term -µ i (x i )(U i -U j ) of the sum over N i models the spatial mobilities of species U i and U j between Ω i and Ω j . We emphasize that if a domain Ω i is connected to two distinct domains Ω j and Ω j (j, j ∈ N i , j = j ), then the mobilities from Ω i to Ω j and from Ω i to Ω j start from the same boundary Γ R i , although they could start from distinct boundaries Γ R i,j , Γ R i,j . We suppose that there exists a positive constant µ 0 such that µ i (x i ) ≥ µ 0 for all x i ∈ Γ R i and all i ∈ {1, . . . , n}. We also suppose that there exists a homeomorphism φ i,j that maps Γ R i onto Γ R j , so that the boundary coupling is defined by the Robin boundary condition [START_REF] Ambrosio | Global attractor of complex networks of reaction-diffusion systems of Fitzhugh-Nagumo type[END_REF], with the convention

U j (x i ) = U j φ i,j (x i ) , x i ∈ Γ R i , 1 ≤ i, j ≤ n. ( 6
)
We emphasize that the latter equation models the conservation of population during the mobility between two domains Ω i and Ω j . In other words, we assume that biological individuals moving from a given patch Ω i necessarily reach some other patch Ω j , hence do not leave the whole habitat. Furthermore, we assume that the mobilities between two domains are instantaneous, although a time delay could be introduced at this stage. However, for simplicity, we do not focus on this point in the present paper. Next, the function f involved in the reaction-diffusion system (1) is a nonlinear operator with values in R m , which models the interactions between the biological species U i,1 , . . . , U i,m living on domain Ω i (1 ≤ i ≤ n); the regularity of the function f will be detailed below. Finally, the functions (U i,0 ) 1≤i≤n are initial conditions defined in Ω 1 , . . . , Ω n , respectively.

Related works.

Although complex networks of dynamical systems determined by ordinary differential equations have been studied for more than two decades (see for instance [START_REF] Ha | Remarks on the complete synchronization of Kuramoto oscillators[END_REF], [START_REF] Lewis | Steady states, limit cycles, and chaos in models of complex biological networks[END_REF], [START_REF] Manrubia | Mutual synchronization and clustering in randomly coupled chaotic dynamical networks[END_REF], [START_REF] Maria Da Conceiçao | Homogeneous three-cell networks[END_REF] or [START_REF] Zhang | Identical synchronization of nonidentical oscillators: When only birds of different feathers flock together[END_REF] and the references therein), the study of complex networks of dynamical systems determined by partial differential equations is very recent. It has encountered a rapidly growing interest, due to the rich dynamics of their trajectories and to the great number of real-world applications. Indeed, various forms of synchronization, such as identical synchronization, have been studied in [START_REF] Ambrosio | Global attractor of complex networks of reaction-diffusion systems of Fitzhugh-Nagumo type[END_REF], [START_REF] Rätz | Symmetry breaking in a bulk-surface reaction-diffusion model for signalling networks[END_REF], [START_REF] Wu | Synchronization of coupled reaction-diffusion stochastic neural networks with time-varying delay via delay-dependent impulsive pinning control algorithm[END_REF] or [START_REF] Yang | Synchronization of coupled reaction-diffusion neural networks with timevarying delays via pinning-impulsive controller[END_REF] for complex networks determined by reaction-diffusion systems. The stability of persistence or extinction equilibria in meta-population models have been studied in [START_REF] Cantin | Large-time dynamics in complex networks of reactiondiffusion systems applied to a panic model[END_REF] for a panic model, in [START_REF] Cantin | Dimension estimate of attractors for complex networks of reaction-diffusion systems applied to an ecological model[END_REF] for a competing species system or in [START_REF] Verma | Network synchronization, stability and rhythmic processes in a diffusive mean-field coupled SEIR model[END_REF] for an epidemiological model. In [START_REF] Cantin | On a generalized diffusion problem: A complex network approach[END_REF], it has been proved that the spatial diffusion of individuals in such meta-population models acts as a combination of short and long range diffusion. In parallel, the dynamics of chemical reactions networks have been studied in [START_REF] Desvillettes | Trend to equilibrium for reaction-diffusion systems arising from complex balanced chemical reaction networks[END_REF], [START_REF] Fellner | The entropy method for reaction-diffusion systems without detailed balance: First order chemical reaction networks[END_REF] via the entropy framework; synchronization of unstable patterns in other chemical reactions networks has been investigated in [START_REF] Miranville | Bifurcations and synchronization in networks of unstable reaction-diffusion systems[END_REF]. It is observed that the main research axis which motivates these works is very often related to the synchronization phenomenon. Although innovative, these works however suffer several limitations. Notably, the modeling of the couplings is often reductive. For instance, displacements of individuals are roughly modeled by point-wise couplings in [START_REF] Cantin | Dimension estimate of attractors for complex networks of reaction-diffusion systems applied to an ecological model[END_REF] and [START_REF] Cantin | Large-time dynamics in complex networks of reactiondiffusion systems applied to a panic model[END_REF], whose adequateness with biological observations can legitimately be criticized (see Figure 1(a)). Furthermore, the local dynamics of the complex networks studied in these papers are often described by identical systems defined in identical domains; the equations of these complex networks can be written in the following form:

Ω Ω Ω Ω x x x x (a) Ω 1 Ω 2 Ω 3 Ω 4 (b) Ω 1 Ω 2 Ω 3 Ω 4 (c)
             ∂U i ∂t = D i ∆U i + f (U i ) -µ j∈Ni (U i -U j ) (t, x) ∈ (0, ∞) × Ω, ( 7
)
∂U i ∂ν i = 0, (t, x) ∈ (0, ∞) × ∂Ω, ( 8
)
U i (0, x) = U i,0 (x), x ∈ Ω, (9)
where Ω denotes a single bounded domain included in R M . Compared with the model determined by ( 1)-( 2)-( 3)-( 4), we observe that the couplings are here integrated into the reaction-diffusion equation [START_REF] Camara | Turing and Hopf patterns formation in a predator-prey model with Leslie-Gower-type functional response[END_REF], whereas they are defined by the Robin boundary condition (3) in our new model. Furthermore, the boundary condition [START_REF] Cantin | Dimension estimate of attractors for complex networks of reaction-diffusion systems applied to an ecological model[END_REF] of Neumann type is seen as incompatible with natural displacements of biological individuals through the boundary of their habitat. Thus it appears highly relevant to study the dynamics of complex networks for which a refined modeling of the couplings is expected.

Contributions.

The model studied in this paper, which results from the latter observations, represents the first contribution of our work: by improving the modeling of the geometry of the domains and of the couplings, this new model better takes into account the heterogeneity of the biological habitat and the complexity of the interactions occurring between distinct patches of the environment. Indeed, the domains Ω 1 , . . . , Ω n involved in the complex network (1)-( 2)-( 3)-( 4) are fully distinct, in the sense that they can admit distinct sizes and distinct shapes. Furthermore, the couplings between those domains are defined along their boundaries, which significantly improves the modeling of individuals displacements, which were so far modeled by point-wise couplings in the interior of the domains. Our second contribution corresponds to the mathematical analysis of this new model. Indeed, our aim is to establish sufficient conditions on the numerous parameters involved in equations ( 1)-( 2)-( 3)-( 4), in order to guaranty the global synchronization of the local dynamics on each domain Ω i , 1 ≤ i ≤ n. Reaching a synchronization state for interacting species systems is relevant, since the dynamics of such systems are usually characterized by the coexistence of extinction and persistence equilibria. Hence, a global synchronization state in the whole system is viewed as a control strategy to avoid extinction, and can be set in place by setting ecological corridors, so as to increase the connectivity of natural habitats and to avoid local extinction of several wildlife species. Indeed, from the ecological point of view, this topic represents a major challenge, since it is observed that anthropic activities exert a high pressure on the habitat of numerous wildlife species (see notably [START_REF] Haddad | Habitat fragmentation and its lasting impact on Earth's ecosystems[END_REF][START_REF] Naiman | The role of riparian corridors in maintaining regional biodiversity[END_REF][START_REF] Rodriguez-Iturbe | River networks as ecological corridors: A complex systems perspective for integrating hydrologic, geomorphologic, and ecologic dynamics[END_REF][START_REF] Salau | Varying effects of connectivity and dispersal on interacting species dynamics[END_REF][START_REF] Tambosi | A framework to optimize biodiversity restoration efforts based on habitat amount and landscape connectivity[END_REF][START_REF] Zhou | Construction of ecological security pattern for Plateau lake based on MSPA-MCR model: A case study of Dianchi lake area[END_REF]). In particular, urban growth and extension of agricultural land profoundly modify the landscape and the geometry of ecological habitats. One relevant strategy for mitigating the effects of anthropization of natural spaces is to implement ecological corridors. We show in Figure 2(a) a picture of such an ecological corridor, that connects both sides of a road. Furthermore, a high density of ecological corridors in a fragmented habitat can help restore wildlife (see Figure 2(b)).

To achieve our study of this innovative complex network model, we shall employ original methods, involving a Poincaré-type inequality for mixed Neumann-Robin boundary conditions. Up to our knowledge, such inequalities have never been used in previous related works. Furthermore, in researching sufficient conditions of synchronization for the complex network model ( 1)-( 2)-( 3)-( 4), we bring out that the diffusion in an isotropic habitat described by the Laplace operator is not compatible with the distinct geometries of the domains Ω 1 , . . . , Ω n . To overcome this delicate issue, we also address in this paper the case of diffusion in an anisotropic habitat, by replacing the standard reaction-diffusion equation (1) by the system of semi-linear equations

∂U i ∂t = D i ∇ Ωi • A i ∇ Ωi U i + f (U i ), (t, x i ) ∈ (0, ∞) × Ω i , ( 10
)
where the matrices A i , 1 ≤ i ≤ n, are square matrices of order M with real coefficients. Afterwards, we apply our theoretical approach to a nonlinear predator-prey model with Leslie-Gower-type functional response, given by

∂u ∂t = d 1 ∆u + u(1 -u) - Quv u + A , ∂v ∂t = d 2 ∆v + Sv 1 - v u + C , ( 11
)
where u, v denote the densities of competing species and d 1 ,d 2 , A, C S and Q denote positive constants. This competing species model has been studied on a single domain in numerous papers (see for instance [START_REF] Abid | Diffusion driven instability and Hopf bifurcation in spatial predator-prey model on a circular domain[END_REF], [START_REF] Aguirre | Three limit cycles in a Leslie-Gower predator-prey model with additive Allee effect[END_REF], [START_REF] Arancibia-Ibarra | Turing patterns in a diffusive Holling-Tanner predator-prey model with an alternative food source for the predator[END_REF], [START_REF] Banerjee | Turing instabilities and spatio-temporal chaos in ratio-dependent Holling-Tanner model[END_REF], [START_REF] Camara | Turing and Hopf patterns formation in a predator-prey model with Leslie-Gower-type functional response[END_REF], [START_REF] Liu | Turing patterns in a predator-prey model on complex networks[END_REF], [START_REF] Liu | Dynamics of a predator-prey system with nonlinear prey-taxis[END_REF], [START_REF] Toaha | Stability analysis of prey predator model with Holling II functional response and threshold harvesting for the predator[END_REF]). The dynamics of such competing species models have also been studied in non-convex domains admitting a "dumbbell" shape [START_REF] Matano | Pattern formation in competition-diffusion systems in nonconvex domains[END_REF][START_REF] Cantrell | Spatial ecology via reaction-diffusion equations[END_REF], which resembles a simple twonodes network. However, to the best of our knowledge, the dynamics of the nonlinear predator-prey model [START_REF] Cantin | On a generalized diffusion problem: A complex network approach[END_REF] have never been studied in a complex network with boundary couplings of the form (3).

In the numerical part of our paper, we focus on the synchronization of heterogeneous Turing patterns, which have been proved to appear in the Leslie-Gower model, in a diffusion-driven instability process. We exhibit various forms of these patterns, mainly labyrinths and spots, and show how to synchronize them. We emphasize that several other applications could be considered (including neural networks, epidemiological networks, behavioral networks as cited before), which shows the wide potential of our approach.

Outline. Our paper is organized as follows. In Section 2, after brief preliminaries for recalling important lemmas, we expose the functional context which guarantees the existence and uniqueness of global solutions to the complex network problem (1)-( 2)-( 3)-( 4). Next, we establish sufficient conditions of synchronization for the case of homothetic domains (Theorem 2), which are applied to complete graph topologies (corollary 1), simple two-nodes networks (corollary 2), cyclic graphs and rings of nearest neighbors topologies. We also investigate the situation of diffusion in an anisotropic habitat, described by equation [START_REF] Cantin | Complex network near-synchronization for non-identical predator-prey systems[END_REF], which covers for instance the case of ellipsoid domains of distinct eccentricities (Theorem 4). In section 3, we present several numerical simulations of the Leslie-Gower predator prey model, in order to supplement our qualitative statements by quantitative experiments, which suggest that synchronization of complex networks with boundary couplings might be more delicate to reach than with point-wise couplings.

Complex networks of reaction-diffusion systems with boundary couplings

In this section, we consider the complex network of reaction-diffusion systems determined by ( 1)-( 2)-( 3)-( 4). As presented in the introduction, each domain Ω i , 1 ≤ i ≤ n, of this complex network admits a finite set N i of neighbors. We denote by G = (V , E ) the underlying graph: the set V of its vertices corresponds to the domains Ω 1 , . . . , Ω n , and the set E of its edges is determined by the sets of neighbors N 1 , . . . , N n . An example of such a graph is illustrated in Figure 1(b).

Preliminaries

Let us here briefly present two important lemmas which shall be useful in the sequel of the paper. We first recall a Poincaré-type inequality for mixed Neumann-Robin boundary conditions (see [START_REF] Smoller | Shock waves and reaction-diffusion equations[END_REF], inequality (11.13) in Theorem 11.11).

Lemma 1.

Let Ω be a bounded domain in R M with regular boundary ∂Ω.

Assume that ∂Ω = Γ 1 ∪ Γ 2 ,
where Γ 2 has a positive measure, with

Γ 1 ∩ Γ 2 = ∅. Let µ ∈ L ∞ (Γ 2 ) be such that µ(x) ≥ µ 0 for all x ∈ Γ 2
, where µ 0 denotes a positive constant. Then there exists a positive constant λ P (Ω, µ) such that

Ω |∇u| 2 dx ≥ λ P (Ω, µ) Ω |u| 2 dx, ( 12
)
for all u ∈ W 1,2 (Ω) satisfying

∂u ∂ν = 0 on Γ 1 , ∂u ∂ν = -µ(x)u on Γ 2 .
In the latter lemma, the positive constant λ P is called the Poincaré constant; it depends on the domain Ω and on the function µ (see Remark 3 below).

Next, recall that a complete graph is a graph in which every pair of distinct vertices is connected by a unique edge. The following lemma establishes a necessary and sufficient condition for a connected graph to be complete. Since its proof is elementary, we may omit it.

Lemma 2. Let G = (V , E ) be a connected graph with a finite set V of n vertices (n ≥ 3) and a finite set E of edges. For each i ∈ V , we denote by N i the subset of vertices which are neighbors of i. Then G is complete if and only if

N i \ {j} = N j \ {i}, (13)
for all i, j ∈ V such that i and j are neighbors.

Global solutions of the complex network problem

Before establishing sufficient conditions of synchronization in the complex network determined by problem ( 1)-( 2)-( 3)-( 4), we intend to prove the existence and uniqueness of global solutions. To that aim, we consider the Banach space

Y = L p (Ω 1 ) × • • • × L p (Ω n ) m ,
equipped with the usual product norm. Following [START_REF] Cantrell | Spatial ecology via reaction-diffusion equations[END_REF] (Theorem 1.11), we require that the Lebesgue exponent p satisfies p > M , where M denotes the dimension space of the domains Ω 1 , . . . , Ω n . We then introduce the linear operator A defined by

A = diag {-D i ∆ Ωi , 1 ≤ i ≤ n} , (14)
with the mixed boundary conditions ( 2) and ( 3). The domain D(A) of the linear operator A is included in the space X defined as

X = U = (U i ) 1≤i≤n ∈ W 2,p (Ω 1 ) × • • • × W 2,p (Ω n ) m ; ∂U i ∂ν i = 0 on Γ N i , ∂U i ∂ν i = -µ i (x i ) j∈Ni (U i -U j ) on Γ R i .
Next, we consider the function F defined for U ∈ X by

F (U ) = f k (U i ) 1≤i≤n,1≤k≤m . ( 15
)
The abstract formulation of the complex network problem (1)-( 2)-( 3)-( 4) is written

   dU dt + AU = F (U ), t > 0, U (0) = U 0 , ( 16
)
with

U 0 = (U i,0 ) 1≤i≤n ∈ X.
Assuming that the function F is locally Lipschitz continuous on any bounded set of R n×m , we can apply Theorem 1.11 in [START_REF] Cantrell | Spatial ecology via reaction-diffusion equations[END_REF] and conclude that the complex network problem (1)-( 2)-( 3)-( 4) admits global solutions.

Theorem 1. Assume that the function F given by ( 15) is locally Lipschitz continuous on any bounded set of R n×m . Then there exists γ ∈ (0, 1) such that the complex network problem (1)-( 2)-( 3)-( 4) admits global solutions in the space X γ ⊂ Y generated by the fractional power A γ , and determines a dynamical system S (t), defined for all U 0 ∈ X γ and all t ≥ 0 by

S (t)U 0 = U (t, U 0 ), ( 17
)
where U (t, U 0 ) denotes the unique solution of the Cauchy problem [START_REF] Ha | Remarks on the complete synchronization of Kuramoto oscillators[END_REF] in X γ .

Here, we assume moreover that the continuous dynamical system S (t) admits a compact phase space Φ ⊂ Y , which is bounded in D(A). We emphasize that this compactness requirement is fulfilled for a wide class of reaction-diffusion systems (see for instance [START_REF] Yagi | Abstract parabolic evolution equations and their applications[END_REF]) and often follows from some dissipative estimations of the solutions. Along with the latter assumptions on the phase space Φ, we easily verify that the functions (f k ) 1≤k≤m satisfy the following globally Lipschitz condition:

f k (v) -f k (w) L 2 (Ω) ≤ L k v -w L 2 (Ω) , 1 ≤ k ≤ m, ( 18
)
for all v, w in Φ, with L k > 0. This globally Lispchitz condition shall be very useful in the rest of the paper, for establishing sufficient conditions of synchronization of the complex network problem ( 1)-( 2)-( 3)-(4).

Sufficient conditions of synchronization

In this section, we investigate sufficient conditions for the complex network model with boundary couplings ( 1)-( 2)-( 3)-( 4) to synchronize. Let us first specify the definition of synchronization for complex networks of reaction-diffusion systems. We distinguish partial synchronization of a subset of nodes in the network, and global synchronization in the whole network. Since the domains Ω 1 , . . . , Ω n underlying the complex network model ( 1)-( 2)-( 3)-( 4) are non-identical, we are led to transport the solutions in a common arbitrary domain, hence the following definition differs from the identical synchronization studied for instance in [START_REF] Cantin | Dimension estimate of attractors for complex networks of reaction-diffusion systems applied to an ecological model[END_REF].

Definition 1.

Let Ω be a bounded domain in R M . We say that two nodes i and j synchronize in Ω if there exist two homeomorphisms φ i and φ j mapping Ω onto Ω i and Ω j , respectively, such that for any U 0 ∈ X, the solution U (t, U 0 ) of the complex network problem ( 1)-( 2)-( 3)-( 4) satisfies

lim t→+∞ Ũi (t, x) -Ũj (t, x) L 2 (Ω) m = 0,
where Ũi (t, x) = U i t, φ i (x) and Ũj (t, x) = U j t, φ j (x) for all x ∈ Ω and all t > 0.

We say that the complex network partially synchronizes in Ω if there exist at least one pair of nodes (i, j) that synchronizes in Ω.

We say that the complex network globally synchronizes in Ω if every pair of nodes synchronizes.

Note that the homeomorphism φ i,j defined in [START_REF] Belykh | Connection graph stability method for synchronized coupled chaotic systems[END_REF], that maps Γ R i onto Γ R j , can be defined as the continuation of the composition φ j • φ -1 i that maps Ω i onto Ω j . We are now ready to investigate sufficient conditions of synchronization in the complex network (1)-( 2)-( 3)-( 4). In the sequel, we focus on two situations: first, we establish a theorem for the synchronization of two nodes in a network of homothetic domains; we deduce several statements of global synchronization for complete graph topologies, cyclic graphs and rings of nearest neighbors graph topologies. Next, we consider the more general case of non-homothetic domains, for which we are led to replace the diffusion operator in an isotropic habitat by a generalized diffusion operator for an anisotropic habitat.

Synchronization conditions for a complex network of homothetic domains

We first investigate the case where the domains Ω 1 , . . . , Ω n of the complex network problem (1)-( 2)-( 3)-( 4) are homothetic, that is, admit the same shape, with possibly distinct sizes. To that aim, we fix

Ω = Ω 1 , Γ R = Γ R 1 , Γ N = Γ N 1 , ( 19
)
and consider n -1 homotheties H 2 , . . . , H n mapping Ω 1 onto Ω 2 , . . . , Ω n , respectively:

Ω i = H i (Ω 1 ), 2 ≤ i ≤ n. ( 20
)
Examples of such homothetic domains can be viewed in Figures 8, 10 below. For each i ∈ {2, . . . , n}, the homothety factor of H i is positive and denoted by h i . We assume that the homotheties preserve the splitting of the boundaries of each domain Ω i , that is

Γ N i = H i (Γ N 1 ), Γ R i = H i (Γ R 1 ), 2 ≤ i ≤ n. ( 21
)
Given an initial condition U 0 ∈ X γ , we denote by

U (t, x) = U 1 (t, x), . . . , U n t, H n (x)
the solution of problem ( 1)-( 2)-( 3)-( 4) starting from U 0 . For all x ∈ Ω and all t > 0, we set

Ũ1 (t, x) = U 1 (t, x), Ũi (t, x) = U i t, H i (x) , 2 ≤ i ≤ n. ( 22
)
In addition, we assume that the functions (f k ) 1≤i≤n,1≤k≤m satisfy the global Lipschitz condition [START_REF] Hecht | New development in FreeFem++[END_REF], with L k > 0. The following theorem establishes sufficient conditions on the parameters involved in the complex network ( 1)-( 2)-( 3)-( 4) for two nodes to synchronize.

Theorem 2. Let assumptions [START_REF] Hecht | New development in FreeFem++[END_REF] and (20) hold. Suppose that two neighbor nodes i and j of the complex network (1)-( 2)-( 3)-( 4) satisfy

D i,k h 2 i = D j,k h 2 j , ( 23
)
D i,k h i µ i H i (x) = D j,k h j µ j H j (x) , ( 24
)
for each k ∈ {1, . . . , m} and for all x ∈ Γ R . Suppose moreover that

N i \ {j} = N j \ {i}. ( 25
)
Finally, suppose that

D i,k λ P (Ω, µ i ) h 2 i > 1 + L k 2 2 , ( 26
)
for each k ∈ {1, . . . , m}, where λ P (Ω, µ i ) is the Poincaré constant determined by Lemma 1.

Then the two nodes i and j of the complex network synchronize in Ω. The proof of Theorem 2 relies of the following technical lemma, which establishes an estimate of the boundary terms modeling the migrations of individuals between two neighbor nodes.

(a) (b) (c) (d) (e)

Lemma 3.

Let assumption [START_REF] Liu | Turing patterns in a predator-prey model on complex networks[END_REF] hold. Suppose that two neighbor nodes i and j of the complex network (1)-( 2)-( 3)-( 4) satisfy [START_REF] Maria Da Conceiçao | Homogeneous three-cell networks[END_REF] and [START_REF] Matano | Pattern formation in competition-diffusion systems in nonconvex domains[END_REF] for each k ∈ {1, . . . , m} and for all x ∈ Γ R . Then we have

D i,k h 2 i ∂ Ũi,k ∂ν (t, x) - ∂ Ũj,k ∂ν (t, x) Ũi,k (t, x) -Ũj,k (t, x) = -(2 + |S i,j |) D i,k h i µ i H i (x) Ũi,k (t, x) -Ũj,k (t, x) 2 - D i,k h i µ i H i (x) l∈Mi,j Ũi,k (t, x) -Ũl,k (t, x) Ũi,k (t, x) -Ũj,k (t, x) + D i,k h i µ i H i (x) l∈Mj,i Ũj,k (t, x) -Ũl,k (t, x) Ũi,k (t, x) -Ũj,k (t, x) ( 27
)
for all t ≥ 0 and all x ∈ Γ R , where M i,j denotes the set of neighbors of i which are not neighbors of j, M j,i denotes the set of neighbors of j which are not neighbors of i, S i,j denotes the set of common neighbors of i and j and |S i,j | its cardinal.

Proof of Lemma 3. Let k ∈ {1, . . . , m} and i, j ∈ {1, . . . , n}. For short, we denote

u(t, x i ) = U i,k (t, x i ), (t, x i ) ∈ (0, +∞) × Ω i , v(t, x j ) = U j,k (t, x j ), (t, x j ) ∈ (0, +∞) × Ω j and ũ(t, x) = Ũi,k (t, x), ṽ(t, x) = Ũj,k (t, x), (t, x) ∈ (0, +∞) × Ω.
We examine the boundary terms of Robin type on Γ R . We have:

∂ ũ ∂ν (t, x) = ν(x) • ∇ Ω ũ(t, x) = h i × ν i H i (x) • ∇ Ωi u t, H i (x) = h i × ∂u ∂ν i t, H i (x) = -h i µ i H i (x) l∈Ni u t, H i (x) -U l,k t, H i (x) , which leads to ∂ ũ ∂ν (t, x) = -h i μi (x) l∈Ni ũ(t, x) -U l,k t, H i (x) , ( 28
)
where μi is defined on Γ R by μi (x) = µ i H i (x) ; similarly, we have

∂ṽ ∂ν (t, x) = -h j μj (x) l∈Nj ṽ(t, x) -U l,k t, H j (x) .
We obtain, applying assumption [START_REF] Maria Da Conceiçao | Homogeneous three-cell networks[END_REF]:

D i,k h 2 i ∂ ũ ∂ν (t, x) - ∂ṽ ∂ν (t, x) = D i,k h 2 i ∂ ũ ∂ν (t, x) - D j,k h 2 j ∂ṽ ∂ν (t, x) = - D i,k h i μi (x) l∈Ni ũ(t, x) -U l,k t, H i (x) + D j,k h j μj (x) l∈Nj ṽ(t, x) -U l,k t, H j (x) = - D i,k h i μi (x) ũ(t, x) -ṽ(t, x) - D i,k h i μi (x) l∈Ni\{j} ũ(t, x) -U l,k t, H i (x) + D j,k h j μj (x) ṽ(t, x) -ũ(t, x) + D j,k h j μj (x) l∈Nj \{i} ṽ(t, x) -U l,k t, H j (x) .
Now, the boundary convention [START_REF] Belykh | Connection graph stability method for synchronized coupled chaotic systems[END_REF] guarantees

that l∈Ni\{j} U k,l t, H i (x) = l∈Ni\{j} U k,l t, H l (x) ,
and analogously

l∈Nj \{i} U k,l t, H j (x) = l∈Nj \{i} U k,l t, H l (x) .
Next, by virtue of assumption [START_REF] Matano | Pattern formation in competition-diffusion systems in nonconvex domains[END_REF], we obtain:

D i,k h 2 i ∂ ũ ∂ν (t, x) - ∂ṽ ∂ν (t, x) = -2 D i,k h i μi (x) ũ(t, x) -ṽ(t, x) - D i,k h i μi (x) l∈Si,j ũ(t, x) -ṽ(t, x) - D i,k h i μi (x) l∈Mi,j ũ(t, x) -Ũl,k (t, x) + D i,k h i μi (x) l∈Mj,i ṽ(t, x) -Ũl,k (t, x) = -(2 + |S i,j |) D i,k h i μi (x) ũ(t, x) -ṽ(t, x) - D i,k h i μi (x) l∈Mi,j ũ(t, x) -Ũl,k (t, x) + D i,k h i μi (x) l∈Mj,i ṽ(t, x) -Ũl,k (t, x) ,
where |S i,j | denotes the cardinal of the set S i,j . Multiplying both sides of the latter equality by ũ(t, x) -ṽ(t, x) leads to [START_REF] Rätz | Symmetry breaking in a bulk-surface reaction-diffusion model for signalling networks[END_REF], which completes the proof of Lemma 3.

We are now ready to present the proof of Theorem 2.

Proof of Theorem 2. We use the same notations as in the proof of Lemma 3 and consider the energy functional defined by

E i,j,k (t) = 1 2 Ω |ũ(t, x) -ṽ(t, x)| 2 dx. ( 29
)
We easily compute its derivative:

dE i,j,k dt (t) = Ω ∂(ũ -ṽ) ∂t (ũ -ṽ)dx = Ω ∂ ũ ∂t - ∂ṽ ∂t (ũ -ṽ)dx,
where we omit the variables t and x under the integral symbol in order to lighten our notations. Elementary computations show that

∇ Ω ũ(t, x) = h i ∇ Ωi u t, H i (x) , (t, x) ∈ (0, +∞) × Ω, (30)
from which it follows that

∆ Ω ũ(t, x) = h 2 i ∆ Ωi u t, H i (x) , (t, x) ∈ (0, +∞) × Ω. ( 31
)
We can deduce that

∂ ũ ∂t (t, x) = ∂u ∂t t, H i (x) = D i,k ∆ Ωi u(t, H i (x) + f k u(t, H i (x)) = D i,k h 2 i ∆ Ω ũ(t, x) + f k ũ(t, x) ,
for all (t, x) ∈ (0, +∞) × Ω. Similarly, we have

∂ṽ ∂t (t, x) = D j,k h 2 j ∆ Ω ṽ(t, x) + f k ṽ(t, x) , (t, x) ∈ (0, +∞) × Ω.
We obtain, by virtue of assumption [START_REF] Maria Da Conceiçao | Homogeneous three-cell networks[END_REF]:

dE i,j,k dt (t) = D i,k h 2 i Ω ∆ Ω (ũ -ṽ)(ũ -ṽ)dx + Ω f k (ũ) -f k (ṽ) (ũ -ṽ)dx. ( 32
)
Since ũ and ṽ satisfy a mixed Neumann-Robin boundary condition, we have, by virtue of Green's formula:

Ω ∆ Ω (ũ -ṽ)(ũ -ṽ)dx = - Ω |∇ Ω (ũ -ṽ)| 2 dx + Γ R ∂(ũ -ṽ) ∂ν (ũ -ṽ)ds = - Ω |∇ Ω (ũ -ṽ)| 2 dx + Γ R ∂ ũ ∂ν - ∂ṽ ∂ν (ũ -ṽ)ds,
where ν = ν 1 denotes the outward unit normal vector of ∂Ω = ∂Ω 1 .

Next, assumption [START_REF] Miranville | Bifurcations and synchronization in networks of unstable reaction-diffusion systems[END_REF] means that the nodes i and j admit only common neighbors, which implies that M i,j = M j,i = ∅ and S i,j = N i \ {j} = N j \ {i}, where the sets M i,j , M j,i and S i,j are defined in Lemma 3. Therefore, we can apply Lemma 3 and we obtain

D i,k h 2 i ∂ ũ ∂ν - ∂ṽ ∂ν (ũ -ṽ) = -(2 + |S i,j |) D i,k h i μi (x) |ũ -ṽ| 2 ≤ 0 on Γ R , from which we deduce Ω ∆ Ω (ũ -ṽ)(ũ -ṽ)dx ≤ - Ω |∇ Ω (ũ -ṽ)| 2 dx.
Now, we use the Poincaré inequality (see Lemma 1) to write

Ω |∇ Ω (ũ -ṽ)| 2 dx ≥ λ P (Ω, µ i ) Ω |ũ -ṽ| 2 dx,
with λ P (Ω, µ i ) > 0, which leads to

D i,k h 2 i Ω ∆ Ω (ũ -ṽ)(ũ -ṽ)dx ≤ - D i,k λ P (Ω, µ i ) h 2 i Ω |ũ -ṽ| 2 dx. ( 33
)
In parallel, we apply the Young inequality ab ≤ a 2 2 + b 2 2 , which is valid for all a, b in R, and assumption [START_REF] Hecht | New development in FreeFem++[END_REF] to obtain

Ω f k (ũ) -f k (ṽ) (ũ -ṽ)dx ≤ 1 2 f k (ũ) -f k (ṽ) 2 L 2 (Ω) + 1 2 ũ -ṽ 2 L 2 (Ω) ,
which leads to

Ω f k (ũ) -f k (ṽ) (ũ -ṽ)dx ≤ (1 + L k 2 )E i,j,k (t). ( 34
)
Combining the inequalities ( 33) and [START_REF] Verma | Network synchronization, stability and rhythmic processes in a diffusive mean-field coupled SEIR model[END_REF], equation [START_REF] Tambosi | A framework to optimize biodiversity restoration efforts based on habitat amount and landscape connectivity[END_REF] becomes

dE i,j,k dt (t) ≤ -2 D i,k λ P (Ω, µ i ) h 2 i E i,j,k (t) + (1 + L 2 k )E i,j,k (t). ( 35
)
Applying Gronwall lemma leads to

E i,j,k (t) ≤ E i,j,k (0) exp -2 D i,k λ P (Ω, µ i ) h 2 i -(1 + L 2 k ) t ,
for each k ∈ {1, . . . , m} and for all t ≥ 0, which guarantees that lim t→+∞ Ũi (t) -Ũj (t) L 2 (Ω) m = 0, provided [START_REF] Naiman | The role of riparian corridors in maintaining regional biodiversity[END_REF] holds. The proof is complete.

Remark 1 (Non-identical internal dynamics and near-synchronization). We emphasize that the function f involved in the reaction-diffusion system (1) is the same on all the domains Ω 1 , . . . , Ω n , which means that the internal dynamics of the complex network are identical on each patch of the fragmented environment. The case of non-identical dynamics can be studied by replacing Equation (1) by

∂U i ∂t = D i ∆ Ωi U i + f i (U i ), (t, x i ) ∈ (0, ∞) × Ω i ,
where the function f i now depends on the index i of the domain

Ω i (1 ≤ i ≤ n).
In particular, this dependence can model local variations of internal parameters, and f i can be written under the form

f i (U i ) = f (U i , λ i ), 1 ≤ i ≤ n,
with λ i ∈ R p . In this case, the Lipschitz condition [START_REF] Hecht | New development in FreeFem++[END_REF] becomes

f i,k (v) -f j,k (w) L 2 (Ω) ≤ L i,j,k v -w L 2 (Ω) + λ i -λ j R p ,
with L i,j,k > 0. Recently, it has been proved in [START_REF] Cantin | Complex network near-synchronization for non-identical predator-prey systems[END_REF] that this situation leads to near-synchronization, which is a relaxed form of synchronization.

Remark 2 (Interpretation of assumption ( 23)). If the domains Ω i and Ω j in Theorem 2 are isometric, then the homothety factors h i and h j are both equal to 1. In that case, condition [START_REF] Maria Da Conceiçao | Homogeneous three-cell networks[END_REF] simply becomes

D i,k = D j,k ,
for all k ∈ {1, . . . , m}. The latter condition means that for each k, the species U i,k and U j,k living in the domains Ω i and Ω j , respectively, should diffuse at the same rate in their respective habitat, in order to synchronize, which is conform to intuition. However, it does not imply that two distinct species U i,k1 and U i,k2 of the same domain Ω i should diffuse at the same rate. Since spatial heterogeneity patterns such as Turing patterns usually appear as the diffusion rates are very distinct, then assumption [START_REF] Maria Da Conceiçao | Homogeneous three-cell networks[END_REF] seems to be compatible with the emergence of Turing patterns. This compatibility will be experimented with a numerical approach in Section 3.

Next, if the domains Ω i and Ω j in Theorem 2 are not isometric, then assumption [START_REF] Maria Da Conceiçao | Homogeneous three-cell networks[END_REF] means that synchronization can occur if the species which live in a large domain diffuse at a greater rate than the species which live in a small domain; this interpretation is also conform to intuition.

Remark 3 (Interpretation of assumption ( 26)). As mentioned previously, we emphasize that the Poincaré constant λ(Ω, µ i ) in the synchronization condition [START_REF] Naiman | The role of riparian corridors in maintaining regional biodiversity[END_REF] depends on the domain Ω and on the boundary coupling function µ i . More specifically, according to Corollaries 2.2 and 2.3 in [START_REF] Cantrell | Spatial ecology via reaction-diffusion equations[END_REF] (see also [START_REF] Daners | A Faber-Krahn inequality for Robin problems in any space dimension[END_REF] or [START_REF] Savo | Optimal eigenvalue estimates for the Robin Laplacian on Riemannian manifolds[END_REF]), λ(Ω, µ i ) decreases if the size of Ω increases, whereas it increases if µ i increases (in the sense that

µ 1 i ≤ µ 2 i if µ 1 i (x) ≤ µ 2 i (x) for all x ∈ Γ R i ).
Therefore, the synchronization condition (26) is more likely to be fulfilled in small domains, connected with strong boundary couplings. It is worth noting that the synchronization of complex networks with point-wise couplings of the form (7)-( 8)-( 9) is not influenced by the size of the domains, as proved in [START_REF] Cantin | Dimension estimate of attractors for complex networks of reaction-diffusion systems applied to an ecological model[END_REF] for instance. Roughly speaking, boundary couplings of the form (3) are able to synchronize the local dynamics in a neighborhood of the boundaries; if the domains are small, this boundary synchronization can extend to the whole domain.

As a direct consequence of Theorem 2 and Lemma 2, we obtain the following corollary, which establishes sufficient conditions for a complex network, whose underlying graph is complete, to globally synchronize. [START_REF] Hecht | New development in FreeFem++[END_REF] and (20) hold. Suppose furthermore that the graph G underlying the complex network problem (1)-( 2)-( 3)-( 4) is a complete graph, and that assumptions (23)-( 24)-( 26) hold for each i, j ∈ {1, . . . , n} and each k ∈ {1, . . . , m}. Then the complex network globally synchronizes in Ω.

Corollary 1. Let assumptions

In the case of two-nodes network, the following corollary is also directly obtained from Theorem 2. [START_REF] Hecht | New development in FreeFem++[END_REF] and (20) hold with n = 2. Suppose furthermore that assumptions (23)-( 24)-( 26) hold for i = 1, j = 2 and for each k ∈ {1, . . . , m}. Then the two-nodes network globally synchronizes in Ω.

Corollary 2. Let assumptions

Synchronization conditions for cyclic graphs and rings of nearest neighbors

Corollaries 1 and 2 describe the global synchronization of a complex network whose underlying graph is complete. However, the complete graph topology requirement is rather restrictive. Indeed, from the ecological point of view, fragmented habitats are rarely densely connected, so that the corresponding graph topology is often far from the complete graph topology. Hence, it is important to investigate the dynamics of a complex network whose underlying graph is not complete. In this section, we analyze two other important configurations. The first configuration corresponds to a cyclic graph. We recall that a cyclic graph is a graph with n vertices, that can be numbered in such a way that node i is connected only to nodes i -1 mod(n) and i + 1 mod(n) (see Figure 3(c)). The second configuration corresponds to a ring of 2K-nearest neighbors [START_REF] Belykh | Connection graph stability method for synchronized coupled chaotic systems[END_REF]. If K is an integer such that 1 ≤ K ≤ n 2 , we say that the graph G is a ring of 2K-nearest neighbors if it admits n vertices that can be numbered in such a way that i is connected only to i ± j mod(n) for each j ∈ {1, . . . , K} (for example, a ring of 4-nearest neighbors and a ring of 6-nearest neighbors are depicted in Figure 3 (d) and (e), respectively). The following Theorem proves that such graph topologies, although weakly densely connected, can globally synchronize.

Theorem 3. Let assumptions [START_REF] Hecht | New development in FreeFem++[END_REF] and (20) hold. Suppose that the graph G underlying the complex network (1)-( 2)-( 3)-( 4) is a cyclic graph or a ring of 2K-nearest neighbors (with K ≥ 1). Finally, assume that each pair (i, j) of nodes satisfies [START_REF] Maria Da Conceiçao | Homogeneous three-cell networks[END_REF], ( 24) and [START_REF] Naiman | The role of riparian corridors in maintaining regional biodiversity[END_REF].

Then the complex network (1)-( 2)-( 3)-( 4) globally synchronizes in Ω.

Proof. Assume for simplicity that the graph G underlying the complex network problem (1)-( 2)-( 3)-( 4) is a cycle. Renumbering the vertices of the graph G if necessary, we can without loss of generality assume that the sets of neighbors are given by

N 1 = {n, 2}, N i = {i -1, i + 1} for 2 ≤ i ≤ n -1, N n = {n -1, 1}. (36)
Now, let k ∈ {1, . . . , n}. We consider the total energy function along the cyclic graph G defined by

E k (t) = n-1 i=1 E i,i+1,k (t) + E n,1,k (t), (37)
where E i,j,k (t) is defined as in [START_REF] Salau | Varying effects of connectivity and dispersal on interacting species dynamics[END_REF]. For 1 ≤ i ≤ n -1, we compute as in the proof of Theorem 2:

dE i,i+1,k dt (t) = Ω ∂ Ũi,k ∂t - ∂ Ũi+1,k ∂t Ũi,k -Ũi+1,k dx = - D i,k h 2 i Ω ∇ Ω Ũi,k -Ũi+1,k 2 dx + D i,k h 2 i Γ R ∂ Ũi,k -Ũi+1,k ∂ν Ũi,k -Ũi+1,k ds + Ω f k ( Ũi,k ) -f k ( Ũi+1,k ) Ũi,k -Ũi+1,k dx.
Similarly, we have

dE n,1,k dt (t) = - D n,k h 2 n Ω ∇ Ω Ũn,k -Ũ1,k 2 dx + D n,k h 2 n Γ R ∂ Ũn,k -Ũ1,k ∂ν Ũn,k -Ũ1,k ds + Ω f k ( Ũn,k ) -f k ( Ũ1,k ) Ũn,k -Ũ1,k dx.
Now, we examine the sum of the integral terms

D i,k h 2 i Γ R ∂ Ũi,k -Ũi+1,k ∂ν Ũi,k -Ũi+1,k ds, 1 ≤ i ≤ n -1, and D n,k h 2 n Γ R ∂ Ũn,k -Ũ1,k ∂ν Ũn,k -Ũ1,k ds.
Since ( 23) and ( 24) are satisfied for each pair (i, j) of nodes, we have:

n-1 i=1 D i,k h 2 i Γ R ∂ Ũi,k -Ũi+1,k ∂ν Ũi,k -Ũi+1,k ds + D n,k h 2 n Γ R ∂ Ũn,k -Ũ1,k ∂ν Ũn,k -Ũ1,k ds = D 1,k h 2 1 Γ R n-1 i=1 ∂ Ũi,k -Ũi+1,k ∂ν Ũi,k -Ũi+1,k + ∂ Ũn,k -Ũ1,k ∂ν Ũn,k -Ũ1,k ds.
Next, since the graph G is cyclic, the neighbors of each vertex are given by [START_REF] Yagi | Abstract parabolic evolution equations and their applications[END_REF]. Hence, the boundary condition (3) leads to:

D 1,k h 2 1 Γ R n-1 i=1 ∂ Ũi,k -Ũi+1,k ∂ν Ũi,k -Ũi+1,k + ∂ Ũn,k -Ũ1,k ∂ν Ũn,k -Ũ1,k ds = D 1,k h 1 Γ R μ1 (x) -2( Ũ1,k -Ũ2,k ) 2 -( Ũ1,k -Ũn,k )( Ũ1,k -Ũ2,k ) + ( Ũ2,k -Ũ3,k )( Ũ1,k -Ũ2,k ) -2( Ũ2,k -Ũ3,k ) 2 -( Ũ2,k -Ũ1,k )( Ũ2,k -Ũ3,k ) + ( Ũ3,k -Ũ4,k )( Ũ2,k -Ũ3,k ) -2( Ũ3,k -Ũ4,k ) 2 -( Ũ3,k -Ũ2,k )( Ũ3,k -Ũ4,k ) + ( Ũ4,k -Ũ5,k )( Ũ3,k -Ũ4,k ) . . . -2( Ũn,k -Ũ1,k ) 2 -( Ũn,k -Ũn-1,k )( Ũn,k -Ũ1,k ) + ( Ũ1,k -Ũ2,k )( Ũn,k -Ũ1,k ) ds.
We can rearrange the terms which are contained in the brackets and write:

D 1,k h 2 1 Γ R n-1 i=1 ∂ Ũi,k -Ũi+1,k ∂ν Ũi,k -Ũi+1,k + ∂ Ũn,k -Ũ1,k ∂ν Ũn,k -Ũ1,k ds = D 1,k h 1 Γ R μ1 (x) -( Ũ1,k -Ũ2,k ) 2 + 2( Ũ1,k -Ũ2,k )( Ũn,k -Ũ1,k ) -( Ũn,k -Ũ1,k ) 2 -( Ũ2,k -Ũ3,k ) 2 + 2( Ũ2,k -Ũ3,k )( Ũ1,k -Ũ2,k ) -( Ũ1,k -Ũ2,k ) 2 -( Ũ3,k -Ũ4,k ) 2 + 2( Ũ3,k -Ũ4,k )( Ũ2,k -Ũ3,k ) -( Ũ2,k -Ũ3,k ) 2 . . . -( Ũn,k -Ũ1,k ) 2 + 2( Ũn,k -Ũ1,k )( Ũn,k -Ũn-1,k ) -( Ũn,k -Ũn-1,k ) 2 ds = D 1,k h 1 Γ R μ1 (x) -( Ũ1,k -Ũ2,k -Ũn,k + Ũ1,k ) 2 -( Ũ2,k -Ũ3,k -Ũ1,k + Ũ2,k ) 2 -( Ũ3,k -Ũ4,k -Ũ2,k + Ũ3,k ) 2 . . . -( Ũn,k -Ũ1,k -Ũn,k + Ũn-1,k ) 2 ds.
Therefore, we have

D 1,k h 2 1 Γ R n-1 i=1 ∂ Ũi,k -Ũi+1,k ∂ν Ũi,k -Ũi+1,k + ∂ Ũn,k -Ũ1,k ∂ν Ũn,k -Ũ1,k ds ≤ 0.
It follows that the derivative of the total energy E k given by ( 37) can be simplified and written

dE k dt (t) = n-1 i=1 - D i,k h 2 i Ω ∇ Ω Ũi,k -Ũi+1,k 2 dx - D n,k h 2 n Ω ∇ Ω Ũn,k -Ũ1,k 2 dx + n-1 i=1 Ω f k ( Ũi,k ) -f k ( Ũi+1,k ) Ũi,k -Ũi+1,k dx + Ω f k ( Ũn,k ) -f k ( Ũ1,k ) Ũn,k -Ũ1,k dx.
Hence, arguing as in the proof of Theorem 2, we can write

dE k dt (t) ≤ n-1 i=1 -2 D i,k λ P (Ω, µ i ) h 2 i -(1 + L 2 k ) E i,i+1,k -2 D n,k λ P (Ω, µ n ) h 2 n -(1 + L 2 k ) E n,1,k , which leads to dE k dt (t) ≤ -C * k E k (t)
, for all t ≥ 0, where C * k is the positive constant defined by

C * k = min min 1≤i≤n-1 2 D i,k λ P (Ω, µ i ) h 2 i -(1 + L 2 k ) , 2 D n,k λ P (Ω, µ n ) h 2 n -(1 + L 2 k ) .
Finally, we apply Gronwall lemma to write

E k (t) ≤ E k (0)e -C * k t ,
for all t ≥ 0 and all k ∈ {1, . . . , m}, which proves that the complex network ( 1)-( 2)-( 3)-( 4) globally synchronizes in Ω. This completes the proof for a cyclic graph. The case of a ring of 2K-nearest neighbors is treated analogously.

Synchronization conditions for non-homothetic domains

We now investigate the case of non-homothetic domains. We observe that the proof of Theorem 2 mostly relies on equation [START_REF] Savo | Optimal eigenvalue estimates for the Robin Laplacian on Riemannian manifolds[END_REF], which expresses how the gradient is modified under the action of a given transformation. However, equation ( 30) is characteristic from homotheties, thus cannot be fulfilled for a general planar transformation (as for instance, a transformation that maps an ellipse onto an another ellipse of distinct eccentricity). Therefore, we are led to consider a modification of the Laplace diffusion operator involved in equation ( 1), so as to guaranty that equations ( 31) and [START_REF] Rodriguez-Iturbe | River networks as ecological corridors: A complex systems perspective for integrating hydrologic, geomorphologic, and ecologic dynamics[END_REF], which are obtained from [START_REF] Savo | Optimal eigenvalue estimates for the Robin Laplacian on Riemannian manifolds[END_REF], are still fulfilled under the action of a non-homothetic transformation.

To that aim, we consider a complex network of semi-linear equations determined by [START_REF] Cantin | Complex network near-synchronization for non-identical predator-prey systems[END_REF] instead of (1); the Neumann and Robin boundary conditions ( 2)-( 3) and the initial condition (4) are conserved, where the derivative with respect to the outward normal unit vector is now defined by

∂U i ∂ν i (t, x i ) = ν i (x i ) • A i ∇ Ωi U i (t, x i ), (38)
with t > 0 and x i ∈ ∂Ω i (see for instance equation (2.23) in [START_REF] Yagi | Abstract parabolic evolution equations and their applications[END_REF]). For simplicity, we assume that the matrices A i involved in equation ( 10) are invertible and diagonal; this assumption guarantees commutativity with the gradient operator, that is

(A i ∇ Ωi )u = ∇ Ωi (A i u), ( 39
)
for all u ∈ W 1,p (Ω i ). Note that the existence and uniqueness of global solutions to the complex network determined by ( 10)-( 2)-( 3)-( 4) can be treated as in Section 2.2. Next, we set as before Ω = Ω 1 , Γ = Γ 1 and suppose that there exists a family (B i ) 2≤i≤n of n -1 planar transformations mapping Ω onto Ω 2 , . . . , Ω n , respectively, preserving the boundaries Γ R i , Γ N i (2 ≤ i ≤ n). We assume that for each i ∈ {2, . . . , n}, the transformation B i satisfies the property

∇ Ω (u • B i ) = B i (∇ Ωi u) • B i , ( 40
)
for all u ∈ W 1,p (Ω i ), where the matrix B i is invertible and diagonal. Unlike the case of homotheties, the planar tranformation is likely to rotate the outward normal vector ν Ω ; thus we denote

ν Ωi • B i (x) = Θ i (x)ν Ω (x), x ∈ Ω, 2 ≤ i ≤ n, ( 41
)
where Θ i (x) is a square matrix of order 2.

The following theorem establishes sufficient conditions of synchronization for two nodes of the complex network determined by problem ( 10)-( 2)-( 3)-( 4). Theorem 4. Let assumptions (18), ( 39) and (40) hold. Suppose moreover that two neighbor nodes i and j of the complex network (10)-( 2)-( 3)-( 4) satisfy equation [START_REF] Miranville | Bifurcations and synchronization in networks of unstable reaction-diffusion systems[END_REF] and the following conditions:

D i,k (B -1 i ) 2 A i = D j,k (B -1 j ) 2 A j = diag {α 1,k , . . . , α M,k } , ( 42
)
D i,k B -1 i Θ i µ i • B i = D j,k B -1 j Θ j µ j • B j , ( 43
)
γ k λ P (Ω, µ i ) > 1 + L 2 k 2 , ( 44
)
for each k ∈ {1, . . . , n}, where diag {α 1,k , . . . , α M,k } is a diagonal matrix of order M with positive coefficients, and

γ k = min {α 1,k , . . . , α M,k }.
Then the two nodes i and j of the complex network determined by the semi-linear equations [START_REF] Cantin | Complex network near-synchronization for non-identical predator-prey systems[END_REF], the boundary conditions (2)-( 3) and the initial conditions (4) synchronize.

Proof. Let k ∈ {1, . . . , m} and i, j ∈ {1, . . . , n}. We use the same notations as before for u, ũ, v and ṽ, and show the steps which differ with the proof of Theorem 2.

We consider again the energy functional defined by

E i,j,k (t) = 1 2 Ω |ũ(t, x) -ṽ(t, x)| 2 dx,
and compute its derivative:

dE i,j,k dt (t) = Ω ∂ ũ ∂t - ∂ṽ ∂t (ũ -ṽ)dx.
Using equations ( 39) and ( 40), we easily compute

∂ ũ ∂t = D i,k ∇ • (B -1 i ) 2 A i ∇ũ + f k (ũ), ∂ṽ ∂t = D j,k ∇ • (B -1 j ) 2 A j ∇ũ + f k (ṽ).
By virtue of (42) and Green's formula, we deduce that

dE i,j,k dt (t) = Ω α 1,k ∂ 2 (ũ -ṽ) ∂x 2 1 (ũ -ṽ) + • • • + α M,k ∂ 2 (ũ -ṽ) ∂x 2 M (ũ -ṽ) dx + Ω f k (ũ) -f k (ṽ) (ũ -ṽ)dx ≤ -γ k Ω |∇(ũ -ṽ)| 2 dx + D i,k Γ R ∂(ũ -ṽ) ∂ν (ũ -ṽ)ds + Ω f k (ũ) -f k (ṽ) (ũ -ṽ)dx.
We use again the Poincaré inequality (1) to write

-γ k Ω |∇(ũ -ṽ)| 2 dx ≤ -2γ k λ P (Ω, µ i )E k .
In parallel, applying again (42), we compute:

D i,k ∂ ũ ∂ν - ∂ṽ ∂ν = D i,k ν • (B -1 i ) 2 A i ∇ũ -D j,k ν • (B -1 j ) 2 A j ∇ṽ = D i,k B -1 i Θ i ν i • A i (∇u) • B i -D j,k B -1 j Θ j ν j • A j (∇v) • B j = D i,k B -1 i Θ i ∂u ∂ν i -D j,k B -1 j Θ j ∂v ∂ν j .
Next, from (3), ( 25) and (43), we deduce as in the proof of Theorem 2 that

D i,k B -1 i Θ i ∂u ∂ν i -D j,k B -1 j Θ j ∂v ∂ν j (ũ -ṽ) ≤ 0.
Finally, using [START_REF] Hecht | New development in FreeFem++[END_REF] and combining the above inequalities leads to

dE i,j,k dt (t) ≤ -2γ k λ P (Ω, µ i )E k (t) + (1 + L 2 k )E k (t).
Applying Gronwall lemma as for Theorem 2 completes the proof.

Remark 4 (Interpretation of assumption (40)). Note that assumption (40) is satisfied for a wide class of planar transformations, such as, for instance, the transformation defined by

(x 1 , x 2 ) ∈ R 2 -→ (ϕx 1 , ψx 2 ),
with ϕ = 0, ψ = 0 and ϕ = ψ, which maps an ellipse onto another ellipse of distinct eccentricity (in that case, we have M = 2). Numerical experiments shown in Figure 9 illustrate the synchronization of such non-homothetic domains.

Remark 5. Assumptions (42), ( 43) and (44) can be viewed as generalizations of assumptions (23), ( 24) and (26), respectively. Therefore, their interpretation is analogous (see Remarks 2 and 3). Furthermore, statements similar to Corollaries 1 and 2 can easily be deduced from Theorem 4, for complete graphs and two-nodes networks.

Finally, it is easily seen that Theorem 3 also extends to the case of anisotropic diffusion for cyclic graphs and rings of nearest neighbors. We obtain the following theorem. Theorem 5. Let assumptions [START_REF] Hecht | New development in FreeFem++[END_REF], ( 39) and (40) hold. Suppose that the graph G underlying the complex network (10)-( 2)-( 3)-( 4) is a cyclic graph or a ring of 2K-nearest neighbors (with K ≥ 1). Finally, assume that each pair (i, j) of nodes satisfies (42), ( 43) and (44).

Then the complex network (10)-( 2)-( 3)-( 4) globally synchronizes in Ω.

Numerical simulations

In this section, our aim is to underpin our theoretical statements by numerical experiments and to explore in a quantitative approach the validity of the assumptions made in Theorem 2 and 4. Of particular interest, as discussed in Remark 3, is the influence of the sizes of the domains Ω 1 , . . . , Ω n and of the coupling strengths µ 1 , . . . , µ n , involved in the complex networks (1)-( 2)-( 3)-( 4) or ( 10)-( 2)-( 3)-( 4), on the Poincaré constants λ P (Ω i , µ i ) appearing in the synchronization assumptions ( 26) and (44).

To that aim, we consider the predator-prey model with Leslie-Gower-type functional response, given by the following system of two reaction-diffusion equations:

∂u ∂t = d 1 ∆u + u(1 -u) - Quv u + A , ∂v ∂t = d 2 ∆v + Sv 1 - v u + C , ( 45
)
in a bounded domain Ω ⊂ R 2 ; u and v denote the densities of preys and predators, respectively. The parameters A, C, Q, S, d 1 and d 2 are positive constants. Various forms of the latter model have been studied in several papers on a single domain (see for instance [START_REF] Abid | Diffusion driven instability and Hopf bifurcation in spatial predator-prey model on a circular domain[END_REF], [START_REF] Aguirre | Three limit cycles in a Leslie-Gower predator-prey model with additive Allee effect[END_REF], [START_REF] Arancibia-Ibarra | Turing patterns in a diffusive Holling-Tanner predator-prey model with an alternative food source for the predator[END_REF], [START_REF] Banerjee | Turing instabilities and spatio-temporal chaos in ratio-dependent Holling-Tanner model[END_REF], [START_REF] Camara | Turing and Hopf patterns formation in a predator-prey model with Leslie-Gower-type functional response[END_REF], [START_REF] Liu | Dynamics of a predator-prey system with nonlinear prey-taxis[END_REF]). In particular, the emergence of Turing patterns in this reaction-diffusion system has been analyzed in [START_REF] Arancibia-Ibarra | Turing patterns in a diffusive Holling-Tanner predator-prey model with an alternative food source for the predator[END_REF].

Here, we consider complex networks with boundary couplings, of the form (1)-( 2)-( 3)-( 4) or ( 10)-( 2)-( 3)-( 4), and with point-wise couplings, of the form ( 7)-( 8)-( 9), built with n instances of the predatorprey model (45). Following [START_REF] Arancibia-Ibarra | Turing patterns in a diffusive Holling-Tanner predator-prey model with an alternative food source for the predator[END_REF], we fix A = 0.15, C = 0.28, Q = 0.575, S = 0.26, d 1 = 1, d 2 = 35 and present several numerical simulations which show how to globally synchronize Turing patterns on each node of the network. Our computations have been performed on the calculation server of the VELO research team (Laboratoire des Sciences du Numérique, Nantes Université, France), in a GNU/LINUX environment, using a finite elements splitting method and the free and open-source software FreeFem++ [START_REF] Hecht | New development in FreeFem++[END_REF]. First, we present a case of synchronization in large domains with pointwise couplings. Next, we experiment the effect of the boundary couplings on small domains and on a complete graph network. We then explore the case of non-homothetic domains with diffusion in an anisotropic habitat. Finally, we test the possibility to globally synchronize a network whose underlying graph is a cyclic graph of four nodes.

Synchronization in large domains with point-wise couplings

In our first numerical simulation, we show how to synchronize two identical domains with point-wise couplings of the form [START_REF] Camara | Turing and Hopf patterns formation in a predator-prey model with Leslie-Gower-type functional response[END_REF]. The domains Ω 1 and Ω 2 are circles of radius 200. The state variables are (u 1 , v 1 ) on domain Ω 1 and (u 2 , v 2 ) on domain Ω 2 . In Figure 4, we illustrate the possible synchronization by showing the pair (u 1 , u 2 ). Note that the dynamics of the pair (v 1 , v 2 ) is similar. The initial conditions are given for (x 1,1 , x 1,2 ) ∈ Ω 1 and (x 2,1 , x 2,2 ) ∈ Ω 2 , respectively by

u 1,0 (x 1,1 , x 1,2 ) = 0.1 + 0.1 1.1 + cos(x 1,1 + x 1,2 ) , v 1,0 (x 1,1 , x 1,2 ) = C + 0.1 1.1 + cos(x 1,1 + x 1,2 ) , u 2,0 (x 2,1 , x 2,2 ) = 0.1 + 0.1 1.1 + sin(0.8x 2,1 + 0.7x 2,2 ) , v 2,0 (x 2,1 , x 2,2 ) = C + 0.1 1.1 + sin(0.8x 2,1 + 0.7x 2,2 ) . ( 46
)
In absence of coupling, that is µ ≡ 0, we observe that the solutions converge towards distinct Turing patterns of labyrinths type (see Figure 4(a)). A careful watch on these Turing patterns convinces on their synchronization when the coupling strength is set to µ ≡ 0.3 (see Figure 4(b)). This shows that, in the case of point-wise couplings, a weak coupling strength is sufficient to synchronize the local dynamics, even in large domains.

Effect of the boundary couplings for small domains

In our second numerical simulation, we show how to synchronize two identical domains with boundary couplings of the form (3). The domains Ω 1 and Ω 2 are now circles of radius 20, thus 10 times smaller than in the first simulation. The Robin type boundaries Γ R 1 and Γ R 2 of Ω 1 and Ω 2 , respectively cover more than 3 quarters of the whole boundaries (from -5π/6 to 5π/6, depicted with black thick lines in Figure 7 (b)). The state variables are again (u 1 , v 1 ) on domain Ω 1 and (u 2 , v 2 ) on domain Ω 2 . In Figures 5-6-7, we illustrate the possible synchronization by showing again the pair (u 1 , u 2 ). The initial conditions are again given by (46) (see Figure 5). In absence of coupling, that is µ 1 = µ 2 ≡ 0, we observe1 that the solutions converge towards distinct Turing patterns of spots type (see Figure 7(a)). When the coupling strengths are set to µ 1 = µ 2 ≡ 30, thus 100 times stronger than in the first simulation, these Turing patterns are significantly modified, but they are synchronized (see Figure 7(b); note that, due to symmetry, the synchronization is viewed with a "mirror" effect). Other numerical experiments show that this synchronization does not occur in larger domains, with a weaker coupling strength, or with smaller Robin type boundaries Γ R 1 and Γ R 2 . These observations confirm the discussion on the Poincaré constants λ P (Ω, µ i ) given in Remark 3.

Four domains complete graph network

Our third simulation explores the case of a small four nodes network, whose underlying graph is complete, so that assumption [START_REF] Miranville | Bifurcations and synchronization in networks of unstable reaction-diffusion systems[END_REF] is fulfilled for each pair of nodes. We vary the geometry of the domains by considering egg-shape domains, as depicted in Figure 8. The state variables are (u 1 , v 1 ) on domain Ω 1 , (u 2 , v 2 ) on domain Ω 2 , (u 3 , v 3 ) on domain Ω 3 and (u 4 , v 4 ) on domain Ω 4 . In Figure 8, we illustrate the possible synchronization by showing the tuple (u 1 , u 2 , u 3 , u 4 ). The initial conditions are given as in (46) in the domains Ω 1 and Ω 2 , and

u 3,0 (x 3,1 , x 3,2 ) = 0.1 + 0.1 1.1 + sin(0.6x 3,1 + 0.1x 3,2 ) , v 3,0 (x 3,1 , x 3,2 ) = C + 0.1 1.1 + sin(0.35x 3,1 + 0.82x 3,2 ) , u 4,0 (x 4,1 , x 4,2 ) = 0.1 + 0.1 1.1 + sin(0.06x 4,1 + 0.91x 4,2 ) , v 4,0 (x 4,1 , x 4,2 ) = C + 0.1 1.1 + sin(0.535x 4,1 + 0.082x 4,2 ) , ( 47
)
in Ω 3 and Ω 4 . The observations are similar to the second simulation: in absence of coupling, the solutions converge to distinct Turing patterns; when the coupling strengths are set to µ 1 = µ 2 = µ 3 = µ 4 ≡ 30, these Turing patterns are modified, but synchronized, according to Theorem 2 (with the same "mirror" effect as in Figure 7(b)). Note that in the case of a four nodes network, the synchronization towards equilibrium occurs faster (t = 5800) than in the case of a two nodes network (t = 7800).

Non-homothetic domains

Our next simulation experiments the synchronization in a two-nodes network of non-homothetic domains, with the semi-linear equation [START_REF] Cantin | Complex network near-synchronization for non-identical predator-prey systems[END_REF]. The matrix A 1 is the identity matrix, whereas the matrix A 2 is given by A 2 = diag {4/9, 9/4}. The elliptic domains Ω 1 and Ω 2 are depicted in Figure 9 and the initial conditions are given by (46). The domain Ω 1 is mapped onto the domain Ω 2 under the action of the transformation given by (x 1 , x 2 ) → (2/3 x 1 , 3/2 x 2 ), which fits with matrix A 2 , so that assumptions (42) and (43) are fulfilled. Roughly speaking, synchronization in such non-homothetic domains is possible only if the diffusion rates vary in space according to the shapes of the domains. Here, the domain Ω 1 is wide along the horizontal axis, whereas the domain Ω 2 is narrow along the same axis. Hence, for reaching a synchronization state, the species u 1 , v 1 living in Ω 1 must move quickly along the horizontal axis while the species u 2 , v 2 must move slowly along the same axis. Otherwise, synchronization cannot occur. As in Subsections 3.1 and 3.2, the state variables are (u 1 , v 1 ) on domain Ω 1 and (u 2 , v 2 ) on domain Ω 2 . In Figure 9, we illustrate the possible synchronization by showing the pair (v 1 , v 2 ), since the forms of the Turing patterns are more precise than for the pair (u 1 , u 2 ). As before, the solutions converge towards Turing patterns in absence of coupling. If the couplings are activated, these Turing patterns are much modified, but they are synchronized, although the domains are not homothetic, which illustrates Theorem 4.

Synchronization in non-convex domains with a cyclic graph

In our last simulation, we test the possibility to synchronize a complex network whose underlying graph topology is not complete. Hence we consider a cyclic graph of four nodes. Moreover, we experiment a non-convex shape, so as to show that our results can be applied to any shape of spatial domain, and not only to disks, elliptic domains or other convex domains. The domains Ω 1 , Ω 2 , Ω 3 , Ω 4 , which are shown in Figure 10, admit a non-convex "C" shape; the initial conditions are again given by ( 46)-(47).

The results of the simulation are depicted in Figure 10. In absence of coupling, the solutions converge towards distinct spot patterns. If the couplings are activated, these patterns are modified, but they are synchronized (according to Theorem 3), although the graph topology underlying the complex network is not complete.

From the ecological point of view, this example shows that a global synchronization state can be reached with only a weakly densely connected network. This suggests that a control strategy for maintaining coexistence and avoiding extinction can be successfully set in place in a fragmented habitat with a relevant distribution of ecological corridors.

Conclusion and perspectives

In this paper, we have studied a new type of complex networks of reaction-diffusion systems, in which the domains are non-identical, and the couplings are defined along the boundaries of the domains. We have proved that it is possible to synchronize the local dynamics of such complex networks, provided the domains are roughly not too large, and the couplings are strong enough, independently of the shapes and sizes of the domains. We have established the following properties:

• synchronization can occur between two nodes, if the two nodes admit only common neighbors;

• synchronization can be global if the graph underlying the complex network is a complete graph, a cyclic graph, or a ring of nearest neighbors;

• compared with complex networks with point-wise couplings, our numerical experiments show that the synchronization with boundary couplings is more delicate to reach.

Our theoretical results have been applied to a complex network of competing species living in a fragmented habitat. The migrations of biological individuals from one patch of the fragmented habitat to another are supported by connections between two nodes of the corresponding network, which model ecological corridors. Our numerical simulations show that the local dynamics of the complex network can be controlled. In particular, spatial instabilities such as Turing patterns are perturbed by the boundary couplings, but they can be synchronized. In a future paper, we aim to deepen our work in several directions. First, it is natural to wonder whether the sufficient conditions considered in Theorems 2, 3 and 4 are necessary or not, so as to complete our synchronization statements with non-synchronization theorems. Next, the migrations of biological individuals between two patches of their fragmented habitat have been here assumed to occur instantaneously, for simplicity; hence, a relevant perspective would be to consider a time delay on these spatial migrations. Finally, the couplings of the complex network (1)-( 2)-( 3)-( 4) have been designed in such a manner that a patch Ω i admitting several neighbors Ω j1 , Ω j2 , . . . , Ω jp (with p = |N i |), is connected to these neighbors by corridors that share the same boundary Γ R i . Hence, it would be relevant to relax this constraint, so that a corridor connecting a pair of patches (Ω i , Ω j ) admits its own starting boundary Γ R i,j , with possibly Γ R i,j = Γ R i,j if j, j ∈ N i are such that j = j . In this way, equation (3) would be rewritten which improves the modeling of the connections in the complex network and can lead to original behaviors that have not been observed in the present work. Overall, we believe that a further study of such complex networks of reaction-diffusion systems will reveal again rich dynamics in a near future.

∂U i ∂ν i = -µ i,j (x i )(U i -U j ), (t, x i ) ∈ (0, ∞) × Γ R i,j , j ∈ N i , (48)

Figure 1 :

 1 Figure 1: Complex networks modeling a fragmented habitat with migrations of individuals between patches. (a) Identical domains and simple point-wise couplings: individuals at point x ∈ Ω migrate towards the same point x of a copy of Ω. (b)-(c) Distinct domains and realistic boundary couplings: individuals migrate from one patch to another through their boundary by crossing a corridor. The graph topology underlying the complex network can be weakly dense (b) or complete (c).

Figure 2 :

 2 Figure 2: The setting of ecological corridors can help restore biodiversity in an ecosystem whose habitat is fragmented. (a) Ecological corridor connecting both sides of a road. Source: https://the-abc-ecologicalcorridor.html. (b) Dense implementation of ecological corridors in a fragmented environment. Source: [39].

Figure 3 :

 3 Figure 3: Particular graph topologies guaranteeing partial synchronization or global synchronization. (a) In an arbitrary graph topology, two nodes admitting the same neighbors (depicted in red) synchronize. (b) Complete graph. (c) Cyclic graph. (d) Ring of 4-nearest neighbors. (e) Ring of 6-nearest neighbors.

4 :

 4 Turing instability with labyrinths patterns in a two-nodes network of predator-prey models with Leslie-Gower-type functional response (illustration of Subsection 3.1). (a) In absence of coupling, the solutions converge towards distinct Turing patterns. (b) With point-wise couplings, the Turing patterns are synchronized.

Figure 5 :

 5 Figure 5: Components u1, u2 of the initial conditions defined by (46) (illustration of Subsection 3.2).

Figure 6 :

 6 Figure 6: Turing instability with spots patterns in a two-nodes network of predator-prey models with Leslie-Gower-type functional response (illustration of Subsection 3.2). In absence of coupling, the solutions converge towards distinct Turing patterns.

Figure 7 :

 7 Figure 7:Turing instability with spots patterns in a two-nodes network of predator-prey models with Leslie-Gower-type functional response (illustration of Subsection 3.2). With boundary couplings (along the black thick frontier of the domains), the Turing patterns are modified but synchronized. We observe that synchronization occurs very quickly, before the convergence to equilibrium.

Figure 8 :

 8 Figure 8: Turing instability with spots patterns in a four-nodes network of predator-prey models with Leslie-Gower-type functional response, whose underlying graph is complete (illustration of Subsection 3.3). (a) In absence of coupling, the solutions converge towards distinct Turing patterns. (b) With boundary couplings (along the black thick frontier of the domains), the Turing patterns are modified but synchronized.

Figure 9 :

 9 Figure 9: Turing instability with spots patterns in a two-nodes network of non-homothetic domains (illustration of Subsection 3.4). (a) In absence of coupling, the solutions converge towards distinct Turing patterns. (b) With boundary couplings (along the black thick frontier of the domains), the Turing patterns are modified but synchronized.

Figure 10 :

 10 Figure 10: Turing instability with spots patterns in a cyclic four-nodes network (illustration of Subsection 3.5). (a) In absence of coupling, the solutions converge towards distinct Turing patterns. (b) With boundary couplings (along the black thick frontier of the domains), the Turing patterns are modified but they are synchronized.

Animations showing the temporal evolution can be freely visualized at https://pagesperso.ls2n.fr/ ˜cantin-g/ turingpatterns.html.

Acknowledgment

The authors wish to express their sincere gratitude to the anonymous reviewers for their valuable comments which greatly improved the presentation of the paper.