# Pattern reconstruction through generalized eigenvectors on defective
networks

Marie Dorchain  Department of Mathematics and naXys, Namur Institute for
Complex Systems, University of Namur, Rue Grafé 2, B5000 Namur, Belgium
Riccardo Muolo  Department of Mathematics and naXys, Namur Institute for
Complex Systems, University of Namur, Rue Grafé 2, B5000 Namur, Belgium
Timoteo Carletti<EMAIL_ADDRESS>Department of Mathematics and
naXys, Namur Institute for Complex Systems, University of Namur, Rue Grafé 2,
B5000 Namur, Belgium

###### Abstract

Self-organization in natural and engineered systems causes the emergence of
ordered spatio-temporal motifs. In presence of diffusive species, Turing
theory has been widely used to understand the formation of such patterns on
continuous domains obtained from a diffusion-driven instability mechanism. The
theory was later extended to networked systems, where the reaction processes
occur locally (in the nodes), while diffusion takes place through the networks
links. The condition for the instability onset relies on the spectral property
of the Laplace matrix, i.e., the diffusive operator, and in particular on the
existence of an eigenbasis. In this work we make one step forward and we prove
the validity of Turing idea also in the case of a network with defective
Laplace matrix. Moreover, by using both eigenvectors and generalized
eigenvectors we show that we can reconstruct the asymptotic pattern with a
relatively small discrepancy. Because a large majority of empirical networks
are non-normal and often defective, our results pave the way for a thorough
understanding of self-organization in real-world systems.

Introduction. \- We are surrounded by patterns. Those spatio-temporal motifs
are the signature of the emergence of order from disorder Anderson (1972)
resulting from the collective behavior of the many nonlinearly interacting
basic units Pastor-Satorras and Vespignani (2010); Nicolis and Prigogine
(1977). In many relevant applications these interactions can be modeled by
using reaction-diffusion equations aiming at describing the behavior of
concentrations in time and space, being the latter a continuous domain Murray
(2001) or a discrete one, e.g., a complex network Nakao and Mikhailov (2010).
Indeed local reactions require, by their very first nature, species to be
spatially close, hence separated from other groups; it is thus natural to
consider species to occupy spatially limited zones, i.e., nodes of a network,
and diffuse across paths connecting different zones, i.e., the links of a
network. This will be the framework we will be interested in this work.

Alan Turing introduced and studied in the 50s a symmetry breaking mechanism
where a spatially homogeneous equilibrium of a reaction-diffusion system loses
its stability once disturbed with an heterogeneous perturbation; eventually
the system achieves a new, generally, patchy stationary or oscillatory
solution Turing (1952). Nowadays, Turing instability finds application beyond
the original framework of morphogenesis or chemical reaction systems Castets
_et al._ (1990); De Kepper  _et al._ (1991); Tompkins  _et al._ (2014) and it
stands for a pillar to explain self-organization in nature Nicolis and
Prigogine (1977); Pismen (2006), having being formalized by the existence of
an interplay between slow diffusing activators and fast diffusing inhibitors
Gierer and Meinhardt (1972). Indeed the latter determines a general mechanism:
a local feedback, i.e., short range production of a given species, which
should be, at the same time, inhibited at distance, by long range interaction.

The onset of Turing instability ultimately relies on the study of the spectral
properties of a suitable operator built by using the Jacobian of the reaction
part and the diffusion term, i.e., the Laplace operator. By assuming the
existence of an eigenbasis for the latter, one can compute the dispersion
relation, that ultimately determines the onset of the instability. Those ideas
have been largely applied to study the emergence of Turing patterns for system
whose underlying network is symmetric Nakao and Mikhailov (2010), directed
Asllani  _et al._ (2014); Carletti and Muolo (2022), but also for multiplex
Busiello  _et al._ (2015) and multilayer networks Asllani  _et al._ (2016),
temporal networks Petit  _et al._ (2017) and even in the novel framework of
higher-order structures, such as hypergraphs Muolo  _et al._ (2023) and
simplicial complexes Giambagli  _et al._ (2022). In particular, it has been
shown that the final pattern can be partially reconstructed by considering the
eigenvectors associated to the unstable modes, namely those for which the
dispersion relation is positive. Indeed in the linear regime, namely close to
the bifurcation, the pattern is completely aligned with those critical
eigenvectors; remarkably enough the nonlinearity of the model only slightly
perturbs this behavior and thus the final pattern can be accurately described
by a linear combination of critical eigenvectors Nakao and Mikhailov (2010).
The agreement is stronger the fewer is the number of unstable modes.

Scholars have recently pointed out that most real-world networks are non-
normal Asllani  _et al._ (2018); Duan  _et al._ (2022), namely their adjacency
matrix $\mathbf{A}$ does satisfy
$\mathbf{A}\mathbf{A}^{\top}\neq\mathbf{A}^{\top}\mathbf{A}$ Trefethen and
Embree (2005), or equivalently $\mathbf{A}$ is not diagonalizable through an
orthonormal transformation. Turing patterns on non-normal networks have been
recently studied with a numerical approach Asllani and Carletti (2018); Muolo
_et al._ (2019). These latter results however still rely on the assumption of
the existence of a basis of eigenvectors for the Laplace matrix.

The goal of this work is thus twofold. We first analytically solve the problem
on defective networks, i.e., networks whose Laplace matrix does not admit an
eigenbasis, and thus the eigenvalues have algebraic multiplicity larger than
one and greater than the geometric multiplicity. Then, we show how the pattern
reconstruction is improved when considering also the generalized eigenvectors
associated to the unstable modes.

Turing theory on defective networks. \- Let us consider two different species
populating a directed network composed by $n$ nodes and let us denote by
$u_{i}(t)$ and $v_{i}(t)$, $i=1,\dots,n$, their respective concentrations on
node $i$ at time $t$. When species happen to share the same node, they
interact via some generic nonlinear functions $f(u_{i},v_{i})$ and
$g(u_{i},v_{i})$. On the other hand, they can diffuse across the available
network links according to Fick’s law by taking into account the link
directionality. The model can hence be mathematically cast in the form

$\begin{dcases}\frac{du_{i}}{dt}&=f(u_{i},v_{i})+D_{u}\sum_{j=1}^{n}L_{ij}u_{j}\\\
\frac{dv_{i}}{dt}&=g(u_{i},v_{i})+D_{v}\sum_{j=1}^{n}L_{ij}v_{j}\end{dcases}\quad\forall
i=1,\dots,n\,,$ (1)

where $D_{u}>0$ (resp. $D_{v}>0$) is the diffusion coefficients of species $u$
(resp. $v$). The Laplace matrix, $L_{ij}=A_{ij}-\delta_{ij}k_{i}^{(in)}$, is
the discrete equivalent of the diffusion operator in the continuous support
case, where $A_{ij}$ is the $(i,j)$ entry of the adjacency matrix that allows
to encode the nodes connections, $A_{ij}=1$ if there is a link pointing from
node $j$ to node $i$, and $k_{i}^{(in)}=\sum_{j}A_{ij}$ is the in-degree of
node $i$.

In the spirit of Turing framework, we assume the existence of a stable
solution of (1) once we silence the diffusive part, namely there exists
$(u_{*},v_{*})$ such that $f(u_{*},v_{*})=g(u_{*},v_{*})=0$ and moreover
$\mathrm{tr}(\mathbf{J}_{0})<0$ and $\mathrm{det}\mathbf{J}_{0}>0$, where
$\mathbf{J}_{0}$ is the Jacobian matrix of the reaction part evaluated at the
equilibrium $(u_{*},v_{*})$

$\mathbf{J}_{0}=\left(\begin{matrix}\partial_{u}f&\partial_{v}f\\\
\partial_{u}g&\partial_{v}g\end{matrix}\right)\Big{\rvert}_{(u_{*},v_{*})}\,,$
(2)

where we denoted by $\partial_{u}f$ the derivative of $f$ with respect to $u$
evaluated on the equilibrium $(u_{*},v_{*})$, and similarly for the other
terms.

We then require such equilibrium to turn out unstable once diffusion is at
play. To verify such condition we perform a linear stability analysis, namely
we introduce a perturbation from the homogeneous solution $\delta
u_{i}(t)=u_{i}(t)-u_{*}$ and $\delta v_{i}(t)=v_{i}(t)-v_{*}$, and expand Eq.
(1), keeping only the first order terms in the perturbation (the latter
assumed to be small). We thus obtain for all $i=1,\dots,n$

$\begin{dcases}\frac{d\delta u_{i}}{dt}&=\partial_{u}f\delta
u_{i}+\partial_{v}f\delta v_{i}+D_{u}\sum_{j=1}^{n}L_{ij}\delta u_{j}\\\
\frac{d\delta v_{i}}{dt}&=\partial_{u}g\delta u_{i}+\partial_{v}g\delta
v_{i}+D_{v}\sum_{j=1}^{n}L_{ij}\delta v_{j}\,,\end{dcases}$ (3)

where we employed the fact that $\sum_{j}L_{ij}=0$ to nullify the terms
$\sum_{j}L_{ij}u_{*}$ and $\sum_{j}L_{ij}v_{*}$.

By introducing the $n\times 2$ vector $\delta\mathbf{x}=(\delta u_{1},\delta
v_{1},\dots,\delta u_{n},\delta v_{n})^{\top}$, we can eventually rewrite the
latter equation in a compact form as:

$\frac{d\delta\mathbf{x}}{dt}=\left[\mathbf{I}_{n}\otimes\mathbf{J}_{0}+\mathbf{L}\otimes\mathbf{D}\right]\delta\mathbf{x}\,,$
(4)

where $\otimes$ denotes the Kronecker product of matrices, $\mathbf{I}_{n}$ is
the $n\times n$ identity matrix and
$\mathbf{D}=\left(\begin{smallmatrix}D_{u}&0\\\
0&D_{v}\end{smallmatrix}\right)$.

