††thanks: These two authors contributed equally††thanks: These two authors
contributed equally

# Turing patterns in systems with high-order interactions

Riccardo Muolo  naXys, Namur Institute for Complex Systems, University of
Namur, Belgium  Department of Mathematics, University of Namur, Belgium
<EMAIL_ADDRESS>Department of Applied Mathematics, Mathematical
Institute, Federal University of Rio de Janeiro, Brazil  Luca Gallo  naXys,
Namur Institute for Complex Systems, University of Namur, Belgium  Department
of Physics and Astronomy & INFN, University of Catania, Italy  Vito Latora
Department of Physics and Astronomy & INFN, University of Catania, Italy
School of Mathematical Sciences, Queen Mary University of London, UK
Complexity Science Hub Vienna, Austria  Mattia Frasca  Department of
Electrical, Electronics and Computer Science Engineering, University of
Catania, Italy  Istituto di Analisi dei Sistemi ed Informatica “A. Ruberti”,
Consiglio Nazionale delle Ricerche (IASI-CNR), Roma, Italy  Timoteo Carletti
naXys, Namur Institute for Complex Systems, University of Namur, Belgium
Department of Mathematics, University of Namur, Belgium

(August 27, 2024)

###### Abstract

Turing theory of pattern formation is among the most popular theoretical means
to account for the variety of spatio-temporal structures observed in Nature
and, for this reason, finds applications in many different fields. While
Turing patterns have been thoroughly investigated on continuous support and on
networks, only a few attempts have been made towards their characterization in
systems with higher-order interactions. In this paper, we propose a way to
include group interactions in reaction-diffusion systems, and we study their
effects on the formation of Turing patterns. To achieve this goal, we rewrite
the problem originally studied by Turing in a general form that accounts for a
microscropic description of interactions of any order in the form of a
hypergraph, and we prove that the interplay between the different orders of
interaction may either enhance or repress the emergence of Turing patterns.
Our results shed light on the mechanisms of pattern-formation in systems with
many-body interactions and pave the way for further extensions of Turing
original framework.

##  I Introduction

Many natural and engineered systems exhibit collective behaviors, manifesting
themselves as spatio-temporal ordered motifs, whose emergence can be explained
utterly by considering the interactions within the system anderson . One of
the most elegant and popular theories for the emergence of self-organized
patterns is due to the British mathematician Alan Turing, who proposed, in the
context of morphogenesis, a mechanism of pattern-formation rooted on a
diffusion-driven instability, which now bears his name Turing . In a nutshell,
the Turing instability results from the combined action of two processes,
(local) reaction and (long-range) diffusion, involving an activator and an
inhibitor species Meinhardt ; Murray2001 . While each process considered
separately would drive the system to a spatially homogeneous state, there can
be conditions on the models and on the interactions such that any
heterogeneous, arbitrarily small, perturbation of the homogeneous state is
amplified and eventually returns a macroscopic patchy (non-homogeneous)
solution, i.e., a Turing pattern. The diffusive terms being the destabilizing
factors, the above mechanism is also known in the literature as a diffusion-
driven instability.

Despite its generality, the framework requires the involved species to diffuse
with quite different rates, a condition that is not often naturally realized
without introducing additional mechanisms, such as convection, electromagnetic
fields, differential adherence. Indeed experimental evidence of the existence
of Turing patterns was obtained almost half a century later in chemical
reactions castets ; DeKepper  and in Cellular Neural Networks chua_Turing1 ;
chua_Turing2 . Without adding unnecessary mechanisms and still looking for
simple models, scholars have proposed several variants to the original Turing
scheme to facilitate the emergence of patterns. For example, noise biancalani
or an upper bound to the signal propagation ZH2016 ; jop_carletti  can be
added to the picture.

In many relevant applications by their very definition, the local reactions
involve very close species, that can thus be considered spatially separated
from other groups. Starting from this observation Turing idea has been
extended to discrete systems: species occupy spatially limited zones, i.e.,
nodes, and diffuse across links connecting different zones. Initially the
theory has been developed as to include regular $1D$ and $2D$ lattices OS1971
, and then extended on complex networks NM2010 ; pastorsatorrasvespignani2010
.

The latter framework has proved particularly fruitful for Turing patterns,
especially for the introduction of asymmetric displacements Asllani2 , as in
the case of directed networks Asllani1  or non-normal ones jtb , where it was
shown that Turing mechanism is enhanced. In particular, such framework is very
general and allows to study the Turing instability on discrete topologies
which are not trivially embedded on continuous domains, as it is the case for
regular lattices. Nevertheless, certain dynamics escape the network framework,
as the interactions between the elementary units are not only pairwise, but
can nvolve several agents at the same time, we are thus in presence of many-
body interactions. Let us here stress that we do not (only) refer to processes
that introduce, a possibly small, corrections to the first order model, i.e.,
the pairwise coupling, but we overcome this setting and consider cases where
high-order interactions are the main driver. Think of social systems where the
group interaction of human being or social animals, determine the system fate
in the case, e.g., of virus spreading iacopini2019simplicial ;
barrat2022social . Or consider processes where crowding is key factor; this
can happen in the case of random walk AsllaniPRL2017 ; AsllaniPRR2020  or
chemical systems FanelliMcKane2010 ; Fanelli_2013 . Let us observe that
nonlinear terms introduced with the high-order interaction can be associated
to nonlinear diffusion processes already studied in the literature, our
results differentiate however from the latter because of the many-body
assumption. One thus needs to resort to more sophisticated mathematical
structures, such as hypergraphs or simplicial complexes, which are an
extension of networks, beyond the framework of pairwise interactions
battiston2020networks ; natphys . Dynamical systems on high-order topologies
have been recently studied carletti2020dynamical ; gambuzza2021stability , but
a general theory of Turing pattern formation on top of such structures is
still lacking. Let us observe that a similar framework has been used to study
synchronization of high-order coupled oscillators maxime2020 . Reaction-
diffusion systems have been recently studied in the frameworks of topological
signals turing_topological , which differs from the problem we are studying,
as the in the former species lie in nodes and links, while in the latter only
in nodes. The aim of this paper is to take one step further by proposing an
extension of Turing instability on high-order topologies. For this scope, we
will employ the formalism developed by Gambuzza et al. gambuzza2021stability
in the context of synchronization.

By performing a linear stability analysis involving a high-order Laplace
matrix and the Jacobian matrices of the reaction and coupling parts, we will
show that the joint action of high-order structures and nonlinear interactions
is the key feature in driving the systems unstable. Then, considering suitable
choices of couplings and high-order topologies, we have analytically
determined the conditions for the instability to occur by focusing on the role
of the diffusion coefficients at every order. As we will show, having
different orders of diffusivity allows either to enhance or to reduce the
formation of Turing patterns with respect to the case of systems with only-
pairwise interactions, by tuning the parameters of the high-order couplings.

The paper is organized as follows. In section II we will set the theoretical
ground, by extending the reaction-diffusion equations originally studied by
Turing to the case of systems coupled through high-order interactions. In
section III, we will choose a particular type of nonlinear coupling, namely
the diffusive-like cubic one, and we will analyze the emergence of patterns in
the case of systems with only-pairwise interactions. In section IV, we will
study the effects of high-order interactions on Turing instability. We will
first restrict to the particular case of the so-called natural coupling, which
allows for a fully analytical treatment, showing how this is a generalization
of some previous attempts to tackle the problem. We will then relax the
natural coupling hypothesis, and examine an intermediate case with some
specific structure, that we call regular topologies, which also allow an
analytical treatment. Lastly, we will consider the case of the most general
high-order coupling topologies, where only a numerical study is possible.
Finally, in the last section, we will discuss further lines of investigation
and possible applications.

##  II Reaction-Diffusion systems on high-order structures

Let us consider a dynamical system composed by $N$ identical units subject to
some (local) nonlinear reaction dynamics. Assume moreover that many-body
interactions, i.e., interactions in groups of more than two units, are
allowed.

Let the state of the $i$–th unit to be described by the vector
$\vec{x}_{i}(t)\in\mathbb{R}^{m}$, then, under the mentioned assumptions, the
evolution rate of the state vector of the $i$–th unit is governed by the
following equation:

$\frac{d{\vec{x}}_{i}}{dt}=\vec{f}(\vec{x_{i}})+\sum_{d=1}^{P}\sigma_{d}\sum_{j_{1}=1}^{N}\dots\sum_{j_{d}=1}^{N}A_{ij_{1}j_{2}\dots
j_{d}}^{(d)}\vec{g}^{(d)}(\vec{x}_{i},\vec{x}_{j_{1}},\vec{x}_{j_{2}},\dots,\vec{x}_{j_{d}})\,$
(1)

with $i=1,\dots,N$, where $\vec{f}:\mathbb{R}^{m}\rightarrow\mathbb{R}^{m}$
describes the local nonlinear dynamics. Let $P+1$ denote the size of the
largest group of interacting units, then
$\vec{g}^{(d)}:\mathbb{R}^{m\times(d+1)}\rightarrow\mathbb{R}^{m}$,
$d\in\\{1,\dots,P\\}$, are the nonlinear coupling functions ruling the
$(d+1)$-body interactions, encoded into the adjacency tensors $A^{(d)}$, with
$A^{(d)}_{ij_{1}\dots j_{d}}=1$ if and only if the $d+1$ units
$\\{i,j_{1},\dots,j_{d}\\}$ interact together, i.e., they are connected by a
hyperedge with the convention that repeated indexes yield a $0$ entry. We
denote by $\sigma_{d}>0$ the coupling strengths. Let us assume the nonlinear
couplings to be diffusive-like, hence for every $d$ there exists a function
$\vec{h}^{(d)}\colon\mathbb{R}^{m\times d}\to\mathbb{R}^{m}$ such that
$\vec{g}^{(d)}$ can be written as:

$\vec{g}^{(d)}(\vec{x}_{i},\vec{x}_{j_{1}},\dots,\vec{x}_{j_{d}})=\vec{h}^{(d)}(\vec{x}_{j_{1}},\dots,\vec{x}_{j_{d}})-\vec{h}^{(d)}(\vec{x}_{i},\dots,\vec{x}_{i})\,$
(2)

