Finite propagation enhances Turing patterns in reaction-diffusion networked systems

We hereby develop the theory of Turing instability for reaction-diffusion systems defined on complex networks assuming finite propagation. Extending to networked systems the framework introduced by Cattaneo in the 40's, we remove the unphysical assumption of infinite propagation velocity holding for reaction-diffusion systems, thus allowing to propose a novel view on the fine tuning issue and on existing experiments. We analytically prove that Turing instability, stationary or wave-like, emerges for a much broader set of conditions, e.g., once the activator diffuses faster than the inhibitor or even in the case of inhibitor-inhibitor systems, overcoming thus the classical Turing framework. Analytical results are compared to direct simulations made on the FitzHugh-Nagumo model, extended to the relativistic reaction-diffusion framework with a complex network as substrate for the dynamics.


Introduction
A blossoming of regular spatio-temporal patterns can be observed in Nature. These are the signature of selforganised processes where ordered structures emerge from disordered ones [1,2]. Very often, the interaction among the microscopic units, by which the system is made of, can be modelled by means of reaction-diffusion equations that describe the deterministic evolution of the concentrations both in time and space, the latter being a regular substrate [2] or a discrete one, e.g., a complex network [3]. Spatially homogeneous equilibria of a reaction-diffusion system may undergo a symmetry breaking instability, when subjected to a heterogeneous perturbation, eventually driving the system towards a patchy, i.e., spatially heterogeneous, stationary or oscillatory solution, as firstly explained by Turing [4] in the 50s and corroborated experimentally almost four decades later [5][6][7].
Nowadays, applications of the Turing instability phenomenon go well beyond the original framework of the morphogenesis or chemical reaction systems and it stands for a pillar to explain self-organisation in Nature [8][9][10]. The conditions for the emergence of Turing patterns have been elegantly grounded on the interplay between slow diffusing activators and fast diffusing inhibitors [11]; indeed this determines a local feedback, short range production of a given species, which should be, at the same time, inhibited at long ranges. Starting from these premises, scholars have been able to extend the original Turing mechanism to non-autonomous systems, e.g., evolving domains [12,13] or time dependent diffusion and reaction rates [14], as well as discrete substrates, e.g., lattices [15] or complex networks [16], and their generalisation, e.g., directed networks [17], multiplex networks [18] and recently to time varying networks [19,20]. The interested reader can consult the recent review [21] for a modern perspective on Turing instability.
As previously observed, at the root of Turing instability there is a reaction-diffusion process which is thus grounded on a (nonlinear) 'heat equation', namely a parabolic partial differential equation (PDE) with a nonlinear source term. The latter PDE is characterised by an infinite fast propagation of the initial datum along the supporting medium and thus it can accurately model the physical phenomenon only in cases of very large diffusivity, D 1. To overcome this drawback, scholars have considered more realistic frameworks. In particular Cattaneo proposed in the 1948 to modify the constitutive equation (Fick's first law) by including a relaxation term with some given characteristic inertial time, τ > 0. Operating in this framework, Fick's second law returns a modified diffusion equation allowing also for a second derivative with respect to time [22][23][24][25]. The resulting equation is nowadays known in the literature as the Cattaneo equation as well as the telegraph equation, the damped nonlinear Klein-Gordon equations or the relativistic heat equation, depending on the research field and on the feature one is interested to emphasise [26]. In any case its main characteristic is to exhibit a finite propagation velocity, v = D/τ , and moreover in the limit of arbitrarily small relaxation time, τ → 0, one recovers Fick's second law and thus a parabolic reaction-diffusion model with an infinite propagation velocity. Our focus being the role of the finite propagation, we will hereby name such framework finite-velocity as done in [27] or sometimes relativistic heat equation although no Lorentz phenomena are at play, to recall the existence of a maximal allowed velocity as for the speed of light in relativity theory; we thus operated a different choice with respect to [28,29], where the name 'hyperbolic reaction-diffusion equations' has been preferred. The aim of this paper is to study the conditions for the onset of Turing instability, being them stationary or oscillatory patterns, for a reaction-diffusion system defined on top of a complex network and modified according to the Cattaneo recipe, to allow for a finite propagation velocity (section 1). We thus consider two different species populating a network composed of n nodes. When species happen to share the same node, they interact via nonlinear functions f(u i , v i ) and g(u i , v i ), being u i and v i the species concentrations at the ith node. On the other hand, they can diffuse across the available network links. The local currents, i.e., associated to each links, are assumed to satisfy a modified constitutive equation, Fick's first law, that includes a relaxation term with a given inertial time. Hence, the continuity equation, Fick's second law, allows to derive a modified local diffusion term, i.e., defined on the node. The latter, together with the reaction part, determines the hyperbolic reaction-diffusion system defined on top of a complex network, we will hereby be interested in.
Our work extends to the network case the study presented in [28], realised under the simplifying hypothesis of equal inertial times for the two species, namely τ u = τ v , and assuming a continuous substrate. Indeed, we hereby assume generic inertial times for each species, τ u = τ v . Let observe that our results established for a discrete substrate, can be straightforwardly extended to the continuous case and thus they complete the work done in [28] to allow for different inertial times. This setting has been recently studied in [29] in the framework of hyperbolic reaction-diffusion models with cross-diffusion defined on a continuous substrate. Differently from our approach of directly adapting the Cattaneo idea to the current flowing on each network link, the authors of [29] have used extended thermodynamics [30]; the resulting characteristic polynomial (see below) we obtained is different from the one given in [29] and it allows us to draw several interesting conclusions.
The Turing mechanism relies on the assumption of the existence of a stable homogeneous equilibrium that loses its stability once subjected to spatially heterogeneous perturbations, in presence of a diffusive term; for this reason, such process is also known as diffusion-driven instability. The system can then exhibit stationary spatially heterogeneous solutions as well as time oscillatory ones. The same mechanism can be proved to hold true in the new proposed framework of the hyperbolic reaction-diffusion systems defined on top of a complex network (section 1). The dispersion relation, which ultimately signals the onset of the instability, is a function of the discrete spectrum of the Laplace matrix, namely the diffusion operator associated to the underlying network. Let us mention that cross-diffusion is excluded. The dispersion relation is obtained from the roots of the fourth order characteristic polynomial. To progress with the analytical understanding of the problem, we resort to the Routh-Hurwitz stability criterion [31][32][33], allowing to prove the (in)stability feature of a real coefficients polynomial. Let us observe that this criterion is a widely used tool in dynamical systems and control theory (see, e.g., [25]).
We have shown that the use of the inertial times strongly enlarges the parameter region for which Turing instability and Turing-waves can emerge, even beyond the classical Turing conditions of fast inhibitor and slow activator. For generic values of the inertial times, τ u = τ v , we have proved that Turing patterns can set up with a fast activator and slow inhibitor, with species exhibiting the same diffusion coefficients and even with an inhibitor-inhibitor system. As in these cases classical Turing instability cannot develop and being the latter solely due to the presence of the inertial times, we propose to call it inertia-driven instability, that can result into stationary or wave-like phenomena. Of course, the proposed framework allows to prove the existence of Turing instability also for an inhibitor diffusing faster than the activator, as for the non-relativistic framework.
In the particular case where both species have the same inertial time, we have shown that the stability of the homogeneous solution is conditional to the inertial time; indeed there exists a threshold, τ max , beyond which the homogeneous equilibrium turns out to be unstable. The system exhibits thus patterns but they cannot be associated to Turing instability, even if they are indistinguishable from the latter. Moreover, the threshold τ max depends on the model parameters and there are combinations of the latter for which it is arbitrary large; stated differently, for such parameters the homogeneous equilibrium is always stable (with respect to the inertial time).
The theoretical framework hereby proposed has been complemented with a dedicated numerical analysis of the FitzHugh-Nagumo model [34][35][36] (see section 4), that is a nonlinear system often used as paradigm for the study of the emergence of Turing patterns [37][38][39][40] as well as for synchronisation phenomena [41,42]. The FitzHugh-Nagumo model has thus been extended to the framework of hyperbolic reaction-diffusion networked systems. We have numerically found stationary patterns as well as synchronised oscillatory ones. We have also found new interesting solutions for which the dispersion relation has limited predictive power; indeed we showed the existence of two sets of model parameters associated to similar dispersion relations, for which the unstable modes have a nontrivial complex component, but in one case the solution oscillates in time whereas in the second it converges to a stationary pattern.
The proposed framework tackles thus the issue of infinite propagation velocity for networked reaction-diffusion systems, and it is general enough to account for novel interesting results, strengthening the importance of self-organisation in nonlinear networked system. In particular the possibility to prove the emergence of Turing patterns, being the latter stationary or wave-like, in the case of activator diffusing faster than the inhibitor but also in the case of inhibitor-inhibitor systems, could provide new insights into the fine tuning problem [21,43] and propose a novel view on experimental results.

Reaction-diffusion system with finite propagation on networks
The aim of this section is to extend Cattaneo's idea to a discrete substrate, i.e. to deal with a networked system. We will briefly show how to modify networked reaction-diffusion systems in order to allow for a finite velocity of propagation.
Let us thus consider a network made of n nodes and connected by a collection of m undirected links allowing form pairwise exchanges among nodes. Such structure can be encoded into the m × n incidence matrix, M. Let e = (i, j) be the link connecting nodes i and j, then M ei = 1, M ej = −1 and M e = 0 for all = i, j. From this matrix we can build the Laplace matrix, L = −M M, where denotes the matrix transpose. The Laplace matrix is symmetric by construction and thus it admits a set of orthonormal eigenvectors, φ (α) , and real non-positive 1 eigenvalues Λ (α) , for α = 1, . . . , n. By construction j L ij = 0, hence the largest eigenvalue is Λ (1) = 0 associated to the eigenvector φ (1) = (1, . . . , 1) / √ n. The diagonal element −L ii defines the nodes degree, say k i , namely the number of incidents links of the ith node; hence we can rewrite L = A − D, where D = diag(k 1 , . . . , k n ) and A is the adjacency matrix, that is A ij = 1 if and only if nodes i and j are connected, encoding thus the coupling network.
Let us now focus on the diffusion of a single species in the network, the generalisation to more species being a direct extension. Let u(t) = (u 1 (t), . . . , u n (t)) to denote the state of the system at time t, where u i (t) is the density of the species in node i at time t. Let e = (i, j) be a link in the network and let χ e (t) be the current flowing through it at time t; then, borrowing the constitutive equation, namely Fick's first law, from the continuous framework, we can state that that is, the current is proportional to the difference of the densities in the nodes forming the link and flowing from higher concentrations to lower ones 2 , being D u the diffusion coefficient of species u. By defining the currents vector χ = (χ e 1 , . . . , χ e m ) , the continuity equation can be written as namely the variation of u i is proportional to the sum of the currents entering and exiting from node i. The classical Fick second law follows by combining the above equations: 1 The matrix L is negative-semidefinite. Indeed let φ (α) be any orthonormal eigenvector, then ( φ (α) , L φ (α) ) = Λ (α) and at the same time 2 Let us observe that, despite the different sign in front of equation (1), the latter is the analogous of the Fick's first law in the continuous setting: the current flows from regions of higher concentration to regions of lower one. Indeed, once we fix the link 'ordering' as e = (i, j), then the current χ e will be positive, i.e., respecting the link ordering if u i > u j , while the current will be negative, i.e., opposite to the link order if u i < u j .
where we can realise [15,16] that L replaces the second order differential operator used in the continuous substrate case and thus the model given by (3) exhibits infinite propagation velocity.
To overcome this problem we modify, as Cattaneo did, the constitutive equation (1) by introducing a relaxation factor with some characteristic inertial time τ u > 0, namely Combining this equation with the continuity equation (2) allows us to obtain eventually providing the generalised Cattaneo equation defined on networks The latter can be seen as a modification of the 'heat equation' defined on network by the inclusion of a second order time derivative, returning thus a relativistic or hyperbolic heat equation. Consider now two different species populating a network composed by n nodes and let us denote by u i and v i , i = 1, . . . , n, their respective concentrations on node i. When species happen to share the same node, they interact via 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 accordingly to the modified Cattaneo equation (6). The model can hence be mathematically cast in the form where D u > 0 (resp. D v > 0) is the diffusion coefficients of species u (reps. v) and τ u > 0 (resp. τ v > 0) the inertial time for species u (resp. v).

Turing instability in networked reaction-diffusion systems with finite propagation
The Turing mechanism is the result of a diffusion-driven instability, namely an homogeneous stable equilibrium of the reaction-diffusion system turns out to be unstable, with respect to inhomogeneous spatial perturbations, once the diffusion is at play. The aim of this section is to determine the conditions for such instability to develop in the relativistic, i.e., in presence of a maximal allowed velocity, reaction-diffusion systems defined on networks given by equation (7). Let us hence assume there exists an homogeneous solution of (7), that is u i (t) = u 0 and v i (t) = v 0 for all i = 1, . . . , n and t > 0. Namely u 0 and v 0 should satisfy f(u 0 , v 0 ) = g(u 0 , v 0 ) = 0. Being the latter equilibrium solely determined by the reaction terms, it happens to be also an equilibrium for the non-relativistic system. Let us denote by δu i (t) = u i (t) − u 0 and δv i (t) = v i (t) − v 0 the perturbations from the homogeneous solution. In order to determine the time evolution of the latter, we use (7), keeping only the first order terms in the perturbation (the latter assumed to be small). We thus obtain where we employed the fact that j L ij = 0 to nullify the terms j L ij u 0 and j L ij v 0 . Let us also stress that throughout the rest of the section the partial derivatives, i.e., ∂ u f ≡ ∂f/∂u and similarly for the other ones, are evaluated at the homogeneous equilibrium (u 0 , v 0 ). To progress with the analytical understanding, we develop the perturbations on the eigenbasis of the Laplace i . Inserting the latter into equation (8) we obtain the equation describing the evolution of the modesû α (t) andv α (t) namely we end up with n linear 2 × 2 systems instead of the initial 2n × 2n one. We further hypothesisê u α (t) ∼ e λ α t andv α (t) ∼ e λ α t , and to ensure the existence of a nontrivial solution we eventually obtain that the linear growth rate λ α should solve where the fourth degree characteristic polynomial is defined by whose coefficients are given by being ⎞ ⎠ the Jacobian of the sole reaction system without diffusive coupling, hence also named aspatial system, evaluated at the homogeneous equilibrium The coefficients a and b are positive and do not depend on the index α.
Turing instability arises if the homogeneous equilibrium (u 0 , v 0 ) is stable, namely if the four roots of the polynomial p 1 (λ) all have negative real part 3 , while there exists at least one α > 1 for which the polynomial p α (λ) does admit at least one root with positive real part. The root with the largest real part seen as a function of Λ (α) , is called in the literature the dispersion relation, λ α := max i=1,...,4 Rλ i (Λ (α) ). Turing instability is thus equivalent to require λ 1 < 0 and λ α > 0 for some α > 1, hereby called critical roots. Indeed, because of the ansatzû α (t) ∼ e λ α t andv α (t) ∼ e λ α t , the former implies an initial exponential divergence from the homogeneous equilibrium. In the following we will also use information about the imaginary part of λ α , we thus define ρ α := max i=1,...,4 { λ i (Λ (α) ) : λ α > 0}, i.e., the largest imaginary part of the critical roots. In the case the imaginary part of the critical root is non-zero, ρ α = 0, we are in presence of a Turing-wave instability, and the perturbation initially exhibits a combination of exponential growth and oscillating behaviour. Eventually the nonlinearities of the model determine the final pattern, that could result to be stationary or wave-like one.
To prove the existence of Turing instability for the system (7) we shall rely on the Routh-Hurwitz criterion [31][32][33], providing necessary and sufficient conditions to prove that p 1 is stable 4 while p α is unstable for some α > 1.

Remark 1 (Connection with the relativistic reaction-diffusion system defined on a continuous substrate).
As already remarked, the Laplace matrix L in equation (7) takes the place of the second order differential operator Δ = i ∂ 2 x i . After linearising the resulting PDE system about the homogeneous equilibrium, the use of the periodic boundary conditions and the Fourier series, is equivalent to project the linear system onto the eigenfunctions of Δ, that is e ikx (for k ∈ Z, in the case of a 1 dimensional spatial domain), whose eigenvalues are −k 2 . Proceeding in this way, one can determine a polynomial similar to the one given in (10), where we have to replace Λ (α) by −k 2 . However, let us observe that now the spectrum of L is discrete and thus the dispersion relation for the networked system will be 'sampled' from the one holding in the continuous case (see red dots on the blue curves in the following figures representing the dispersion relations). This may introduce finite size effects, as the continuous support case is capable to exhibit Turing patterns, while the networked one cannot because the Laplace spectrum has a gap that exactly avoids the region of positive dispersion relation. To control for this phenomenon, one should be able to relate topological features of the network to the Laplace spectrum [44,45].

Conditions for the stability of p 1
The aim of this section is to introduce the conditions for the linear stability of the homogeneous solution of (7). As already noticed, the coefficients a = τ u τ v and b = τ u + τ v are positive, hence the necessary and sufficient conditions (see appendix A) to ensure the stability of p 1 are given by: Before proceeding with the analysis in the general setting, let us consider a special but relevant case, namely τ u = τ v = τ . Assuming equation (16) to hold true, then equations (15) and (18) easily follow. Moreover, if 4 det(J 0 ) < (tr(J 0 )) 2 , then equation (19) is always satisfied, while if 4 det(J 0 ) > (tr(J 0 )) 2 , the following upper bound for τ is obtained to satisfy (19): This last results will be important in the following, because it states that the stability of the homogeneous equilibrium depends on τ (see panel (b) in figure 1). More importantly, if τ is large enough, the system (7) exhibits patterns; they are not emerging from a Turing mechanism but instead from the instability of the homogeneous equilibrium, observe that one cannot discriminate them with respect to Turing patterns by simple visual inspection.

Conditions for the instability of p α
Using again the Routh-Hurwitz criterion we can prove the existence of (at least) an α > 1 for which p α is unstable conditioned on the stability of p 1 . Observe again that the coefficients a and b are positive. Moreover, by assuming equations (15) and (16) to hold true and by recalling that −Λ (α) > 0 for all α > 1, then c α > 0 and d α > 0 (see equations (12) and (13)). In conclusion the unique coefficient of p α that can be negative is e α . Hence (see appendix A) the instability can arise if one of the following couples of conditions is verified: or Let us observe that equations (21a) and (21b) do not depend on τ u and τ v and are indeed the same conditions one imposes to obtain the Turing instability in the classical, i.e., non-relativistic setting [16]. In particular, they require D v > D u . However, equations (22a) and (22b) do not ask for such condition on the diffusivities, implying that the hypothesis of a finite propagation velocity allows to enlarge the parameter region for which Turing instability arises, in particular allowing for D v D u .
Based on the above one can conclude that if the patterns with positive inertial times are due to equations (22a) and (22b), then they persist also in the non-relativistic limit, τ u → 0 and τ v → 0. On the other hand, if the instability has been initiated by conditions equations (21a) and (21b), we can show (see appendix B) that in the non-relativistic limit, the patterns fade out and disappear for positive and sufficiently small inertial times.
To start our analysis, let us thus consider the case D u = D v = D. As already observed equation (22a) cannot be satisfied having by equation (16) the fact that tr(J 0 ) < 0, thus this cannot be a path towards Turing . . , n, as a function of β and μ: the black regions denote stability while white ones instability. Panel (a) corresponds to the classical setting, i.e., τ u = τ v = 0, the remaining panels are associated to positive values of the inertial times, τ u = τ v = 1 (panel (b)), τ u = 5 and τ v = 1 (panel (c)) and τ u = 1 and τ v = 5 (panel (d)). In all the panels the red line denotes the condition tr(J 0 ) = 0, while det(J 0 ) = 0 is represented by the yellow one; these two lines determine the boundary of the stability region in the classical setting. Such region is shrunk in the case of positive inertial times because of the additional constraints, equation (18) (green line) and equation (19) (blue one). The grey shaded region in panel b), coloured according to ln τ max , is associated to a stability of the homogenous equilibrium constrained to a bound on τ , see equation (20), while in the black region any positive value of τ is admissible.
instability. On the other hand, let us reorganise terms and rewrite condition (21a) as follows Indeed, ∂ u f > 0, being u the activator species, and tr(J 0 ) < 0 by the stability assumption on p 1 (λ); hence the term in brackets on the right-hand side is the sum of three positive terms, from which the claim follows. On the contrary, if Finally, condition (21b) no longer depends on D and can thus be verified by a suitable choice of the remaining parameters.
In conclusion, we can have Turing instability also in the case of equal diffusivities, D u = D v , provided the inhibitor has a larger inertial time than the activator, τ v > τ u .
Let us conclude this section by considering again the case τ u = τ v = τ . Because of the previous analysis we have to assume D u = D v , otherwise no Turing instability can develop. Then equation (21a) simplifies into Being v the inhibitor species we have ∂ v g < 0, hence, if D v > D u , we can conclude that the term in brackets on the right-hand side is the sum of positive terms and thus Finally, the remaining condition (21b) can be rewritten as and a straightforward computation allows to show that it is satisfied if .
Let us stress that in this setting, D v < D u , the conditions (21a) and (21b) cannot be satisfied, hence the emergence of Turing instability is solely due to the finite propagation velocity and imposes a lower bound on the inertial time. Before introducing the model we will use to present our results, let us emphasise two more relevant results. First, the proposed framework allows to prove the emergence of Turing instability also in an inhibitor-inhibitor system, that is ∂ u f < 0 and ∂ v g < 0; indeed, while equation (22a) cannot hold true, equations (21a) and (21b) can be satisfied for suitable choice of the parameters, as we will show below (see figure 8 and the associated discussion). Second, inertia-driven Turing instability cannot manifest in suitable m-species linear kinetic models, as described in the following remark. [46] an m-species non-relativistic kinetic system can be expressed as a polynomial differential equation assuming mass-action for the reaction rates; however not all polynomials can be considered models of chemical reactions [47] because of the possible presence of negative cross terms, i.e., the abundance of a species decreases in a process where it is not involved. The absence of negative cross effects in the case of first order kinetic differential systems, prevents the Turing instability [46]. We can prove a similar result to hold true in the relativistic framework, provided all the species have the same inertial time and the ratios of the diffusion coefficients divided by the inertial time are large enough (see appendix C).