To make some analytical progress, the standard step is thus to simplify the
previous $2n\times 2n$ system into $n$ systems $2\times 2$ by assuming the
existence of an eigenbasis for the Laplace matrix and projecting the
perturbations $\delta u_{i}$ and $\delta v_{i}$ upon such basis. Our goal is
to show that one can obtain a similar understanding of the onset of Turing
instability also in the case the Laplace matrix is defective. Such framework
has been studied in Nishikawa and Motter (2006) in the study of
synchronization of coupled oscillators.

To achieve this goal, we can invoke the Jordan canonical form to determine an
invertible $n\times n$ matrix $\mathbf{P}$ such that

$\mathbf{P}^{-1}\mathbf{L}\mathbf{P}=\mathbf{B}=\mathrm{diag}(\mathbf{B}_{1},\dots,\mathbf{B}_{\ell})\,,$

where the $\mathbf{B}_{j}$ is the $m_{j}\times m_{j}$ Jordan block, being
$m_{j}$ the algebraic multiplicity of the eigenvalue $\Lambda^{(j)}$,
$m_{1}+\dots+m_{\ell}=n$ and

$\mathbf{B}_{j}=\left(\begin{matrix}\Lambda^{(j)}&&&\\\ 1&\Lambda^{(j)}&&\\\
&\ddots&\ddots&\\\ &&1&\Lambda^{(j)}\end{matrix}\right)\,.$ (5)

Let us consider again Eq. (4). By defining
$\mathbf{Q}=\mathbf{P}\otimes\mathbf{I}_{2}$ and
${\color[rgb]{0,0,0}{\delta\mathbf{y}=\mathbf{Q}^{-1}\delta\mathbf{x}}}$ we
get

$\displaystyle\frac{d\delta\mathbf{y}}{dt}$ $\displaystyle=$
$\displaystyle\mathbf{Q}^{-1}\frac{d\delta\mathbf{x}}{dt}$ (6)
$\displaystyle=$
$\displaystyle(\mathbf{P}^{-1}\otimes\mathbf{I}_{2})\left[\mathbf{I}_{n}\otimes\mathbf{J}_{0}+\mathbf{L}\otimes\mathbf{D}\right](\mathbf{P}\otimes\mathbf{I}_{2})\mathbf{Q}^{-1}\delta\mathbf{x}$
$\displaystyle=$
$\displaystyle\left[\mathbf{I}_{n}\otimes\mathbf{J}_{0}+\mathbf{B}\otimes\mathbf{D}\right]\delta\mathbf{y}\,.$

The vector $\delta\mathbf{y}$ inherits the Jordan decomposition, hence we can
write
$\delta\mathbf{y}=((\delta\mathbf{y}^{(1)})^{\top},\dots,(\delta\mathbf{y}^{(\ell)})^{\top})^{\top}$,
where $\delta\mathbf{y}^{(j)}$ is a $(2\times m_{j})$-dimensional vector.

The stability properties of $\delta\mathbf{y}$ will thus be determined by
analyzing the behavior of $\delta\mathbf{y}^{(j)}$, $j=1,\dots,\ell$. Let us
consider separately the case $\Lambda^{(1)}=0$ and
${\color[rgb]{0,0,0}{\Re\Lambda^{(j)}<0}}$ for $j\geq 2$. Assume thus
$\Lambda^{(1)}=0$ to be degenerate and be $m_{1}$ its multiplicity. If
$m_{1}=1$ then $\delta\mathbf{y}^{(1)}$ evolves accordingly to
$\mathbf{I}_{n}\otimes\mathbf{J}_{0}$ and thus $\delta\mathbf{y}^{(1)}$ is
stable because of the condition imposed on the homogeneous equilibrium
$(u_{*},v_{*})$ 111Let us observe that a similar result can be obtained once
the degeneracy is a consequence of the presence of nodes without incoming
links, called leader nodes O’Brien  _et al._ (2021), as the corresponding
eigenvectors are the canonical ones, and the degeneration is algebraic and not
geometric.. Otherwise $\mathbf{B}_{1}$ is a $m_{1}\times m_{1}$ matrix of the
form (5) with $0$ on the diagonal.

The part of Eq. (6) relative to $\Lambda^{(1)}$ can thus be rewritten as

$\displaystyle\frac{d\delta\mathbf{y}^{(1)}}{dt}$ $\displaystyle=$
$\displaystyle\left[\left(\begin{matrix}\mathbf{J}_{0}&&\\\ &\ddots&\\\
&&\mathbf{J}_{0}\end{matrix}\right)+\left(\begin{matrix}0&&&\\\
\mathbf{D}&0&&\\\ &\ddots&\ddots&\\\
&&\mathbf{D}&0\end{matrix}\right)\right]\delta\mathbf{y}^{(1)}\,,$

the matrix on the right hand side has the same eigenvalues of $\mathbf{J}_{0}$
from which we can conclude that $\delta\mathbf{y}^{(1)}$ is stable.

We can now analyze the remaining cases $j=2,\dots,\ell$. The part of Eq. (6)
involving $\delta\mathbf{y}^{(j)}$ is thus

$\displaystyle\frac{d\delta\mathbf{y}^{(j)}}{dt}$ $\displaystyle=$
$\displaystyle\left(\begin{matrix}\mathbf{J}_{0}&&\\\ &\ddots&\\\
&&\mathbf{J}_{0}\end{matrix}\right)\delta\mathbf{y}^{(j)}+$
$\displaystyle+\left(\begin{matrix}\Lambda^{(j)}\mathbf{D}&&&\\\
\mathbf{D}&\Lambda^{(j)}\mathbf{D}&&\\\ &\ddots&\ddots&\\\
&&\mathbf{D}&\Lambda^{(j)}\mathbf{D}\end{matrix}\right)\delta\mathbf{y}^{(j)}\,.$

We can reformulate the previous equation by writing
$\delta\mathbf{y}^{(j)}=(\xi^{(j)}_{1},\eta^{(j)}_{1},\dots,\xi^{(j)}_{m_{j}},\eta^{(j)}_{m_{j}})^{\top}$,
where
$\delta\mathbf{y}_{i}^{(j)}$$=(\xi^{(j)}_{i},\eta^{(j)}_{i})\in\mathbb{R}^{2}$
for all $i=1,\dots,m_{j}$, thus obtaining

$\displaystyle\frac{d}{dt}\binom{\xi^{(j)}_{1}}{\eta^{(j)}_{1}}$
$\displaystyle=$
$\displaystyle\mathbf{J}_{0}\binom{\xi^{(j)}_{1}}{\eta^{(j)}_{1}}+\Lambda^{(j)}\mathbf{D}\binom{\xi^{(j)}_{1}}{\eta^{(j)}_{1}}\,,$
(8) $\displaystyle\frac{d}{dt}\binom{\xi^{(j)}_{2}}{\eta^{(j)}_{2}}$
$\displaystyle=$
$\displaystyle\mathbf{J}_{0}\binom{\xi^{(j)}_{2}}{\eta^{(j)}_{2}}+\Lambda^{(j)}\mathbf{D}\binom{\xi^{(j)}_{2}}{\eta^{(j)}_{2}}+\mathbf{D}\binom{\xi^{(j)}_{1}}{\eta^{(j)}_{1}}\,,$
(9) $\displaystyle\vdots$
$\displaystyle\frac{d}{dt}\binom{\xi^{(j)}_{m_{j}}}{\eta^{(j)}_{m_{j}}}$
$\displaystyle=$
$\displaystyle\mathbf{J}_{0}\binom{\xi^{(j)}_{m_{j}}}{\eta^{(j)}_{m_{j}}}+\Lambda^{(j)}\mathbf{D}\binom{\xi^{(j)}_{m_{j}}}{\eta^{(j)}_{m_{j}}}+\mathbf{D}\binom{\xi^{(j)}_{m_{j}-1}}{\eta^{(j)}_{m_{j}-1}}\,.$

The first Eq. (8) is the same equation one would get once the Laplace matrix
admits an eigenbasis, thus one can determine a condition on $\Lambda^{(j)}$ to
make the projection
$\delta\mathbf{y}_{1}^{(j)}=(\xi^{(j)}_{1},\eta^{(j)}_{1})^{\top}$ unstable
Asllani  _et al._ (2014), namely to compute the eigenvalue with the largest
real part of the matrix
$\mathbf{M}_{j}=\mathbf{J}_{0}+\Lambda^{(j)}\mathbf{D}$ (see also Appendix A).

Let us now consider Eqs. (9) and observe that each of them is composed by two
terms, the first one involving the same matrix of the former equation,
$\mathbf{M}_{j}$, while the second one depends on the projection
$\delta\mathbf{y}^{(j)}(t)$. If, for the choice of $\Lambda^{(j)}$ the matrix
$\mathbf{M}_{j}$ is unstable, namely its spectrum contains eigenvalues with
positive real part, and thus $\delta\mathbf{y}_{1}^{(j)}(t)$ has an
exponential growth, then the same is true for
$\delta\mathbf{y}_{2}^{(j)}(t)=(\xi^{(j)}_{2},\eta^{(j)}_{2})^{\top}$. By
considering the remaining equations and by exploiting the peculiar lower
triangular shape of the system, we can prove that if Eq. (8) returns an
unstable solution, then all the solutions
$\delta\mathbf{y}_{i}^{(j)}(t)=(\xi^{(j)}_{i},\eta^{(j)}_{i})^{\top}$ are
unstable as well.