This definition of $\vec{g}^{(d)}$ guarantees that the coupling vanishes when
all units have the same state vector, namely
$\vec{g}^{(d)}(\vec{x}_{i},\vec{x}_{i},\dots,\vec{x}_{i})=0$. This could be
intuitively considered as a generalization of Fickean diffusion on networks,
diffusion tends to homogenize local differences, hence to vanish in the case
of equal system states. If, in addition, we assume the existence of
$\vec{x}^{*}\in\mathbb{R}^{m}$ such that $\vec{f}(\vec{x}^{*})=0$, meaning
that a solution of the $N$ isolated systems exists; then by setting
$\vec{x}_{i}=\vec{x}^{*}$, for all $i=1,\dots,N$, and by assuming Eq. (2) to
hold true, then
$(\vec{x}_{1},\dots,\vec{x}_{N})^{\top}=(\vec{x}^{*},\dots,\vec{x}^{*})^{\top}$
results to be a spatially independent solution of Eq. (1). This condition is a
prerequisite for Turing instability that, as mentioned above, is based on the
existence of a homogeneous stable equilibrium eventually destabilized by the
diffusion. Notice that the diffusion terms are determined by the topology of
the connections as well as by the parameters entering in the diffusive
couplings, so that all these factors concurrently contribute to the mechanism
of Turing pattern formation in the general model (1).

The standard framework to study Turing instability consists of two species
reaction-diffusion systems, thus even if the proposed framework is quite
general, we preferred in the following to limit our analysis to a
$2$-dimensional case. Let us thus set $m=2$ and denote the two components of
the state vector $\vec{x}_{i}$ by $(u_{i},v_{i})$. Then by defining
$\vec{f}(\vec{x})=(f_{1}(u,v),f_{2}(u,v))$ and
$\vec{h}^{(d)}(\vec{x}_{1},\dots,\vec{x}_{d})=({h}^{(d)}_{1}(u_{1},\dots,u_{d},v_{1},\dots,v_{d}),{h}^{(d)}_{2}(u_{1},\dots,u_{d},v_{1},\dots,v_{d}))$,
for all $d=1,\dots,P$, we can rewrite Eq. (1) as

$\begin{cases}\displaystyle\frac{d{u}_{i}}{dt}=f_{1}(u_{i},v_{i})+\sum_{d=1}^{P}\sigma_{d}\sum_{j_{1}=1}^{N}\dots\sum_{j_{d}=1}^{N}A_{i,j_{1},\dots,j_{d}}^{(d)}\bigg{[}h_{1}^{(d)}(u_{j_{1}},\dots,u_{j_{d}},v_{j_{1}},\dots,v_{j_{d}})\\\
\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad-
h_{1}^{(d)}(u_{i},\dots,u_{i},v_{i},\dots,v_{i})\bigg{]}\\\
\displaystyle\frac{d{v}_{i}}{dt}=f_{2}(u_{i},v_{i})+\sum_{d=1}^{P}\sigma_{d}\sum_{j_{1}=1}^{N}\dots\sum_{j_{d}=1}^{N}A_{i,j_{1},\dots,j_{d}}^{(d)}\bigg{[}h_{2}^{(d)}(u_{j_{1}},\dots,u_{j_{d}},v_{j_{1}},\dots,v_{j_{d}})\\\
\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad-
h_{2}^{(d)}(u_{i},\dots,u_{i},v_{i},\dots,v_{i})\bigg{]}\end{cases}$ (3)

where we have taken into account condition (2). To focus on the role of high-
order interactions, we can further simplify the model by assuming that the
nonlinear diffusion does not contain any cross-diffusion term, namely for all
$d\in\\{1,\dots,P\\}$ the function $h_{1}^{(d)}$ (resp. $h_{2}^{(d)}$) depends
only on $\\{u_{j_{1}},\dots,u_{j_{d}}\\}$ (resp.
$\\{v_{j_{1}},\dots,v_{j_{d}}\\}$). A throughout analysis of the general case
goes beyond the scope of this work and it could be consider in a forthcoming
study. Under this hypothesis, the examined system is for all $i=1,\dots,N$

$\begin{cases}\displaystyle\frac{d{u}_{i}}{dt}=f_{1}(u_{i},v_{i})+\sum_{d=1}^{P}\sigma_{d}\sum_{j_{1}=1}^{N}\dots\sum_{j_{d}=1}^{N}A_{i,j_{1},\dots,j_{d}}^{(d)}\left[h_{1}^{(d)}(u_{j_{1}},\dots,u_{j_{d}})-h_{1}^{(d)}(u_{i},\dots,u_{i})\right]\\\
\\\
\displaystyle\frac{d{v}_{i}}{dt}=f_{2}(u_{i},v_{i})+\sum_{d=1}^{P}\sigma_{d}\sum_{j_{1}=1}^{N}\dots\sum_{j_{d}=1}^{N}A_{i,j_{1},\dots,j_{d}}^{(d)}\left[h_{2}^{(d)}(v_{j_{1}},\dots,v_{j_{d}})-h_{2}^{(d)}(v_{i},\dots,v_{i})\right]\,\end{cases}$
(4)

In the next sections, we will study the conditions for the emergence of Turing
patterns in systems of the form (4) focusing, one-by-one, on the two novel
aspects of our model. First, we will deal with the study of Turing patterns in
the standard case when only two-body interactions are present (i.e., the model
with $P=1$), but the diffusive coupling is nonlinear (Sec. III). Then, we will
investigate the model in presence of multi-body nonlinear diffusive coupling
(Sec. IV).

##  III Turing theory with nonlinear diffusive-like coupling

Let us start with the analysis of the conditions on the fixed point for the
isolated system, that we here indicate as $(u^{*},v^{*})$. This means to only
consider local reaction and silence the interactions among the different
units. The fixed point satisfies the equations:

$f_{1}(u^{*},v^{*})=f_{2}(u^{*},v^{*})=0\,$

and should be stable, a condition that can be obtained by imposing

$\mathrm{tr}\leavevmode\nobreak\ \mathbf{J}_{0}<0\text{ and
}\det\mathbf{J}_{0}>0\,$ (5)

where
$\mathbf{J}_{0}=\left(\begin{matrix}\partial_{u}f_{1}&\partial_{v}f_{1}\\\
\partial_{u}f_{2}&\partial_{v}f_{2}\end{matrix}\right)$ is the Jacobian matrix
of the reaction function, being $\partial_{a}f_{\ell}=\frac{\partial
f_{\ell}}{\partial a}(u^{*},v^{*})$, with $\ell\in\\{1,2\\}$ and
$a\in\\{u,v\\}$. We remark that with the notation
$\partial_{a}f_{\ell}=\frac{\partial f_{\ell}}{\partial a}(u^{*},v^{*})$ we
indicate that all the derivatives in the Jacobian matrix are evaluated at the
equilibrium point $(u^{*},v^{*})$.

Before studying the effects of high-order terms on the emergence of Turing
patterns in the most general system in Eq. (4), let us first analyze the
effect of a nonlinear diffusive coupling in a system with pairwise
interactions only, i.e., when $P=1$. This corresponds thus to consider a
reactive system where nonlinear diffusion is present and interactions occur
among units that are mapped as the nodes of a complex network. Let us observe
that the problem of Turing instability with nonlinear diffusion has already
been studied, but on a continuous support for the dynamics nonlin_Tur ;
nonlin_Tur2 ; nonlin_Tur3 . The main purpose of the following sections is to
extend this analysis to the case of networked systems and to introduce the
reader to the framework of Turing theory on networks.

For sake of definitiveness, we will present our results by using a cubic
diffusion term. We adopt the same assumption even in the case of high-order
interactions, so that it will be easier to examine the effects of the latter
on the dynamical behavior of the system.

###  III.1 Turing patterns in networked systems with nonlinear diffusion

By setting $P=1$ in Eq. (4) we obtain

$\begin{cases}\displaystyle\frac{d{u}_{i}}{dt}=f_{1}(u_{i},v_{i})+\sigma_{1}\sum_{j=1}^{N}A_{ij}^{(1)}\left(h^{(1)}_{1}(u_{j})-h^{(1)}_{1}(u_{i})\right)\\\
\displaystyle\frac{d{v}_{i}}{dt}=f_{2}(u_{i},v_{i})+\sigma_{1}\sum_{j=1}^{N}A_{ij}^{(1)}\left(h^{(1)}_{2}(v_{j})-h^{(1)}_{2}(v_{i})\right)\end{cases}\forall
i=1,\dots,N\,.$ (6)

These equations can be linearized around the the equilibrium point
$(u^{*},v^{*})$, by giving:

$\begin{cases}\displaystyle\frac{d{\delta
u}_{i}}{dt}=\partial_{u}f_{1}(u^{*},v^{*})\delta
u_{i}+\partial_{v}f_{1}(u^{*},v^{*})\delta
v_{i}+\sigma_{1}\sum_{j=1}^{N}A_{ij}^{(1)}\partial_{u}h^{(1)}_{1}(u^{*})(\delta
u_{j}-\delta u_{i})\\\ \displaystyle\frac{d{\delta
v}_{i}}{dt}=\partial_{u}f_{2}(u^{*},v^{*})\delta
u_{i}+\partial_{v}f_{2}(u^{*},v^{*})\delta
v_{i}+\sigma_{1}\sum_{j=1}^{N}A_{ij}^{(1)}\partial_{v}h^{(1)}_{2}(v^{*})(\delta
v_{j}-\delta v_{i})\,\end{cases}$ (7)

where $\delta u_{i}=u_{i}-u^{*}$, $\delta v_{i}=v_{i}-v^{*}$ and
$i=1,2,\dots,N$. Let $L^{(1)}_{ij}=A_{ij}^{(1)}-k^{(1)}_{i}\delta_{ij}$ be the
$i,j$ element of the network Laplacian matrix111Let us observe that such
matrix has a non positive spectrum and it is the discrete analogous of the
continuous diffusion operator $\nabla^{2}$ once the underlying network is a
regular lattice; for this reason, the network Laplacian matrix plays a
relevant role in the context of Turing pattern formation. However, we would
like to point out that in the literature, e.g., consensus and synchronization,
it is common to find the Laplacian matrix defined as
$L^{(1)}_{ij}=k^{(1)}_{i}\delta_{ij}-A_{ij}^{(1)}$ (e.g., newmanbook ;
gambuzza2021stability ), having thus a non negative spectrum. Then in such
cases the coupling term will exhibit a negative sign., where
$k^{(1)}_{i}=\sum_{j}A_{ij}^{(1)}$ is the node degree. We can rewrite the
latter equation in compact form to emphasize the $2$-dimensional nature of the
problem as