The FitzHugh-Nagumo model
The aim of this section is to present an application of the theory hereby developed. For the sake of definitiveness, we decided to use the FitzHugh-Nagumo model [34][35][36], but our results go beyond the chosen model. The FitzHugh-Nagumo model (for short FHN) is a paradigmatic nonlinear system already used in the literature to study the emergence of Turing patterns [37][38][39][40] as well as synchronisation phenomena [41,42]. Let us observe that this model is not a kinetic one, since the −v term appearing in the rate evolution for u, expresses a negative cross-effect [46], at the same time this supports our statement that Turing instability finds applications beyond the morphogenesis and chemical frameworks. Our choice relies also on the observation that such model has been conceived in the framework of neuroscience as a schematisation of an electric impulse propagating through an axon. For this reason, we believe that it would make a suitable setting to account for a finite velocity propagation of signals and it could be interesting for future applications.
The FHN model can be described by the system of ODE where the parameters μ, γ and β are assumed to be positive. We will hereby focus on its behaviour close to the fixed point (u 0 , v 0 ) = (0, 0). The linear stability analysis ensures stability of the latter under the conditions μ < γβ and μβ < 1 (see panel (a) in figure 1). Let us observe that, once such conditions are not met, the system undergoes a supercritical Hopf-Andronov bifurcation [48]: the equilibrium point becomes unstable giving birth to a limit cycle solution. In this study we will limit ourselves to the former case, leaving the oscillating case for a future work. Consider now n identical copies of the FitzHugh-Nagumo model (23) interacting with each other through a diffusive-like coupling and assume to work in the Cattaneo framework presented in section 1. The resulting model can thus be written as where D u (resp. D v ) is the diffusion coefficients of species u (reps. v) and τ u (resp. τ v ) the inertial time for species u (resp. v). The matrix L is the Laplace matrix describing the diffusive coupling among the FitzHugh-Nagumo systems.
Remark 3 (About the network substrate). The possible onset of Turing instability depends both on the dynamical system as well as the network substrate via the eigenvalues Λ (α) of the associated Laplace matrix, L. As previously stated, a discrete topology may affect the dynamics due to finite size effects. Because a comprehensive study of such impact on Turing instability goes beyond the scope of this paper and for the sake of definitiveness, we decided to use an Erdös-Rényi random graph [49] made of n nodes and each couple of nodes having a probability p ∈ (0, 1) to be linked. In the following we fixed n = 30 and p = 0.1 and we also checked that the resulting network is connected.
In the rest of this section we will present our analysis about the emergence of Turing instability in the relativistic FHN defined on networks (24). Let us start by determining the parameter region associated to a stable homogeneous equilibrium, (u 0 , v 0 ) = (0, 0). The panel (a) of figure 1 represents the classical case where we assume an infinite propagation velocity, namely τ u = τ v = 0. The stability region (black) is delimited by the conditions μ < γβ (red line) and μβ < 1 (yellow line). In panel (b) we report the case of equal inertial times, τ u = τ v = 1; we can observe that the stability region (black and shades of grey) is contained in the previous one, being delimited by the same conditions and in addition by equation (19) (blue line). As previously observed, the stability of the homogeneous solution depends on the value of τ u = τ v = τ , meaning that the equilibrium loses its stability if the inertial time is too large, as shown by equation (20). The grey shaded region in panel (b) has thus been coloured according to ln τ max : smaller values are associated to lighter shades of grey. On the contrary, in the black region any positive value of τ returns a stable homogeneous equilibrium (being τ max = ∞). In the remaining panels of figure 1 we considered different inertial times, τ u = 5 and τ v = 1 in panel (c), and τ u = 1 and τ v = 5 in panel (d). The stability region (black) is delimited by the same lines as before, with the exception of the case τ u < τ v , where an extra condition needs to be considered, i.e., equation (18) (green line).
We are now able to study the emergence of Turing instability under the assumption τ u = τ v = τ . In the panel (a) of figure 2 we report the region (black) in the parameter space allowing for classical Turing instability to arise for a choice of the diffusivities D u < D v . Such region is contained in the one associated to a stable homogeneous solution and it is now also bounded by the conditions The same values of the parameters are used in panel (b) assuming now the inertial times to be positive, τ u = τ v = 1; the Turing region (black and shades of grey) is smaller, as it is also delimited by the condition equation (19) (blue line). Once again, the shades of grey represent the values of ln τ max to ensure the stability of the homogeneous solution (see equation (20)). Finally in panel (c) we report the analysis of a setting for which classical Turing instability can never emerge because the inhibitor diffuses slower than the activator, D u = 2.2 > D v = 0.2. The instability being determined by the inertial time τ > 0, we named it inertia-driven instability. The Turing region (black and shades of grey) is now delimited also by the condition equation (21b), where again the shades of grey represent the values of ln τ max .
The impact of τ max can be appreciated in figure 3 where we report for few generic sets of parameters the dispersion relation, λ α , as a function of Λ (α) , the eigenvalues of the Laplace matrix, L. In panel (a) we show the dispersion relation for the choice τ u = τ v = 1 and (β, μ) = (0.8, 1.0), lying in the Turing instability region (yellow star in the panel (c) of figure 2). We can observe that the homogeneous equilibrium is stable (the dispersion relation is negative for Λ (1) = 0) but it turns out to be unstable under heterogeneous perturbations (there exist α > 1 (red dots) for which λ α > 0) and synchronised oscillatory patterns emerge (see inset where  (22b)). Together with the blue line in panel (b) corresponding to equation (19), these lines delimitate the parameter region allowing for Turing instability in the case D u < D v . In panel (c), corresponding to D u > D v , a similar parameter region is bounded by the same blue line but also by the dashed black line, namely equation (21b). The grey shaded region in panels (b) and (c), coloured according to ln τ max , is associated to a stability of the homogeneous equilibrium constrained to a bound on τ , see equation (20), while in the black region any positive value of τ is admissible. . The homogeneous equilibrium is stable (the dispersion relation is negative for Λ (1) = 0), but it turns out to be unstable under heterogeneous perturbations and synchronised oscillatory patterns emerge (inset), indeed the critical root has positive imaginary part, ρ α > 0 (conditions (21a) and (21b) are satisfied). In panel (b) we fix τ u = τ v = 2.2 and (β, μ) = (0.7, 1.0) (red triangle in the panel (c) of figure 2), still in the Turing region but conditioned to the value of τ max . The behaviour is similar to the one reported in panel (a) but now the homogeneous equilibrium is weakly stable, the dispersion relation is negative but very close to 0 for Λ (1) = 0, indeed for these values of the parameters we have τ max ∼ 2.31. Again, an oscillatory behaviour is obtained (inset) associated to ρ α > 0. In panel (c) we used the same parameters (β, μ) but we increased τ u = τ v = 3.5 > τ max and indeed the homogeneous equilibrium is unstable, the dispersion relation is positive for Λ (1) = 0. Again synchronised oscillatory patterns emerge (inset), they are indistinguishable from the ones one could obtain with the parameters used in panels (a) and (b) but they are not the result of Turing instability.
we report u i (t) vs t 5 ), indeed, we are in presence of an oscillatory Turing instability because ρ α > 0 (data not shown). Panel (b) (τ u = τ v = 2.2 and (β, μ) = (0.7, 1.0), red triangle in the panel (c) of figure 2) corresponds to a similar behaviour, being the parameters still in the Turing region but conditioned to the value of τ max ; the dispersion relation assumes positive values but the homogeneous equilibrium is weakly stable, the dispersion 5 Throughout this work the numerical simulations have been performed using a 4th order Runge-Kutta scheme implemented in Matlab [50]; the code is available upon request to the corresponding author. The initial conditions have been realised by drawing uniformly random perturbations δ-close to the homogeneous equilibrium and the simulation time has been taken of the order of −log δ/max α Rλ α . Indeed according to the ansatzû α ∼ e λα t andû α ∼ e λαt , this is the time necessary to (possibly) increase the δ-perturbation up to a macroscopic size. In the rest of the work we set δ = 10 −2 , small enough to discriminate between the onset of the instability using a reasonable simulation time for the values of max α Rλ α we are dealing with.   figure 4) and the conditions (21a) and (21b) are satisfied. In both cases the aspatial equilibrium is stable (λ 1 < 0), but it turns out to be unstable under heterogeneous perturbations and synchronised oscillatory patterns emerge. Being ρ α > 0 we are in presence of a Turing-wave instability driven by the inertial times. relation is negative but very close to 0 for Λ (1) = 0, being τ u = τ v = 2.2 close to τ max ∼ 2.31; again, an oscillating synchronous behaviour is observed (inset). In panel (c) we used the same parameters (β, μ) but we increased the inertial times beyond the critical values, τ u = τ v = 3.5 > τ max , and indeed the homogeneous equilibrium is unstable, the dispersion relation is positive for Λ (1) = 0. Again synchronised oscillatory patterns emerge (inset), they are indistinguishable from the ones one can obtain from the setting presented in panels (a) and (b) but they are not the result of Turing instability.
We can now consider the more general case of different inertial times and show the onset of Turing instability for a choice of the diffusivities that cannot allow for the classical Turing phenomenon, notably because   μ) and (τ u , τ v ) we show the dispersion relation (panels (a) and (d)), λ α as a function Λ (α) , the imaginary part of the root with the largest real part, ρ α (panels (b) and (e)), and the time evolution of the solutions u i (t) (panels (c) and (f)). The top three panels correspond to the choice τ u = 1, τ v = 2, β = 0.9, μ = 1.0, D u = 0.2 and D v = 2.2, associated to Turing instability (see panel (b) figure 6); the conditions (22a) and (22b) hold true. While in the bottom three panels we invert the sizes of the inertial times, τ u = 5 and τ v = 1, and the remaining parameters are (β, μ) = (0.7, 1.15), D u = 0.2 and D v = 2.2, still in the Turing region (see panel (a) figure 6) and the conditions (22a) and (22b) are satisfied. In both cases the aspatial equilibrium is stable (λ 1 < 0), but it turns out to be unstable under heterogeneous perturbations and synchronised oscillatory patterns can emerge. Let us observe that ρ α vanished on an interval containing all the unstable modes, −Λ (α) , and it is positive elsewhere, i.e. in correspondence to decaying modes, the obtained instability possesses thus both the characteristic of a Turing instability and a Turing-wave. Indeed the pattern associated to the first set of parameters (top panels) keeps oscillating after a transient time (panel (c)) while the pattern resulting from the second set of parameters (bottom panels) settle onto a stationary solution (panel (f)). the activator can diffuse faster than the inhibitor. For this reason we hereby stress again that inertia-driven instability should be a suitable name for such phenomena.
In figure 4 we report the region (black) in the parameter space (β, μ) allowing for Turing instability under the assumptions τ u = τ v and D u D v . Such region is contained in the region associated to a stable homogeneous solution (see panels (c) and (d) of figure 1) and delimited in addition by the conditions (18) (green line), (19) (magenta line) and (21b) (dashed black line). We can observe that for all the choices of the inertial times and diffusion constants (τ u = 5, τ v = 1, D u = 2.2 and D v = 0.2 panel (a), τ u = 1, τ v = 5, D u = 2.2 and D v = 0.2 panel (b) and τ u = 1, τ v = 5, D u = D v = 2.2 panel (c)) there are always parameters (β, μ) allowing Turing instability to occur. One can show that in all the presented cases we are dealing with a Turing oscillatory instability, driven by the inertial times in the panels (b) and (c). We report in figure 5 two generic dispersion relations for the latter settings. In both cases we can appreciate the fact that the aspatial solution is stable, indeed λ 1 < 0, while there are α > 1 (red dots) for which λ α > 0, testifying the instability of such equilibrium once subjected to heterogeneous perturbations and resulting into synchronised oscillatory patterns (see panels (c) and (f)) associated with a positive ρ α (see panels (b) and (e)).
The bifurcation regions reported in figure 6 correspond to a parameters setting for which Turing instability could emerge because the inhibitor diffuses faster than the activator, D v > D u , even without the presence of positive inertial times, indeed the conditions (22a) and (22b) are satisfied. However, the resulting patterns and dispersion relations (see figure 7) are quite different in the relativistic case, τ u > 0 and τ v > 0. The top three panels refer to a generic point in the Turing region (see panel (b) figure 6), here γ = 4, β = 0.9, μ = 1.0, τ u = 1, τ v = 2, D u = 0.2 and D v = 2.2; we can clearly appreciate that Turing instability is at play. Indeed, the aspatial equilibrium is stable, λ 1 < 0, and there are modes α > 1 for which the dispersion relation is positive, λ α > 0 (panel (a)); moreover, such unstable modes are real, being their imaginary part zero, ρ α = 0 (panel (b)). One should thus expect the system to settle into stationary patterns, but that is not the case (panel (c)); in fact, the solution departs from the homogeneous equilibrium and it spends a transient time (much longer that the initial period needed to depart from the equilibrium) oscillating with very small amplitudes around different values, only after this phase the amplitudes increase and a wave develops. Observe also that each node oscillates about a different average value, which is not the case in the oscillating patterns shown before. The bottom three panels correspond to a generic point still in the Turing region (see panel (a) figure 6), with γ = 4, β = 0.7, μ = 1.15, τ u = 5, τ v = 1, D u = 0.2 and D v = 2.2 (conditions (22a) and (22b) hold true). The dispersion relation and its imaginary part behave similarly to the previous case, however now the solution diverging from the unstable homogeneous equilibrium settles onto a stationary heterogeneous equilibrium (panel (f)), as one should expect from the classical, i.e., non-relativistic, Turing instability.
To the best of our knowledge, this is a remarkable phenomenon that should be taken into account in the problem of patterns prediction [51]. Indeed, since the seminal paper by A. Turing [4], scholars are aware of the existence of stationary Turing patterns, often associated to a real dispersion relation, and of oscillatory Turing patters resulting from a Turing wave instability. The use of discrete substrates such as networks questioned this dichotomy and a rule of thumb seems to apply [52]: oscillatory patterns develop if the most unstable mode has a large imaginary part, ρ α λ α . The last example goes in the opposite direction because here λ α > ρ α = 0 and the system can exhibit stationary patterns as well as waves, recalling that the final patterns are initiated by the linear behaviour, but rather shaped by the nonlinear character of the system.