In conclusion one can compute the dispersion relation, $\lambda(\zeta)$,
namely to determine the largest real part of the eigenvalues of the
$1$-(complex) parameter family system $\mathbf{J}_{0}+\zeta\mathbf{D}$; if the
Laplace matrix $\mathbf{L}$ is defective one can check the instability
condition on the available eigenvalues. This accounts to study the sign of
$\lambda(\Lambda^{(j)})$, and conclude about the emergence of patterns solely
based on this information. Let us observe that this is a sufficient condition,
indeed it can happen that the matrix $\mathbf{M}_{j}$ is stable for all $j$,
i.e., all its eigenvalues have negative real part, but the presence of Jordan
blocks introduces a transient (polynomial) growth in the linear regime that
results strong enough to limit the validity of the linear approximation.
Stated differently, the size of the basin of attraction of the stable fixed
point considerably shrinks because of this transient growth. Thus the
nonlinear system could exhibit orbits departing from the homogeneous reference
solution; only infinitesimal perturbations will be attracted to the latter,
the solution is thus stable but finite perturbations can be amplified. Hence
the latter result extends and completes the numerical analysis performed in
the case of diagonalizable non-normal networks Muolo  _et al._ (2019, 2021).

A case study: the Brusselator model. \- Let us present the described theory by
considering the Brusselator model Prigogine and Nicolis (1967); Prigogine and
Lefever (1968); Boland  _et al._ (2008), often invoked in the literature as a
paradigm nonlinear reaction scheme for studying self-organized phenomena such
as synchronization Muolo  _et al._ (2021), Turing patterns Asllani  _et al._
(2014) and oscillation death Koseska  _et al._ (2013); Lucas  _et al._ (2018).
The key feature of the model is the presence of two species, reacting via a
cubic nonlinearity

$\begin{dcases}\frac{du}{dt}&=1-(b+1)u+cu^{2}v\\\ \frac{dv}{dt}&=bu-
cu^{2}v\,,\end{dcases}$ (10)

where $b>0$ and $c>0$ act as tunable model parameters. One can easily realize
the existence of a unique equilibrium $u_{*}=1$ and $v_{*}=b/c$, that results
stable if the Jacobian of the reaction part evaluated on it,
$\mathbf{J}_{0}=\left(\begin{smallmatrix}b-1&c\\\
-b&-c\end{smallmatrix}\right)$, has a negative trace,
$\mathrm{tr}(\mathbf{J}_{0})=b-c<1$, and a positive determinant
$\det\mathbf{J}_{0}=c>0$.

By considering $n$ identical copies of the Brusselator model, each one
anchored on a node of a network and interacting with the first neighbors, we
obtain $\forall i=1,\dots,n$

$\begin{dcases}\frac{du_{i}}{dt}&=1-(b+1)u_{i}+cu_{i}^{2}v_{i}+D_{u}\sum_{j=1}^{n}{L}_{ij}u_{j}\\\
\frac{dv_{i}}{dt}&=bu_{i}-cu_{i}^{2}v_{i}+D_{v}\sum_{j=1}^{n}{L}_{ij}v_{j}\,.\end{dcases}$
(11)

By linearizing the above equation about the homogeneous equilibrium and by
using the Jordan blocks we can obtain the analogous of Eqs. (8) - (9). Then,
we can determine the region in the complex plane (see Appendix A) associated
to the Turing instability if at least one eigenvalue of the Laplace matrix
falls into this region. Let us remark that in the following we always deal
with an unstable dispersion relation, i.e., there exists at least one
eigenvalue $\Lambda^{(j)}$ such that Eq. (Pattern reconstruction through
generalized eigenvectors on defective networks) is unstable. In the top panel
of Fig. 1 we show the region of instability for the Brusselator model defined
on a defective network composed by $N=10$ nodes (see Fig. 2), where the stable
eigenvalues are $\Lambda^{(1)}=0$, with multiplicity $2$, $\Lambda^{(2)}=-1$
with multiplicity $5$ and $\Lambda^{(3)}=-2$ with multiplicity $1$. There is
only one unstable eigenvalue $\Lambda^{(4)}=-4$ with multiplicity $2$.

Figure 1: Region of the complex plane associated to Turing instability (top
panel) and dispersion relation (bottom panel) computed for the Brusselator
model with parameters $b=3.92$, $c=3$, $D_{u}=0.2$ and $D_{v}=0.8$. Top panel:
The black dots denote the eigenvalues of the Laplace matrix, $\Lambda^{(1)}=0$
(multiplicity $2$), $\Lambda^{(2)}=-1$ (multiplicity $5$), $\Lambda^{(3)}=-2$
(multiplicity $1$), $\Lambda^{(4)}=-4$ (multiplicity $2$), the green region is
associated to a positive dispersion relation, while the white one to the
negative case. Bottom panel: the largest real part of the spectrum of the
matrix $\mathbf{J}_{0}+\zeta\mathbf{D}$ is shown in blue as a function of
$\zeta$, the dispersion relation evaluated on the Laplace spectrum is reported
by using red triangles.

In the bottom panel of Fig. 1 we report the dispersion relation computed for
the Brusselator defined on the same defective network. The Laplace eigenvalues
are represented by symbols (red triangles) while the continuous curve is the
dispersion relation computed for the $1$-(complex) parameter family of linear
systems introduced above with the matrix $\mathbf{M}_{j}$. One can observe
that $\lambda(\Lambda^{(j)})$ is positive for $\Lambda^{(4)}$ and thus the
equilibrium $(u_{*},v_{*})\sim(1,1.3067)$ of the coupled system is unstable,
as we can appreciate by inspecting Fig. 3 where we report the time evolution
of the concentrations $u_{i}(t)$ vs. time. The same conclusion can be obtained
by observing the top panel of Fig. 1 where we can realize that $\Lambda^{(4)}$
lies inside the instability region (green area).

Figure 2: Random non-normal defective network composed by $n=10$ nodes, built
by using a directed Erdős-Rényi algorithm where the probability to create a
bidirectional link is $0.2$ and the probability to transform it into a
directed one is $0.6$. Nodes have been colored according to the value of
species $u$ at time $\hat{t}=200$ (see colorbar). Figure 3: Evolution of the
concentration of the specie $u_{i}$ over time for the Brusselator model with
parameters $b=3.92$, $c=3$, $D_{u}=0.2$ and $D_{v}=0.8$. The underlying
network is the one shown in Fig. 2. The nodes concentrations have been
initialized to the homogeneous equilibrium, $u_{*}=1$, upon which a random,
node dependent, perturbation of order $0.01$ has been added. The resulting
orbits have been obtained by using a Runge-Kutta 4th method with time step
$0.01$. Each trajectory has been represented by using the same color of the
corresponding node in Fig. 2, namely the value at time $\hat{t}=200$.

Pattern reconstruction through generalized eigenvectors. \- When there exists
an eigenbasis for the Laplace matrix, we can show that pattern can be
described by using the eigenvectors related to the unstable eigenvalues, i.e.,
the eigenvalues of the Laplace matrix that return a positive dispersion
relation or equivalently they lie in the instability region as shown above.
Our goal is to show that a similar result can be obtained in the case of
defective Laplace matrix by recurring to generalized eigenvectors to
reconstruct the pattern. To achieve such goal, we use again the Brusselator
model (11) defined on top of the defective random non-normal network presented
above (see Fig. 2). Let us remember that the Laplace matrix of the network has
an unstable eigenvalue, $\Lambda^{(4)}$ with multiplicity $2$, to which we
associate the eigenvector $\phi^{(4)}$ and a generalized eigenvector
$v^{(4)}$.

To support our claim, we considered the solution $u_{i}(t)$ of model (11) up
to certain (large) time, $\hat{t}=250$, and we thus obtain the vector
$\hat{u}=(u_{1}(\hat{t}),\dots,u_{10}(\hat{t}))^{\top}$, i.e., the asymptotic
pattern. We then proceed by reconstructing (see Appendix B) such pattern by
using the eigenvector $\phi^{(4)}$ (see upper panel of Fig. 4) or the above
eigenvector together with the generalized eigenvector $v^{(4)}$ (see lower
panel of Fig. 4). We can observe that in both cases the reconstructed pattern
is very close to the original one and moreover the one obtained by using also
the generalized eigenvector is noticeably improved, having a smaller error.
Let us observe that as the reconstruction would obviously be better with two
vectors than with one, a weighted absolute error was used to compare the
accuracy of the reconstruction. Indeed Each absolute error is weighted by the
number of (eigen)vectors used, divided by the total number of available
unstable (eigen)vectors, namely the geometric multiplicity. In this case we
hence obtain $\nicefrac{{1}}{{2}}$ in the first case and $\nicefrac{{2}}{{2}}$
in the second one (see Appendix B for a rigorous definition of the used
weighted error).

Figure 4: Pattern vs. reconstructed pattern for the Brusselator model with
parameters $b=3.92$, $c=3$, $D_{u}=0.2$ and $D_{v}=0.8$. Upper panel: the
reconstruction is obtained with only the eigenvectors (EV). Lower panel: the
eigenvectors and generalized eigenvectors (GEV) are used for the
reconstruction. In the first case, the reconstruction error is
$\varepsilon=0.088$ while in the second one we have $\varepsilon=0.015$.

Let us notice that our results are robust with respect to the number of
(generalized) unstable eigenvectors and their multiplicity (see Appendix D and
Fig. 16), the chosen value of $\hat{t}$ or, in the case of oscillatory
pattern, to the fact of considering the time-average of $u_{i}(t)$, i.e.,
$\langle u\rangle=(\langle u_{1}\rangle,\dots,\langle u_{10}\rangle)^{\top}$
where $\langle u\rangle_{i}$ is the time-average of the orbit of the $i$-th
node (see Appendix D and Fig. 11). Let us conclude by observing that the
proposed pattern reconstruction method scales well with the increasing size of
the network. Moreover, the results we obtain by relying on a family of
directed defective networks, support the claim that the use of generalized
eigenvectors always returns a better estimate of the pattern, than the
eigenvectors alone (see Appendix D and Fig. 17).

Conclusions. \- In this paper we have further extended the Turing theory of
pattern formation by studying the case of defective networks. After solving
the problem analytically and showing the effects on the instability mechanism
given by the presence of Jordan blocks, we have shown the pivotal role of
generalized eigenvectors in the reconstruction of the asymptotic pattern.
Considering that most real-world networks are non-normal, our results become
particularly relevant even looking at other nonlinear phenomena beyond Turing
pattern formation and may help in filling the gap between theory and
observations. Moreover, the proposed pattern reconstruction method further
improves the understanding of the interplay between the dynamics and the
underlying topology, paving the way for finer methods of network
reconstruction from observational data.