$\frac{d}{dt}\binom{\delta u_{i}}{\delta v_{i}}=\mathbf{J}_{0}\binom{\delta
u_{i}}{\delta
v_{i}}+\sigma_{1}\sum_{j=1}^{N}L_{ij}^{(1)}\mathbf{J}_{H^{(1)}}\binom{\delta
u_{j}}{\delta v_{j}}\quad\forall i=1,\dots,N\,$ (8)

where
$\mathbf{J}_{H^{(1)}}=\left(\begin{smallmatrix}\partial_{u}h^{(1)}_{1}(u^{*})&0\\\
0&\partial_{v}h^{(1)}_{2}(v^{*})\end{smallmatrix}\right)$. Let us observe that
matrix $\mathbf{J}_{H^{(1)}}$ is diagonal, as there are no off-diagonal terms,
due to the assumption of no cross-diffusion terms. Let $\vec{\xi}=(\delta
u_{1},\delta v_{1},\dots,\delta u_{N},\delta v_{N})^{\top}$ we can eventually
rewrite the last equation in a compact form

$\frac{d\vec{\xi}}{dt}=\left(\mathbb{I}_{N}\otimes\mathbf{J}_{0}+\sigma_{1}\mathbf{L}^{(1)}\otimes\mathbf{J}_{H^{(1)}}\right)\vec{\xi}\,$
(9)

where $\mathbb{I}_{N}$ is the $N\times N$ identity matrix and $\otimes$ is the
Kronecker product.

The eigenvalues of the $2N\times 2N$ linear system (9) determine the stability
of the solution $\vec{\xi^{*}}=0$, which corresponds to
$(u_{i},v_{i})=(u^{*},v^{*})$ for all $i$. Being $\mathbf{L}^{(1)}$ a
symmetric matrix, one can find a set of orthonormal eigenvectors
$\phi_{\alpha}^{(1)}$ associated to eigenvalues $\Lambda^{(1)}_{\alpha}$,
$\alpha=1,\dots,N$. We can then make one step further by projecting (9) onto
this basis, and thus obtaining $N$ linear and decoupled $2\times 2$ systems,
each one depending on a single eigenvalue, namely

$\frac{d}{dt}\binom{\delta\hat{u}_{\alpha}}{\delta\hat{v}_{\alpha}}=\left[\mathbf{J}_{0}+\sigma_{1}\Lambda^{(1)}_{\alpha}\mathbf{J}_{H^{(1)}}\right]\binom{\delta\hat{u}_{\alpha}}{\delta\hat{v}_{\alpha}}:=\mathbf{J}^{(\alpha)}\binom{\delta\hat{u}_{\alpha}}{\delta\hat{v}_{\alpha}}\quad\forall\alpha=1,\dots,N\,$
(10)

Here $\delta\hat{u}_{\alpha}=\sum_{i}\delta u_{i}\phi^{(1)}_{\alpha,i}$ and
$\delta\hat{v}_{\alpha}=\sum_{i}\delta v_{i}\phi^{(1)}_{\alpha,i}$, are the
projection respectively of $\delta u_{i}$ and $\delta v_{i}$ on such
eigenbasis, and the last equality defines the matrix $\mathbf{J}^{(\alpha)}$.
Its eigenvalues can be obtained by solving

$\det\left(\mathbf{J}^{(\alpha)}-\lambda\mathbb{I}_{2}\right)=0\,$

that is

$\lambda^{2}-2\mathrm{tr}\leavevmode\nobreak\
\mathbf{J}^{(\alpha)}\lambda+\det\mathbf{J}^{(\alpha)}=0\,$

The root with the largest real part, considered as a function of
$\Lambda^{(1)}_{\alpha}$, is named dispersion relation,
$\lambda_{\alpha}=\max\lambda\left(\Lambda^{(1)}_{\alpha}\right)$. If there
exists $\hat{\alpha}$ such that $\lambda_{\hat{\alpha}}>0$, then the
equilibrium solution $(u_{i},v_{i})=(u^{*},v^{*})$ is unstable and a Turing
instability is observed, the system drives away from the homogeneous
equilibrium to eventually reach a new, possibly heterogeneous, solution.

Let us observe that $\mathrm{tr}\leavevmode\nobreak\
\mathbf{J}^{(\alpha)}=\mathrm{tr}\leavevmode\nobreak\
\mathbf{J}_{0}+\Lambda^{(1)}_{\alpha}\mathrm{tr}\leavevmode\nobreak\
\mathbf{J}_{H}$. Because $h_{1}$ and $h_{2}$ encode a diffusive coupling, one
can safely assume that $\partial_{u}h_{1}(u^{*})>0$ and
$\partial_{v}h_{2}(v^{*})>0$, moreover being $\Lambda^{(1)}_{\alpha}\leq 0$
and $\mathrm{tr}\mathbf{J}_{0}<0$, one can conclude that
$\mathrm{tr}\mathbf{J}^{(\alpha)}<0$ for all $\alpha$. In conclusion, a
sufficient condition to have $\lambda_{\hat{\alpha}}>0$ results to be

$\det\mathbf{J}^{(\hat{\alpha})}<0$ (11)

for some $\hat{\alpha}>1$.

###  III.2 The Brusselator model with cubic diffusion

For the sake of simplicity and without loss of generality, let us illustrate
the above results using as a reaction part the Brusselator model
PrigogineNicolis1967 ; galla , an extensively adopted dynamical system, when
it comes to study the emergence of self-organized patterns. In addition, let
us consider a cubic diffusion term. This accounts to set

$\displaystyle f_{1}(u,v)$ $\displaystyle=$ $\displaystyle 1-(b+1)u+cu^{2}v$
$\displaystyle f_{2}(u,v)$ $\displaystyle=$ $\displaystyle bu-cu^{2}v\,$ (12)

and

$\displaystyle h^{(1)}_{1}(u)$ $\displaystyle=$ $\displaystyle
D_{u}^{(1)}u^{3}$ $\displaystyle h^{(1)}_{2}(u)$ $\displaystyle=$
$\displaystyle D_{v}^{(1)}v^{3}\,$ (13)

where $b$ and $c>0$ are model parameters, and $D_{u}^{(1)}>0$ and
$D_{v}^{(1)}>0$ are generalized diffusion coefficients. Hence, the system
under investigation reads:

$\begin{cases}\displaystyle\frac{d{u}_{i}}{dt}=1-(b+1)u_{i}+cu_{i}^{2}v_{i}+\sigma_{1}D_{u}^{(1)}\sum_{j_{1}=1}^{N}A_{ij_{1}}^{(1)}(u_{j_{1}}^{3}-u_{i}^{3})\\\
\\\
\displaystyle\frac{d{v}_{i}}{dt}=bu_{i}-cu_{i}^{2}v_{i}+\sigma_{1}D_{v}^{(1)}\sum_{j_{1}=1}^{N}A_{ij_{1}}^{(1)}(v_{j_{1}}^{3}-v_{i}^{3})\,\end{cases}$
(14)

for which it is straightforward to show the existence of a unique fixed point
$(u^{*},v^{*})=(1,b/c)$. Correspondingly, the matrices $\mathbf{J}_{0}$ and
$\mathbf{J}_{H^{(1)}}$ are:

$\mathbf{J}_{0}=\begin{bmatrix}b-1&c\\\
-b&-c\end{bmatrix}\quad\text{and}\quad\mathbf{J}_{H^{(1)}}=3\begin{bmatrix}(u^{*})^{2}D_{u}^{(1)}&0\\\
0&(v^{*})^{2}D_{v}^{(1)}\end{bmatrix}\,$ (15)

Thus, the equilibrium $(u^{*},v^{*})$ is stable provided that
$\mathrm{tr}\mathbf{J}_{0}=b-c-1<0$ and $\det\mathbf{J}_{0}=c>0$. The
conditions (11) for the onset of the instability are given by:

$\begin{cases}-c^{3}D_{u}^{(1)}+(b-1){b}^{2}D_{v}^{(1)}>0\\\
4b^{2}c^{3}D_{u}^{(1)}D_{v}^{(1)}-{\left(-c^{3}D_{u}^{(1)}+(b-1)b^{2}D_{v}^{(1)}\right)^{2}}<0\,\end{cases}$
(16)

Let us observe that the above conditions do not depend on $\sigma_{1}$. As in
the case of Turing instability resulting from linear diffusion, the relevant
parameter is the ratio between the diffusive coefficients.

In Fig. 1 we compare the Turing instability regions (i.e., the set of
parameters for which Turing patterns can emerge) obtained in the case of a
linear diffusion term and in the case of a cubic one. By fixing in both cases
the diffusion coefficients $D^{(1)}_{u}=0.1$ and $D^{(1)}_{v}=1$, and the
coupling strength $\sigma_{1}=1$, we can observe the following behavior in the
plane $(b,c)$. The instability region is wider in the case of linear coupling
than in the cubic one for large values of the parameters, while the
instability region shrinks to zero more slowly in the case of cubic coupling
for small values of the parameters. This implies that there is a region for
large enough values of $b$ and $c$ where the linear diffusion allows the
emergence of Turing patterns, while in the case of cubic diffusion any initial
perturbation about the homogeneous equilibrium fades out. On the contrary, if
$b$ and $c$ are small enough, Turing patterns can emerge under the assumption
of cubic diffusion but not if diffusion is linear. These considerations are
confirmed by the dispersion relation as reported in the left panels of Fig. 1:
the reader can clearly identify the existence of eigenmodes associated to a
positive dispersion relation for the cubic diffusion case (upper panel), while
under the assumption of linear diffusion the dispersion relation is always
negative, or the opposite case (bottom panel) where is the dispersion curve
associated to nonlinear diffusion to be always negative. Let us remark that
the spectrum of the Laplacian $\Lambda^{(\alpha)}$ is discrete, however in
order to help the reader to better visualize the results the dispersion
relations are plotted also as continuous curves.

