Analyzing Diffusion and Flow-driven Instability using Semidefinite Programming

Analyzing Diffusion and Flo w-driv en Instabilit y using Semidefinite Programming Y utak a Hori ∗ , Hiroki Miy azako ∗ Abstract Diffusion and flow-driv en instability , or transp ort-driv en instability , is one of the cen tral mec hanisms to generate inhomogeneous gradient of concentrations in spatially distributed c hemical systems. Ho wev er, verifying the transp ort-driv en instability of reaction-diffusion-adv ection systems requires chec king the Jacobian eigenv alues of in- finitely man y F ourier mo des, whic h is computationally intractable. T o ov ercome this limitation, this paper proposes mathematical optimization algorithms that determine the stability/instabilit y of reaction-diffusion-advection systems b y finite steps of alge- braic calculations. Specifically , the stability/instabilit y analysis of F ourier mo des is form ulated as a sum-of-squares (SOS) optimization program, which is a class of conv ex optimization whose solvers are widely av ailable as softw are pack ages. The optimiza- tion program is further extended for facile computation of the destabilizing spatial mo des. This extension allo ws for predicting and designing the shap e of concentra- tion gradient without sim ulating the gov erning equations. The streamlined analysis pro cess of self-organized pattern formation is demonstrated with a simple illustrativ e reaction mo del with diffusion and adv ection. Keyw ords: reaction-diffusion-adv ection model, transp ort-driv en instabilit y , con vex op- timization, semidefinite programming, self-organized pattern formation * The authors contributed equally to this work. Y utak a Hori is with Departmen t of Applied Ph ysics and Ph ysico-Informatics, Keio Universit y . 3-14-1 Hiy oshi, Kohoku-ku, Y okohama, Kanagaw a 223-8522, Japan. Hiroki Miyazak o is with Department of Information Physics and Computing, The Universit y of T okyo. 7-3-1 Hongo, Bunky o-ku, T oky o 113-8656, Japan. Corresp ondence should b e addressed to Y utak a Hori (yhori@appi.keio.ac.jp). This work should b e cited as Y. Hori, H. Miyazak o, “Analysing diffusion and flow-driv en instability using semidefini te programming,” Journal of the R oyal So ciety Interfac e , V ol. 16, No. 150, 20180586, 2019. https://doi.org/10.1098/ rsif.2018.0586 1 1 In tro duction Molecular transportation is a fundamen tal mec hanism that couples spatially distributed c hemical reaction net works and enables self-organized pattern formation in biology and c hemistry . F or example, in developmen tal biology , passiv e diffusion and advectiv e flow of signaling molecules within and b et ween cells are known to play a cen tral role in identifying and regulating positions and shap es during em bryo-genesis [1 – 3] and cell division [4 – 8]. In c hemistry , diffusion-driv en self-organized patterns were reconstituted in crafted reactors with commo dit y c hemicals [9, 10]. More recen tly , synthetic biologists hav e b een attempting to utilize quorum sensing, a mechanism of diffusion-based cell-to-cell signaling of E. coli, to regulate a p opulation of synthetic biomolecular reaction systems that comm unicate [11 – 15] and synchronize with each other [16 – 19]. The theoretical foundation of transp ort-driv en self-organization was established b y T ur- ing [20], where he sho wed using a reaction-diffusion equation that molecular concentrations can form a spatially p erio dic gradien t due to the in teraction of reaction and diffusion de- spite the av eraging nature of diffusion. Later, Rovinsky and Menzinger [21, 22] verified b oth mathematically and exp erimen tally that differen tial flo w, or advection, of molecules can also induce an oscillatory spatio-temporal concentration gradien t. In these w orks, math- ematical analyses rev ealed that the pattern formation was due to the destabilization of spatial oscillation mo des, and the destabilization w as induced by the difference of diffusion and flow rates b et ween reactiv e molecules. Mathematically , this can b e in terpreted that molecular transp ortation caused sp on taneous growth of some spatial F ourier mo des and led to spatially p eriodic pattern formation at steady state. Curren tly , a widely used approac h to analyzing the transp ort-driv en pattern formation is the linear stability/instabilit y analysis of spatial F ourier comp onen ts. Specifically , the Jacobian eigenv alues of the gov erning reaction-diffusion-adv ection equations are computed for different spatial F ourier mo des to explore the existence of unstable harmonic com- p onen ts. Then, spatially p erio dic pattern formation is exp ected if an eigen v alue of the Jacobian resides in the right-half complex plane for some non-zero frequency comp onen ts. F or relatively simple reaction systems, analytic instability conditions we re obtained using this approac h, and the parameter space for pattern formation w as thoroughly character- ized [1, 23, 24]. F or large-scale and complex reaction systems, on the other hand, insta- bilit y analysis requires substanti al computational efforts to iteratively compute Jacobian eigen v alues for eac h F ourier component. How ever, this approac h is essentially incapable of 2 dra wing mathematically rigorous conclusion on the stabilit y of reaction-diffusion-adv ection systems since the spatial gradient of chemical concen trations is represented with infinitely man y F ourier mo des. Another approach to analyzing the stabilit y/instabilit y of reaction-diffusion equations is to use Ly apunov’s method, where the core idea is to guarantee the energy dissipation of certain F ourier comp onen ts b y constructing a Ly apunov function. The exploration of a Ly a- puno v function typically reduces to an algebraic optimization problem called semidefinite programming (SDP), whic h is an efficien tly solv able class of conv ex optimization [25, 26]. Existing w orks presented sufficient conditions for the stabilit y of reaction-diffusion equa- tions, certifying non-existence of spatially inhomogeneous solutions [27 – 29]. A more general SDP based framework was also developed for the stability analysis of univ ariate nonlin- ear partial differen tial equations by exploring p olynomial integral Lyapuno v functions [30]. Ho wev er, the focus of these works is spatially homogeneous behavior, and the conditions are not necessarily suitable for studying diffusion-driv en instability that leads to spatially p eriodic pattern formation. This paper presents computationally tractable necessary and sufficient conditions for the lo cal stability/instabilit y of reaction-diffusion-adv ection systems. The prop osed approach is amenable to computational implementation and is particularly con venien t for analyzing transp ort-driv en instabilit y . Sp ecifically , w e derive linear matrix inequality (LMI) con- ditions that certify the stabilit y/instability of infinitely many F ourier comp onen ts. This algebraic condition can b e efficien tly verified with semidefinite programming (SDP) [25, 26]. The deriv ation of these conditions is based on the linear stability analysis of F ourier com- p onen ts. Although the analysis essen tially requires chec king the ro ots of infinitely many c haracteristic p olynomials with complex co efficien ts, we sho w that the lo cal stabilit y is equiv alent to an existence of sum-of-squares (SOS) decomp osition of certain Hurwitz p oly- nomials, which can b e form ulated as a semidefinite optimization problem [31]. Based on this result, we further deriv e conditions for certifying the stability/instabilit y of reaction- diffusion-adv ection systems for a presp ecified set of spatial frequency . This extension allows for facile computation of the destabilizing spatial mo des, enabling the prediction of spatial oscillations without sim ulating the gov erning equation. The prop osed instability analy- sis helps not only b etter understanding of reaction-diffusion-advection kinetics but also engineering of chemical reactions in syn thetic biology and chemistry . The follo wing notations are used in this pap er. N is a set of p ositiv e in tegers. Z is a 3 set of all integers. R is a set of real n umbers. R n × n is a set of n b y n matrices with real en tries. | A | is a determinant of the matrix A . A  O means that A is p ositiv e semidefinite. deg( p ( x )) is the degree of a p olynomial p ( x ) . d x e is the ceiling function, that is, the smallest in teger that is greater than or equal to x . 2 Mo del of reaction-diffusion-advection system W e consider a reaction-diffusion-advection process of n molecular sp ecies defined in a finite spatial domain Ω . F or notational simplicit y , w e consider only one dimensional space Ω := [0 , L ] , whose co ordinate is sp ecified b y the sym b ol x , but all theoretical results sho wn in this paper can be generalized to higher dimensions. Let C i ( x, t ) denote the concen tration of the i -th molecule at p osition x ∈ Ω and time t and C ( x, t ) b e the vector of the molecular concentrations C ( x, t ) := [ C 1 ( x, t ) , C 2 ( x, t ) , · · · , C n ( x, t )] T . The spatio- temp oral dynamics of the molecular concentrations are then describ ed by the reaction- diffusion-adv ection equation ∂ C ( x, t ) ∂ t = f ( C ( x, t )) + D ∂ 2 C ( x, t ) ∂ x 2 + V ∂ C ( x, t ) ∂ x , (1) where the function f ( · ) is a C 1 v ector-v alued function go verning lo cal reactions, and D := diag( d 1 , d 2 , · · · , d n ) and V := diag( v 1 , v 2 , · · · , v n ) are the co efficien ts of diffusion and flow v elo cit y , resp ectiv ely . W e assume d i > 0 and v i ≥ 0 in the follo wing theoretical dev elopment ( i = 1 , 2 , · · · , n ) . The reaction-diffusion-advection system sho ws a v ariety of spatio-temporal dynamics ranging from spatially uniform steady state to spatio-temporal oscillations dep ending on the parameters and the stoichiometry of the reactions. As a motiv ating example, we consider the following set of reactions that consists of molecules P , Q , and a pro duct R : P + 2 Q − − → 3 Q , Q − − → R , (2) where the molecule Q catalyzes its o wn pro duction using the substrate P , and the product R is inert to the reactions [32, 33]. The substrate P is constan tly supplied at a constant rate and all molecules are drained at the same rate as illustrated in Fig. 1(A). W e assume that all molecules are spatially distributed in one dimensional space Ω := [0 , 30 π ] with p eriodic b oundary conditions, and they are transported by constant flo w and passiv e diffusion. The 4 Position x 0 Ti me t 20 40 60 80 0 1000 2000 3000 4000 5000 6000 Position x 0 Ti me t 20 40 60 80 0 1000 2000 3000 4000 5000 6000 Position x 0 Ti me t 20 40 60 80 01 0 0 0 2 0 0 0 3 0 0 0 4 0 0 0 5 0 0 0 6 0 0 0 P Q Q Q Q Q P P P P P R R R a a a aa b Const ant flo w Const ant flo w A BC D 0.5 0.4 0.2 0.3 0.8 0.6 0.2 0.4 0.5 0.4 0.3 0.2 0.1 Figure 1: A reaction-diffusion-adv ection system inducing flo w-driv en instabilit y (A) Sc hematics of the reaction model with diffusion and adv ection. (B-D) Spatio-temp oral dynamics of concentration gradien t for different v alues of parameters. spatio-temp oral dynamics of the molecular concentrations are then mo deled by ∂ C 1 ∂ t = − C 1 C 2 2 + a (1 − C 1 ) + d ∂ 2 C 1 ∂ x 2 + v 1 ∂ C 1 ∂ x , ∂ C 2 ∂ t = C 1 C 2 2 − ( a + b ) C 2 + ∂ 2 C 2 ∂ x 2 + v 2 ∂ C 2 ∂ x , (3) where C 1 ( x, t ) and C 2 ( x, t ) denote the concen trations of P and Q , resp ectively . In the mo del (3), the reaction rates are normalized by that of the auto catalytic reaction in (2). The constan t a and b represent the normalized supply rate of the substrate P and the pro duction rate of R (see Fig. 1(A)). The spatial co ordinate x is defined so that the diffusion co efficien t of the molecule Q b ecomes one. Fig. 1(B)-(D) illustrate qualitatively differen t spatio-temp oral dynamics for the reaction- diffusion-adv ection system (3) for different choices of b , v 1 and v 2 . Sp ecifically , ( b, v 1 , v 2 ) = (0 . 040 , 0 , 0) for Fig. 1(B), ( b, v 1 , v 2 ) = (0 . 055 , 0 , 0) for Fig. 1(C), and ( b, v 1 , v 2 ) = (0 . 040 , 1 . 897 , 0 . 3162) for Fig. 1(D) w ere used (see Metho d section for the other param- 5 eters and initial v alues). The concen tration of P forms spatially p eriodic oscillations as b increases from 0 . 040 to 0 . 055 despite the a veraging effect of the passive diffusion (Fig. 1(B), (C)). This pattern formation, which is widely known as T uring pattern formation, is caused by the destabilization of certain spatial oscillation mo des due to the difference of the diffusion rates (diffusion instability) [2]. On the other hand, the spatio-temp oral oscillations in Fig. 1(D) is induced by a differen t destabilization mechanism based on the adv ective transp ortation of the molecules. In the follo wing sections, we first review that these destabilizing effects can b e explained b y probing the lo cal instability of F ourier modes of the reaction-diffusion-adv ection system. W e then present no vel algebraic stability conditions for the stability/instabilit y analysis of infinitely many F ourier mo des with semidefinite programming. 3 Stabilit y analysis of spatial F ourier comp onen ts T o analyze instability of spatial mo des asso ciated with spatial pattern formation, w e lin- earize the equation (1) around a spatially homogeneous equilibrium p oin t C ( x, t ) = ¯ C . In other words, ¯ C is an equilibrium of lo cal reactions satisfying f ( ¯ C ) = 0 . Assuming the exis- tence of suc h equilibrium p oin t, we can write the ev olution of the molecular concen trations around ¯ C b y ∂ c ( x, t ) ∂ t = A c ( x, t ) + D ∂ 2 c ( x, t ) ∂ x 2 + V ∂ c ( x, t ) ∂ x , (4) where c ( x, t ) := C ( x, t ) − ¯ C is the v ector of relativ e concen trations, and A is the Jacobian of f ( · ) ev aluated at ¯ C . The lo cal b eha vior of c ( x, t ) can b e characterized b y its spatial F ourier comp onen ts ˜ c ( ζ , t ) when the boundaries of the spatial domain ∂ Ω satisfy certain conditions. Sp ecifically , let ˜ c ( ζ , t ) b e the F ourier co efficien ts of c ( x, t ) satisfying c ( x, t ) = X ζ ∈Z ˜ c ( ζ , t ) e j ζ x , (5) where Z is a set of discrete frequency v ariables that dep ends on the b oundary conditions 6 (see Remark 1). Multiplying by e − j ζ x and taking the in tegral of b oth sides of (4), w e ha ve d ˜ c ( ζ , t ) dt = ( A − ζ 2 D + j ζ V ) ˜ c ( ζ , t ) . (6) This equation represents the dynamics of frequency comp onen t ˜ c ( ζ , t ) . It should b e no- ticed that, for eac h fixed ζ , (6) is an n -th order linear time-in v ariant ordinary differential equation (ODE) with complex co efficien ts. Th us, the reaction-diffusion-advection equation is decomp osed into a set of infinitely many ODEs. Remark 1. The frequency v ariable ζ in (5) takes discrete v alues that dep end on the b oundary conditions. Sp ecifically , let the set of all frequency v ariables for a given b oundary condition b e denoted by Z . Then, Z = { 2 k π /L } k ∈ Z for p eriodic b oundary conditions. When V = O , Z is defined by Z = { k π /L } k ∈ Z for the Neumann b oundary condition and Z = { k π /L } k ∈ Z \{ 0 } for the Dirichlet b oundary condition, respectively [34]. The reader is referred to [34] for other boundary conditions. It should also b e noted that Im[ ˜ c ( ζ , t )] = 0 for the Neumann b oundary with V = 0 and Re[ ˜ c ( ζ , t )] = 0 for the Diric hlet b oundary with V = 0 , implying that the F ourier cosine and sine transforms are used to obtain ˜ c ( ζ , t ) , resp ectiv ely . Since the equation (6) represents the dynamics of F ourier components, c ( x, t ) asymptot- ically conv erges to zero if the growth rate of ˜ c ( ζ , t ) is negativ e for all frequency comp o- nen ts ζ . On the other hand, if there exists a frequency comp onen t ˜ c ( ζ , t ) with non-zero frequency ζ whose growth rate is positive, the corresp onding non-zero spatial mo de is amplified around the spatially homogeneous equilibrium ¯ C . In the nonlinear reaction- diffusion-adv ection system (1), the unstable frequency component p oten tially induces pe- rio dic pattern formation in space. In other words, the existence of a growing frequency comp onen t with non-zero frequency is a necessary condition for the generation of spatially p eriodic pattern in the system (1). More formally , this can b e stated as lo cal stabilit y/instability of the nonlinear reaction- diffusion-adv ection system (1). The nonlinear system (1) is said to b e lo cally stable around ¯ C if the linear system (6) is asymptotically stable. F or linear systems, asymptotic stability is determined from the eigen v alues of the matrix A − ζ 2 D + j ζ V . More specifically , the follo wing lemma holds (see Supplementary Information for the pro of ). Lemma 1. Consider the linear reaction-diffusion-advection system (4). The system (4) is 7 asymptotically stable if and only if all ro ots of the c haracteristic p olynomial ϕ ( ζ , s ) := | sI − A + D ζ 2 − j ζ V | = 0 (7) are in the op en left-half complex plane { s ∈ C | Re[ s ] < 0 } for all ζ ∈ Z , where Z is defined in Remark 1. Th us, the stabilit y analysis of the reaction-diffusion-adv ection system (4) reduces to finding the ro ots of the p olynomial ϕ ( ζ , s ) = 0 for ζ ∈ Z . It should, how ever, b e noted that verifying the stabilit y condition requires solving infinitely man y p olynomial equations with complex coefficients (7), and analytic conditions are hardly obtained except for some sp ecial cases [21, 35]. Consequen tly , most of the existing works resort to approximate stabilit y analysis b y solving (7) for a finite range of quan tized ζ . In the next section, we prop ose a nov el computational approach that ov ercomes this limitation. Sp ecifically , we in tro duce a computationally tractable optimization problem that certifies the stability of the reaction-diffusion-advection system. 4 Linear matrix inequalit y condition for stabilit y/instabilit y analysis In this section, w e sho w matrix inequality conditions for computational verification of the stabilit y/instability of the system (4). F or theoretical dev elopment, we here consider a sligh tly mo dified condition of Lemma 1 that all ro ots of the p olynomial (7) lie in the op en left-half complex plane for ζ ∈ R instead of ζ ∈ Z . Mathematically , this b ecomes only a sufficien t condition for the asymptotic stabilit y of the system (4). How ever, conv en tionally , this condition has been used for the stabilit y/instability analysis of reaction-diffusion- adv ection systems [20 – 22] since the elements of Z , whic h lie on the real axis ( −∞ , ∞ ) tend to b e close to each other as the length of the domain L becomes large, and the gap from the necessary condition b ecomes small enough in practical applications. T o this end, we first review an algebraic stabilit y condition due to Hurwitz (Theorem 3.4.68 of [36]). Let ϕ Re ( ζ , s ) and ϕ Im ( ζ , s ) b e real and imaginary parts of the complex p olynomial ϕ ( ζ , j s ) , resp ectiv ely . Sp ecifically , we denote b y ϕ ( ζ , j s ) = ϕ Re ( ζ , s ) + j ϕ Im ( ζ , s ) (8) 8 with ϕ Re ( ζ , s ) = p n ( ζ ) s n + p n − 1 ( ζ ) s n − 1 + p n − 2 ( ζ ) s n − 2 + · · · + p 0 ( ζ ) , ϕ Im ( ζ , s ) = q n ( ζ ) s n + q n − 1 ( ζ ) s n − 1 + q n − 2 ( ζ ) s n − 2 + · · · + q 0 ( ζ ) , where p i ( ζ ) and q i ( ζ ) are polynomials of ζ with real co efficients. Note that ϕ Re ( ζ , s ) and ϕ Im ( ζ , s ) are not defined for ϕ ( ζ , s ) but for ϕ ( ζ , j s ) in (8). Using these p olynomials, we define the following 2 n × 2 n Sylvester matrix S :=                q n q n − 1 q n − 2 · · · · · · q 0 0 · · · 0 p n p n − 1 p n − 2 · · · · · · p 0 0 · · · 0 0 q n q n − 1 q n − 2 · · · · · · q 0 · · · 0 0 p n p n − 1 p n − 2 · · · · · · p 0 · · · 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . 0 · · · · · · · · · q n q n − 1 q n − 2 · · · q 0 0 · · · · · · · · · p n p n − 1 p n − 2 · · · p 0                (9) and define 2 i -th order principal minors of the matrix b y ∆ i ( ζ ) ( i = 1 , 2 , · · · , n ) . F or example, ∆ 1 ( ζ ) =      q n q n − 1 p n p n − 1      , ∆ 2 ( ζ ) =            q n q n − 1 q n − 2 q n − 3 p n p n − 1 p n − 2 p n − 3 0 q n q n − 1 q n − 2 0 p n p n − 1 p n − 2            and ∆ n = | S | . A ccording to the Hurwitz criterion for complex p olynomials (Theorem 3.4.68 of [36]), all roots of ϕ ( ζ , s ) = 0 lie in the op en left-half complex plane for a giv en ζ if and only if the co efficien ts satisfy the p olynomial inequalities ∆ i ( ζ ) > 0 ( i = 1 , 2 , · · · , n ) . Thus, for eac h fixed v alue of ζ , w e can c heck the stability of the corresp onding F ourier mo de using computationally tractable algebraic conditions. T o analyze the stabilit y of the reaction- diffusion-adv ection system (4), how ev er, w e need to guaran tee the stabilit y for all frequency ζ , which is computationally intractable. As a result, the stability analysis of the reaction- diffusion-adv ection system often resorts to approximation by iterativ ely chec king the sign of ∆ i ( ζ ) for a finite range of discretized v alues of ζ . T o ov ercome this issue, w e in tro duce 9 a computationally tractable condition for certifying non-negativit y of ∆ i ( ζ ) . Prop osition 1. Let M i ∈ R ( ` i +1) × ( ` i +1) b e an y one of the matrices satisfying ∆ i ( ζ ) = z T i M i z i , (10) where ` i := d deg(∆ i ( ζ )) / 2 e and z i := [1 , ζ , ζ 2 , · · · , ζ ` i ] T ∈ R ` i +1 ( i = 1 , 2 , · · · , n ) . Then, the following (i) and (ii) are equiv alent. (i) ∆ i ( ζ ) ≥ 0 for all ζ ∈ R and i = 1 , 2 , · · · , n . (ii) There exists a symmetric matrix N i ∈ R ( ` i +1) × ( ` i +1) suc h that M i + N i  0 and X ( j,k ) ∈ Θ ` ν ( i ) j k = 0 (11) for ` = 2 , 3 , · · · , 2 ` i + 2 and i = 1 , 2 , · · · , n , where ν ( i ) j k is the ( j, k ) -th entry of the matrix N i , and Θ ` := { ( j, k ) ∈ N × N | j + k = `, 1 ≤ j ≤ d `/ 2 e , 1 ≤ k ≤ d `/ 2 e} . (12) The core idea of this prop osition is to find a sum of squares (SOS) decomp osition of ∆ i ( ζ ) (see Supplementary Information for the pro of ). That is, w e certify non-negativit y of ∆ i ( ζ ) b y showing that ∆ i ( ζ ) can b e represented as a sum of non-negative terms. In general, w e can alwa ys find a constan t matrix M i that makes the quadratic form (10) for a giv en p olynomial ∆ i ( ζ ) . Th us, a sufficient condition for the non-negativity of the p olynomial ∆ i ( ζ ) is M i  O . How ever, this do es not constitute a necessary condition since the c hoice of M i is not unique. In other words, w e need to explore all p ossible quadratic expressions of ∆ i ( ζ ) to sho w the necessary condition. The matrix N i satisfying z T i N i z i = 0 is added to M i in the condition (ii) for this purp ose. The condition (ii) is amenable to computational implemen tation since it requires finding a single set of matrices N i ( i = 1 , 2 , · · · , n ) that satisfies (11) instead of verifying the non- negativit y of ∆ i ( ζ ) for all ζ . In fact, the problem of finding N i in (11) can b e reduced to semidefinite programming (SDP) [25, 26], whic h is a class of conv ex optimization program with a linear ob jective function and semidefinite constrain ts. Thus, we can utilize existing SDP solvers such as SeDuMi [37] and SDPT3 [38], which implement in terior p oin t metho ds 10 [39] to efficiently search for the matrix N i . It should b e noted that the linear equality constrain ts in (11) can b e equiv alently transformed to semidefinite constrain ts L  O and L  O with a diagonal matrix L whose entries are the left-hand side of the equalit y constrain ts, i.e., L := diag ( P ( j,k ) ∈ Θ 2 ν (1) j k , P ( j,k ) ∈ Θ 3 ν (1) j k , · · · , P ( j,k ) ∈ Θ 2 ` n ν ( n ) j k ) . Th us, all of the constraints in the condition (ii) can b e reduced to semidefinite conditions. Remark 2. Prop osition 1 can b e also used for the stability analysis of reaction-diffusion- adv ection pro cess in multi-dimensional space Ω . T o see this, we consider m dimensional space Ω , and define a vector of spatial frequency ζ := [ ζ 1 , ζ 2 , · · · , ζ m ] T . Then, the dynamics of multi-dimensional F ourier comp onen ts ˜ c ( ζ , t ) is obtained as d ˜ c ( ζ , t ) dt = ( A − ¯ ζ 2 D + j ¯ ζ V ) ˜ c ( ζ , t ) , (13) where ¯ ζ := P m i =1 ζ i . This leads to the c haracteristic p olynomial that determines the stabilit y of the reaction-diffusion-advection system as ϕ ( ζ , s ) = | sI − A + ¯ ζ 2 D − j ¯ ζ V | . Here, an important implication is that { A − ¯ ζ 2 D + j ¯ ζ V | ¯ ζ := P m i =1 ζ i , ζ i ∈ R } = { A − ζ 2 D + j ζ V | ζ ∈ R } . Thus, ϕ ( ζ , s ) = 0 has all ro ots in the open left-half complex plane for all ζ ∈ R m if and only if ϕ ( ζ , s ) = 0 do es so for all ζ ∈ R . The latter condition b oils down to chec king the non-negativit y of ∆ i ( ζ ) as discussed in this section. Therefore, the matrix inequality conditions in Prop osition 1 can b e used for the stabilit y analysis of m ulti-dimensional reaction-diffusion-adv ection systems. Example. W e demonstrate the SDP based stability test using the reaction-diffusion- adv ection system (3). Let [ C ∗ 1 , C ∗ 2 ] T b e a spatially homogeneous equilibrium p oin t of the system and consider the linearized system (4). It follows from (3) that C ∗ 1 = 1 2 (1 − √ w ) , C ∗ 2 = a 2( a + b ) (1 + √ w ) (14) is an equilibrium p oin t, where w := 1 − 4( a + b ) 2 /a . The Jacobian linearization leads to ∂ ∂ t " c 1 c 2 # = " − a − C ∗ 2 2 − 2 C ∗ 1 C ∗ 2 C ∗ 2 2 − ( a + b ) + 2 C ∗ 1 C ∗ 2 # " c 1 c 2 # + " d 0 0 1 # " ∇ 2 c 1 ∇ 2 c 2 # + " v 1 0 0 v 2 # " ∇ c 1 ∇ c 2 # . (15) T o analyze the stability of the equilibrium p oin t, we substitute the parameters in to (7) and compute the p olynomials ϕ Re ( ζ , s ) and ϕ Im ( ζ , s ) , and define the Sylvester matrix S 11 0.03 0.04 0.05 0.06 0.07 0 0.1 0.2 0.3 0.4 v 0 0.02 0.04 0.06 0.05 0.1 0.15 0.2 0.25 A B C D E F No equilibrium exists Stable Unstable Unstable No equilibrium exists Stable U nstable Stable D E F Position x 0 Time t 20 40 60 80 0 1000 2000 3000 4000 5000 6000 0.5 0.4 0.2 0.3 Position x 0 Time t 20 40 60 80 0 1000 2000 3000 4000 5000 6000 0.5 0.4 0.3 0.2 0.1 Position x 0 Time t 20 40 60 80 0 1000 2000 3000 4000 5000 6000 0.5 0.4 0.3 0.2 0 0.02 0.04 0.06 0.05 0.1 0.15 0.2 0.25 No equilibrium exists Stable Unstable Figure 2: P arameter space analysis for reaction-diffusion-adv ection system (3) . (A, B) P arameter maps sho wing stable and unstable regions without advection ( v = 0 ) and with advection ( v = 0 . 3162 ). (C) Parameter map sho wing the robustification of transp ort- driv en instability . (D-F) Spatio-temp oral profiles of the molecule P for different v alues of parameters (see Fig. 2 (C)). The panels D and F are rep eated from Fig. 1 (B) and Fig. 1 (D), resp ectiv ely . defined in (9). Here we first consider the parameter v alues used in Fig. 1(B), where ( b, v 1 , v 2 ) = (0 . 040 , 0 , 0) (see Metho d section for the other parameters). The corresp onding p olyno- mials ∆ i ( ζ ) ( i = 1 , 2) are then obtained as ∆ 1 ( ζ ) = 7 ζ 2 + 0 . 184 , (16) ∆ 2 ( ζ ) = 294 ζ 8 − 0 . 0382 ζ 6 + 0 . 192 ζ 4 + 0 . 0315 ζ 2 + 0 . 000555 . (17) Hurwitz’s stabilit y condition implies that the equilibrium [ C ∗ 1 , C ∗ 2 ] is stable if and only if ∆ 1 ( ζ ) > 0 and ∆ 2 ( ζ ) > 0 for all ζ ∈ R . It is obvious that ∆ 1 ( ζ ) > 0 , but the sign of ∆ 2 ( ζ ) needs further examination. Thus, we use the condition (ii) of Proposition 1 to derive conditions for non-negativity of ∆ 1 ( ζ ) and ∆ 2 ( ζ ) . The p ositivit y of ∆ 1 ( ζ ) can b e easily confirmed since ∆ 1 ( ζ ) = z T 1 M 1 z 1 with M 1 = diag (0 . 184 , 7)  O and z 1 = [1 , ζ ] T . On the other hand, ∆ 2 ( ζ ) is expressed as ∆ 2 ( ζ ) = z T 2 M 2 z 2 with monomials z 2 := [1 , ζ , ζ 2 , ζ 3 , ζ 4 ] T and M 2 = diag (0 . 000555 , 0 . 0315 , 0 . 192 , − 0 . 0382 , 294) . Thus, we examine the existence of 12 a matrix N 2 suc h that M 2 + N 2  O and P ( j,k ) ∈ Θ ` ν (2) j k = 0 ( ` = 2 , 3 , · · · , 10) b y solving the feasibility problem of the optimization program. The optimization solver yields N 2 =          0 0 − 0 . 0264 0 0 . 0710 0 0 . 0528 0 − 1 . 3655 0 . 0012 − 0 . 0264 0 2 . 5889 − 0 . 0012 − 16 . 797 0 − 1 . 3655 − 0 . 0012 33 . 595 0 0 . 0710 0 . 0012 − 16 . 797 0 0          (18) whic h indeed satisfies M 2 + N 2  O and P ( j,k ) ∈ Θ ` ν (2) j k = 0 ( ` = 2 , 3 , · · · , 10) . This implies that ∆ 2 ( ζ ) = z T 2 ( M 2 + N 2 ) z 2 = z T 2 M 2 z ≥ 0 for all ζ ∈ R . Th us, the equilibrium p oin t is lo cally stable, and the concentrations [ C 1 , C 2 ] T con verge to the equilibrium when they are p erturbed in the vicinity of [ C ∗ 1 , C ∗ 2 ] (Fig. 1(B)). T o see an example of an unstable case, w e next consider the parameter sets used in Fig. 1(C), that is, ( b, v 1 , v 2 ) = (0 . 055 , 0 , 0) (see Metho d section for the other parameters). In this case, ∆ 1 ( ζ ) = 7 ζ 2 + 0 . 0679 , (19) ∆ 2 ( ζ ) = 294 ζ 8 − 19 . 141 ζ 6 − 0 . 0999 ζ 4 + 0 . 0045 ζ 2 + 0 . 000033 . (20) Th us, the sign of ∆ 2 ( ζ ) determines the stability of the equilibrium p oin t. The SDP solver returns “infeasible”, implying that there is no N 2 satisfying the conditions (ii) in Prop osition 1. This means that there exists a spatial frequency ζ for which the characteristic polynomial ϕ ( ζ , s ) = 0 has a ro ot in the right-half complex plane. Thus, the system (6) is unstable for the frequency ζ . In particular, the destabilization o ccurs for some non-zero spatial frequency ζ since ∆ 1 ( ζ ) > 0 and ∆ 2 ( ζ ) > 0 when ζ = 0 . This observ ation is consistent with the sim ulation result in Fig. 1(C) in that the steady state solution con verges to the spatial p erio dic oscillations with some non-negativ e frequency . Using the SDP based stabilit y test, w e further inv estigate ho w the adv ective flo w affects the formation of spatial patterns in the reaction system (3). Specifically , we first char- acterize the parameter space that induces transp ort-driv en instability with and without flo w (Fig. 2 (A) and (B), resp ectiv ely). In Fig. 2(B), the flow rates of P and Q w ere set ( v 1 , v 2 ) = (0 . 3162 d, 0 . 3162) = (1 . 897 , 0 . 3162) , resp ectiv ely . The difference of the flow rate is due to the Stokes-Einstein relation, where the difference of the molecular size results 13 in the diffusion and advection rates. Figures 2(A) and (B) illustrate that the advectiv e flo w broadens the parameter region for instability , implying that the formation of p eriodic patterns b ecomes more robust when there is adv ective flo w in the system. T o further analyze the relation b et w een the velocity of the flo w and the parameter region for transp ort-driv en instability , we define the flo w rate as ( v 1 , v 2 ) = ( dv , v ) and v ary v b et ween 0 and 0.3162. It should b e noted that v is the flow rate of Q and is prop ortional to the velocity of the constan t flo w in the reactor (see Fig. 1(A)). Fig. 2(C) illustrates the parameter region for transp ort-driv en instabilit y for differen t v alues of v , where v = 0 means that the transp ortation of the molecules is exclusiv ely due to passiv e diffusion. The result suggests that the transp ort-driv en instabilit y is induced b y adv ective flo w once the sp eed of the flo w reaches a threshold. In particular, the instability region increases almost linearly with the sp eed of the flow and becomes almost t wice as large as the diffusion only case when v = 0 . 40 . The spatio-temp oral sim ulations of the reaction- diffusion-adv ection system show that p erio dic pattern formation do es o ccur by increasing the sp eed of adv ection (Fig. 2(D-F)). Although the link b et w een diffusion-driv en instability and biological pattern formation is often criticized due to the lack of robustness, these results suggest that the advectiv e transp ortation of molecules could comp ensate for the fragilit y of the diffusion-driven pattern formation. 5 P attern formation with sp ecified spatial profiles In the design and analysis of reaction-diffusion-advection systems, we are often interested in the shap e of spatial patterns, whic h is roughly determined by the spatial frequency , or the characteristic wa v elength, of the spatial oscillations. This requires identification of the destabilizing frequency , that is, the frequency ζ that mak es the p olynomials ∆ i ( ζ ) negativ e. In what follows, we extend Prop osition 1 to enable lo cal stability analysis of reaction-diffusion-adv ection systems for a sp ecified range of spatial frequency I ⊂ R . Using the extended theorem, we can specify a range of spatial frequency I that destabilizes the reaction-diffusion-adv ection system. This facilitates the analysis and design of the spatial profiles of transp ort-driv en patterns. A ccording to the Hurwitz condition presented in Section 4, the reaction-diffusion-advection system (4) is stable for a given range of spatial frequency I ∈ R if and only if the p oly- nomials ∆ i ( ζ ) ( i = 1 , 2 , · · · , n ) is p ositiv e for all ζ ∈ I . T o chec k the sign of ∆ i ( ζ ) for a given rage I , w e again in tro duce a set of linear matrix inequalities that is amenable to 14 semidefinite programming. Prop osition 2. Let I denote an open interv al on real num b ers, i.e., I ⊂ R . The following (i)–(iii) are equiv alent. (i) ∆ i ( ζ ) ≥ 0 for all ζ ∈ I and i = 1 , 2 , · · · , n . (ii) There exist p olynomials f ( ζ ) and g ( ζ ) such that f ( ζ ) ≥ 0 and g ( ζ ) ≥ 0 for all ζ ∈ R , and ∆ i ( ζ ) = f ( ζ ) + h ( ζ ) g ( ζ ) , where h ( ζ ) :=          − ( ζ − ζ )( ζ − ζ ) when I = [ ζ , ζ ] ζ − ζ when I := [ ζ , ∞ ) − ζ + ζ when I := ( −∞ , ζ ] . (21) (iii) The following conditions are satisfied. • In the case of a finite interv al I := [ ζ , ζ ] , define δ ( i ) ` as the co efficients of the p olynomial ∆ i ( ζ − ζ ) ζ + ( ζ + ζ ) 2 ! = X ` δ ( i ) ` ζ ` . ( i = 1 , 2 , · · · , n ) (22) Then, there exists K i ∈ R ( ` i +1) × ( ` i +1) and L i ∈ R ` i × ` i suc h that K i  O , L i  O , (23) δ ( i ) ` = X ( j,k ) ∈ Θ ` +2 ( κ ( i ) j k + λ ( i ) j k ) − X ( j,k ) ∈ Θ ` λ ( i ) j k , (24) where ` i := d deg(∆ i ( ζ )) / 2 e , and κ ( i ) j k and λ ( i ) j k represen t the ( j, k ) -th entry of the matrices K i and L i , resp ectiv ely . Θ ` is defined in Prop osition 1. • In the case of a semi-infinite interv al I := [ ζ , ∞ ) , define δ ( i ) ` as the coefficients of the p olynomial ∆ i  ζ + ζ  = X ` δ ( i ) ` ζ ` . ( i = 1 , 2 , · · · , n ) , (25) 15 or in the case of I := ( −∞ , ζ ] , define δ ( i ) ` as ∆ i  − ζ + ζ  = X ` δ ( i ) ` ζ ` . ( i = 1 , 2 , · · · , n ) . (26) Then, there exists K i ∈ R ( ` i +1) × ( ` i +1) and a matrix L i suc h that K i  O , L i  O , (27) δ ( i ) ` = X ( j,k ) ∈ Θ ` +2 κ ( i ) j k + X ( j,k ) ∈ Θ ` +1 λ ( i ) j k , (28) where ` i := d deg(∆ i ( ζ )) / 2 e , and κ ( i ) j k and λ ( i ) j k represen t the ( j, k ) -th entry of the matrices K i and L i , resp ectiv ely . The size of L i is ` i + 1 by ` i + 1 when deg(∆ i ( ζ )) is o dd, and ` i b y ` i when deg (∆ i ( ζ )) is even. Θ ` is defined in Prop osition 1. It should b e noted that Θ ` is an empty set for ` = 0 , 1 . Similar to Proposition 1, the basic idea is to explore the SOS decomp osition of ∆ i ( ζ ) to guaran tee its non-negativit y on the in terv al I (see Supplementary Information for the pro of ). This problem essen tially b oils do wn to finding f ( ζ ) and g ( ζ ) in the condition (ii). In the condition (iii) of Prop osition 2, w e con vert the condition (ii) into the form of matrix inequalities by using the change of v ariable. That is, the in terv al I is con verted into [ − 1 , 1] and [0 , ∞ ) . This leads to the simple LMI based form ulation (23), (24), (27) and (28), where the existence of the matrices K i and L i is equiv alent to the existence of non-negative p olynomials f ( ζ ) := z T i K i z and g ( ζ ) := z T i L i z such that ∆ i ( ζ ) = f ( ζ ) + h ( ζ ) g ( ζ ) with h ( ζ ) = − ( ζ − 1)( ζ + 1) for the finite in terv al and h ( ζ ) = ζ for the infinite in terv al case. Since the existence problem of the matrices K i and L i can b e implemen ted as a feasibilit y problem of a semidefinite program, it is p ossible to efficien tly explore the matrices satisfying the constraints using existing solv ers. The condition (ii) of Prop osition 2 enables iden tification of the stable and unstable spatial frequency of the reaction-diffusion-adv ection system. This is particularly helpful to find the parameter space of the reaction that leads to the formation of presp ecified spatial patterns as demonstrated in the next example. Remark 3. When m ulti-dimensional space Ω is considered, the v ariable ζ in Prop osition 2 is equal to the sum of the spatial frequency ¯ ζ = P m i =1 ζ i as sho wn in Remark 2. Thus, there can b e m ultiple differen t combinations of frequency ζ i ( i = 1 , 2 , · · · , m ) for a giv en 16 0 0.05 0.1 0.15 0.2 0.25 -1 0 1 2 3 4 2 ( ) 10 -3 Figure 3: V alue of ∆ 2 ( ζ ) . The reaction-diffusion-adv ection system is destabilized at the spatial frequency ζ satisfying ∆ 2 ( ζ ) < 0 . destabilizing ζ . This means that the gro wth rate of frequency comp onen ts ˜ c ( ζ , t ) are equal b et w een all combinations of ζ i that sums to ζ . Therefore, an y of these frequency compo- nen ts could appear as the resulting spatial pattern dep ending on the initial perturbation to the homogeneous steady state ¯ C . Example. W e consider the reaction-diffusion-advection system (3) and its linearized mo del (15). Let the parameters a, b and d b e set as ( a, b, d ) = (0 . 06 , 0 . 04 , 6 . 0) and v = 0 . 3162 , with which the equilibrium is unstable as illustrated in Fig. 2C (mark ed as “F”). As ha ve already seen in the previous section, the matrix inequalities in Prop osi- tion 1(ii) are infeasible since ∆ 2 ( ζ ) in (17) is negative for some ζ ∈ R . This implies that the equilibrium p oin t of the reaction-diffusion-advection system is lo cally unstable, lead- 17 ing to the formation of the spatial pattern in Fig. 2F. T o further analyze, w e in vestigate the v alue of ∆ 2 ( ζ ) as sho wn in Fig. 3. The figure implies that the reaction-diffusion- adv ection system is destabilized around 0 . 12 ≤ ζ ≤ 0 . 2 , whose wa v elength corresp onds to 31 . 4 ≤ 2 π /ζ ≤ 48 . 3 . W e observ e that the w av elength of the p erio dic spatial pattern in Fig. 2F agrees with the destabilizing frequency . Prop osition 2 enables optimization-based analysis of the stabilizing/destabilizing fre- quency ζ without actually computing ∆ i ( ζ ) unlike Fig. 3. As an example, we solv e SDP with the condition (ii) in Prop osition 2 by setting I = [0 , 0 . 1] . Since ∆ 2 ( ζ ) is positive on I , the SDP is feasible and the solver finds K 2 = 10 − 4 ×          2 . 6498 0 . 5060 − 11 . 327 − 0 . 7936 9 . 7421 0 . 5060 9 . 7828 − 8 . 9253 − 12 . 081 10 . 089 − 11 . 327 − 8 . 9253 79 . 539 11 . 580 − 76 . 375 − 0 . 7936 − 12 . 081 11 . 580 15 . 730 − 13 . 547 9 . 7421 10 . 089 − 76 . 375 − 13 . 547 75 . 743           O , L 2 = 10 − 4 ×        1 . 8741 − 1 . 5595 − 4 . 2902 3 . 3575 − 1 . 5595 22 . 153 4 . 7537 − 38 . 461 − 4 . 2902 4 . 7537 15 . 638 − 13 . 547 3 . 3575 − 38 . 461 − 13 . 547 75 . 743         O , whic h verifies non-negativit y of ∆ 2 ( ζ ) for 0 ≤ ζ ≤ 0 . 1 . Similarly , for I = [0 . 25 , ∞ ) , the solv er finds K 2 =          0 . 0039 − 0 . 1259 0 . 1373 0 . 1207 0 . 0830 − 0 . 1259 8 . 2473 − 22 . 293 − 12 . 915 − 10 . 193 0 . 1373 − 22 . 293 139 . 67 52 . 697 75 . 458 0 . 1207 − 12 . 915 52 . 697 264 . 83 166 . 16 0 . 0830 − 10 . 193 75 . 458 166 . 16 294 . 00           O , L 2 =        0 . 4259 − 2 . 9496 1 . 0891 0 . 9705 − 2 . 9496 62 . 322 − 11 . 063 12 . 481 1 . 0891 − 11 . 063 169 . 72 56 . 854 0 . 9705 12 . 481 56 . 854 255 . 68         O . On the other hand, the solver certifies that there is no K i and L i that satisfies conditions 18 (ii) in Prop osition 2 when 0 . 1 < ζ < 0 . 25 . These results imply that ∆ 2 ( ζ ) is p ositiv e for ζ ∈ [0 , 0 . 1] and ζ ∈ [0 . 25 , ∞ ) and is negative for some ζ ∈ (0 . 1 , 0 . 25) . Thus, the destabilizing frequency is identified as 0 . 1 < ζ < 0 . 25 without explicitly computing the v alue of ∆ 2 ( ζ ) . The SDP feasibilit y test in tro duced ab ov e allo ws the identification of the interv al of destabilizing frequency . Using this feature, it is p ossible to further narro w do wn the parameter sets for diffusion/advection-driv en instability (Fig. 2(A-C)) with the additional constraints of the wa velength of spatial patterns. When there is a single connected interv al of destabilizing frequency I , w e can also use a bisection-t yp e algorithm to precisely estimate the destabilizing frequency . Sp ecifically , w e first solve the SOS program in Prop osition 2 for some initial in terv al [0 , ζ ] with sufficiently large ζ , and iterativ ely divide ζ in to a half un til the SOS program in Proposition 2 returns a feasible solution. Once we find a frequency that gives a feasible solution, say ζ ∗ , the lo wer b ound of the destabilizing frequency is b et ween ζ ∗ and 2 ζ ∗ . Th us, w e further solve the SOS program by setting ζ = 1 . 5 ζ ∗ and divide the domain into a half until it b ecomes small enough, whic h conv erges to the low er b ound of the destabilizing frequency . The upp er b ound of the destabilizing frequency can b e found b y essentially the same idea, where w e initially solve the SOS program in Prop osition 2 for [ ζ , ∞ ) with sufficien tly small ζ , and iteratively narro w the domain by [2 ζ , ∞ ) . It should, ho wev er, b e noted that this bisection-t yp e algorithm giv es only the upp er and low er b ounds of destabilizing frequency when there are m ultiple destabilizing frequency in terv als. F or example, ζ 1 and ζ 4 are obtained when there are tw o disconnected interv als of destabilizing frequency I 1 = ( ζ 1 , ζ 2 ) and I 2 = ( ζ 3 , ζ 4 ) with ζ 2 < ζ 3 . 6 Conclusion and Discussion Self-organizing phenomena in spatially in teracting chemical systems hav e b een studied in biology and chemistry for a long time due to its link with biological developmen tal pro cess. In microbiology , it is kno wn that culture conditions may affect the spatial profile of bacterial colonies, and its link with reaction-diffusion systems hav e b een extensiv ely studied (see Chapter 11 of [2] for example). In developmen tal biology , exp erimen tal works revealed that spatio-temp oral molecular patterning in em bryos affects morphology , the structure of organisms [3, 40, 41]. This motiv ates spatio-temp oral control of morphogen concentrations in tissue engineering [42]. In addition to the biological applications, there are p otential engineering applications of molecular based comm unications, where spatially dispersed 19 micro and nano-scale systems comm unicate with eac h other by small diffusible molecules [43]. Therefore, understanding the dynamics of transp ort-driv en pattern formation will not only contribute to biological science but also helps op en up new engineering applications. This pap er has prop osed a no vel computational approac h to analyzing the lo cal stabil- it y/instability of reaction-diffusion-advection systems. The prop osed algebraic conditions in Prop ositions 1 and 2 bypass the iterative stabilit y analysis of individual F ourier com- p onen ts and provide a route to ward the direct c haracterization of the lo cal b ehavior by a single run of mathematical optimization, sp eeding the analysis and syn thesis of self- organizing chemical systems in biology and chemistry . The condition in Prop osition 2 further allows the computation of the destabilizing spatial frequency . This helps detailed quan titative analysis of the spatio-temp oral profile of the self-organizing chemical concen- trations. W e hav e n umerically illustrated the optimization-based analysis pro cess using the auto-catalytic reaction mo del with diffusion and adv ection. F rom a theoretical viewp oint, our developmen t hinges up on the sum-of-squares (SOS) optimization technique [31] to pro ve non-negativity of ∆ i ( ζ ) , whic h implies stability of the reaction-diffusion-advection systems. In the last decade, the SOS optimization was extensiv ely studied in con trol engineering communit y to analyze the stabilit y of nonlinear dynamical systems. Softw are pack ages were developed to facilitate the implemen tation of SOS optimization programs [44]. In general, the existence of SOS decomp osition of a p olynomial implies non-negativit y , but the conv erse is not necessarily the case. In other w ords, not all non-negativ e p olynomials can b e represented as a sum of squares (see [31] for example). In our theoretical developmen t, how ever, we hav e used the fact that the p olynomial inequalities asso ciated with the stabilit y analysis are univ ariate, in which case there alw ays exists an SOS decomp osition if the p olynomial is non-negative, leading to not only the sufficien t but also the necessary condition for the stability of reaction-diffusion- adv ection systems. This is particularly useful for the study of transp ort-driv en pattern formation since the condition can provide rigorous instability certificate of the reaction- diffusion-adv ection systems. It is also w orth mentioning that the p olynomial inequality conditions ∆ i ( ζ ) ≥ 0 ( i = 1 , 2 , · · · , n ) are univ ariate even for m ulti-dimensional reaction- diffusion-adv ection systems as discussed in Remark 2. This is attractive from a viewp oin t of computational burden since, in general, the size of the matrices in volv ed in the SOS programs suc h as M i , N i , K i and L i gro w in a combinatorial fashion with the n umber of v ariables and the degrees of p olynomials, whic h is one of actively studied issues in re- 20 cen t y ears [45 – 47]. A limitation of the prop osed approach is that it can only v erify lo cal prop erties around a given equilibrium p oin t due to the linearization, and thus, it cannot rig- orously guarantee the emergence of inhomogeneous spatial patterns in nonlinear reaction- diffusion-adv ection systems. T o o vercome this issue, an alternative approac h w ould b e to use Lyapuno v’s approach to obtain sufficient conditions for global stabilit y [27 – 30]. A ma jor criticism ab out T uring’s diffusion-driv en instability is that the parameter space for the instability is to o small to achiev e even in engineered chemical systems although spatial pattern formation is quite a common phenomenon in biology . Recent theoretical w orks, on the other hand, presented that sto c hastic intrinsic noise of c hemical reactions in biological cells enhances diffusion-driv en instabilit y , resulting in the increase of the instabilit y parameter regime [48 – 51]. Our transp ort-driv en instability analysis in Fig. 2(A)-(C) suggests that adv ective flow is another factor that broadens the parameter space for transp ort-driven instabilit y . In fact, Fig. 2F sho ws that the originally stable reaction- diffusion system (Fig. 2D) exhibits spatially inhomogeneous patterns b y adding the flo w. These results imply that the robustness of self-organization in nature could b e guaranteed b y circulatory systems in addition to intrinsic noise. The idea of flo w-driven instability was previously presented in [21, 52 – 54], where p eriodic spatial patterns w ere observ ed due to the destabilizing effect of advection. The theoret- ical prediction was also verified with engineered c hemical and biological systems [22, 55]. Unlik e the mo del (3) in this pap er, the previous analyses w ere limited for systems where only one molecule could diffuse and flow while other molecules were immobilized. This assumption can b e easily relaxed with the proposed stabilit y analysis method as it is ca- pable of dealing with arbitrarily many diffusing and flowing molecules in principle despite the complex eigenv alues of the advection op erator. Although this pap er has only shown a simple reaction example to fo cus on the v erification and demonstration of the theoretical dev elopment, it will b e helpful to examine more biologically relev an t mo dels, in future, to understand the underlying mechanisms of the flow-driv en robustification of the spatial pattern formation. The authors also envision that there are p otential extensions of the pro- p osed SOS programs to robust stabilit y analysis and (reaction) design problems, given that these topics are activ ely studied using SOS programs in control systems communit y [56]. 21 Metho d The spatio-temp oral dynamics in Fig. 1(B)-(D) and Fig. 2(D)-(F) w ere simulated with the p eriodic boundary condition using W olfram Mathematica 11.0.1.0 on macOS 10.12.6. These colormaps illustrate the concen tration of P , or C 1 ( x, t ) . In all sim ulations, a = 0 . 06 and d = 6 . 0 were used. The spatial length was L = 30 π . Figure 2(D) and (F) are rep eated from Fig. 1(B) and (D), resp ectiv ely , for the clarity of presentation. F or Fig. 2(E), the parameters were set a = 0 . 06 , b = 0 . 04 , d = 6 . 0 , v 1 = 1 . 518 , v 2 = 0 . 2530 . F or all spatio- temp oral simulations, the initial v alues w ere set C 1 ( x, 0) = 0 . 5+0 . 0025 P 20 k =1 { cos(2 k π x/L )+ sin(2 k π x/L ) } and C 2 ( x, 0) = 0 . 2 + 0 . 0025 P 20 k =1 { cos(2 k π x/L ) + sin(2 k π x/L ) } . The semidefinite programs were run with SeDuMi 1.3 [57] and Y ALMIP to olbox [58] on MA TLAB 2016b. Data, co des and materials The program co des used in this study are a v ailable at GitHub ( https://github.com/ hori- group/SDP_for_flow- driven_stability_analysis ) Comp eting in terests The authors hav e no comp eting in terests. Authors’ con tributions Y.H. conceiv ed of and designed the study . Both authors performed the mathematical deriv ation and implemented the computational co des for optimization and simulations. Both authors drafted the manuscript and gav e final approv al for publication. F unding This w ork was supp orted in part by JSPS KAKENHI Grant Num b er JP18H01464 and 15J09841. 22 References [1] K o c h AJ, Meinhardt H. Biological pattern formation: from basic mechanisms to complex structures. Reviews of Mo dern Physics. 1994;66(4):1481–1507. [2] Murra y JD. Mathematical Biology I I. 3rd ed. Springer; 2003. [3] K ondo S, Miura T. Reaction-diffusion mo del as a framew ork for understanding bio- logical pattern formation. Science. 2010;329(5999):1616–1620. [4] Hu Z, Lutk enhaus J. T op ological regulation of cell division in Esc herichia coli inv olves rapid p ole to p ole oscillation of the division inhibitor MinC under the con trol of MinD and MinE. Molecular Microbiology . 1999;34(1):82–90. [5] Meinhardt H, de Boer P AJ. Pattern formation in Esc heric hia coli: A mo del for the p ole-to-p ole oscillations of Min proteins and the lo calization of the division site. Pro ceedings of National Academ y of Sciences of the United States of Amer- ica. 2001;98(25):14202–14207. [6] Lutk enhaus J. Assem bly Dynamics of the Bacterial MinCDE System and Spatial Regulation of the Z Ring. Annual Review of Bio c hemistry . 2007;76:539–562. [7] Rudner DZ, Losick R. Protein sub cellular lo calization in bacteria. Cold Spring Harb or P ersp ectiv es in Biology . 2010;2:a000307. [8] V ecchiarelli A G, Li M, Mizuuc hi M, Hwang LC, Seol Y, Neuman KC, et al. Mem brane- b ound MinDE complex acts as a toggle switch that driv es Min oscillation coupled to cytoplasmic depletion of MinD. Pro ceedings of National Academ y of Sciences of the United States of America. 2016;113(11):E1479–E1488. [9] Castets VV, Dulos E, Boissonade J, Kepp er PD. Experimental evidence of a sus- tained standing T uring-type nonequilibrium c hemical pattern. Physical Review Let- ters. 1990;64(24):2953–2956. [10] Ouy ang Q, Swinney HL. T ransition from a uniform state to hexagonal and strip ed T uring patterns. Nature. 1991;352:610–612. [11] Y ou L, I I I RSC, W eiss R, Arnold FH. Programmed p opulation control by cell-cell comm unication and regulated killing. Nature. 2004;428(6985):868–871. 23 [12] Basu S, Gerchman Y, Collins CH, Arnold FH, W eiss R. A synthetic multicellular system for programmed pattern formation. Nature. 2005;434(7037):1130–1134. [13] Liu C, F u X, Liu L, Ren X, Chau CKL, Li S, et al. Sequential establishment of stripe patterns in an expanding cell p opulation. Science. 2011;334(6053):238–241. [14] Mo on TS, Lou C, T amsir A, Stan ton BC, V oigt CA. Genetic programs constructed from lay ered logic gates in single cells. Nature. 2012;491(7423):249–253. [15] T amsir A, T ab or JJ, V oigt CA. Robust m ulticellular computing using genetically enco ded NOR gates and c hemical âĂ Ÿ wiresâĂŹ. Nature. 2011;469(7329):212–215. [16] Danino T, Mondragón-Palomino O, T simring L, Hasty J. A synchronized quorum of genetic clo c ks. Nature. 2010;463(7279):326–330. [17] Din MO, Danino T, Prindle A, Sk alak M, Selimkhano v J, Allen K, et al. Synchronized cycles of bacterial lysis for in viv o delivery . 536. 2016;536:81–85. [18] Scott SR, Din MO, Bittihn P , Xiong L, T simring LS, Hasty J. A stabilized microbial ecosystem of self-limiting bacteria using syn thetic quorum-regulated lysis. Nature Microbiology . 2017;2(12):17083. [19] Baumgart L, Mather W, Hasty J. Synchronized DNA cycling across a bacterial p op- ulation. Nature genetics. 2017;49(8):1282–1285. [20] T uring AM. The chemical basis of morphogenesis. Philosophicals T ransactions of the Ro yal So ciet y of London B. 1952;237(641):37–72. [21] Ro vinsky AB, Menzinger M. Chemical instabilit y induced by a differential flow. Phys- ical Review Letters. 1992;69(8):1193–1196. [22] Ro vinsky AB, Menzinger M. Self-organization induced b y the differen tial flow of activ ator and inhibitor. Physical Review Letters. 1993;70(6):778–781. [23] Gierer A, Meinhardt H. A theory of biological pattern formation. Kyb ernetik. 1972;12(1):30–39. [24] Hori Y, Miyazak o H, Kumagai S, Hara S. Co ordinated Spatial Pattern F ormation in Biomolecular Comm unication Netw orks. IEEE T ransactions on Molecular, Biological, and Multi-Scale Communications. 2015;1(2):111–121. 24 [25] Bo yd S, Ghaoui LE, F eron E, Balakrishnan V. Linear matrix inequalities in system and control theory . So ciet y for Industrial and Applied Mathematics; 1994. [26] Bo yd S, V andenberghe L. Conv ex Optimization. Cambridge Univ ersity Press; 2008. [27] Jo v ano vić MR, Arcak M, Sontag ED. A Passivit y-Based Approach to Stabilit y of Spatially Distributed Systems With a Cyclic Interconnection Structure. IEEE T rans- actions on Automatic Con trol. 2008;53:75–86. (sp ecial issue). [28] Arcak M. Certifying spatially uniform b ehavior in reaction-diffusion P DE and com- partmen tal ODE systems. Automatica. 2011;47(6):1219–1229. [29] Shafi SY, Aminzare Z, Arcak M, Sontag ED. Spatial uniformity in diffusively-coupled systems using weigh ted L2 norm contractions. In: Pro ceedings of American Con trol Conference; 2013. p. 5639–5644. [30] V almorbida G, Ahmadi M, P apachristodoulou A. Stabilit y analysis for a class of partial differential equations via semidefinite programming. IEEE T ransactions on Automatic Control. 2016;61(6):1649–1654. [31] P arrilo P A. Semidefinite programming relaxations for semialgebraic problems. Math- ematical Programming. 2003;96(2):293–320. [32] Gra y P , Scott SK. Auto catalytic reactions in the isothermal contin uous stirred tank reactor. Chemical Engineering Science. 1984;39(6):1087–1097. [33] P earson JE. Complex patterns in a simple system. Science. 1993;261(5118):189–192. [34] Strauss W A. P artial differential equations: an in tro duction. 2nd ed. John Wiley & Sons, Inc.; 2007. [35] Bamforth JR, Merkin JH, Scott SK, Tóth R, Gáspár V. Flow-distributed oscil- lation patterns in the Oregonator mo del. Ph ysical Chemistry Chemical Physics. 2001;3(8):1435–1438. [36] Hinric hsen D, Pritchard AJ. Mathematical Systems Theory I. 3rd ed. Springer; 2011. [37] Sturm JF. Using SeDuMi 1.02, a MA TLAB to olb o x for optimization ov er symmetric cones. Optimization Metho ds and Soft ware. 1999;11–12:625–653. 25 [38] T oh KC, T odd MJ, T utuncu RH. SDPT3 — a Matlab softw are pac k age for semidefinite programming. Optimization Metho ds and Soft ware. 1999;11:545–581. [39] F orsgren A, Gill PE, W right MH. Interior Metho ds for Nonlinear Optimization. SIAM review. 2002;44(4):525–597. [40] P almeirim I, Henrique D, Ish-Horo wicz D, Pourquie O. A vian hairy Gene Expression Iden tifies a Molecular Clo c k Linked to V ertebrate Segmentation and Somitogenesis. Cell. 1997;91(5):639–648. [41] Bessho Y, Sak ata R, Komatsu S, Shiota K, Y amada S, Kageyama R. Dynamic expres- sion and essen tial functions of Hes7 in somite segmen tation. Genes and Dev elopment. 2001;15(20):2642–2647. [42] Biondi M, Ungaro F, Quaglia F, Netti P A. Con trolled drug delivery in tissue engi- neering. Adv anced Drug Delivery Reviews. 2008;60(2):229–242. [43] Nak ano T, Suda T, Ok aie Y, Mo ore MJ, V asilakos A V. Molecular comm unication among biological nanomac hines: a la yered architecture and researc h issue. IEEE T ransactions on NanoBioscience. 2014;13(3):169–197. [44] P apachristodoulou A, Anderson J, V almorbida G, Pra jna S, Seiler P , Parrilo P A. SOSTOOLS: Sum of squares optimiza- tion to olbox for MA TLAB. ; 2013. A v ailable from http://www.eng.ox.ac.uk/control/sostools , http://www.cds.caltech.edu/sostools and http://www.mit.edu/˜parrilo/sostools . [45] W aki H, Kim S, K o jima M, Muramatsu M. Sums of squares and semidefinite program relaxations for p olynomial optimization problems with structured sparsity . SIAM Journal on Control and Optimization. 2006;17(1):218–242. [46] Ahmadi AA, Ma jumdar A. DSOS and SDSOS optimization: more tractable alterna- tiv es to sum-of-squares and semidefinite optimization; 2017. [47] W eisser T, Lasserre JB, T oh KC. Sparse-BSOS: a b ounded degree sos hierarc hy for large scale p olynomial optimization with sparsity . Mathematical Programming Computation. 2018;10(1):1–32. 26 [48] Biancalani T, F anelli D, P atti FD. Sto c hastic T uring patterns in the Brusselator mo del. Ph ysical Review E. 2010;81(4):046215. [49] Butler T, Goldenfeld N. Fluctuation-driven T uring patterns. Ph ysical Review E. 2011;84(1):011112. [50] Scott M, Poulin FJ, T ang H. Approximating intrinsic noise in contin uous m ultisp ecies mo dels. Pro ceedings of the Roy al So ciet y A. 2011;467(2127):718–737. [51] Hori Y, Hara S. Noise-induced spatial pattern formation in sto c hastic reaction- diffusion systems. In: Pro ceedings of the IEEE Conference on Decision and Con trol; 2012. p. 1053–1058. [52] Malc how H. Motional Instabilities in Prey-Predator Systems. Journal of Theoretical Biology . 2000;204(4):639–647. [53] Gholami A, Steinbo ck O, Zyk ov V, Bo densc hatz E. Flow-driv en instabilities during pattern formation of Dict yostelium discoideum. New Journal of Physics. 2015;17(6):063007. [54] Vidal-Henriquez E, Zyk ov V, Bo densc hatz E, Gholami A. Conv ective instabilit y and b oundary driv en oscillations in a reaction-diffusion- adv ection mo del. Chaos. 2017;27:103110. [55] Ec kstein T, Vidal-Henriquez E, Bae A, Zyko v V, Bo densc hatz E, Gholami A. Influence of fast advectiv e flows on pattern formation of Dicty ostelium discoideum. PLoS ONE. 2018;13(3):e0194859. [56] P apachristodoulou A, Pra jna S. A tutorial on sum of squares techniques for systems analysis. In: Pro ceedings of American Control Conference; 2005. p. 2686–2700. [57] Sturm JF. Using SeDuMi 1.02, a MA TLAB to olb o x for optimization ov er symmetric cones. Optimization Metho ds and Soft ware. 1999;11–12:625–653. [58] Löfb erg J. Y ALMIP : A T o olbox for Modeling and Optimization in MA TLAB. In: In Pro ceedings of the CACSD Conference. T aip ei, T aiw an; 2004. . 27