## Acknowledgements

R.M. acknowledges funding from the FNRS, Grant FC 33443, funded by the Walloon
region.

## References

  * Anderson (1972) P. Anderson, Phys. Rev. Lett. 177 (1972).
  * Pastor-Satorras and Vespignani (2010) R. Pastor-Satorras and A. Vespignani, Nat. Phys. 6, 480 (2010).
  * Nicolis and Prigogine (1977) G. Nicolis and I. Prigogine,  _Self-organization in nonequiibrium systems: From dissipative structures to order through fluctuations_ (J. Wiley and Sons, 1977).
  * Murray (2001) J. Murray,  _Mathematical biology II: Spatial models and biomedical applications_ (Springer-Verlag, 2001).
  * Nakao and Mikhailov (2010) H. Nakao and A. Mikhailov, Nat. Phys. 6, 544 (2010).
  * Turing (1952) A. Turing, Phil. Trans. R. Soc. Lond. B 237, 37 (1952).
  * Castets  _et al._ (1990) V. Castets, E. Dulos, J. Boissonade, and P. De Kepper, Phys. Rev. Lett. 64, 2953 (1990).
  * De Kepper  _et al._ (1991) P. De Kepper, V. Castets, E. Dulos, and J. Boissonade, Physica D 49 (1991).
  * Tompkins  _et al._ (2014) N. Tompkins, N. Li, C. Girabawe, M. Heymann, G. Ermentrout, I. Epstein, and S. Fraden, Proc. Natl. Acad. Sci. U.S.A. 111, 4397 (2014).
  * Pismen (2006) L. Pismen,  _Patterns and interfaces in dissipative dynamics_ (Springer Science & Business Media, 2006).
  * Gierer and Meinhardt (1972) A. Gierer and H. Meinhardt, Kybernetik 12, 30 (1972).
  * Asllani  _et al._ (2014) M. Asllani, J. Challenger, F. Pavone, L. Sacconi, and D. Fanelli, Nature Communication 5 (2014).
  * Carletti and Muolo (2022) T. Carletti and R. Muolo, Chaos Solit. Fractals 164, 112638 (2022).
  * Busiello  _et al._ (2015) D. Busiello, G. Planchon, M. Asllani, T. Carletti, and D. Fanelli, Eur. Phys. J. B 88, 222 (2015).
  * Asllani  _et al._ (2016) M. Asllani, T. Carletti, and D. Fanelli, Eur. Phys. J. B , 89 (2016).
  * Petit  _et al._ (2017) J. Petit, B. Lauwens, D. Fanelli, and T. Carletti, Phys. Rev. Lett. 119, 148301 (2017).
  * Muolo  _et al._ (2023) R. Muolo, L. Gallo, V. Latora, M. Frasca, and T. Carletti, Chaos Solit. Fractals 166, 112912 (2023).
  * Giambagli  _et al._ (2022) L. Giambagli, M. Calmon, R. Muolo, T. Carletti, and G. Bianconi, Phys. Rev. E 106 (2022).
  * Asllani  _et al._ (2018) M. Asllani, R. Lambiotte, and T. Carletti, Sci. Adv. 4, Eaau9403 (2018).
  * Duan  _et al._ (2022) C. Duan, T. Nishikawa, D. Eroglu, and A. E. Motter, Science Advances 8, eabm8310 (2022).
  * Trefethen and Embree (2005) L. Trefethen and M. Embree,  _Spectra and Pseudospectra: The Behavior of Nonnormal Matrices and Operators_ (Princeton Univ. Press, 2005).
  * Asllani and Carletti (2018) M. Asllani and T. Carletti, Phys. Rev. E 97 (2018).
  * Muolo  _et al._ (2019) R. Muolo, M. Asllani, D. Fanelli, P. Maini, and T. Carletti, J. Theor. Biol. 480, 81 (2019).
  * Nishikawa and Motter (2006) T. Nishikawa and A. Motter, Phys. Rev. E 73, 065106(R) (2006).
  * Muolo  _et al._ (2021) R. Muolo, T. Carletti, J. Gleeson, and M. Asllani, Entropy 23, 36 (2021).
  * Prigogine and Nicolis (1967) I. Prigogine and G. Nicolis, J. Chem. Phys. 46, 3542 (1967).
  * Prigogine and Lefever (1968) I. Prigogine and R. Lefever, J. Chem. Phys. 48, 1695 (1968).
  * Boland  _et al._ (2008) R. Boland, T. Galla, and A. McKane, J. Stat. Mech. , P09001 (2008).
  * Koseska  _et al._ (2013) A. Koseska, E. Volkov, and J. Kurths, Phys. Rep. 531, 173 (2013).
  * Lucas  _et al._ (2018) M. Lucas, D. Fanelli, T. Carletti, and J. Petit, EPL 121, 58008 (2018).
  * O’Brien  _et al._ (2021) J. O’Brien, K. Oliveira, J. Gleeson, and M. Asllani, Phys. Rev. Res. 3, 023117 (2021).
  * Forrow  _et al._ (2018) A. Forrow, F. Woodhouse, and J. Dunkel, Phys. Rev. X 8, 041043 (2018).
  * Nicoletti  _et al._ (2021) S. Nicoletti, T. Carletti, D. Fanelli, G. Battistelli, and L. Chisci, J. Phys. Complex. 2 (2021).

##  Appendix A Condition for the Turing instability.

In the main text we showed (see Eqs. (8)-(9)) that Turing instability can be
determined by studying the spectrum of the matrix
$\mathbf{M}_{j}=\mathbf{J}_{0}+\Lambda^{(j)}\mathbf{D}$, where
$\mathbf{J}_{0}$ is the Jacobian of the reaction part (2) evaluated at the
homogeneous equilibrium $(u_{*},v_{*})$, $\mathbf{D}$ is the diagonal matrix
of the diffusion coefficients and $\Lambda^{(j)}$ is the $j$-th eigenvalue of
the Laplace matrix $\mathbf{L}$. Moreover because the homogeneous solution is
stable we also have $\mathrm{tr}(\mathbf{J}_{0})<0$ and
$\det\mathbf{J}_{0}>0$.

The eigenvalues of $\mathbf{M}_{j}$ can be straightforwardly obtained by using
the formula

$\mu_{j}=\frac{\mathrm{tr}(\mathbf{M}_{j})+\sqrt{(\mathrm{tr}(\mathbf{M}_{j}))^{2}-4\det\mathbf{M}_{j}}}{2}\,;$

because the network is directed the Laplace matrix is generally asymmetric,
thus its eigenvalues are complex numbers and so does $\mu_{j}$. By writing
$\Lambda^{(j)}=\Re\Lambda^{(j)}+i\Im\Lambda^{(j)}$ we can obtain

$\displaystyle\Re\mathrm{tr}(\mathbf{M}_{j})$ $\displaystyle=$
$\displaystyle\mathrm{tr}(\mathbf{J}_{0})+(D_{u}+D_{v})\Re\Lambda^{(j)}$
$\displaystyle\Im\mathrm{tr}(\mathbf{M}_{j})$ $\displaystyle=$
$\displaystyle(D_{u}+D_{v})\Im\Lambda^{(j)}$
$\displaystyle\Re\det(\mathbf{M}_{j})$ $\displaystyle=$
$\displaystyle\det(\mathbf{J}_{0})+\left(({J}_{0})_{11}D_{v}+({J}_{0})_{22}D_{u}\right)\Re\Lambda^{(j)}+D_{u}D_{v}\left[(\Re\Lambda^{(j)})^{2}-(\Im\Lambda^{(j)})^{2}\right]$
$\displaystyle\Im\det(\mathbf{M}_{j})$ $\displaystyle=$
$\displaystyle\left(({J}_{0})_{11}D_{v}+({J}_{0})_{22}D_{u}\right)\Im\Lambda^{(j)}+2D_{u}D_{v}\Re\Lambda^{(j)}\Im\Lambda^{(j)}\,,$

and eventually express the real part of $\mu_{j}$ as follows

$\Re\mu_{j}=\frac{\Re\mathrm{tr}(\mathbf{M}_{j})+\gamma}{2}\,,$

where

$\gamma=\sqrt{\frac{A+\sqrt{A^{2}+B^{2}}}{2}}\,,A=\left(\Re\mathrm{tr}(\mathbf{M}_{j})\right)^{2}-\left(\Im\mathrm{tr}(\mathbf{M}_{j})\right)^{2}-4\Re\det(\mathbf{M}_{j})\text{
and
}B=2\Re\mathrm{tr}(\mathbf{M}_{j})\Im\mathrm{tr}(\mathbf{M}_{j})-4\Re\det(\mathbf{M}_{j})\,.$

The condition for instability is $\Re\mu_{j}>0$ for some $j$, namely

$\Re\mathrm{tr}(\mathbf{M}_{j})>-\gamma\,,$

that can be rewritten as

$S_{2}\left(\Re\Lambda^{(j)}\right)\left(\Im\Lambda^{(j)}\right)^{2}<-S_{1}\left(\Re\Lambda^{(j)}\right)\,,$
(12)

where $S_{1}$ (resp. $S_{2}$) is a polynomial of fourth (resp. second) degree
in $\Re\Lambda^{(j)}$. More precisely

$\displaystyle S_{1}\left(x\right)$ $\displaystyle=$ $\displaystyle
C_{14}x^{4}+C_{13}x^{3}+C_{12}x^{2}+C_{11}x+C_{10}$ $\displaystyle
S_{2}\left(x\right)$ $\displaystyle=$ $\displaystyle
C_{22}x^{2}+C_{21}x+C_{20}\,,$

where the coefficients $C_{ij}$ are explicitly given by