Figure 1: Turing instability in networks with nonlinear diffusion. Right
panel: Turing instability regions for the Brusselator model with
$D_{u}^{(1)}=0.1$, $D_{v}^{(1)}=1$ and $\sigma_{1}=1$. The analytical red
curves indicate the instability regions for the linear coupling, while the
blue region is for the cubic coupling, as obtained numerically. We can observe
that, for lower values of the parameters $(b,c)$ the cubic coupling allows for
pattern formation, where the linear one does not, while for greater values we
find the opposite situation. This can be visualized through the dispersion
relations reported in the left panels: on the bottom left for lower values of
the parameter (only the system subject to cubic coupling can go unstable),
while on the upper left for greater values of the parameter (only the linear
case yields patterns).

Let us conclude this section by observing that the standard Turing framework
of linear Fickean diffusion, operating both with networks or continuous
supports, requires the activator species to diffuse faster than the inhibitor
one Meinhardt , i.e., $D_{v}^{(1)}>D_{u}^{(1)}$. With the introduction of
nonlinear diffusion, such condition can be relaxed, as the system can yield
patterns even with equal diffusivities or with faster activator nonlin_Tur .
This claim follows directly from conditions (16) for the specific case of the
Brusselator model. However, it is important to note that if one considers the
entries of $J_{H^{(1)}}$ as the ”effective” diffusion coefficients, one finds
a generalization of the Turing condition $D_{v}^{(1)}>D_{u}^{(1)}$, which is

$(v^{*})^{2}D_{v}^{(1)}>(u^{*})^{2}D_{u}^{(1)}$ (17)

for the cubic case. In fact, conditions (5) and (11) imply that
$\partial_{u}f_{1}<|\partial_{v}f_{2}|$, from which we obtain

$\partial_{u}f_{1}(v^{*})^{2}D_{v}^{(1)}>\partial_{v}f_{2}(u^{*})^{2}D_{u}^{(1)}\Rightarrow\partial_{u}f_{1}(v^{*})^{2}D_{v}^{(1)}>|\partial_{v}f_{2}|(u^{*})^{2}D_{u}^{(1)}>\partial_{u}f_{1}(u^{*})^{2}D_{u}^{(1)}$

hence the new condition $(v^{*})^{2}D_{v}^{(1)}>(u^{*})^{2}D_{u}^{(1)}$.

Let us remark that the above condition is necessary but sufficient for Turing
patterns. This means that if conditions (16) are verified, then we have
condition (17), but the vice versa is not true.

In the following, we will make use of the above results and focus our
attention on the impact of high-order terms on the onset of the instability.
The nonlinear coupling is necessary to make such many-body interactions
meaningful, as pointed out in neuhauser2020multibody : in fact, if the high-
order coupling is linear, high-order interactions are nothing more than the
sum of pairwise ones. Moreover, nonlinear couplings are of particular physical
interest, as they are associated to anomalous diffusion sadighi2007exact .

##  IV High-order interactions

Let us now consider the most general case of Eq. (4) with high-order
interactions. For simplicity and without loss of generality, we will limit
ourselves to consider only first and second order interactions, i.e., $P=1$
and $P=2$. Nevertheless, the theory goes beyond the exposed examples and can
be straightforwardly generalized as to include higher $P$. The system we are
thus interested in is given by

$\begin{cases}\displaystyle\frac{d{u}_{i}}{dt}=f_{1}(u_{i},v_{i})+\sigma_{1}D_{u}^{(1)}\sum_{j_{1}=1}^{N}A_{ij_{1}}^{(1)}(h^{(1)}_{1}(u_{j_{1}})-h^{(1)}_{1}(u_{i}))\\\
\qquad\qquad\qquad\qquad\qquad\qquad+\sigma_{2}D_{u}^{(2)}\sum\limits_{j_{1}=1}^{N}\sum\limits_{j_{2}=1}^{N}A_{ij_{1}j_{2}}^{(2)}(h^{(2)}_{1}(u_{j_{1}},u_{j_{2}})-h^{(2)}_{1}(u_{i},u_{i}))\\\
\displaystyle\frac{d{v}_{i}}{dt}=f_{2}(u_{i},v_{i})+\sigma_{1}D_{v}^{(1)}\sum_{j_{1}=1}^{N}A_{ij_{1}}^{(1)}(h^{(1)}_{2}(v_{j_{1}})-h^{(1)}_{2}(v_{i}))\\\
\qquad\qquad\qquad\qquad\qquad\qquad+\sigma_{2}D_{v}^{(2)}\sum\limits_{j_{1}=1}^{N}\sum\limits_{j_{2}=1}^{N}A_{ij_{1}j_{2}}^{(2)}(h^{(2)}_{2}(v_{j_{1}},v_{j_{2}})-h^{(2)}_{2}(v_{i},v_{i}))\,\end{cases}$
(18)

where $h^{(1)}_{1}(u)$ and $h^{(1)}_{2}(v)$ encode the first order
interaction, while $h^{(2)}_{1}(u_{1},u_{2})$ and $h^{(2)}_{2}(v_{1},v_{2})$
model the second order coupling. We again assume the existence of a
homogeneous solution $(u_{i},v_{i})=(u^{*},v^{*})$ and assume it to be stable
once we silence both the pairwise coupling, $P=1$, and the high-order one,
$P>1$. To study its stability under spatially dependent perturbations, we can
follow the derivation presented in the previous section (see also
gambuzza2021stability ), that ultimately relies on the computation of a Master
Stability Function pecora1998master  in a setting involving a stationary
equilibrium and high-order coupling.

Let us introduce again the perturbation vector $\vec{\xi}=(\delta u_{1},\delta
v_{1},\dots,\delta u_{N},\delta v_{N})^{\top}$, where $\delta
u_{i}=u_{i}-u^{*}$ and $\delta v_{i}=v_{i}-v^{*}$, then one can
straightforwardly show that it evolves according to

$\frac{d\vec{\xi}}{dt}=\left(\mathbb{I}_{N}\otimes\mathbf{J}_{0}+\sigma_{1}\mathbf{L}^{(1)}\otimes\mathbf{J}_{H^{(1)}}+\sigma_{2}\mathbf{L}^{(2)}\otimes\mathbf{J}_{H^{(2)}}\right)\vec{\xi}\,$
(19)

where
$\mathbf{J}_{H^{(2)}}=\left(\begin{smallmatrix}\partial_{u_{1}}h^{(2)}_{1}(u^{*},u^{*})+\partial_{u_{2}}h^{(2)}_{1}(u^{*},u^{*})&0\\\
0&\partial_{v_{1}}h^{(2)}_{2}(v^{*},v^{*})+\partial_{v_{2}}h^{(2)}_{2}(v^{*},v^{*})\end{smallmatrix}\right)$
and $\mathbf{L}^{(2)}$ represents a generalized Laplacian matrix
gambuzza2021stability  accounting for the three-body interactions, whose
elements are given by

$\mathbf{L}^{(2)}_{ij}=\begin{cases}-\sum_{j,k=1}^{N}A^{(2)}_{ijk}&\mathrm{for}\
i=j\\\ \\\ \phantom{+}\sum_{k=1}^{N}A^{(2)}_{ijk}&\mathrm{for}\ i\neq j\\\
\end{cases}$ (20)

In the following, we will examine different cases of coupling: first, in
section IV.1, the so called natural coupling gambuzza2021stability ; carletti
, which allows an analytical treatment of the problem, but does not permit to
fully unveil the potential of high-order terms. The latter will be fully
exploited in section IV.2, where we will present two special cases of high-
order structures in which, despite the general form of the coupling functions,
a complete analytical study can be performed. Finally, in section IV.3 the
scenario where both the coupling and the topology are general will be dealt
with.

###  IV.1 Natural coupling

A largely used assumption on the coupling terms $\vec{h}^{(d)}$, is that once
they are evaluated on the synchronous manifold, they return the same value
given by $\vec{h}^{(1)}$, more precisely

$\vec{h}^{(d)}(\vec{x},\dots,\vec{x})=...=\vec{h}^{(2)}(\vec{x},\vec{x})=\vec{h}^{(1)}(\vec{x})\,$

Such condition is known in the literature as natural coupling
gambuzza2021stability . In the present framework this implies that

${h}_{1}^{(2)}(u,u)={h}_{1}^{(1)}(u)\text{ and
}{h}_{2}^{(2)}(v,v)={h}_{2}^{(1)}(v)\,$

To compare the case $P>1$ with the one presented for $P=1$, we consider
nonlinearities based on polynomials and product of variables with total power
equal to three, both in the first and in the second order terms. In
particular, we set $h^{(1)}_{1}(u)=D^{(1)}_{u}u^{3}$ and
$h^{(1)}_{2}(v)=D^{(1)}_{v}v^{3}$ for the two-body coupling, while we take
$h^{(2)}_{1}(u_{1},u_{2})=D^{(2)}_{u}u_{1}^{2}u_{2}$ and
$h^{(2)}_{2}(v_{1},v_{2})=D^{(2)}_{v}v_{1}^{2}v_{2}$, for the three-body
interaction. Note that this choice for the coupling functions satisfies the
condition of natural coupling as long as

$D_{u}^{(1)}=D_{u}^{(2)}\text{ and }D_{v}^{(1)}=D_{v}^{(2)}\,$ (21)

Given these coupling functions, we can rewrite Eq. (18) as