Discussion
In this work we have improved the Cattaneo framework of relativistic reaction-diffusion systems to allow for complex network substrates. We have thus analytically studied the conditions for the emergence of Turing instability, stationary or wave-like, for hyperbolic reaction-diffusion networked systems. The introduction of the inertial times removes the unphysical assumption of infinite propagation velocity and, more importantly this new framework allows for Turing patterns to emerge also for parameter values for which classical, i.e., non-relativistic, Turing instability cannot arise, e.g., once the activator diffuses faster than the inhibitor or even in the case of inhibitor-inhibitor systems.
To support the last claim, let us consider a generic quadratic Lotka-Volterra system [53,54] involving two species, namely f(u, v) = u(a 1 where all the parameters are positive numbers, and let us consider its relativistic networked extension: We can observe that the homogeneous equilibrium is stable (λ 1 < 0), but it turns out to be unstable under heterogeneous perturbations (red dots in the panel (a)) and oscillatory patterns can emerge (panel (b)). The imaginary part of the largest roots is positive (data not shown) and we are thus dealing with an inertia-driven wave-instability.
The homogenous nontrivial equilibrium is , and one can easily show that ∂ u f = −b 1 u 0 < 0 and ∂ v g = −b 2 v 0 < 0, provided u 0 > 0 and v 0 > 0 as we will hereby assume; we are hence dealing with an inhibitor-inhibitor system. Contrary to the classical setting where Turing pattern are not allowed for, in the relativistic framework parameters can be chosen in such a way that the above system exhibits an inertia-driven instability resulting in an oscillatory behaviour (see figure 8).
We have shown that the stability of the homogeneous solution is conditional to the inertial time common to both species. There exists a threshold, τ max , that, if exceeded, returns an unstable homogeneous solution: the system exhibits patterns but they are not ascribed to a Turing instability; let us observe that the latter are indistinguishable from the ones emerging following the Turing mechanism. Interestingly enough, such threshold depends on the model parameters and it can become arbitrarily large for a specific range of the latter; in such case, the homogeneous equilibrium is always stable (with respect to the inertial time). For generic values of the inertial times, τ u = τ v , we have proved that Turing instability can set up both for the inhibitor diffusing faster than the activator, D v > D u , as it occurs in the classical setting, but also in the complementary regime, i.e., D v D u , which is forbidden in the absence of inertial time. Even more striking, Turing patterns can emerge also in the case of inhibitor-inhibitor systems. The framework we propose allows to relax the severe parameters conditions for the patterns onset and thus provide new insights into the fine tuning problem [21,43]. Hence, existing experiments could also be read with this novel perspective and analysed in the proposed framework.
We have complemented our general analytical results with a numerical study of the FitzHugh-Nagumo model extended to the framework of hyperbolic reaction-diffusion networked systems. We have found stationary patterns as well as synchronised oscillatory ones; we have also found an interesting class of solutions where the system spends a transient time into a stationary-like regime but then it evolves into an oscillatory one. This example raises relevant questions about the prediction of the patterns following a Turing instability, which is, up to now, an open problem [51].
The investigation discussed in this paper could be further extended in several directions. Previous studies have shown that different kinds of networks, such as directed [17] or non-normal ones [55,56], extend the conditions for the emergence of patterns and allow for a richer spectrum of instabilities. Moreover, it has been shown that an instability similar to the Turing mechanism can be obtained by perturbing a stable limit cycle [57] and that non-normal networks further enhance such instability [58]. Given the oscillatory behaviours of neurones [34,35] and the non-normal nature of neural networks [55], an extension towards such direction would open the way to interesting new results and applications.
As already observed, the coefficients a and b are positive. Moreover, because of the assumption on the stability of p 1 (λ), we also have c 1 > 0 and d 1 > 0. Finally, observing that −Λ (α) 0 for all α we can conclude that c α > 0 and d α > 0 ∀ α.
The Routh-Hurwitz criterion ensures that p 1 (λ) is unstable if at least one of the following conditions is met: Let us first show that condition (a) is never met under the assumption of stability of p 1 (λ). From the definitions of the coefficients a, b, c α and d α we obtain We can now conclude that bc α − d α a > 0; indeed because of the stability of p 1 (λ), bc 1 − d 1 a > 0, and being −Λ (α) > 0 for all α > 1, the claim easily follows.
Let us now consider condition (b) and look for the existence of α > 1 such that We firstly rewrite this equation by using the definition of the involved coefficients and then we reorganise the terms in the latter, to write it as a second order polynomial in the variable Λ (α) , hence: where (after some algebraic manipulation): The coefficient A is positive, as well as the coefficient C, under the assumption of stability for p 1 (λ). Then the second order polynomial in Λ (α) can exhibit negative values if and only if that are exactly the conditions (21a) and (21b). Finally, an eigenvalue Λ (ᾱ) ,ᾱ > 1, must exist such that where x 1 and x 2 are the two real and negative roots of second order polynomial in Λ (α) . Let us finally consider condition (c) and look for the existence of α > 1 such that This is a second order polynomial in the variable Λ (α) whose leading coefficient, D u D v , is positive as well as the constant term, det(J 0 ), because of the stability of p 1 (λ). The polynomial can thus assume negative values if and only if Namely, the conditions (22a) and (22b). Let us observe that the latter do not depend on τ u and τ v and indeed they are the classical conditions required for the Turing instability to arise [16]: an eigenvalue Λ (ᾱ) ,ᾱ > 1, must exist such that where η 1 and η 2 are the two real and negative roots of e α = 0.

Appendix B. Non-relativistic limit of inertia-driven instability
In the main text we have proved that a Turing instability sets up driven by the inertial times τ u and τ v if any couple of conditions equations (21a) and (21b) or equations (22a) and (22b) are satisfied. Let us observe that the latter do not depend on the inertial times and thus they can be satisfied for a suitable choice of the model parameters, also for τ u = τ v = 0. The same could not hold true for the former one, explicitly dependent on the inertial times. The aim of this section is thus to study the non-relativistic limit of the inertia-driven instability.
For a sake of definitiveness let us assume τ v = θτ u for some θ > 0, that is the inertial times approach zero with the same rate. The characteristic polynomial given by (10) can thus be rewritten as where we have used (11) to rewrite the coefficients of λ 4 and λ 3 , and we have definedĉ (12)). Let us observe that d α and e α do not depend on the inertial times (see (13) and (14)).
In the limit τ u → 0 the latter results to be a singular polynomial, indeed its degree jumps from 4 once τ u > 0 to 2 for τ u = 0. Mathematically this means that two of the four roots of p(λ) should diverge to infinity, by determining which ones and the followed path will allow to conclude about the non-relativistic limit of the inertia-driven instability.
Let us start by looking at the roots that remain in a bounded domain. To do this let us set 6 λ = λ 0 + τ u λ 1 + . . . , impose p(λ) = 0 and by reordering the involved terms (see equation (B1)) according to the powers of τ u , we eventually get We can thus conclude that λ 0 is a solution of the second degree equation while λ 1 is obtained by solving (1 + θ)λ 3 0 + 2λ 0 λ 1 −ĉ α λ 0 + dλ 1 = 0.
In conclusion we get for the two roots: where we denoted by λ 0,± the two roots of (B3). Let us observe that the latter is the same second order equation one will obtain in the classical Turing framework; we have thus shown that in the non-relativistic limit two roots of the fourth degree characteristic polynomial p(λ) converge to the roots of the second order polynomial one should deal with in the classical Turing case. Let us now study the remaining two roots and determine their path towards infinity. As already stated the characteristic polynomial is singular, one should thus resort to the singular perturbation theory [60]. Let us set λ = ω/τ u and evaluate p α (λ) on λ = ω/τ u , then we get where the fourth degree polynomial q α (ω) has been defined by the last equality. Let us observe that p α (λ) vanishes if and only if q α (ω) does.  (19)) does not touch the boundary τ u = 0 or τ v = 0. Panel (b) refers to a choice of the parameters for which the classical Turing instability can emerge, indeed the black region (bounded above by equation (19)) reaches the boundary τ u = 0 or τ v = 0.
In conclusion if τ u and τ v are positive and sufficiently small, then the onset of Turing instability is ruled out by the roots (B5), i.e., those associated to the ones arising in the classical setting. Stated differently, if for τ u > 0 and τ v > 0, the instability can be initiated by conditions equations (21a) and (21b), then by decreasing the inertial times the patterns fade out and disappear before reaching the limit and thus the transition is not abrupt.
In figure 9 we report numerical results to complement the analytical findings described above. We selected two generic sets of parameter values γ = 4.0, β = 0.6, μ = 1.0, D u = 2.2 and D v = 0.2 (left panel) and γ = 4.0, β = 0.7, μ = 1.15, D u = 0.2 and D v = 2.2 (right panel), and we study the emergence of an inertiadriven instability as a function of τ u and τ v , a black dot corresponds to the presence of the instability while a white one to its absence. We can observe that the first set of parameters does not allow the onset of the instability for small enough values of the inertial times, indeed the black region, bounded below by condition (19) and above by condition (21b), does not intersect the axes τ u = 0 or τ v = 0. On the other hand the second set of parameters allows the existence of patterns for τ 0 = 0 or τ v = 0, the black region (bounded above by condition (21b)) reaches the axes. This means that in this case also classical Turing patterns are allowed.
Let ρ (α) j , j = 1, . . . , m, be the eigenvalues of A − C (α) = T −1 a + Λ (α) T −1 D. Then the above quoted result [46] ensures that Rρ (α) j < 0 for all j = 1, . . . , m and α = 1, . . . , n. Because T −1 D is a diagonal matrix, the diagonal elements of A − C (α) are translated to the left, while the off-diagonal elements do not vary. Invoking the Gershgorin circle theorem [61], we can take the elements of T −1 D sufficiently large, i.e., D j 1 or τ 1, such that τ |Iρ (α) j | < −Rρ (α) j for all j = 1, . . . , m, α = 1, . . . , n. We can thus conclude that the relation dispersion associated to equation (C3) is always negative and thus the Turing instability cannot develop. In conclusion linear kinetic systems without cross inhibition and equal (small) inertial times, or large diffusion coefficients, cannot exhibit Turing instability in the relativistic framework. The necessity of the latter assumptions remains an open question that we believe goes beyond the scope of this work.