$\displaystyle C_{14}$ $\displaystyle=$ $\displaystyle
D_{u}D_{v}(D_{u}+D_{v})^{2}$ $\displaystyle C_{13}$ $\displaystyle=$
$\displaystyle(D_{u}+D_{v})^{2}\left((J_{0})_{11}D_{v}+(J_{0})_{22}D_{u}\right)+2\mathrm{tr}(\mathbf{J}_{0})D_{u}D_{v}(D_{u}+D_{v})$
$\displaystyle C_{12}$ $\displaystyle=$
$\displaystyle\det(\mathbf{J}_{0})(D_{u}+D_{v})^{2}+\left(\mathrm{tr}(\mathbf{J}_{0})\right)^{2}D_{u}D_{v}+2\mathrm{tr}(\mathbf{J}_{0})(D_{u}+D_{v})\left((J_{0})_{11}D_{v}+(J_{0})_{22}D_{u}\right)$
$\displaystyle C_{11}$ $\displaystyle=$ $\displaystyle
2\mathrm{tr}(\mathbf{J}_{0})(D_{u}+D_{v})\det(\mathbf{J}_{0})+\left(\mathrm{tr}(\mathbf{J}_{0})\right)^{2}\left((J_{0})_{11}D_{v}+(J_{0})_{22}D_{u}\right)$
$\displaystyle C_{10}$ $\displaystyle=$
$\displaystyle\det(\mathbf{J}_{0})\left(\mathrm{tr}(\mathbf{J}_{0})\right)^{2}\,,$

and

$\displaystyle C_{22}$ $\displaystyle=$ $\displaystyle
D_{u}D_{v}(D_{u}-D_{v})^{2}$ $\displaystyle C_{21}$ $\displaystyle=$
$\displaystyle\left((J_{0})_{11}D_{v}+(J_{0})_{22}D_{u}\right)(D_{u}-D_{v})^{2}$
$\displaystyle C_{20}$ $\displaystyle=$
$\displaystyle(J_{0})_{11}(J_{0})_{22}(D_{u}-D_{v})^{2}\,.$

In conclusion, whenever the Laplace matrix admits at least one eigenvalue such
that condition (12) is satisfied, then the homogeneous equilibrium turns out
to be unstable once submitted to heterogeneous perturbations. The region of
the complex plane $(\Re\Lambda,\Im\Lambda)$ where this condition holds true is
the green region shown in Fig. 1 in the main text and Fig. 8 in Appendix D.

##  Appendix B Pattern reconstruction

To reconstruct the pattern, $u(\hat{t})$, we first centered it by subtracting
the homogeneous equilibrium $u_{*}$. The goal is thus to project the obtained
vector, ${p}=u(\hat{t})-u_{*}\in\mathbb{R}^{n}$, on the subspace $V$ generated
by the $d\geq 1$ linearly independent vectors considered for the
reconstruction, i.e., unstable eigenvectors and generalized ones. Let us
denote by $\pi_{V}({p})$ the projection of ${p}$ onto $V$. This vector can
thus be expressed by a linear combination of those generating vectors

$\pi_{V}({p})=\mathbf{B}{a}\,,$ (13)

where ${a}$ is the vector containing the coefficients of the linear
combination and $\mathbf{B}$ is the matrix with the $d$ linearly independent
vectors as columns. By using basic algebra, one can express the centered
pattern as ${p}=\pi_{V}({p})+{z}$, with $\pi_{V}({p})\in V$ and ${z}\in
V^{\perp}$. By invoking the image-kernel theorem one can write
$V^{\perp}=\Im(\mathbf{B})^{\perp}=\ker(\mathbf{B}^{\top})$,
$z={p}-\pi_{V}({p})$ lies therefore in $\ker(\mathbf{B}^{\top})$, namely
$\mathbf{B}^{\top}({p}-\pi_{V}({p}))={0}$. By using the expression (13) we can
rewrite the latter as

$\mathbf{B}^{\top}({p}-\mathbf{B}{a})={0}\,.$

By developing the computations we obtain

${a}=(\mathbf{B}^{\top}\mathbf{B})^{-1}\mathbf{B}^{\top}{p}\,,$

and eventually the following expression for the projected centered pattern
$\pi_{V}({p})=\mathbf{B}(\mathbf{B}^{\top}\mathbf{B})^{-1}\mathbf{B}^{\top}{p}$.
Let us observe that if the (generalized) eigenvectors are complex, and so does
the matrix $\mathbf{B}$, we decided to replace each complex vector with two
real ones obtained by taking the real and the imaginary part of the former
complex vector.

For sake of clarity in the Figures presenting the pattern, the equilibrium has
been added back to the projection to better compare $\pi_{V}({p})+u_{*}$ with
the pattern $u(\hat{t})$.

We can now focus more on the vectors we use for the reconstruction. As said in
the main text, when the Laplace matrix is diagonalizable we just compute the
eigenvectors and use those associated with unstable eigenvalues. When we work
with defective networks, the Laplace matrix does not have a linearly
independent set of eigenvectors. We therefore compute the Jordan canonical
form of the Laplace matrix, as seen in the main text,
$\mathbf{P}^{-1}\mathbf{L}\mathbf{P}=\mathbf{B}=\mathrm{diag}(\mathbf{B}_{1},\dots,\mathbf{B}_{\ell})\,$
where the $\mathbf{B}_{j}$ is the $m_{j}\times m_{j}$ Jordan block,
$m_{1}+\dots+m_{\ell}=n$. The matrix $\mathbf{P}$ has the (generalized)
eigenvectors as columns, which are linearly independent. We then use the
(generalized) eigenvectors corresponding to unstable eigenvalues in the
reconstruction.

In order to compare the reconstructions with and without the generalized
eigenvectors, we used a weighted absolute error. If we note by
$\tilde{u}=(\tilde{u}_{1},\dots,\tilde{u}_{n})^{\top}=\pi_{V}({p})+u_{*}$ the
reconstructed pattern starting from the real one
$\hat{u}=(u_{1}(\hat{t}),\dots,u_{n}(\hat{t}))^{\top}$ for the fixed time
$\hat{t}$, the absolute error is $\|\hat{u}-\tilde{u}\|_{1}/n$, i.e., the
$1$-norm of the $n$-dimensional vectors divided by the network size. The
latter error is then weighted by the number of (generalized) eigenvectors used
in the reconstruction, i.e., the dimension $d$ of the subspace $V$, divided by
the total number of (generalized) unstable available eigenvectors, $N_{u}$. In
conclusion the error used to evaluate the goodness of the reconstructed
pattern is

$\varepsilon=\frac{d}{N_{u}}\frac{\|\hat{u}-\tilde{u}\|_{1}}{n}\,.$ (14)

##  Appendix C Generating directed networks with prescribed defective
Laplacian spectra

The goal of this section is to present a novel algorithm allowing to determine
a directed network with a prescribed defective Laplacian spectrum. To the best
of our knowledge this is an open problem and in the literature few results are
available, the interested reader can consult Forrow  _et al._ (2018) for the
case of symmetric networks or Nicoletti  _et al._ (2021) in the case of
directed ones, and the references therein. Let us observe that none of the
previous works can be applied to the present case because both assume the
existence of a basis for the Laplace matrix. For the scope of this work we
thus decided to develop the algorithm under the simplifying assumption of a
real non-positive spectrum.

Let us thus consider a collection of $s+1\geq 2$ real eigenvalues,
$\Lambda^{(0)}=0$, $\Lambda^{(j)}<0$ and let us assume moreover that each
eigenvalue has algebraic multiplicity $m_{j}\geq 1$, for $j=1,\dots,s$,
strictly larger than the geometric multiplicity. Let us initially impose the
null eigenvalue to be simple, i.e., $m_{0}=1$; we will relax this assumption
in the following. We thus trivially have $\sum_{j=1}^{s}m_{j}=n$ and the
network will thus have $n+1$ nodes.

To each eigenvalue $\Lambda^{(j)}$, $j=1,\dots,s$, we associate a Jordan block
of dimension $m_{j}\times m_{j}$

$\mathbf{B}^{(j)}=\Lambda^{(j)}\mathbf{I}_{m_{j}}+\mathbf{N}_{m_{j}}\,,$ (15)

where $\mathbf{I}_{m_{j}}$ is the $m_{j}\times m_{j}$ identity matrix and
$\mathbf{N}_{m_{j}}$ the nilpotent matrix

$\mathbf{N}_{m_{j}}=\left(\begin{matrix}0&1&0&\dots&0\\\ 0&0&1&\dots&0\\\
\vdots&\vdots&\ddots&\ddots&\vdots\\\ 0&0&0&\ddots&1\\\
0&0&0&0&0\end{matrix}\right)\,.$

Let us introduce the $(n+1)\times(n+1)$ block diagonal matrix $\mathbf{J}_{L}$

$\mathbf{J}_{L}=\left(\begin{matrix}0&\vec{0}_{n}^{\top}\\\
\vec{0}_{n}&\mathbf{B}\end{matrix}\right)\,,$ (16)

where $\vec{0}_{n}=(0,\dots,0)^{\top}$ is the $n$-dimensional null vector and
$\mathbf{B}$ the $n\times n$ block diagonal matrix build with the Jordan
block, namely

$\mathbf{B}=\mathrm{diag}(\mathbf{B}^{(1)},\dots,\mathbf{B}^{(s)})\,.$ (17)

Let us observe that we do not consider in $\mathbf{B}$ the null eigenvalue
that has been already set in the entry $(1,1)$ of $\mathbf{J}_{L}$. The latter
matrix will be the Jordan Canonical Form of the Laplace matrix we are looking
for.

Let us define the $(n+1)\times(n+1)$ non-singular matrix $\mathbf{S}$ given by

$\mathbf{S}=\left(\begin{matrix}1&\vec{0}_{n}^{\top}\\\
\mathbf{u}_{n}&\mathbf{I}_{n}\end{matrix}\right)\,,$ (18)