$\begin{cases}\displaystyle\frac{d{u}_{i}}{dt}=f_{1}(u_{i},v_{i})+\sigma_{1}D_{u}^{(1)}\sum_{j_{1}=1}^{N}A_{ij_{1}}^{(1)}(u_{j_{1}}^{3}-u_{i}^{3})+\sigma_{2}D_{u}^{(2)}\sum\limits_{j_{1}=1}^{N}\sum\limits_{j_{2}=1}^{N}A_{ij_{1}j_{2}}^{(2)}(u_{j_{1}}^{2}u_{j_{2}}-u_{i}^{3})\\\
\displaystyle\frac{d{v}_{i}}{dt}=f_{2}(u_{i},v_{i})+\sigma_{1}D_{v}^{(1)}\sum_{j_{1}=1}^{N}A_{ij_{1}}^{(1)}(v_{j_{1}}^{3}-v_{i}^{3})+\sigma_{2}D_{v}^{(2)}\sum\limits_{j_{1}=1}^{N}\sum\limits_{j_{2}=1}^{N}A_{ij_{1}j_{2}}^{(2)}(v_{j_{1}}^{2}v_{j_{2}}-v_{i}^{3})\,\end{cases}$
(22)

Again, to study the stability of the homogeneous solution
$(u_{i},v_{i})=(u^{*},v^{*})$ we linearize Eq. (22). By resorting to the
natural coupling assumption we can conclude that
$\mathbf{J}_{H^{(1)}}=\mathbf{J}_{H^{(2)}}$ and thus rewrite Eq. (19) as
follows:

$\frac{d\vec{\xi}}{dt}=\left[\mathbb{I}_{N}\otimes\mathbf{J}_{0}+\left(\sigma_{1}\mathbf{L}^{(1)}+\sigma_{2}\mathbf{L}^{(2)}\right)\otimes\mathbf{J}_{H^{(1)}}\right]\vec{\xi}\,$
(23)

where

$\mathbf{J}_{H^{(1)}}=3\begin{bmatrix}(u^{*})^{2}D_{u}^{(1)}&0\\\
0&(v^{*})^{2}D_{v}^{(1)}\end{bmatrix}\,$

and $\vec{\xi}=(\delta u_{1},\delta v_{1},\dots,\delta u_{N},\delta
v_{N})^{\top}$, with $\delta u_{i}=u_{i}-u^{*}$ and $\delta
v_{i}=v_{i}-v^{*}$. Let us observe that Eq. (23) recalls the similar result
obtained in carletti2020dynamical . However, the mathematical framework here
proposed is more general, as it allows to study a broader class of high-order
couplings (see Appendix A for more details).

Figure 2: Turing instability with nonlinear diffusion in high-order
structures. Brusselator model with $b=3$, $c=3.5$, $D_{u}^{(1,2)}=0.1$,
$D_{v}^{(1,2)}=2$, $\sigma_{1,2}=1$; the initial perturbation is $\sim
10^{-2}$. On the upper left panel, the dispersion law: the blue line is the
continuous one, while the cyan dots indicate the discrete counterpart. On the
upper right panel, Turing patterns for the $u$ species. On the bottom panel
the hypergraph on which the system has been simulated: $11$ nodes, $11$ links
and $3$ $2$-hyperedges (i.e., triangles).

The natural coupling assumption allowed us to introduce an “effective”
Laplacian that encodes the high-order structure into a weighted network, whose
weights are self-consistently determined,
$\mathbf{M}=\sigma_{1}\mathbf{L}^{(1)}+\sigma_{2}\mathbf{L}^{(2)}$. Resorting
to the eigenbasis for $\mathbf{M}$ one can make one step further and project
the $2N\times 2N$ linear system onto $N$ linear systems of size $2\times 2$
depending each one on a single eigenvalue, $\Lambda_{\alpha}$, of
$\mathbf{M}$, more precisely

$\frac{d}{dt}\binom{\delta\hat{u}_{\alpha}}{\delta\hat{v}_{\alpha}}=\left[\mathbf{J}_{0}+\Lambda_{\alpha}\mathbf{J}_{H^{(1)}}\right]\binom{\delta\hat{u}_{\alpha}}{\delta\hat{v}_{\alpha}}:=\mathbf{J}^{(\alpha)}\binom{\delta\hat{u}_{\alpha}}{\delta\hat{v}_{\alpha}}\quad\forall\alpha=1,\dots,N\,$
(24)

where again $\delta\hat{u}_{\alpha}=\sum_{i}\delta u_{i}\phi_{\alpha,i}$,
resp. $\delta\hat{v}_{\alpha}=\sum_{i}\delta v_{i}\phi_{\alpha,i}$, is the
projection of $\delta u_{i}$, resp. $\delta v_{i}$, on the eigenvector
$\phi_{\alpha}$.

Let us observe that the latter equation is formally the same as the one
obtained in the case of pairwise interactions (10) and thus the same analysis
follows, except that now the eigenvalues depend on the coupling strengths. In
Fig. 2, we show the dispersion relation and an example of Turing patterns for
the Brusselator model (Eqs. (III.2)). Let us notice that the blue curve (left
panel Fig. 2) has been computed by considering $\Lambda_{\alpha}$ to be a
continuous variable, corresponding thus to the dispersion relation for a
system defined on a continuous support with periodic boundary conditions.
Having fixed the topology of binary interactions and the high-order ones,
hence the matrices $\mathbf{L}^{(1)}$ and $\mathbf{L}^{(2)}$, we can vary the
coupling strengths $\sigma_{1}$ and $\sigma_{2}$, and thus letting
$\Lambda_{\alpha}$ “to slide” along this curve (cyan dots). Stated
differently, the Turing instability can be obtained (or repressed) by simply
changing the coupling and keeping the remaining part of the model unchanged.

###  IV.2 Regular topologies

Let us now relax the natural coupling assumption by setting different
diffusion coefficients at the different considered order, i.e.,
$D_{u}^{(1)}\neq D_{u}^{(2)}$ and $D_{v}^{(1)}\neq D_{v}^{(2)}$, but keeping
the same form for the coupling functions. Namely we take
$h^{(1)}_{1}(u)=D^{(1)}_{u}u^{3}$ and $h^{(1)}_{2}(v)=D^{(1)}_{v}v^{3}$ for
the two-body coupling, while we take
$h^{(2)}_{1}(u_{1},u_{2})=D^{(2)}_{u}u_{1}^{2}u_{2}$ and
$h^{(2)}_{2}(v_{1},v_{2})=D^{(2)}_{v}v_{1}^{2}v_{2}$ for the three-body
coupling. It would also be possible to adopt different coupling functional
forms for each order, but, as it will be clear in the following, it is more
interesting to focus on the diffusion coefficients, the latter having a key
role in the Turing mechanism of pattern formation. Nevertheless, the following
analysis could be easily extended to a setting in which coupling functions
have different form, as long as they remain diffusive-like.

The starting point for the linear stability analysis is again Eq. (19);
however, because $\mathbf{J}_{H^{(1)}}\neq\mathbf{J}_{H^{(2)}}$ one cannot
simplify the latter equation to obtain the analogous of (23). In the same
spirit of simplifying the stability equation, we can observe that in general
the two Laplacians do not commute and thus they cannot be simultaneously
diagonalized. In conclusion one cannot determine a single equation, depending
on the spectrum of the involved Laplacian matrices, and proceed with the
analysis as done in the case of the Master Stability Function. The numerical
computation of the eigenvalues of the $2N\times 2N$ linear system would not
allow a clear understanding of the role of the involved parameters.

There are however some high-order structures allowing for a complete
analytical description: this is the case whenever the high-order Laplacian
matrix is a multiple of the network Laplacian one. We will call topologies
with the above property regular topologies. Tetrahedra or icosaedra are
examples of regular topologies as well as the triangular lattice with periodic
boundary conditions, i.e., a $2$-torus paved with triangles, that we hereby
analyze in details. Let us observe that we assume the nodes forming a triangle
(of first order interactions) to be also part of a three-body interaction; for
this reason we will call it triangular $2$-lattice. In this case, each node
interacts with its six neighbors through six two-body interactions and six
three-body interactions

Let us consider a triangular $2$-lattice of $N$ nodes, where each node $i$ has
$6$ incident links. Its first order Laplacian will then be

$\textbf{L}_{ij}^{(1)}=\begin{cases}-6\leavevmode\nobreak\
\leavevmode\nobreak\ \mbox{ if }\leavevmode\nobreak\ i=j\\\
A_{ij}^{(1)}\leavevmode\nobreak\ \leavevmode\nobreak\ \mbox{
otherwise}\end{cases}$

If now each node will also have $6$ incident triangles, we have that

$\sum_{j,k=1}^{N}A^{(2)}_{ijk}=2\cdot 6=12\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \mbox{ and }\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
\sum_{k=1}^{N}A^{(2)}_{ijk}=2A^{(1)}_{ij}$

From the definition of the second order Laplacian, Eq. (20), we obtain the
relation

$\mathbf{L}^{(2)}=2\mathbf{L}^{(1)}$ (25)

hence the two Laplacians can be both diagonalized through the eigenvectors of
$\mathbf{L}^{(1)}$.
Eq. (19) now takes the form

$\frac{d\vec{\xi}}{dt}=\left(\mathbb{I}_{N}\otimes\mathbf{J}_{0}+\mathbf{L}^{(1)}\otimes(\sigma_{1}\mathbf{J}_{H^{(1)}}+2\sigma_{2}\mathbf{J}_{H^{(2)}})\right)\vec{\xi}\,$
(26)

We can then proceed in projecting on the eigenvectors of $\mathbf{L}^{(1)}$,
by obtaining

$\frac{d}{dt}\binom{\delta\hat{u}_{\alpha}}{\delta\hat{v}_{\alpha}}=\left[\mathbf{J}_{0}+\Lambda_{\alpha}(\sigma_{1}\mathbf{J}_{H^{(1)}}+2\sigma_{2}\mathbf{J}_{H^{(2)}})\right]\binom{\delta\hat{u}_{\alpha}}{\delta\hat{v}_{\alpha}}\quad\forall\alpha=1,\dots,N\,$
(27)

which is analogous to Eq. (24) due to the properties of the topology. The
Laplacian eigenvalues $\Lambda_{\alpha}$ can be effectively replaced by a
continuous parameter, which is the continuous spectrum of the corresponding
diffusion operator, allowing us to have an analytical expression of the
dispersion law. The above result can be straightforwardly generalized to all
regular topologies..

