Stabilization of hyperbolic reaction-diffusion systems on directed networks through the generalized Routh-Hurwitz criterion for complex polynomials
Riccardo Muolo Anthony Hastir Hiroya Nakao


Heavy Dark Matter in White Dwarfs: Multiple-Scattering Capture and Thermalization
    [
    April 24, 2024
=================================================================================







§ INTRODUCTION



Studying stability of a dynamical system is a paramount step before designing appropriate control laws in engineering. For linear differential equations, this task may be carried out by looking at the eigenvalues of the matrix describing the dynamics. When moving to the nonlinear setting, the eigenvalues of the Jacobian matrix of the vector field also give information, but only on the local stability of a particular equilibrium, as shown by the Hartman-Grobman Theorem for the case of fixed point equilibria, or the Andronov-Witt theorem for limit cycle equilibria. Similar results exist when the state-space is infinite-dimensional, typically a Banach or a Hilbert space. This setting is particularly useful to study partial differential equations (PDEs) in an abstract way. In this context, appropriate definitions of linearization have to be given in order to link the stability of a linearized system and its nonlinear version. These questions have been recently analyzed in <cit.> and <cit.> in order to study the stability of the equilibria of a nonisothermal axial dispersion tubular reactor based on a new version of the so-called "Fréchet differentiability". Therein, it was notably shown that the multiplicity and the stability of the equilibria are "diffusion driven" in the sense that the diffusion coefficient in the process determines whether one or three equilibria are exhibited. Moreover, in the case of only one equilibrium, the latter has been shown locally exponentially stable for the nonlinear system while bistability is depicted in the case of three equilibria, that is, stable – unstable – stable.


In such context, the goal is to annihilate the formation of instabilities, which disrupt the homogeneous equilibrium of our system. However, when discussing reaction-diffusion system in a stationary state, the objective may as well be antagonistic, as scholars often look for ways to obtain Turing instabilities, which are consequence of a diffusion-driven mechanism, i.e., caused by a difference in the diffusivities of the two interacting species <cit.>, or other kinds of instabilities. Turing theory was conceived in the framework of (parabolic) PDEs, but the setting has then been extended to lattices by Othmer and Scriven <cit.> and, successively, on complex networks by Nakao and Mikhailov <cit.>. The formalism developed in the latter work enabled the study of asymmetric couplings, where it was shown that such asymmetry enhanced the diffusion-driven mechanism of instability <cit.>. In fact, it was then generally proved that non-reciprocal interactions enhance inhomogeneities <cit.>. 

When dealing with a more realistic setting, an upper bound to the signal propagation speed is considered and one obtains hyperbolic reaction-diffusion system. Now, besides the "classical" diffusion-driven instability, there is a whole new class of patterns caused by the inertia in the signal propagation, as shown by Zemskov and Horsthemke <cit.>. The latter study has then been extended to networked systems <cit.>, but with symmetric coupling. In this work, we make a step further by studying the case of hyperbolic reaction-diffusion systems on directed networks. From the above discussion, we can already understand that in such setting it is very difficult to achieve stability: in fact, diffusion, inertia and the asymmetry of the coupling are all factors conspiring against the stability of the homogeneous equilibrium. Moreover, the fact of dealing with an asymmetric coupling has the consequence that the coupling operator has a complex spectrum, making the usual stability analysis techniques, such as the Routh-Hurwitz criterion <cit.>, not applicable. Hence, we put to use a recent generalization of such techniques, developed to deal with complex polynomials <cit.>. This allows us to visualize the stability regions, which can be fundamental in the design of systems where stability is fundamental, such as chemical reactors arrays. 


The paper is structured as follows: in the next section we will discuss the linear stability analysis for hyperbolic reaction-diffusion systems on networks; then, in Sec. <ref>, we will put to use the extended Routh-Hurwitz criterion for complex polynomials and analyze the stability conditions. Finally, in Sec. <ref>, we will aid the visualization of the stability regions through numerical simulations for the case of a chemical model, the Brusselator, and, in Sec. <ref>, we conclude by discussing some caveats and limitations of the linear stability analysis.




§ STABILITY OF HYPERBOLIC REACTION-DIFFUSION SYSTEMS ON DIRECTED NETWORKS


Let us consider a network on n nodes, defined by its adjacency matrix A, whose elements A_ij take value 1 if there is a link between nodes i and j, and 0 otherwise. A diffusion process on such network is modeled through the discrete diffusion equation. For instance, if we consider the diffusion of a chemical species u, we have, for each node i 
    u̇_i=D_u∑_i=1^nL_iju_j,

where D_u is the diffusion coefficient of species u, whose dimensions are [m^2/s] and L_ij=A_ij-δ_ijk_i is the discrete Laplacian matrix, a negative semi-definite matrix playing the role of the diffusion operator on networks <cit.>. k_i is the node degree, i.e., to how many other nodes it is connected. The spectrum of the discrete Laplacian is denoted by {Λ^α}_α∈{1,...,n}, where Λ^1=0 is the Laplacian's spectral abscissa, i.e., the maximum real part of the spectrum. If the network is connected, then we have that Λ^1>Λ^α ∀α∈{2,...,n}. Let us remark that, the network being symmetric, the spectrum of the Laplacian is real.

The above description works well when the domain where the dynamics takes place (e.g., the network) is small or the diffusivity D_u is large, i.e., one can consider the signal propagation as instantaneous. However, in many practical applications this approximation does not work, hence the diffusion equation needs to be corrected by accounting for an inertia (hence, with dimension [1/s]) in the signal propagation, generally indicated with τ, as done by Cattaneo in the 50s <cit.>. Eq. (<ref>) then becomes

    τ_uü_i+u̇_i=D_u∑_i=1^nL_iju_j.

The above equation takes several names (Cattaneo equation, telegrapher's equation, relativistic heat equation, to name a few) and it is usually studied in the framework of PDEs, where it is a hyperbolic equation. It has been rigorously derived on networks in <cit.>, where the name hyperbolic has been kept despite being an ODE in such framework.


Let us now consider a system of two species, u and v, which interact by means of two nonlinear functions f and g. We consider n identical copies of such system, whose dynamics is described by the following 
    
        u̇_i=f(u_i,v_i), 
 v̇_i=g(u_i,v_i).
    
 The above system could describe, for instance, the dynamics into a well-stirred environment, such as a chemical reactor. Such dynamics could be, for instance, a stationary state or an oscillating one. Let us observe that, since we are considering a homogeneous domain, there is no propagation, hence the second time derivatives do not appear. Let us know study the case of n copies of the above system and then couple them, by letting the chemical concentrations u_i and v_i diffuse between the reactor in a Fickean (linear) fashion. We are, hence, assuming that the dynamics take place within the nodes, while the links are the means by which the particles can move between adjacent systems. Such setting is schematized in Fig. <ref>.



If we also account for a finite propagation, we obtain the following hyperbolic reaction-diffusion networked system: 


    τ_uü_i+u̇_i=f(u_i,v_i)+D_u∑_i=1^nL_iju_j, 
    
    τ_vv̈_i+v̇_i=g(u_i,v_i)+D_v∑_i=1^nL_ijv_j.

As explained in the introduction, let us stress that the described framework is general, and can describe different systems, from chemical systems to ecological neural ones, to name a few. 


Let us now discuss the stability of system (3) . For sake of simplicity, let us assume that each isolated system lies in a stationary state (u^*,v^*). Such equilibrium point can be found by solving the system f=0, g=0. The case of oscillating, or even chaotic, dynamics can be treated with the same techniques, but one cannot obtain analytical results and one must resort to numerical tools, such as the Master Stability Function <cit.>. The first step is to impose the stability of the isolated system, meaning that the Jacobian matrix computed at the stationary state needs to be stable, i.e., given 

    J_0=[ 
f_u     f_v 
 g_u     g_v
],
 where f_u is the derivative of the function f with respect to the variable u computed in (u^*,v^*), we need the trace to be negative and the determinant to be positive, i.e, 

    f_u+g_v<0 ,
 f_ug_v -f_vg_u>0.
    
 
Let us point out that, in general, the inertia in the signal propagation can affect the stability of the homogeneous equilibrium <cit.>. In what follows, however, we will consider settings in which the stability is not affected by the inertia. One can then study whether an inhomogeneous perturbation about such stable homogeneous equilibrium disrupts it, giving rise to an instability, or goes back to the stationary state. The full discussion of the case in which the network is symmetric can be found in <cit.>. In this work, we want to extend such framework by considering asymmetric couplings, which are more common in practical applications. In this case, the equation describing the dynamical process are the same as Eq. (<ref>), but with the important difference that now the Laplacian is a directed matrix, defined as L_ij=A_ij-δ_ijk_i^in, where k_i^in is the in-degree, i.e., how many directed links are incoming to node i. Hence, its spectrum will be, in general, complex, meaning that now Λ^α=Λ_Re^α+iΛ_Im^α ∀α∈{2,...,n}, while Λ^1=0 remains the largest eigenvalue, i.e., Λ_Re^α<0=Λ^1 for any α∈{2,...,n}.


The linear stability analysis breaks down to compute 
    [ 
    J_11    J_12
     J_21    J_22]=0,
 
with


J_11= τ_uλ^2+λ-f_u-D_u Λ_Re^α-i D_u Λ_Im^α,

J_12=-f_v, 

J_21= -g_u, 

J_22 = τ_vλ^2+λ-g_v-D_v Λ_Re^α-i D_v Λ_Im^α,


where λ is the growth rate of the perturbation in the linear regime, i.e., if it is negative, the system will go back to its stable state, while if it is positive, there will be an exponential instability driving the system away from the equilibrium.


From the above expressions, (<ref>) is equivalent to the following equation  

    λ^4+∑_j=1^4 (a_j+ib_j)λ^4-j=0,
 
whose coefficients a_j, b_j, j=1,2,3,4, are given by


a_1=ε(τ_u+τ_v),


b_1=0,

a_2=ε(1-g_vτ_u-f_uτ_v-(τ_u D_v+τ_v D_u)Λ_Re^α),

 b_2=ε(-Λ_Im^α(τ_u D_v+τ_v D_u)),

a_3=ε(-g_v+f_u+(D_u+D_v)Λ_Re^α),

b_3=ε(-Λ_Im^α(D_u+D_v)),

a_4=ε[(f_uD_v +g_vD_u)Λ_Re^α+f_u g_v-f_v g_u 

                             +D_u D_v((Λ_Re^α)^2+(Λ_Im^α)^2)],

b_4=ε((g_vD_u+f_uD_v)Λ_Im^α+D_uD_vΛ_Re^αΛ_Im^α),


with ε := (τ_uτ_v)^-1. The stability of the above polynomial cannot be studied through the Routh-Hurwitz criterion <cit.>, which can be applied only to polynomials with real coefficients. However, we can rely on a recent generalization to the complex case, as we discuss in the next section.



§ THE GENERALIZED ROUTH-HURWITZ CRITERION


This section is devoted to the characterization of stability/instability of the polynomial (<ref>). To this end, one shall rely on a generalization of the classical Routh-Hurwitz criterion for polynomials with complex coefficients. This method has been initiated in <cit.> and recently presented in an algorithmic way in <cit.>. The main advantage of <cit.> is the ease of implementation of the method as a pedagogical tool. From an algorithmic point of view, the method presented as in <cit.> does not involve more than determinants of 2× 2-matrices at each step, no matter the degree of the polynomial.  The idea of the proposed algorithm consists in building a table of coefficients similar to the one introduced by Hurwitz in <cit.> and then to exploit its first column to characterize stability/instability. For a 4-th order polynomial of the form (<ref>), the table constructed via <cit.> is given in Tab. <ref>.


As a consequence, the stability/instability of (<ref>) is characterized in the following proposition.


The polynomial (<ref>) is stable if and only if the following three conditions are satisfied, namely

    υ := ε(τ_u + τ_v)[τ_u+τ_v - g_vτ_u^2 - f_uτ_v^2 - τ_u^2 D_vΛ_Re^α.
       .- τ_v^2 D_uΛ_Re^α] - Λ_Im^α(D_u+D_v)^2 > 0
       γ := α(a_3α - a_1^3 a_4 + a_1 b_2 b_4) + [α b_2 - a_1β]β > 0
       (γβ-α)[α^2 b_4 - (a_1^2 b_4 - b_2 b_4)(α b_2 - a_1β)]
       + (a_1^2 a_4 - b_2 b_4)γ^2 > 0,

where α := ε^2υ and β := a_1(a_1 b_3 - b_4) + a_3 b_2.


First, observe that, thanks to <cit.>, the polynomial (<ref>) is stable if and only if the coefficients a_k^(k) in Tab. <ref> are all positive.
Then, note that the inertia constants τ_u and τ_v are positive, which makes condition a_1^(1)>0 trivially satisfied. By looking at the constant a_2^(2), there holds

    a_2^(2) = a_1(a_1 a_2-a_3)-b_2^2
       = (1/τ_u+1/τ_v)^2ε[1-g_vτ_u-f_uτ_v.
       .-(τ_u D_v+τ_v D_u)Λ_Re^α] + (1/τ_u+1/τ_v)ε[g_v+f_u.
       .+ (D_u+D_v)Λ_Re^α] - ε^2[Λ_Im^α(D_u+D_v)]^2
       =(1/τ_u+1/τ_v)ε[(1/τ_u+1/τ_v) - g_vτ_u/τ_v - f_uτ_v/τ_u.
       .- τ_u/τ_vD_vΛ_Re^α - τ_v/τ_uD_uΛ_Re^α] - ε^2Λ_Im^α(D_u+D_v)^2
       =ε^2υ,

where υ is defined in (<ref>). Then, a_2^(2)> 0 if and only if υ >0, which shows (<ref>).  Let us now move to the coefficient a_3^(3). One has that 

    a_3^(3) = α (α a_3 - a_1 a_4^(2)) + (α b_2 - a_1 b_3^(2))b_3^(2)
       = a_3α^2 - a_1α(a_1 a_4^(1) - b_2 b_4) + b_2α(a_1 b_3^(1)+b_2 a_3)
        - a_1(a_1 b_3^(1) + b_2 a_3)^2.

The condition (<ref>) then follows from the equalities a_4^(1) = a_1a_4 and b_3^(1) = a_1b_3-b_4. The condition (<ref>) may be obtained in a similar way, by expressing a_4^(4) as a function of a_3^(3) = γ and a_2^(2) = ε^2υ.






§ NUMERICAL RESULTS



As stated in the introduction, our goal is to annihilate the formation of instabilities, which disrupt the homogeneous equilibrium of our system. In what follows, we will focus our attention to instability arising from the inertia in the signal propagation, or inertia-driven instabilities. In fact, the case of diffusion-driven instabilities has been thoroughly studied in past works, even though not specifically in the context of asymmetrically coupled hyperbolic reaction-diffusion systems, while inertia-driven instabilities constitute a rather unexplored setting; a gap which we aim to contribute in filling with the present work. 


As we have shown in the previous section, we can find explicit expressions for the stability conditions thanks to the generalized Routh-Hurwitz criterion. Nonetheless, such expressions remain cumbersome and difficult to handle analytically and the stability region is difficult to visualize. Hence, we resort to numerical simulation to validate the developed theory. For this, we will make use of the Brusselator, a model describing the dynamics of an autocatalytic reaction, paradigmatic in the study of chemical instabilities and synchronization in reaction diffusion-systems <cit.>. The version of the model that we consider <cit.> models the nonlinear interactions between the two chemical species u and v by mean of two positive parameters b and c, namely 
    f(u,v)=1-(b-1)u+cu^2v, 
    
        g(u,v)=bu-cu^2v.


The above system has a unique fixed point in 
(u^*=1,v^*=b/c),  hence, the Jacobian matrix becomes


    J_0=[ b-1     c 
 -b     -c
].


Now, if we consider the above model on an asymmetric networks of n nodes, diffusively coupled, and with a finite signal propagation, we obtain the following hyperbolic reaction-diffusion networked system

    τ_uü_i+u̇_i=1-(b-1)u_i+cu_i^2v_i
    +D_u∑_i=1^nL_iju_j, 
    τ_vv̈_i+v̇_i=bu_i-cu_i^2v_i+D_v∑_i=1^nL_ijv_j.



If we perform a linear stability analysis, as in Sec. <ref>, and then apply the generalized Routh-Hurwitz criterion, as in Sec. <ref>, we obtain the stability condition for system (<ref>). As already mentioned, we want to focus on the instability driven by the inertia in the signal propagation, hence let us fix the diffusion coefficients and the model parameters and focus on the inertia terms τ_u and τ_v. The other free parameters will be the real and imaginary parts spectrum of the Laplacian, i.e., Λ_Re^α and Λ_Im^α, which will be considered as continuous parameters. In fact, in principle, the spectrum of the coupling operator could be anywhere in ℂ^-, but then there will be only n eigenvalues, i.e., the size of the network. Except for Λ^1=0, all the other eigenvalues are potential instability modes and need to be controlled onto the stability region, as can be visualized in Fig. <ref>.



The numerical results can be visualized in Fig. <ref>, where we plot the stability region for system (<ref>), i.e., the region where the conditions (<ref>), (<ref>) and (<ref>) are satisfied. Let us observe that, for this choice of parameters in which the two diffusion coefficients are equivalent, i.e., D_u≡ D_v, the only factors causing the instability are the presence of the inertia terms τ_u and τ_v and the imaginary part of the spectrum Λ_Im^α, induced by the asymmetry of the coupling. In fact, diffusion driven instabilities arise only when D_u<D_v <cit.>.



In Fig. <ref>a), we can see that, for that choice of parameters, the system becomes unstable also when the coupling is symmetric for a sufficiently large τ_v, while in the other three panels (b-d), only through the imaginary part of the coupling operator instability can be reached. This is very important, for instance, in the design of couplings between reactors, as a certain degree of asymmetry may trigger unwanted instabilities.






§ CONCLUSIONS



In this work we have discussed the stability for hyperbolic reaction-diffusion networked systems with asymmetric coupling, showing that achieving stability is a difficult task due to several factors known to trigger instability. Since an asymmetric coupling yields a complex spectrum, i.e., a complex polynomial whose stability needs to be determined, we have carried out the analysis through a recent generalization of the Routh-Hurwitz criterion for complex polynomials. We have analytically studied the stability conditions and we have numerically visualized the stability region for the Brusselator model, a toy model of an autocatalytic chemical reaction, often adopted in the study of reaction-diffusion systems. Nonetheless, let us remark that our results are general, and the analysis can be carried out for any model of the form of Eq. (3). Moreover, the generalization of the Routh-Hurwitz criterion can be used for any polynomial with complex coefficients. 

Lastly, let us point out that, in many real-world systems, the coupling is not only asymmetric, but also non-normal, meaning that the adjacency matrix A satisfies the condition AA^*≠ A^*A <cit.>. In that case, the prediction given by a linear stability analysis may be valid only when the systems is subjected to infinitesimal perturbations. However, in applications the perturbations are always finite and instabilities can arise even if the Routh-Hurwitz conditions are satisfied <cit.>. In such cases, the stability analysis needs to be complemented with simulations, in order to be sure that the system will remain stable against small perturbations. 



The authors are grateful to Timoteo Carletti for interesting discussions. R.M. and H.N. acknowledge JSPS KAKENHI JP22K11919, JP22H00516, and JST CREST JP-MJCR1913 for financial support. A.H. is supported by a FNRS Postdoctoral Fellowship, Grant CR 40010909.


9
Hastir_TAC
Hastir, A., Winkin, J. & Dochain, D. On Exponential Bistability of Equilibrium Profiles of Nonisothermal Axial Dispersion Tubular Reactors. IEEE Transactions On Automatic Control. 66, 3235-3242 (2021)
Hastir_Automatica
Hastir, A., Winkin, J. & Dochain, D. Exponential stability of nonlinear infinite-dimensional systems: Application to nonisothermal axial dispersion tubular reactors. Automatica. 121 pp. 109201 (2020)
Turing
Turing, A. The chemical basis of morphogenesis. Phil. Trans. R. Soc. Lond. B. 237 pp. 37 (1952)
OS1971Othmer, H.G. & Scriven, L.E. Instability and dynamic pattern in cellular networks. J. Theor. Biol. 32 pp. 507 (1977)
NM2010Nakao, H. & Mikhailov, A. Turing patterns in network-organized activator-inhibitor systems. Nat. Phys. 6 pp. 544 (2010)
Asllani1Asllani, M., Challenger, J., Pavone, F., Sacconi, L. & Fanelli, D. The theory of pattern formation on directed networks. Nature Communication. 5 (2014)
carlettiCarletti, T. & Muolo, R. Non-reciprocal interactions enhance heterogeneity. Chaos Solit. Fractals. 164 pp. 112638 (2022)
ZH2016Zemskov, E. & Horsthemke, W. Diffusive instabilities in hyperbolic reaction-diffusion equations. Phys. Rev. E. 93 pp. 032211 (2016)
jop_carletti
Carletti, T. & Muolo, R. Finite propagation enhances Turing patterns in reaction–diffusion networked systems. J. Phys. Complex. 2, 045004 (2021)
Routh1877Routh, E. Stability of a given state of motion. (MacMillan,1877)
Hurwitz1895Hurwitz, A. Über die Bedingungen, unter welchen eine Gleichung nur Wurzeln mit negativen reellen Theilen besitzt. Math. Ann. 46, 273 (1895)
HastirMuolo_RH_GenHastir, A. & Muolo, R. A generalized Routh-Hurwitz criterion for the stability analysis of polynomials with complex coefficients: application to the PI-control of vibrating structures. IFAC J. Syst. Control. 26 pp. 100235 (2023)
latora_nicosia_russo_2017Latora, V., Nicosia, V. & Russo, G. Complex Networks: Principles, Methods and Applications. (Cambridge University Press,2017)
Cattaneo1958
Cattaneo, C. Sur une forme de l'équation de la chaleur éliminant le paradoxe d'une propagation instantanée. C. R. Acad. Sci..247, 431 (1958)
pecora1998master
Pecora, L. & Carroll, T. Master stability functions for synchronized coupled systems. Phys. Rev. Lett.. 80, 2109 (1998)
Frank_1946Frank, E. On the zeros of polynomials with complex coefficients. Bulletin Of The American Mathematical Society. 52, 144-157 (1946)
PrigogineNicolis1967Prigogine, I. & Nicolis, G. Symmetry breaking instabilities in dissipative systems. J. Chem. Phys. 46 pp. 3542 (1967)
gallaBoland, R., Galla, T. & McKane, A. How limit cycles and quasi-cycles are related in a system with intrinsic noise. J. Stat. Mech. pp. P09001 (2008)
malbor_teoAsllani, M., Lambiotte, R. & Carletti, T. Structure and dynamical behavior of non-normal networks. Sci. Adv.. 4, Eaau9403 (2018)
jtbMuolo, R., Asllani, M., Fanelli, D., Maini, P. & Carletti, T. Patterns of non-normality in networked systems. J. Theor. Biol. 480 pp. 81 (2019)