where we have introduced the $n$-dimensional vector
$\mathbf{u}_{n}=(1,\dots,1)^{\top}$. One can easily prove that the inverse of
$\mathbf{S}$ exists and it is given by

$\mathbf{S}^{-1}=\left(\begin{matrix}1&\vec{0}_{n}^{\top}\\\
-\mathbf{u}_{n}&\mathbf{I}_{n}\end{matrix}\right)\,.$ (19)

By using $\mathbf{S}$ and $\mathbf{S}^{-1}$, we define the matrix

$\mathbf{L}=\mathbf{S}\mathbf{J}_{L}\mathbf{S}^{-1}\,.$ (20)

Clearly this matrix has the same spectrum and same eigenvalues multiplicity of
$\mathbf{J}_{L}$. It remains to prove that $\mathbf{L}$ is indeed a Laplace
matrix of a suitable network whose adjacency matrix $\mathbf{A}$ is given by
$A_{ij}=L_{ij}$ for $i\neq j$ and $A_{ii}=0$. In this way we will have built a
network with a prescribed spectrum and eigenvalues multiplicities, and
positively weighted adjacency matrix.

Let us introduce the $(n+1)$-dimensional vector
$\mathbf{e}_{1}=(1,0,\dots,0)^{\top}$, then we trivially get

$\mathbf{S}\mathbf{e}_{1}=\mathbf{u}_{n+1}\,,$

and thus

$\mathbf{S}^{-1}\mathbf{u}_{n+1}=\mathbf{e}_{1}\,,$

hence

$\mathbf{L}\mathbf{u}_{n+1}=\mathbf{S}\mathbf{J}_{L}\mathbf{S}^{-1}\mathbf{u}_{n+1}=\mathbf{S}\mathbf{J}_{L}\mathbf{e}_{1}=0\,,$
(21)

where in the last step we have used the block structure of the matrix
$\mathbf{J}_{L}$ given by (16). We have thus proved that the rows of
$\mathbf{L}$ sum to zero or equivalently that the constant vector
$\mathbf{u}_{n+1}$ is an eigenvector associated to the eigenvalue
$\Lambda^{(0)}=0$.

It remains to prove that $\mathbf{L}$ has non-positive diagonal and non-
negative out of diagonal entries. To do this, let us compute explicitly
$\mathbf{S}\mathbf{J}_{L}\mathbf{S}^{-1}$, moreover let us rewrite
$\mathbf{S}$ and $\mathbf{S}^{-1}$ by using the block structure induced by the
one of $\mathbf{B}$, namely

$\mathbf{L}=\mathbf{S}\mathbf{J}_{L}\mathbf{S}^{-1}=\left(\begin{matrix}1&\vec{0}_{m_{1}}^{\top}&\vec{0}_{m_{2}}^{\top}&\dots&\vec{0}_{m_{s}}^{\top}\\\
\mathbf{u}_{m_{1}}&\mathbf{I}_{m_{1}}&\mathbf{O}_{m_{1},m_{2}}&\dots&\mathbf{O}_{m_{1},m_{s}}\\\
\mathbf{u}_{m_{2}}&\mathbf{O}_{m_{2},m_{1}}&\mathbf{I}_{m_{2}}&\dots&\mathbf{O}_{m_{2},m_{s}}\\\
\vdots&\vdots&\ddots&\ddots&\vdots\\\
\mathbf{u}_{m_{s}}&\mathbf{O}_{m_{s},m_{1}}&\dots&\mathbf{O}_{m_{s},m_{s-1}}&\mathbf{I}_{m_{s}}\end{matrix}\right)\left(\begin{matrix}0&\vec{0}_{N}^{\top}\\\
\vec{0}_{N}&\mathbf{B}\end{matrix}\right)\left(\begin{matrix}1&\vec{0}_{m_{1}}^{\top}&\vec{0}_{m_{2}}^{\top}&\dots&\vec{0}_{m_{s}}^{\top}\\\
-\mathbf{u}_{m_{1}}&\mathbf{I}_{m_{1}}&\mathbf{O}_{m_{1},m_{2}}&\dots&\mathbf{O}_{m_{1},m_{s}}\\\
-\mathbf{u}_{m_{2}}&\mathbf{O}_{m_{2},m_{1}}&\mathbf{I}_{m_{2}}&\dots&\mathbf{O}_{m_{2},m_{s}}\\\
\vdots&\vdots&\ddots&\ddots&\vdots\\\
-\mathbf{u}_{m_{s}}&\mathbf{O}_{m_{s},m_{1}}&\dots&\mathbf{O}_{m_{s},m_{s-1}}&\mathbf{I}_{m_{s}}\end{matrix}\right)\,,$

where we have introduced the null $m\times m^{\prime}$ matrix denoted by
$\mathbf{O}_{m,m^{\prime}}$.

By performing the computation, we obtain

$\displaystyle\mathbf{L}$ $\displaystyle=$
$\displaystyle\left(\begin{matrix}0&\vec{0}_{m_{1}}^{\top}&\vec{0}_{m_{2}}^{\top}&\dots&\vec{0}_{m_{s}}^{\top}\\\
\vec{0}_{m_{1}}&\mathbf{B}_{m_{1}}&\mathbf{O}_{m_{1},m_{2}}&\dots&\mathbf{O}_{m_{1},m_{s}}\\\
\vec{0}_{m_{2}}&\mathbf{O}_{m_{2},m_{1}}&\mathbf{B}_{m_{2}}&\dots&\mathbf{O}_{m_{2},m_{s}}\\\
\vdots&\vdots&\ddots&\ddots&\vdots\\\
\vec{0}_{m_{s}}&\mathbf{O}_{m_{s},m_{1}}&\dots&\mathbf{O}_{m_{s},m_{s-1}}&\mathbf{B}_{m_{s}}\end{matrix}\right)\left(\begin{matrix}1&\vec{0}_{m_{1}}^{\top}&\vec{0}_{m_{2}}^{\top}&\dots&\vec{0}_{m_{s}}^{\top}\\\
-\mathbf{u}_{m_{1}}&\mathbf{I}_{m_{1}}&\mathbf{O}_{m_{1},m_{2}}&\dots&\mathbf{O}_{m_{1},m_{s}}\\\
-\mathbf{u}_{m_{2}}&\mathbf{O}_{m_{2},m_{1}}&\mathbf{I}_{m_{2}}&\dots&\mathbf{O}_{m_{2},m_{s}}\\\
\vdots&\vdots&\ddots&\ddots&\vdots\\\
-\mathbf{u}_{m_{s}}&\mathbf{O}_{m_{s},m_{1}}&\dots&\mathbf{O}_{m_{s},m_{s-1}}&\mathbf{I}_{m_{s}}\end{matrix}\right)$
$\displaystyle=$
$\displaystyle\left(\begin{matrix}0&\vec{0}_{m_{1}}^{\top}&\vec{0}_{m_{2}}^{\top}&\dots&\vec{0}_{m_{s}}^{\top}\\\
-\mathbf{B}_{m_{1}}\mathbf{u}_{m_{1}}&\mathbf{B}_{m_{1}}&\mathbf{O}_{m_{1},m_{2}}&\dots&\mathbf{O}_{m_{1},m_{s}}\\\
-\mathbf{B}_{m_{2}}\mathbf{u}_{m_{2}}&\mathbf{O}_{m_{2},m_{1}}&\mathbf{B}_{m_{2}}&\dots&\mathbf{O}_{m_{2},m_{s}}\\\
\vdots&\vdots&\ddots&\ddots&\vdots\\\
-\mathbf{B}_{m_{s}}\mathbf{u}_{m_{s}}&\mathbf{O}_{m_{s},m_{1}}&\dots&\mathbf{O}_{m_{s},m_{s-1}}&\mathbf{B}_{m_{s}}\end{matrix}\right)\,.$

The diagonal elements are clearly non-positive being $0$ and the diagonal
elements of the matrices $\mathbf{B}_{m_{j}}$, $j=1,\dots,n$, namely

$\mathrm{diag}(\mathbf{L})=(0,\Lambda^{(1)},\dots,\Lambda^{(1)},\dots,\Lambda^{(s)},\dots,\Lambda^{(s)})\,.$
(22)

Let us now consider the out of diagonal elements. It is straightforward to
realize that for any $i\geq 2$ we have $L_{ij}\geq 0$, indeed those terms are
associated to the null matrices $\mathbf{O}_{m,m^{\prime}}$ or to the
nilpotent ones, $\mathbf{N}_{m}$, hence those entries are $0$ or $1$. It
remains to check the sign of the first column, $L_{1j}$. Those elements are of
the form $-\mathbf{B}_{m_{k}}\mathbf{u}_{m_{k}}$ for some $k=1,\dots,s$.
Because of the definition of the Jordan block (15) we have

$-\mathbf{B}_{m_{k}}\mathbf{u}_{m_{k}}=-\left(\Lambda^{(k)}\mathbf{I}_{m_{k}}+\mathbf{N}_{m_{k}}\right)\mathbf{u}_{m_{k}}=-\Lambda^{(k)}\mathbf{u}_{m_{k}}-\mathbf{N}_{m_{k}}\mathbf{u}_{m_{k}}=(-\Lambda^{(k)}-1,\dots,-\Lambda^{(k)}-1,-\Lambda^{(k)})^{\top}\,,$

that are positive by assuming $\Lambda^{(k)}<-1$.

Let us observe that this last constraint can be relaxed by considering the
nilpotent matrix

$\mathbf{N}_{m_{j}}=\left(\begin{matrix}0&a&0&\dots&0\\\ 0&0&a&\dots&0\\\
\vdots&\vdots&\ddots&\ddots&\vdots\\\ 0&0&0&\ddots&a\\\
0&0&0&0&0\end{matrix}\right)\,,$

for some $a>0$; indeed the condition for non-negativity of the entries becomes
$\Lambda^{(k)}<-a$, and taking $a$ close to $0$ we can relax the latter
constraint.