Let us point out that the assumption of $\mathbf{L}^{(2)}$ to be proportional
to $\mathbf{L}^{(1)}$ allows to emphasize the effects of high-order
interactions which would not be possible to fully appreciate under the natural
coupling assumption. In fact, one can deal with sets of parameters which do
not allow the formation of patterns by considering the sole pairwise
interactions, i.e., condition (17) would not be satisfied,
$(v^{*})^{2}D_{v}^{(1)}<(u^{*})^{2}D_{u}^{(1)}$, and thus patterns could not
develop, but the appropriate values of the second order couplings would make
possible the Turing instability. This is the case shown in Fig. 3, where we
can observe the dispersion law and an example of a pattern obtained with the
Brusselator model on a triangular $2$-lattice of $16$ nodes; we notice that,
while the pairwise case (purple dashed line) never allows instability, the
high-order one yields Turing patterns. Let us remark that the dispersion law
is a function of the spectrum of the Laplacian $\textbf{L}^{(1)}$.

Figure 3: Triangular $2$-lattice of $16$ nodes with periodic boundary
conditions; Brusselator model with $b=5.5$, $c=7$, $D_{u}^{(1)}=1$,
$D_{v}^{(1)}=0.5$, $D_{u}^{(2)}=0.1$, $D_{v}^{(2)}=1$, $\sigma_{1}=0.01$ and
$\sigma_{2}=1$; the initial perturbation is $\sim 10^{-2}$. In the left panel,
the dispersion law for the high-order case (blue line and magenta dots)
compared with the case where only pairwise interactions are present (purple
dashed line). In the right panel, an example of a Turing pattern for the $u$
species.

Let us stress that it is also possible to find the opposite case, in which
pairwise interactions would normally give rise to Turing patterns, but the
presence of high-order ones annihilates them and stabilizes the system, as
exemplified in Fig. 4.

The interplay between nonlinear diffusion and regular topologies could be
better highlighted by rewriting (27). Under our working assumption we have

$\mathbf{J}_{H^{(1)}}=3\begin{bmatrix}D_{u}^{(1)}&0\\\
0&\Big{(}\frac{b}{c}\Big{)}^{2}D_{v}^{(1)}\end{bmatrix}\quad\text{ and
}\quad\mathbf{J}_{H^{(2)}}=3\begin{bmatrix}D_{u}^{(2)}&0\\\
0&\Big{(}\frac{b}{c}\Big{)}^{2}D_{v}^{(2)}\end{bmatrix}\,.$ (28)

We can hence define the effective diffusion coefficients

$\mathcal{D}_{u}^{eff}=\sigma_{1}D_{u}^{(1)}+\sigma_{2}D_{u}^{(2)}\quad\text{
and }\quad\mathcal{D}_{v}^{eff}=\sigma_{1}D_{v}^{(1)}+\sigma_{2}D_{v}^{(2)}$
(29)

so that we can cast (26) in the following way

$\frac{d\vec{\xi}}{dt}=\left(\mathbb{I}_{N}\otimes\mathbf{J}_{0}+\mathbf{L}^{(1)}\otimes\mathbf{J}_{\mathcal{H}}^{eff}\right)\vec{\xi}\,,$
(30)

where

$\mathbf{J}_{\mathcal{H}}^{eff}=3\begin{bmatrix}\mathcal{D}_{u}^{eff}&0\\\
0&\Big{(}\frac{b}{c}\Big{)}^{2}\mathcal{D}_{v}^{eff}\end{bmatrix}\,.$ (31)

By projecting on the eigenvectors of $\mathbf{L}^{(1)}$, we obtain a new form
of Eq. (27), namely

$\frac{d}{dt}\binom{\delta\hat{u}_{\alpha}}{\delta\hat{v}_{\alpha}}=\left[\mathbf{J}_{0}+\Lambda_{\alpha}\mathbf{J}_{\mathcal{H}}^{eff}\right]\binom{\delta\hat{u}_{\alpha}}{\delta\hat{v}_{\alpha}}\quad\forall\alpha=1,\dots,N\,.$
(32)

Following the steps presented above to determine the onset of Turing
instability, we can now obtain

$\begin{cases}-c^{3}\mathcal{D}_{u}^{eff}+(b-1){b}^{2}\mathcal{D}_{v}^{eff}>0\\\
4b^{2}c^{3}\mathcal{D}_{u}^{eff}\mathcal{D}_{v}^{eff}-{\left(-c^{3}\mathcal{D}_{u}^{eff}+(b-1)b^{2}\mathcal{D}_{v}^{eff}\right)^{2}}<0\,.\end{cases}$
(33)

From the latter relations it emerges that Turing patterns are the result of
the interplay not only between the model parameters and the diffusion
coefficients, but also with the strength of the interactions at every order.

Figure 4: Triangular $2$-lattice of $16$ nodes with periodic boundary
conditions; Brusselator model with $b=5.5$, $c=7$, $D_{u}^{(1)}=0.1$,
$D_{v}^{(1)}=1.5$, $D_{u}^{(2)}=1$, $D_{v}^{(2)}=0.5$, $\sigma_{1}=0.7$ and
$\sigma_{2}=0.2$; the initial perturbation is $\sim 10^{-2}$. In the left
panel, the dispersion law for the high-order case (blue line and magenta dots)
compared with the case where only pairwise interactions are present (purple
dashed line). As the high-order dispersion law is always negative, there is no
emergence of Turing patterns, as also shown in the right panel for the $u$
species.

In the next section, we will study the general case where, in absence of
specific assumptions on the coupling or on the structure, an analytical study
cannot be performed, and we thus have to resort to numerical simulations. But
before that, let us examine a particular case of regular topology, were the
Laplacians of all orders can be diagonalized simultaneously: the all-to-all
coupling, where every possible $d$-body interaction is active. When dealing
with such coupling, the analysis is not restricted to some orders of
interactions but can be extended to all of them. Nonetheless, without further
ado, let us again restrict the analysis to first- and second-order
interactions. For $N$ interconnected systems we have
$\mathbf{L}^{(2)}=(N-2)\mathbf{L}^{(1)}$ (see gambuzza2021stability  for a
detailed analysis); let us observe that this is again a case where
$\mathbf{L}^{(2)}$ is proportional to $\mathbf{L}^{(1))}$, however here the
proportionality constant depends on the system size. Based on the above, Eq.
(19) can be rewritten as

$\frac{d}{dt}\vec{\xi}=\left(\mathbb{I}_{N}\otimes\mathbf{J}_{0}+\sigma_{1}\mathbf{L}^{(1)}\otimes\left(\mathbf{J}_{H^{(1)}}+\frac{\sigma_{2}}{\sigma_{1}}(N-2)\mathbf{J}_{H^{(2)}}\right)\right)\vec{\xi}\,$
(34)

Projecting again on the eigenvectors of $\mathbf{L}^{(1)}$ we obtain

$\frac{d}{dt}\binom{\hat{\delta u}_{\alpha}}{\hat{\delta
v}_{\alpha}}=\left[\mathbf{J}_{0}+\sigma_{1}\Lambda_{\alpha}\left(\mathbf{J}_{H^{(1)}}+\frac{\sigma_{2}}{\sigma_{1}}(N-2)\mathbf{J}_{H^{(2)}}\right)\right]\binom{\hat{\delta
u}_{\alpha}}{\hat{\delta v}_{\alpha}}\quad\forall\alpha=1,\dots,N\,$

Let us observe that $\Lambda_{\alpha}\in\\{0,N\\}$, hence the dispersion law
depends on the number of nodes. The above equation indicates that it is
possible to decompose the instability modes.

For a given value of $N$, having fixed the ratio $\sigma_{2}/\sigma_{1}$, we
can vary the value of $\sigma_{1}$, thus generating a continuous curve.
Observe that such a curve allows us to check for the onset of Turing
instability on the “all-to-all” configuration. Indeed, when the curve is
negative Turing instability cannot occur for any value of $\sigma_{1}$, while
Turing patterns can emerge in the opposite case. Hence, the continuous curve,
as function of a parameter $\gamma=\sigma_{1}\Lambda_{\alpha}$, can be
subsumed in the framework of the Master Stability Function approach
pecora1998master .

Figure 5: All-to-all 2-hypergraph of $5$ nodes; Brusselator model with
$b=5.5$, $c=7$, $D_{u}^{(1)}=1$, $D_{v}^{(1)}=0.1$, $D_{u}^{(2)}=0.07$,
$D_{v}^{(2)}=1$, $\sigma_{1}=0.01$ and $\sigma_{2}=1$; the initial
perturbation is $\sim 10^{-2}$. On the left panel, the discrete dispersion
law; note that the continuous curve is a fictitious dispersion law, but it is
a Master Stability Function, as explained in the text; for this reason the
continuous curve can be thought as a function of a parameter $\gamma$. On the
right panel, an example of a Turing pattern for the $u$ species.

To provide a numerical example, let us consider an ensemble of $N$ Brusselator
systems interacting in the all-to-all configuration, with
$h^{(1)}_{1}=D^{(1)}_{u}u^{3}$, $h^{(1)}_{2}=D^{(1)}_{v}v^{3}$,
$h^{(2)}_{1}(u_{1},u_{2})=D^{(2)}_{u}u_{1}^{2}u_{2}$ and
$h^{(2)}_{2}(v_{1},v_{2})=D^{(2)}_{v}v_{1}^{2}v_{2}$, and with the values of
the diffusion coefficients set such that the coupling functions do not
satisfies the natural coupling condition. In Fig. 5, we report a case
analogous to that displayed in Fig. 3, namely patterns emerge due to high-
order interaction. Also in this case, with the appropriate coupling functions
we can obtain the opposite situation, i.e., high-order interactions annihilate
patterns, which would otherwise emerge with the sole pairwise ones.

###  IV.3 General topologies

Let us now focus on the most general case, in which arbitrary high-order
topologies are considered and the coupling is not restricted to the natural
coupling form. Following the same reasoning of the previous section,we use
different diffusion coefficients, while the coupling functions remain of the
same form.

In such setting, we cannot have a semi-analytical form of the dispersion law,
because neither the Laplacians nor the matrices $J_{H^{(d)}}$ can be
diagonalized simultaneously. Hence, we will resort to numerical simulations to
determine the onset of the Turing instability.

Figure 6: Left panel: Brusselator model with $D_{u}^{(1)}=1$,
$D_{v}^{(1)}=0.5$, $D_{u}^{(2)}=0.1$, $D_{v}^{(2)}=2$, $\sigma_{1}=0.1$ and
$\sigma_{2}=1$. The setting does not allow for Turing patterns when only
pairwise interaction are considered, hence to obtain the instability region we
need to include also high-order (three-body) interactions (the region of
parameters where patterns are obtained is shown in orange). Numerical
simulations of the system with $b=5.5$ and $c=7$ without (resp. with) high-
order terms show that Turing patterns are not obtained (resp. are obtained) as
confirmed by the total pattern amplitude NM2010  in blue (resp. orange) on the
lower (resp. upper) inset. Right panel: Brusselator model with
$D_{u}^{(1)}=0.1$, $D_{v}^{(1)}=0.5$, $D_{u}^{(2)}=0.01$, $D_{u}^{(2)}=1$,
$\sigma_{1}=0.2$ and $\sigma_{2}=1$. The instability region when only pairwise
interactions are considered is depicted in blue, while in orange the one
obtained when also high-order ones are active; note that the blue region is a
subset of the orange one. For both panels, the initial perturbation is $\sim
10^{-2}$ and the hypergraph is that of in Fig. 2.

Based on the above discussion, we thus performed a dedicated numerical
analysis of this scenario that allows us to fully exploit the presence of
high-order terms. The results are reported in Fig. 6 where we show the
parameters region associated to a Turing instability in the case of the
Brusselator model in Eq. (22) for two different sets of values of the the
weights $\sigma_{i}$. On the left panel, the first order diffusion
coefficients (i.e., pairwise) do not allow Turing instability (since
$(v^{*})^{2}D_{v}^{(1)}>(u^{*})^{2}D_{u}^{(1)}$), on the other hand the
presence of high-order diffusion with a suitable choice of the weights
$\sigma_{i}$ allows the formation of patterns. On the right panel, instead,
Turing instability arises even in the pairwise setting, but many-body
interactions are still beneficial as they yield a larger region of parameters
for which patterns can be obtained.

Let us conclude by remarking that, in analogy with the setting studied in the
previous section (i.e., regular topologies), it is also possible to find
couplings such that high-order interaction hamper the formation of Turing
patterns, which would normally arise when the interactions are limited to be
pairwise.

##  V Discussion

In this work we have formulated a general theory to study the emergence of
Turing patterns for dynamical systems where multi-body interactions are taken
into account, modeled by using high-order interactions. Our framework goes
beyond the one recently proposed carletti2020dynamical  that, as shown in the
Appendix, can be recovered as a special case of the theory hereby developed.
Our framework is inspired by the work of gambuzza2021stability  dealing with
synchronization in high-order structures and allows to obtain Turing patterns
in settings where it would be otherwise impossible, e.g., by restricting to
pairwise interactions. At the same time it also permits to suppress the
instability which would occur with only pairwise interactions, depending on
the desired applications. We have shown that including diffusion terms from
high-order interactions may either widen the region of parameters where
patterns occur, or, on the contrary, reduce its extension. This can be
achieved acting solely on the diffusion coefficients but also using different
coupling functions. The detailed analytical study we performed provided us a
clear understanding of the role of the model parameters as well as of the
high-order topology in the pattern emergence. This further flexibility in the
choice of the parameters can shed light on the fine tuning problem HG2021  and
benefit the field of optimal control of patterns control_hata ; control_frasca
; control_holme . In fact, one could choose the appropriate coupling functions
and/or diffusion coefficients in order to obtain the desired inhomogeneous
state.

A natural extension of our framework would be to study Turing (wave)
instability emerging from a limit cycle: in fact, it has been shown that a
similar instability mechanism can occur also when the homogeneous stable state
is a limit cycle chall ; entropy . In many applications where it is important
to control the emergence of synchronization, such as neuroscience, our theory
could provide to be extremely useful. Moreover, exploring the onset of Turing
patterns in the newly developed framework of $M$-directed hypergraphs
gallo2022synchronization  could provide another interesting direction where to
focus the attention in future works.

#### Acknowledgements

The authors would like to thank Lorenzo Giambagli and Lucia Valentina Gambuzza
for useful discussions and feedback. R.M. is supported by a FRIA-FNRS PhD
fellowship, Grant FC 33443, funded by the Walloon region. R.M. and L.G.
acknowledge funds from the Erasmus+ program.

##  Appendix A Formalism mapping

In this appendix, we show how the formalism introduced in
carletti2020dynamical  to study the formation of Turing patterns on
hypergraphs can be recovered from the mathematical framework developed in
gambuzza2021stability  for the analysis of synchronization of chaotic
oscillators in simplicial complexes. Note that the latter also applies to the
case of hypergraphs, as discussed in gallo2022synchronization .

We start from the dynamical system (1). First, we assume the coupling
functions $\vec{g}^{(d)}$ to be of the form

$\begin{array}[]{lll}\vec{g}^{(d)}(\vec{x}_{i},\vec{x}_{j_{1}},\vec{x}_{j_{2}},\dots,\vec{x}_{j_{d}})&=&\varphi(d+1)[\gamma(\vec{x}_{j_{1}})-\gamma(\vec{x}_{i})+\dots+\gamma(\vec{x}_{j_{d}})-\gamma(\vec{x}_{i})]\\\
&=&\varphi(d+1)\left[\sum\limits_{n=1}^{d}\gamma(\vec{x}_{j_{n}})-d\gamma(\vec{x}_{i})\right]\end{array}$
(A1)

where $\varphi:\mathbb{R}\longrightarrow\mathbb{R}$ is a generic function of
the number of nodes involved in the high-order interaction, i.e., the size of
the hyperedge or of the simplex, while
$\gamma:\mathbb{R}^{m}\longrightarrow\mathbb{R}^{m}$ is a generic function
encoding the contribution of the node state vectors to the coupling. Note
that, given the form in Eq. (A1), all the coupling functions are non-invasive.
Eq. (1) can be hence rewritten as

$\begin{array}[]{lll}\dot{\vec{x}}_{i}&=&\vec{f}(\vec{x_{i}})+\sigma_{1}\sum\limits_{j_{1}=1}^{N}A_{ij_{1}}^{(1)}\varphi(2)[\gamma(\vec{x}_{j_{1}})-\gamma(\vec{x}_{i})]+\\\\[10.0pt]
&&+\sigma_{2}\sum\limits_{j_{1}=1}^{N}\sum\limits_{j_{2}=1}^{N}A_{ij_{1}j_{2}}^{(2)}\varphi(3)[\gamma(\vec{x}_{j_{1}})+\gamma(\vec{x}_{j_{2}})-2\gamma(\vec{x}_{i})]+\dots\\\\[10.0pt]
&&+\sigma_{P}\sum\limits_{j_{1}=1}^{N}\dots\sum\limits_{j_{P}=1}^{N}A_{ij_{1}j_{2}\dots
j_{P}}^{(P)}\varphi(P+1)\left[\sum\limits_{n=1}^{P}\gamma(\vec{x}_{j_{n}})-P\gamma(\vec{x}_{i})\right]\end{array}$
(A2)

Noticing that each adjacency tensor $A^{(d)}$ is symmetric with respect to its
$d+1$ indices, i.e., $A_{ij_{1}j_{2}\dots j_{d}}^{(d)}=A_{\pi(ij_{1}j_{2}\dots
j_{d})}^{(d)}$, where $\pi$ is a generic permutation of the indices, we can
simplify the coupling terms, thus writing Eq. (A2) as

$\begin{array}[]{lll}\dot{\vec{x}}_{i}&=&\vec{f}(\vec{x_{i}})+\sigma_{1}\sum\limits_{j_{1}=1}^{N}A_{ij_{1}}^{(1)}\varphi(2)[\gamma(\vec{x}_{j_{1}})-\gamma(\vec{x}_{i})]+\\\\[10.0pt]
&&+2\sigma_{2}\sum\limits_{j_{1}=1}^{N}\varphi(3)[\gamma(\vec{x}_{j_{1}})-\gamma(\vec{x}_{i})]\sum\limits_{j_{1}=1}^{N}A_{ij_{1}j_{2}}^{(2)}+\dots\\\\[10.0pt]
&&+P\sigma_{P}\sum\limits_{j_{1}=1}^{N}\varphi(P+1)[\gamma(\vec{x}_{j_{1}})-\gamma(\vec{x}_{i})]\sum\limits_{j_{2}=1}^{N}\dots\sum\limits_{j_{P}=1}^{N}A_{ij_{1}j_{2}\dots
j_{P}}^{(P)}\end{array}$ (A3)

By recalling the definition gambuzza2021stability  of the generalized
$d$-degree, $k^{(d)}_{ij}$, of a link $(i,j)$, which represents the number of
$d$-order structures the link $(i,j)$ is part of

$(d-1)!k^{(d)}_{ij}=\sum\limits_{l_{2}=1}^{N}\dots\sum\limits_{l_{d}=1}^{N}A_{il_{1}l_{2}\dots
l_{d}}^{(d)}$ (A4)

and by noticing that $k^{(1)}_{ij}=A_{ij}^{(1)}$, we can write the equations
governing the dynamics of the system as

$\begin{array}[]{lll}\dot{\vec{x}}_{i}&=&\vec{f}(\vec{x_{i}})+\sigma_{1}\sum\limits_{j=1}^{N}k_{ij}^{(1)}\varphi(2)\Delta_{j}+2\sigma_{2}\sum\limits_{j=1}^{N}k_{ij}^{(2)}\varphi(3)\Delta_{j}+\dots\\\
&&+P(P-1)!\sigma_{P}\sum\limits_{j=1}^{N}k_{ij}^{(P)}\varphi(P+1)\Delta_{j}\end{array}$
(A5)

where we have relabeled index $j_{1}$ as $j$, and where we have defined
$\Delta_{j}=\gamma(\vec{x}_{j})-\gamma(\vec{x}_{i})$.

By considering the high-order incidence matrix carletti2020dynamical ,
$e_{i\alpha}$, defined as