Let us conclude this section by discussing the initial assumption about
$\Lambda^{(0)}=0$ to be simple. The proposed method allows to create networks
with a prescribed real spectrum and multiplicity, but with one node, say the
number $1$, with zero in-degree and maximal out-degree, that is
$k^{(in)}_{1}=0$ and $k^{(out)}_{1}=n$ (see Fig. 5 for an example with
$m_{0}=1$, $m_{1}=2$, $m_{2}=3$ and $m_{3}=1$). From the latter figure, one
can also observe the interesting structures induced by the multiplicity, for
each $m_{j}>1$ there is a sort of “folding fan” with $m_{j}$ directed “sticks”
pointing from the node number $1$ to other nodes connected among them with a
directed path. For instance one can observe on the right part the “folding
fan” $1\rightarrow 3\rightarrow 2\leftarrow 1$ associated to $m_{1}=2$, on the
bottom part the “folding fan” $1\rightarrow 6\rightarrow 5,1\rightarrow
5\rightarrow 4\leftarrow 1$ associated to $m_{2}=3$. Notice that for $m_{j}=1$
the folding fan “collapses” into a directed link, here $1\rightarrow 7$
associated to $m_{3}=1$. To the best of our knowledge, the relations between
the algebraic multiplicity of the Laplace eigenvalues and the above
topological network motifs is new and deserves to be further investigated in
the future.

Figure 5: An example of directed defective network obtained with the previous
algorithm by using the multiplicity : $m_{0}=1$, $m_{1}=2$, $m_{2}=3$ and
$m_{3}=1$.

On the other hand, one can obtain a larger network by “gluing” together
several networks built with the proposed method, seen thus as basic building
blocks. Moreover by carefully choosing where to add the new glue-links we can
create a network with the eigenvalue $0$ with higher algebraic multiplicity
(see Fig. 6). In the left panel of the Figure we show two basic building
blocks, the rightmost one, $\mathcal{N}_{1}$ made by $5$ nodes, with
$(m_{0},m_{1},m_{2},m_{3})=(1,1,2,1)$ and
$\Lambda_{\mathcal{N}_{1}}\sim\left(0,-5.90,-7.25,-8.27\right)^{\top}$, and
the rightmost one, $\mathcal{N}_{2}$ containing $8$ nodes, with
$(m_{0},m_{1},m_{2},m_{3})=(1,1,1,2,3)$ and
$\Lambda_{\mathcal{N}_{2}}\sim\left(0,-1.36,-3.36,-6.74,-9.33\right)^{\top}$.
On the right panel we show the network, $\mathcal{N}$ with $16$ nodes,
obtained by adding the glue-link $7\rightarrow 4$. One can observe that the
latter network contains two nodes, $1$ and $6$, with $k_{i}^{(in)}=0$ and can
thus show that this implies that $\Lambda^{(0)}=0$ has now algebraic
multiplicity $2$, indeed we can compute
$\Lambda_{\mathcal{N}}\sim\left(0,-1.36,-3.36,-6.74,-6.90,-7.25,-8.27,-9.33\right)^{\top}$
with multiplicity $m_{\mathcal{N}}=\left(2,1,1,2,1,2,1,3\right)^{\top}$. Let
us finally observe that the spectrum of $\Lambda_{\mathcal{N}}$ is almost the
union of the spectra of $\Lambda_{\mathcal{N}_{1}}$ and
$\Lambda_{\mathcal{N}_{2}}$, the only difference being the eigenvalue $-6.90$
replacing $-5.90$, this is because the added glue-link modified the weighted
in-degree of node $4$ by adding a new weight $1$.

Figure 6: Creating a larger network by gluing basic building blocks. The
obtained network exhibits a $0$ eigenvalue with algebraic multiplicity $2$.

Let us observe that by using the idea of gluing together basic building
blocks, we can obtain networks that preserve the degenerate eigenvalues, their
multiplicity and only slightly change the remaining eigenvalues, hence without
modifying the multiplicity of $0$, an example is shown in Fig. 7. On the left
panel we present two basic building blocks similar to the ones used in the
construction presented in Fig. 6, i.e., same multiplicity but different
eigenvalues, while on the right part we show the network obtained by adding
the glue-link $7\rightarrow 1$, say $\tilde{\mathcal{N}}$.

The new network has a Laplace matrix whose spectrum is the union of the
spectra of the Laplace matrices for the two smaller networks together with a
new eigenvalue $-1$, in particular the degenerate eigenvalues did not change
their values neither their multiplicity. The reason being that the glue-link
$7\rightarrow 1$ will not modify the in-degree of node $7$, or of any other
node, but the one of node $1$, the latter initially was $0$ and now becomes
$1$. Algebraically the adjacency matrix of $\tilde{\mathcal{N}}$ is given by

$\tilde{\mathbf{A}}=\left(\begin{matrix}\mathbf{A}^{(1)}&\mathbf{E}^{(17)}\\\
\mathbf{O}_{5,8}&\mathbf{A}^{(2)}\end{matrix}\right)$

where the matrix $\mathbf{E}^{(17)}$ has all zero entries but
$E^{(17)}_{17}=1$, and $\mathbf{A}^{(j)}$ denotes the adjacency matrix of the
network $\tilde{\mathcal{N}}_{j}$, $j=1,2$. Hence the Laplace matrix of
$\tilde{\mathcal{N}}$ results to be

$\tilde{\mathbf{L}}=\left(\begin{matrix}\hat{\mathbf{L}}^{(1)}&\mathbf{E}^{(17)}\\\
\mathbf{O}_{5,8}&\tilde{\mathbf{L}}^{(2)}\end{matrix}\right)$

where $\hat{\mathbf{L}}^{(1)}=\tilde{\mathbf{L}}^{(1)}-\mathbf{E}^{(11)}$,
namely there is a “$-1$” in position $11$, while $\tilde{L}^{(1)}_{11}=0$
(node $1$ has in-degree $0$ in the network $\mathcal{\mathcal{N}}_{1}$).
Because the spectra of $\hat{\mathbf{L}}^{(1)}$ and $\tilde{\mathbf{L}}^{(2)}$
are disjoint, the Jordan Canonical Form of $\tilde{\mathbf{L}}$ is the direct
sum of the Jordan Canonical Form of $\hat{\mathbf{L}}^{(1)}$ and
$\tilde{\mathbf{L}}^{(2)}$, from which the conclusion about the spectra and
the multiplicity easily follows.

Figure 7: Creating a larger network by gluing basic building blocks. The
spectrum of the obtained network exhibits the same degenerate eigenvalues,
multiplicity and only slightly change the remaining eigenvalues with respect
to the spectra of the basic building blocks.

##  Appendix D Robustness of pattern reconstruction

The aim of this Section is to present some results supporting the claim that
the proposed method for pattern reconstruction based on the use of generalized
eigenvectors, is robust with respect to the number of involved (generalized)
unstable eigenvectors and the time at which the pattern is considered. In a
second moment we will also study the impact of the network size on the
reconstruction error.

Let us first consider a case where there are two unstable eigenvectors and two
unstable generalized eigenvectors. In the left panel of Fig. 8 we show the
region of instability for the Brusselator model defined on a defective network
composed by $N=10$ nodes (see Fig. 9), where $\Lambda^{(1)}=0$ and
$\Lambda^{(2)}=-1$ both have multiplicity $3$ and are stable, thus they do not
intervene in the Turing instability, while $\Lambda^{(3)}=-3$ has multiplicity
$2$ and is unstable as well as the complex eigenvalue
$\Lambda^{(4)}=-1.5+\mathrm{i}0.866$ and its complex conjugated
$\Lambda^{(5)}=-1.5-\mathrm{i}0.866$, both with multiplicity one.

Figure 8: Region of the complex plane associated to Turing instability (panel
a)) and dispersion relation (panel b)) computed for the Brusselator model with
parameters $b=3.92$, $c=3$, $D_{u}=0.2$ and $D_{v}=0.8$. The black dots denote
the eigenvalues of the Laplace matrix, $\Lambda^{(1)}=0$ (multiplicity $3$),
$\Lambda^{(2)}=-1$ (multiplicity $3$), $\Lambda^{(3)}=-3$ (multiplicity $2$),
$\Lambda^{(4)}=-1.5+i0.866$ (multiplicity $1$) and
$\Lambda^{(5)}=-1.5-\mathrm{i}0.866$ (multiplicity $1$).

In the right panel of Fig. 8 we report the dispersion relation computed for
the Brusselator defined on the defective network shown in Fig. 9. The Laplace
eigenvalues are represented by symbols (red triangles) while the continuous
curve is the dispersion relation computed for the $1$-(complex) parameter
family of linear systems introduced above with the matrix $\mathbf{M}_{j}$.
One can observe that $\lambda(\Lambda^{(j)})$ is positive for $\Lambda^{(3)}$
and $\Lambda^{(4)}$ and thus the equilibrium $(u_{*},v_{*})\sim(1,1.3067)$ is
unstable, as confirmed from the results shown in panel a) of Fig. 10 where we
report $u_{i}(t)$ vs. time. Observe that the same conclusion can be obtained
by looking at the left panel of Fig. 8 where we can realize that
$\Lambda^{(3)}$ and $\Lambda^{(4)}$ lie inside the instability region.

Figure 9: Random non-normal defective network composed by $n=10$ nodes, built
by using a directed Erdős-Rényi where the probability to create a
bidirectional link is $0.2$ and the probability to transform it into a
directed one is $0.6$. Each node has been shown by using a color code
corresponding to the concentration of species $u$ at time $\hat{t}=200$ (see
color map).

In the right panel of Fig. 10 we represent the pattern reconstruction at a
given fixed $\hat{t}$ both with only unstable eigenvectors and with
generalized ones. In the first case the reconstruction error is given by
$\varepsilon=0.0623$ while in the second we have $\varepsilon=0.0012$, let us
observe that the former is larger even if the factor accounting for the number
of used vectors is smaller than one; indeed $d/N_{u}=3/4$, because we used the
eigenvector associated to $\Lambda^{(3)}=-3$ and the two real vectors obtained
from the complex eigenvectors associated to $\Lambda^{(4)}=-1.5+i0.866$. We
can then again conclude that the reconstructed pattern obtained by using both
eigenvectors and generalized ones has a smaller error than the pattern
obtained with only the eigenvectors.