$\displaystyle e_{i\alpha}$ $\displaystyle=\begin{cases}1&v_{i}\in
E_{\alpha}\\\ 0&\mathrm{otherwise}\end{cases}$ (A6)

where $E_{\alpha}$ can represent either a hyperedge or a simplex, we can now
write Eq. (A5) as

$\begin{array}[]{lll}\dot{\vec{x}}_{i}&=&\vec{f}(\vec{x_{i}})+\sigma_{1}\sum\limits_{j=1}^{N}\sum\limits_{\alpha:|E_{\alpha}|=2}e_{i\alpha}e_{j\alpha}\varphi(2)\Delta_{j}+2\sigma_{2}\sum\limits_{j=1}^{N}\sum\limits_{\alpha:|E_{\alpha}|=3}e_{i\alpha}e_{j\alpha}\varphi(3)\Delta_{j}+\dots\\\
&&+P!\sigma_{P}\sum\limits_{j=1}^{N}\sum\limits_{\alpha:|E_{\alpha}|=P+1}e_{i\alpha}e_{j\alpha}\varphi(P+1)\Delta_{j}\end{array}$
(A7)

In fact, as for each order $d$ of the high-order interactions, the factor
$e_{i\alpha}e_{j\alpha}$ is equal to one when both nodes $i$ and $j$ belong to
$E_{\alpha}$, by summing over all $j\neq i$, each summation has exactly
$k_{ij}^{(d)}$ non-zero terms.

Lastly, by assuming the $d$-th coupling strength to be
$\sigma_{d}=\varepsilon/d!$, we finally obtain

$\dot{\vec{x}}_{i}=\vec{f}(\vec{x_{i}})+\varepsilon\sum\limits_{j=1}^{N}\sum\limits_{\alpha}e_{i\alpha}e_{j\alpha}\varphi(|E_{\alpha}|)[\gamma(\vec{x}_{j})-\gamma(\vec{x}_{i})]$
(A8)

which recovers the dynamical system analyzed in carletti2020dynamical  as we
set $\varphi(|E_{\alpha}|)=|E_{\alpha}|-1$. Hence, the mathematical framework
here considered permits to extend the analysis of Turing pattern formation in
hypergraphs, as it allows to consider more general coupling functions.

## References

  * (1) Anderson PW. 1972 More is different. Broken symmetries and the nature of hierarchical structure of science. Physical Review Letters 177.
  * (2) Turing AM. 1952 The chemical basis of morphogenesis. Phil. Trans. R. Soc. Lond. B 237, 37.
  * (3) Gierer A, Meinhardt H. 1972 A theory of biological pattern formation. Kybernetik 12, 30–39.
  * (4) Murray JD. 2001 Mathematical biology II: Spatial models and biomedical applications. Springer-Verlag.
  * (5) Castets V, Dulos E, Boissonade J, De Kepper P. 1990 Experimental evidence of a sustained standing Turing-type nonequilibrium chemical pattern. Phys. Rev. Lett. 64, 2953.
  * (6) De Kepper P, Castets V, Dulos E, Boissonade J. 1991 Turing-type chemical patterns in the chlorite–iodide–malonic acid reaction. Physica D 49.
  * (7) Goras L, Chua L. 1995 Turing patterns in CNNs. II. Equations and behaviors. IEEE Transactions on Circuits and Systems I: Fundamental Theory and Applications 42, 612–626.
  * (8) Buscarino A, Corradino L, Fortuna L, Frasca M, Chua L. 2016 Turing patterns in memristive cellular nonlinear networks. IEEE Transactions on Circuits and Systems I: Regular Papers 63, 1222–1230.
  * (9) Biancalani T, Di Patti F, Fanelli D. 2010 Stochastic Turing patterns in the Brusselator model. Phys. Rev. E 4, 046215.
  * (10) Zemskov EP, Horsthemke W. 2016 Diffusive instabilities in hyperbolic reaction-diffusion equations. Phys. Rev. E 93, 032211.
  * (11) Carletti T, Muolo R. 2021 Finite propagation enhances Turing patterns in reaction–diffusion networked systems. Journal of Physics: Complexity 2, 045004.
  * (12) Othmer HG, Scriven LE. 1971 Instability and dynamic pattern in cellular networks. J. Theor. Biol. 32, 507.
  * (13) Nakao H, Mikhailov AS. 2010 Turing patterns in network-organized activator-inhibitor systems. Nature Physics 6, 544.
  * (14) Pastor-Satorras R, Vespignani A. 2010 Patterns of complexity. Nature Physics 6, 480.
  * (15) Asllani M, Carletti T, Fanelli D, Maini PK. 2020 A universal route to pattern formation in multicellular systems. The European Physics Journal B 93.
  * (16) Asllani M, Challenger JD, Pavone FS, Sacconi L, Fanelli D. 2014 The theory of pattern formation on directed networks. Nature Communication 5.
  * (17) Muolo R, Asllani M, Fanelli D, Maini PK, Carletti T. 2019 Patterns of non-normality in networked systems. Journal of Theoretical Biology 480, 81.
  * (18) Iacopini I, Petri G, Barrat A, Latora V. 2019 Simplicial models of social contagion. Nature Communications 10, 2485.
  * (19) Barrat A, Ferraz de Arruda G, Iacopini I, Moreno Y. 2022 Social contagion on higher-order structures. In Higher-Order Systems pp. 329–346. Springer.
  * (20) Asllani M, Carletti T, Di Patti F, Fanelli D, Piazza F. 2018 Hopping in the Crowd to Unveil Network Topology. Phys. Rev. Lett. 120.
  * (21) Carletti T, Asllani M, Fanelli D, Latora V. 2020 Nonlinear walkers and efficient exploration of congested networks. Phys. Rev. Res. 2.
  * (22) Fanelli D, McKane A. 2010 Diffusion in a crowded environment. Phys. Rev. E 82, 021113.
  * (23) Fanelli D, McKane AJ, Pompili G, Tiribilli B, Vassalli M, Biancalani T. 2013 Diffusion of two molecular species in a crowded environment: theory and experiments. Physical Biology 10, 045008.
  * (24) Battiston F, Cencetti G, Iacopini I, Latora V, Lucas M, Patania A, Young JG, Petri G. 2020 Networks beyond pairwise interactions: structure and dynamics. Physics Reports.
  * (25) Battiston F, et al.. 2021 The physics of higher-order interactions in complex systems. Nature Physics 17, 1093–1098.
  * (26) Carletti T, Fanelli D, Nicoletti S. 2020 Dynamical systems on hypergraphs. Journal of Physics: Complexity 1, 035006.
  * (27) Gambuzza LV, Di Patti F, Gallo L, Lepri S, Romance M, Criado R, Frasca M, Latora V, Boccaletti S. 2021 Stability of synchronization in simplicial complexes. Nature communications 12, 1–13.
  * (28) Lucas M, Cencetti G, Battiston F. 2020 A multi-order Laplacian framework for the stability of higher-order synchronization. Physical Review Research 2, 033410.
  * (29) Giambagli L, Calmon L, Muolo R, Carletti T, Bianconi G. 2022 Diffusion-driven instability of topological signals coupled by the Dirac operator. arXiv preprint https://arxiv.org/abs/2207.07787.
  * (30) Kumar N, Horsthemke W. 2010 Turing bifurcation in a reaction-diffusion system with density-dependent dispersal. Physica A 389, 1812–1818.
  * (31) Gambino F, Lombardo M, Sammartino M, Sciacca V. 2013 Turing pattern formation in the Brusselator system with nonlinear diffusion. Phys. Rev. E 88, 042925.
  * (32) Sadighi A, Ganji D. 2007 Exact solutions of nonlinear diffusion equations by variational iteration method. Comp. Math. Appl. 54, 1112–1121.
  * (33) Newman ME. 2010 Networks: An Introduction. Oxford: Oxford University Press.
  * (34) Prigogine I, Nicolis G. 1967 Symmetry breaking instabilities in dissipative systems. J. Chem. Phys. 46, 3542.
  * (35) Boland RP, Galla T, McKane AJ. 2008 How limit cycles and quasi-cycles are related in a system with intrinsic noise. J. Stat. Mech. p. P09001.
  * (36) Neuhäuser L, Mellor A, Lambiotte R. 2020 Multibody interactions and nonlinear consensus dynamics on networked systems. Physical Review E 101, 032310.
  * (37) Sadighi A, Ganji D. 2007 Exact solutions of nonlinear diffusion equations by variational iteration method. Computers & Mathematics with Applications 54, 1112–1121.
  * (38) Pecora LM, Carroll TL. 1998 Master stability functions for synchronized coupled systems. Physical Review Letters 80, 2109.
  * (39) Carletti T, Muolo R. 2022 Non-reciprocal interactions enhance heterogeneity. Chaos, Solitons & Fractals 164, 112638.
  * (40) Haas PA, Goldstein RE. 2021 Turing’s Diffusive Threshold in Random Reaction-Diffusion Systems. Phys. Rev. Letters 126, 238101.
  * (41) Hata S, Nakao H, Mikhailov AS. 2012 Global feedback control of Turing patterns in network-organized activator-inhibitor systems. Eur. Phys. Lett. 98, 64004.
  * (42) Buscarino A, Corradino C, Fortuna L, Frasca M. 2019 Turing patterns via pinning control in the simplest memristive cellular nonlinear networks. Chaos 29, 103145.
  * (43) Gao S, Chang L, Romić I, Wang Z, Jusup M, Holme P. 2022 Optimal control of networked reaction–diffusion systems. Journal of the Royal Society Interface 19, 20210739.
  * (44) Challenger J, Fanelli D, Burioni R. 2015 Turing-like instabilities from a limit cycle. Phys. Rev. E 92, 022818.
  * (45) Muolo R, Carletti T, Gleeson JP, Asllani M. 2021 Synchronization dynamics in non-normal networks: the trade-off for optimality. Entropy 23, 36\.
  * (46) Gallo L, Muolo R, Gambuzza LV, Latora V, Frasca M, Carletti T. 2022 Synchronization induced by directed higher-order interactions. arXiv preprint arXiv:2202.08707.