Figure 10: Evolution of the concentration of specie $u_{i}$ over time (panel
a)) for the Brusselator model and pattern reconstruction by using only
unstable eigenvectors (panel b)) and unstable eigenvectors and generalized
ones (panel c)). In the former case, the reconstruction error is
$\varepsilon=0.0623$ while in the latter we have $\varepsilon=0.0012$. The
model parameters are given by $b=3.92$, $c=3$, $D_{u}=0.2$ and $D_{v}=0.8$.
The underlying network is the one shown in Fig. 9. Each orbit has been shown
by using a color code corresponding to the concentration of species $u$ at
time $\hat{t}=200$, the same color used in Fig. 9.

In the last example, the concentrations $u_{i}(t)$ oscillate in time, it would
thus be interesting to reconstruct the time average pattern

$\langle u_{i}\rangle=\frac{1}{T}\int_{t_{0}}^{t_{0}+T}u_{i}(s)\,ds\,,$

where $t_{0}>0$ is a sufficiently large lag of time needed to remove a
transient phase in the orbit behavior, and $T$ is a (multiple of the) orbit
period. Indeed we can show that the proposed scheme works equally well if we
want to reconstruct the time average pattern, as shown in Fig. 11, where we
report the time average patterns and their reconstruction by using the
unstable eigenvectors (top panel) and the unstable eigenvectors and
generalized ones for the same networks and model parameters used for Fig. 10.
In the former case the reconstruction error is $\varepsilon=0.0620$ while in
the second one we obtain $\varepsilon=8.5\,10^{-4}$.

Figure 11: Time average pattern $\langle u\rangle$ vs. reconstructed pattern
for the Brusselator model with parameters $b=3.92$, $c=3$, $D_{u}=0.2$ and
$D_{v}=0.8$. Panel a): the reconstruction is obtained with only the
eigenvectors (EV). Panel b): the eigenvectors and generalized eigenvectors
(GEV) are used for the reconstruction. In the former case the error is
$\varepsilon=0.0620$ while in the latter we get $\varepsilon=8.5\,10^{-4}$.

To conclude this section let us consider a larger network, shown in Fig. 12,
made by $20$ nodes and exhibiting three unstable eigenvalues,
$\Lambda^{(4)}=-4$ with multiplicity four, $\Lambda^{(5)}=-5$ and
$\Lambda^{(6)}=-7$ each with multiplicity one, and three stable ones,
$\Lambda^{(3)}=-2$ with multiplicity seven, $\Lambda^{(2)}=-1$ with
multiplicity three and $\Lambda^{(1)}=0$ with multiplicity four.

Figure 12: Random non-normal defective network composed by $n=20$ nodes, built
by using a directed Erdős-Rényi where the probability to create a
bidirectional link is $0.2$ and the probability to transform it into a
directed one is $0.6$. Each node has been shown by using a color code
corresponding to the concentration of species $u$ at time $\hat{t}=200$ (see
color map).

The region of instability shown in panel a) of Fig. 13 or equivalently the
dispersion relation, panel b) of the same figure, testify the existence of
three unstable modes and indeed one can observe the emergence of patterns (see
panel a) of Fig. 14). Looking at the dispersion relation one can observe that
the latter assumes values very close once evaluated at $\Lambda^{(4)}=-4$ and
$\Lambda^{(5)}=-5$, indeed $\lambda(-4)=0.315$ and $\lambda(-5)=0.312$.

Figure 13: Region of the complex plane associated to Turing instability (panel
a)) and dispersion relation (panel b)) computed for the Brusselator model with
parameters $b=3.92$, $c=3$, $D_{u}=0.2$ and $D_{v}=0.8$. The black dots denote
the eigenvalues of the Laplace matrix, $\Lambda^{(1)}=0$ (multiplicity $4$),
$\Lambda^{(2)}=-1$ (multiplicity $3$), $\Lambda^{(3)}=-2$ (multiplicity $7$),
$\Lambda^{(4)}=-4$ (multiplicity $4$), $\Lambda^{(5)}=-5$ and
$\Lambda^{(6)}=-7$ each with multiplicity one. The used network is shown in
Fig. 12.

The most unstable mode drives the onset of the instability, however it is not
clear a priori if the pattern would be better reconstructed by using the most
unstable mode alone, $\Lambda^{(4)}=-4$, together with the associated
generalized eigenvectors or the mode $\Lambda^{(5)}=-5$. By eyeball analysis
of panels b)-e) of Fig. 14 one would conclude that the strategy relying on the
use of $\Lambda^{(4)}=-4$ and the generalized eigenvectors provides the better
results. Let us however observe that the smallest reconstruction error is
found by using $\Lambda^{(4)}=-4$ alone; this is due to the dimensionality
rescaling factor, here $1/4$, corresponding thus to the use of one eigenvector
out of four possible (generalized) eigenvectors. This analysis should thus be
considered as a preliminary step toward a deeper understanding of the problem,
that will be addressed in a forthcoming study.

Figure 14: Evolution of the concentration of specie $u_{i}$ over time (panel
a)) for the Brusselator model and pattern reconstruction by using only the
unstable eigenvector associated to $\lambda^{(4)}=-4$ (panel b)) and the
unstable eigenvector associated to $\lambda^{(5)}=-5$ (panel d)). In panel c)
we report the pattern reconstructed by using the unstable eigenvector for
$\lambda^{(4)}=-4$ and its associated generalized eigenvectors. Finally in
panel e) we report the patterns obtained by using the unstable eigenvectors
associated to $\lambda^{(4)}=-4$ and $\lambda^{(5)}=-5$ and the generalized
eigenvectors associated to the former vector. The model parameters are given
by $b=3.92$, $c=3$, $D_{u}=0.2$ and $D_{v}=0.8$. The underlying network is the
one show in Fig. 12. Orbits have been reported by using the same color of the
corresponding node in Fig. 12, which corresponds to the value of $u$ at
$\hat{t}=200$.

We conclude this section by studying the impact of the network size on the
error reconstruction. We will use the directed defective networks obtained by
using the algorithm presented in Section C as support for the Brusselator
model to perform this analysis. Moreover, to simplify the analysis and
removing for possible confounding factors, we assumed the existence of a
unique unstable degenerate eigenvalue while all the remaining stable
eigenvalues can be degenerate or not. We then fixed some generic values for
the Brusselator model and computed the dispersion relation, namely the largest
real part of the eigenvalues of the $1$-(complex) parameter family of the
matrix $\mathbf{J}_{0}+\zeta\mathbf{D}$, as defined in the main text (see Fig.
15).

Figure 15: Dispersion relation. The largest real part of the spectrum of the
matrix $\mathbf{J}_{0}+\zeta\mathbf{D}$ is shown in blue as a function of
$\zeta$, the dispersion relation evaluated on stable Laplace eigenvalues is
reported by using red triangles, while the unique unstable eigenvalue is shown
with a green square.

The size of the network and the multiplicity are related by
$n=\sum_{j=1}^{s}m_{j}$. A first study concerns thus the dependence of the
pattern reconstruction error (14) as a function of the multiplicity, $m$, of
the unique unstable eigenvalues for a fixed network size. The results reported
in Fig. 16 concern a network whose size is $n=1000$ and we let the
multiplicity $m$ to vary from $2$ to $20$; for each value of $m$ we compute
the pattern reconstruction error by using the unique eigenvector (EV - blue
points) and the latter eigenvector together with the $m-1$ generalized
eigenvectors (GEV - black squares), let us stress that in the former case the
pre-factor ${d}/{N_{u}}$ in Eq. (14) is equal to $1/m$, while it is the unity
in the case of generalized eigenvectors. Each point is the average over $100$
independent replicas of the construction, i.e., different networks with a
different spectrum but with the same multiplicity. Two main messages can be
drawn from these results: first of all, the use of the generalized unstable
eigenvectors provides an error or several order of magnitude smaller than the
use of the unstable eigenvector. Second, the dependence of the error on the
value $m$ is relatively small.

Figure 16: The pattern reconstruction error as a function of the multiplicity
$m$ of the unique unstable eigenvalues for a network of size $n=1000$.

Because the multiplicity of the unstable eigenvalues does not play a relevant
role, we decided to study the error as a function of the network size for
randomly generated networks of increasing sizes obtained by applying the
gluing construction presented in Section C, with the constraint of having a
unique unstable eigenvalue, whose geometric multiplicity is a random number
uniformly drawn in $\\{2,3,4,5\\}$. We show in Fig. 17 the results for
networks whose sizes ranges from $100$ to $1000$ and we can conclude that the
error is always smaller in the case of the generalized eigenvectors (GEV (dir)
- black squares) are used to reconstruct the pattern with respect to the case
where only the unstable eigenvector is used (EV (dir) - blue circles). In the
same Figure we report the error computed by using a symmetric network with a
single unstable eigenvector built by using the algorithm presented in Forrow
_et al._ (2018) (EV (symm) - red diamonds) and we can observe that the error
lies in between the previous two cases.

Figure 17: The pattern reconstruction error as a function of the network size.
In the case of symmetric networks (EV symmetric, red diamonds) only the
unstable eigenvector has been used to reconstruct the pattern; in the case of
directed defective networks we can use the unstable eigenvector alone (EV dir,
blue circles) or together with the generalized eigenvectors (GEV dir, black
squares). Each symbols is the result of the average over $100$ independent
networks reconstructions.

We conclude this analysis with some cautionary remarks: the family of directed
defective networks we obtain with the above presented algorithm could not be
the most general one, in particular because we allow only for real spectra and
well localized eigenvectors as in the algorithm proposed in Forrow  _et al._
(2018), and those facts could induce an unwanted and uncontrolled bias.