Complex oscillatory motion of multiple spikes in a three-component Schnakenberg system
CComplex oscillatory motion of multiple spikes in athree-component Schnakenberg system
Shuangquan Xie ∗ , Theodore Kolokolnikov † and Yasumasa Nishiura ‡ June 11, 2020
Abstract
In this paper we introduce a three-component Schnakenberg model. Its key feature is that ithas a solution consisting of N spikes that undergoes a Hopf bifurcation with respect to N distinctmodes nearly simultaneously. This results in complex oscillatory dynamics of the spikes, not seen intypical two-component models. For parameter values above the Hopf bifurcations, we derive reducedequations of motion which consist of coupled ordinary differential equations (ODEs) of order 2 N forspike positions and their velocities. These ODEs fully describe the slow-time evolution of the spikesnear the Hopf bifurcations. We then apply the method of multiple scales to the resulting ODEs toderive long-time dynamics. For a single spike, we find that long-time motion consist of oscillationsnear the steady state whose amplitude can be computed explicitly. For two spikes, the long-timebehaviour can be either in-phase or out-of-phase oscillations. Both in and out of phase orbits arestable, coexist for the same parameter values, and the fate of orbit depends solely on the initialconditions. Further away from the Hopf bifurcation point, we offer numerical experiments indicatingexistence of highly complex and chaotic oscillations. Keywords—
Activate-substrate-inhibitor system, Coexistence of multiple oscillatory spikes, matched asymp-totics, reduction method, chaotic behaviour.
Nonlinear reaction-diffusion (RD) systems are commonly used to model self-organized phenomena in nature,such as vegetation patterns [1], species invasion phenomena [2], chemical reactions [3, 4], animal skin patterns[5, 6, 7, 8], and morphogenesis [9]. One of the simplest RD models is the Schnakenberg system [10, 11, 12, 13].It is a two-component model of a simplified activator-substrate reaction, and is often used as one of the simplestmodels for studying spike dynamics in reaction-diffusion systems; it is also a limiting case of both the Gray-Scott model [4] and the Klausmeyer model for vegetation [1]. It describes the space-time dependence of theconcentrations of the intermediate products u (the activator) and v (the substrate) in a sequence of reactions,and has the form (cid:26) u t = D u u xx − u + u vv t = D v v xx + A − u v . (1)Here, u and v represent the concentrations of activator and substrate respectively, and D { u,v } is their respectivediffusion coefficients. The activator decays exponentially whereas the substrate v is fed into the system with aconstant rate A . The nonlinear term corresponds to the reaction P + 2 Q → Q. In this paper, motivated by the three-component model of gas-discharge phenomena [14, 15, 16, 17], we introducean “extended” Schankenberg model having the following form: u t = D u u xx − k u − k w + u vθv t = D v v xx + A − u vτ w t = D w w xx + u − k w . (2)The system (2) falls into the category of three-component activator-substrate-inhibitor systems. It has an ad-ditional reactant w acting as an inhibitor to u , which interacts with u and v indirectly through the intake of ∗ WPI-AIMR center, Tohoku university, Japan † Dalhousie University, Candada ‡ WPI-AIMR center,Tohoku university. Research Institute for Electronic Sciences, Hokkaido University, MathAM-OIL,Tohoku University and AIST, Japan a r X i v : . [ n li n . PS ] J un ctivator u . The reactant w can be treated as a out-product that is being removed continuously during thereaction process, see [18]. The “classical” Schnakenberg model corresponds to choosing k = 0 , so that w playsno role in the reaction for u and v. Figure 1: (Color online) Simulation of (3) with ε = 0 . , D = 0 . , κ = 0 . u ( x,
0) = 0 . − x − x ) ) + 0 . − x − x ) ) , v ( x,
0) = 1 , w ( x,
0) = u ( x, . In (a-b), x = x = − . τ is as indicated. There is a Hopf bifurcation at τ ∼ /κ + O ( ε ) . Decaying oscillations in spike position are observed in (a) with sustained oscillations in (b). (c): with τ = 1 . /α, even oscillations are unstable with solution eventually settling into odd oscillation (in-phaseoscillation) after (possibly long) initial transient. (d-e): For τ = 1 . /κ, the system converges to eitherodd or even oscillation depending only on initial conditions. Both even (out-of-phase) and odd (in-phase)oscillations are stable in this case. We are interested in the new dynamics that the extra variable w introduces, as compared to the previous studiesof dynamics in RD systems. We will use τ as a bifurcation parameter, which controls the reaction ratio of u and w . It has been shown in many RD system [19, 20, 21, 22] that oscillatory behaviors of fronts and localizedpatterns are observed when τ varies. We expect a similar behavior in this three-component system. To simplifythe analysis, we will also assume that θ and D w is sufficiently small and can be set to zero. In addition, we willwrite ε = D u and we will assume that ε (cid:28) O ( D v ). Hereafter we use D instead of D v for simplicity. By furtherrescaling, u = ε u, v = εv, w = ε w , normalizing the coefficients k , k , k and A , and dropping the hat, we willconsider the following system as our starting point, u t = ε u xx − (1 − κ ) u − κw + u v Dv xx + − u v/ετ w t = u − w , x ∈ ( − , , t ≥ x = ± . (3)This scaling simplifies the calculation; in particular the Schnakenberg model is a special case correspondingto κ = 0. In this case, the solution is well known to consist of N spikes whose stability and dynamics havebeen extensively studied [23, 10, 12, 24, 25, 26, 13]. The basic steady state consisting of stationary N spikespersists even when τ > . On the other hand, we will show using a simple argument in § τ is increased slightly past τ = 1 /κ, leading to an oscillatory behaviour inspike positions. Moreover, N small eigenvalues (controlling the motion of N spikes) undergo a Hopf bifurcationnearly simultaneously. As a result, a complex interaction between the different modes can be observed, leadingto a coexistence of multiple possible stable orbits. The main goal of this paper is to shed light on this complexbehaviour using finite-dimensional reduction and multiple scales techniques.Figure 1 illustrates this phenomenon. For a single spike, there is a single small eigenvalue that undergoes a Hopfbifurcation for values of τ slightly above 1 /κ . It causes the spike to oscillate periodically. In § mplitude of this oscillation as a function τ. For two spikes, the long-time dynamics are even more interesting.We show § τ slightly above the Hopf bifurcation, the dynamics settle to one of the two possible orbits,corresponding to either odd (in-phase) or even (out-of-phase) oscillations in spike positions (see Figure 1(d) and(e), respectively). Which one is chosen depends both on τ and the initial conditions. When τ = . κ , only evenoscillation is stable (Figure 1(c)). On the other hand, when τ = . κ both even and odd oscillations coexist forthe same parameter values and orbit selection mechanism depends only on initial conditions.The main result of this paper is as follows Principal Result 1.
Let τ = 1 κ + ε ˆ τ and assume that ˆ τ = O (1) as ε → . Then there exists a solution consisting of N spikes nearly-uniformly spaced,but whose centers evolve near the symmetric configurations on a slow time-scale according to the following. Let x k be the center of the k -th spike. Then x k ∼ − k − N + εp k where p k = N (cid:88) i =1 Q ki B i ( ε t ) cos ( εω k t + θ k ) (4) and B k ( s ) solves the ODE system Eq. (96), θ k are constants depending on the initial condition, ω k = (cid:113) − κλ k N ,and λ k is defined by Eq. (15). The article is structured as follows. In §
2, we study the instabilities of multiple spikes pattern triggered byincreasing reaction time ratio τ . It turns out that in the regime τ ∼ κ + ˆ τ ε , only eigenvalues which areasymptotically small as ε → § τ near κ . The reduced system contains 2 N variables correspondingto both positions and velocities of N spikes. In § § In this section, we study the stability of N-spike profiles for system (3). This analysis is an extension of thework done in [10]. In the limit ε →
0, there exist large eigenvalues which tend to a non-zero limit and the smalleigenvalues which tend to zero. We begin by formulating the linear stability problem of steady state. Since theequilibrium solution are exactly the same as two-component system, we follow the conclusion in [10].
Lemma 1. As ε → , the system (3) admit a N-spike solution, whose leading order is given by w = u = 1 ξ N (cid:88) j =1 ρ (cid:16) x − x j ε (cid:17) (5) where x j = − j − N , j = 1 ..N , ρ ( y ) = sech ( y/ is the unique positive solution to ρ yy − ρ + ρ = 0 , ρ (cid:48) (0) = 0 , ρ → as y → ±∞ ; (6) ξ = N (cid:90) ∞−∞ ρ ( y ) dy, and v satisfies Dv xx + − N (cid:80) Nj =1 δ ( x − x j ) = 0 v ( x j ) = ξ v (cid:48) ( −
1) = v (cid:48) (1) = 0 (7)To analyze the stability of the equilibrium solution, we introduce small perturbations u = u s + e λt φ ( x ); v = v s + e λt η ( x ); w = w s + e λt ψ ( x ) . (8)Substituting (8) into the system (3) gives the following eigenvalue problem: λφ = ε φ xx − (1 − κ ) φ + u s η + 2 u s v s φ − κψ (9a)0 = Dη xx − ε − (cid:0) u s η + 2 u s v s φ (cid:1) (9b) τ λψ = φ − ψ (9c) ith Neumann boundary conditions.From (9c), we obtain ψ = φ τ λ . (10)Substituting (10) back to (9a), we obtain: λ (1 − τ κ τ λ ) φ = ε φ xx − φ + u s η + 2 u s v s φ (11a)0 = Dη xx − ε − (cid:0) u s η + 2 u s v s φ (cid:1) (11b)Let γ j and Φ j be the eigenvalues and eigenfunctions of the eigenvalue problem Eqs. (11) at τ = 0. It is ready tosee that the root of λ (1 − τ κ τ λ ) = γ j (12)are the eigenvalues of Eqs. (11). Since we are interested in the role of parameter τ , we required that N-spikessolution is stable with respect to D . Let’s recall the following lemma from [10], Lemma 2.
For N ≥ , let D N := 12 (cid:82) ρ dy N (13) and suppose that O (cid:0) ε (cid:1) (cid:28) D ≤ O (1) . Then for ε (cid:28) , • one-spike solution is stable. • for D < D N , N-spike solution is stable while for
D > D N , N spikes solution is unstable.
Note that the equation (12) is a quadratic equation and γ j is negative when D < D N , a simple calculation showsthat the equation only admit solution with negative real part when τ < | κ + γ j | ≤ κ + γ . (14)We arrive at the following proposition: Proposition 1.
Assume
D < D N for N ≥ , the eigenvalue problem (9) has only eigenvalues with negative realpart when τ < | κ + γ | and N-spikes solutions are stable. When τ = τ := κ + γ , the root of Eq (12) is pure imaginary, indicating that the system is undergo a Hopfbifurcation as τ pass through τ . It has been shown in [10] that, the first N eigenvalues with largest real part isof order ε , whose corresponding eigenfunctions are the translation modes.Define λ := − D and for i = 2 , · · · , N , λ i := − D − N D π ( i − N (cid:32) − DN sin π ( i − N (cid:33) − (15)such that λ i > λ j for i < j . Then γ i = λ i ε N for i = 1 ..N .Define τ N := 1 κ + γ N ∼ κ − λ N ε Nκ + · · · , (16)then τ i > τ j for i > j . When τ > τ N , there exist N unstable modes. It is natural to further investigate whathappens when τ > τ N , especially when all of those modes become unstable. Remark.
Note that all the τ N are of order O ( ε ) distance to the value τ c := κ , if we confine our analysis closeto τ c , the large eigenvalues will remain strictly negative for all small ε , which allows us to focus only on howdifferent translation modes interact. We will carry out multiple time scale analysis to study the dynamics beyond Hopf bifurcations in the next section. Dynamics of multiple spikes near the threshold
In this section, we study dynamical solutions that appear through a Hopf bifurcation at τ c := κ . Near the Hopfbifurcation point, we rewrite the bifurcation parameter τ as τ = τ c + ε ˆ τ . In this regime, we derive a reduceddynamical system describing the spike locations { p k , k = 1 ..N } . Since the velocity of the spike is small close tothe bifurcation point, this can be done by means of a multiple time scale perturbation approach as described inthe framework of front bifurcations [15]. For convenience, we rescale the time ˆ t = εt , after dropping the hat, thesystem (3) becomes: εu t = ε u xx − (1 − κ ) u + u v − κw Dv xx + − u v/ετ εw t = u − w (17)We then derive the reduced system for the motion of spikes in the extended Schnakenberg model. The resultsare compared with the numerical simulation obtained from the full system in § Inner solution : near k -th spike, we expand: u ( x ) = U ( y ) , v ( x ) = V ( y ) , w ( x ) = W ( y ) , y = x − x k − εp k ( t ) ε then Eqs. (17) becomes: − U y εp (cid:48) k = U yy − (1 − κ ) U + U V − κW DV yy + ε − εU V − W y εp (cid:48) k + εW t = τ ( U − W ) (18)We introduce slow time scales: T = εt, T = ε t T = ε t and the following anasatz to facilitate the analysis: UVW = U V W + ε U V W + α k ( t, T , T ) U y + ε r u r v r w + ε R u R v R w + h.o.t (19)To make the expansion of U and W unique, we also demand that (cid:90) ( r u − r w ) U y dy = 0 (20) (cid:90) r u U y dy = 0 (21)Note that the anasatz (19) include a term [0 , , α k U y ] T , it is attributed to an important feature of the linearizedoperator (27) around the steady state. We will explain it later in the order of ε .We expand the system in power of ε and collect terms with equal powers of ε .To leading order, we obtain U yy − (1 − κ ) U + U V − κU V yy U − U (22)It follows that V ( y ) = C ,k ; U ( y ) = ρ ( y ) /C ,k where ρ satisfies (6) and C ,k is a constant to be determined by matching with the outer solution near x k .In the order of O ( ε ) : − ∂p k ∂t U y = U yy − (1 − κ ) U + U V + 2 U V U − κ ( α k U y + W ( y ))0 = DV yy − U V − ∂p k ∂t U y = τ c ( U − α k U y − W ) (23)From (23) we obtain: V = 1 D (cid:90) y −∞ (cid:90) z U V dz + B ,k y + C ,k (24)The constant B ,k , C ,k is determined as follows: B ,k = ∂V ∂y (+ ∞ ) + ∂V ∂y ( −∞ )2 , C ,k = lim y →∞ ( V ( y ) − B ,k y ) (25) e rewrite Eqs. (23) into matrix form: (cid:16) − ∂p k ∂t + κα k (cid:17) U y U V (cid:16) − ∂p k ∂t + α k τ c (cid:17) U y = L U V W (26)where L = ∂ ∂y − (1 − κ ) + 2 U V U − κ ∂ ∂y τ c − τ c (27)The linear operator L is degenerate and its zero eigenvalue has a algebraic multiplicity of two. The system ofeigenfunctions then is incomplete and has to be supplemented with generalized eigenfunctions of the eigenvaluezero. The corresponding eigenfunction and generalized eigenfuncion are G := [ U y , V y , W y ] and P := [0 , , U y ]respectively. Thus for convenience, we include the generalized eigenfunction [0 , , α k U y ] T into our expansion(19).The adjoint operator of L is L † = ∂ ∂y − (1 − κ ) + 2 U V U
20 1 τ c ∂ ∂y − κ − τ c (28)One can directly check that G † = ( U y , , − U y ) T is its eigenfunction and P † = ( κ U y , , T is its correspondinggeneralized eigenfunction such that L † P = G . The conditions (20) and (21) correspond to the orthogonalityconditions [ r u , r v , r w ] · G † = 0 and [ r u , r v , r w ] · P † = 0.We then using the Fredholm alternative to obtain the solvability condition. By projection of O ( ε ) Eqs. (23) onto G † = ( U y , , − U y ), we obtain the first solvability condition:( 1 τ c − κ ) α k (cid:90) U y dy = − (cid:90) U ( y ) ∂V ∂y dy (29)Note that (cid:82) ∞−∞ U (cid:82) y U V dy = 0, Eq. (29) becomes:( 1 τ c − κ ) α k (cid:90) U y dy = − B ,k (cid:90) U ( y ) dy (30)Projection of O ( ε ) Eqs. (23) onto P † = ( κ U y , ,
0) yields a second condition: − ∂p k ∂t (cid:90) U y dy = − B ,k (cid:90) U ( y ) dy − κα k (cid:90) U y dy (31)Note that B ,k = 0 by Eq. (72), Eq. (30) is an identity and Eq. (31) becomes ∂p k ∂t = κα k (32)In the order of O ( ε ), we obtain: − ∂p k ∂T U y = ∂ r u ∂y − (1 − κ ) r u + 2 U V r u + U r v − κr w + U V + 2 U U V D ∂ r v ∂y − U V U − U V (cid:16) − ∂p k ∂T + ∂α k ∂t (cid:17) U y − α k ∂p k ∂t ∂ U ∂y = τ c ( r u − r w ) (33)This system can be solved for r u and r w separately, but we skip this step since only the last equation will beneed in the following. From (33) we obtain: ∂r v ∂y = 1 D (cid:90) y (2 U V U + U V ) dy + C ,k (34)The constant C ,k is determined as follows: C ,k = ∂r v ∂y (+ ∞ ) + ∂r v ∂y ( −∞ )2 (35)Note that (cid:90) ∞−∞ ( U V + 2 U U V ) U y dy = 0 (36) y projection of O ( ε ) Eqs (33) onto G † = ( U y , , − U y ), we obtain: ∂α k ∂t = 13 (cid:90) ∞−∞ U ( y ) ∂r v ∂y dy = C ,k (cid:90) ∞−∞ U ( y ) dy. (37)Projection of O ( ε ) Eqs (33) onto P † = ( κ U y , ,
0) yields: ∂p k ∂T = 13 (cid:90) ∞−∞ U ( y ) ∂r v ∂y dy = C ,k (cid:90) ∞−∞ U ( y ) dy (38)In the order of O ( ε ), we obtain: − ∂p k ∂T U y + ∂r u ∂t = ∂ R u ∂y − (1 − κ ) R u + 2 U V R u + U R v − κR w + U V + 2 U V r u + 2 U V r u + 2 U U r v D ∂ R v ∂y − U r v − U V r u − U V − U U V − ˆ ττ c ∂p k ∂t U y − (cid:16) ∂p k ∂T − ∂α k ∂T (cid:17) U y − α k ∂p k ∂T ∂ U ∂y + ∂r w ∂t = τ c ( R u − R w ) (39)It is easy to solve R v by marching with outer solution. ∂R v ∂y = 1 D (cid:90) y ( U r v + 2 U V r u + U V + 2 U U V ) dy + C ,k (40)The constant C ,k is determined as follows: C ,k = 12 (cid:18) ∂R v ∂y (+ ∞ ) + ∂R v ∂y ( −∞ ) (cid:19) (41)Note that (cid:90) ∞−∞ U y ( U V + 2 U V r u + 2 U V r u + 2 U U r v ) dy = 0 (42)By projection of O ( ε ) Eqs. (39) onto G † = ( U y , , − U y ), we obtain: (cid:90) ˆ ττ c U y ∂p k ∂t − ∂α k ∂T U y + ∂p k ∂T ∂ U ∂y U y + ( ∂r u ∂t − ∂r w ∂t ) U y dy = − (cid:90) U ( y ) ∂R v ∂y dy (43)Projection of O ( ε ) Eqs. (39) onto P † = ( κ U y , ,
0) yields: − ∂p k ∂T (cid:90) U y dy + (cid:90) ∂r u ∂t U y dy = − (cid:90) U ( y ) ∂R v ∂y dy (44)We use Eq. (20) to find0 = ∂∂t (cid:90) ( r u − r w ) U y dy = (cid:90) ∂ ( r u − r w ) ∂t U y − ∂p k ∂t ( r u − r w ) U yy dy (45)and, finally, using the last equation of system (33) , and Eq. (31) to obtain: (cid:90) ∂ ( r u − r w ) ∂t U y dy = − τ c κ α k (cid:90) U yy dy (46)Similarly, using orthogonal condition (21), we conclude (cid:90) ∂r u ∂t U y dy = (cid:90) r u ∂p k ∂t U yy dy = 0 (47)Making use of Eq. (46) to eliminate r u and r w in Eq. (43), we obtain ∂α k ∂T = ˆ ττ c κα k + C ,k (cid:90) U ( y ) dy − τ c κ α k (cid:82) U yy dy (cid:82) U y dy (48) ∂p k ∂T = C ,k (cid:90) U ( y ) dy (49)Hence, from Eq. (32), Eq. (37) and Eq. (48), we obtain our main result for the leading order dynamics of the k -th spike. (cid:40) ∂p k ∂t = κα k∂α k ∂t = C ,k (cid:82) ∞−∞ U ( y ) dy (50) ∂p k ∂T = C ,k (cid:82) ∞−∞ U ( y ) dy ∂α k ∂T = ˆ ττ c κα k + C ,k (cid:82) U ( y ) dy − τ c κ α k (cid:82) U yy dy (cid:82) U y dy (51) Outer solution . We determine the constant in the inner region B ,k , C i,k by matching the inner solution andouter solution. Away from spike centers, u ( x ) is assumed to be exponentially small so that Dv xx + = 0 for x (cid:54) = x k . Near x k + √ εp k the term u vε in (17) acts like a delta function, we then obtain: Dv xx + 12 = N (cid:88) j =1 s j δ ( x − x j − εp j ) (52)Here, the weights s j are defined by s k = (cid:90) ∞−∞ U V dy = (cid:90) U V dy + ε (cid:90) ( U V + 2 U V U ) dy + · · · = (cid:82) ρ dyC ,k − C ,k (cid:82) ρ dyC ,k + · · · = 6 C ,k − ε (cid:32) C ,k C ,k + E (cid:33) + · · · (53)where U and V is the inner solution near the k -th spike, E is a constant independent of B ,k and C ,k and E = (cid:90) ∞−∞ U (cid:90) y −∞ (cid:90) z U V dxdzdy + (cid:90) ∞−∞ U V U dy. (54)Integrating (52), we obtain (cid:90) − dx = N (cid:88) j =1 s j (55)The solution of (52) is then given by v ( x ) = N (cid:88) j =1 s j G ( x, x j + εp j ) + ¯ v (56)Where G is the green’s function satisfies: DG xx + 12 = δ ( x − z ) (57) G x ( −
1) = G x (1) = 0 , (cid:90) − Gdx = 0 (58)We can decompose G ( x, z ) as follows G ( x, z ) = 12 D | x − z | + H ( x, z ) (59)where H is the regular part of G . We define matrix G as G = ( G ( x i , x j )) , (60)Let’s denote the ∂∂x j as ∇ x i , When i (cid:54) = j , we can define ∇ x i G ( x i , x j ) and ∇ x j G ( x i , x j ) in the classical way. When i = j , we define ∇ x i G ( x i , x i ) := ∂∂x (cid:12)(cid:12) x = x i H ( x, x i ) . (61)Then we define the derivative of matrix G as follows. ∇G = ( ∇ x i G ( x i , x j )) . (62)Near the the k th spike x = x k + ε p k + εy , we have v ( x ) = N (cid:88) j =1 s j G ( x k + εp k + εy, x j + εp j ) + ¯ v (63) atching: we match Eq. (63) with the inner region. In the order of O (1), we obtain N (cid:88) j =1 C ,k G ( x k , x j ) + ¯ v = C ,k (64)On the other hand, the leading order of Eq. (55) implies N (cid:88) k =1 C ,k = 1 (65)we solve Eq. (64) and Eq. (65) to obtain C ,k = 16 N (66)Matching the order of O ( ε ), we obtain N (cid:88) j =1 (cid:32) − C ,j C ,k + E (cid:33) G ( x k , x j ) + N (cid:88) j =1 C ,j (cid:0) ∇ x k G ( x k , x j ) p k + ∇ x j G ( x k , x j ) p j (cid:1) + ¯ v = C ,k (67)The order O ( ε ) of Eq. (55) reads, N (cid:88) j =1 (cid:32) − C ,j C ,k + E (cid:33) = 0 (68)Note that by a direct computation we have N (cid:88) j =1 ∇ x k G ( x k , x j ) = 0 , N (cid:88) k =1 ∇ x k G ( x k , x j ) = 0 , ∇ x k G ( x k , x j ) = ∇ x k G ( x j , x k ) (69)Eq. (67), Eq. (68) and Eq. (69) implies ¯ v = − N (cid:88) j =1 EG ( x k , x j ) (70)Hence C = C , · · · C ,N = 1 N ( I + 16 N G ) − ( ∇G ) T p · · · p N (71)Another two constants B ,k , C ,k and C ,k depend on the derivative of v ( x ) near ( x k + εp k ) as follows: B ,k = 12 (cid:18) ∂v k, (0 + ) ∂y + ∂v k, (0 − ) ∂y (cid:19) = k (cid:88) j =1 C ,j ∇ x k G ( x k , x j ) = 0 (72) C ,k = 12 (cid:18) ∂v k, (0 + ) ∂y + ∂v k, (0 − ) ∂y (cid:19) = N (cid:88) j =1 (cid:20) C ,k (cid:0) ∇ x k ∇ x k G ( x k , x j ) p k + ∇ x j ∇ x k G ( x k , x j ) p j (cid:1)(cid:21) − C ,j N ∇ x k G ( x k , x j )= − D p k − N (cid:88) j =1 C ,j N ∇ x k G ( x k , x j ) (73)Note that the third derivatives of v in outer region is 0, we obtain C ,k = 12 (cid:18) ∂v k, (0 + ) ∂y + ∂v k, (0 − ) ∂y (cid:19) = 0 (74)Hence C = C , · · · C ,N = − D I − N ∇G ( I + 16 N G ) − ( ∇G ) T p · · · p N (75)Define M = − D I − N ∇G ( I + 16 N G ) − ( ∇G ) T (76) ubstituting all the constants back into the Eqs. (50) and recall that (cid:82) ρ dz (cid:82) ρ z dz = 6 , (cid:82) U yy dy (cid:82) U y dy = 57 , we obtain the reduced dynamic system for spike locations, as formally stated in the following proposition: Proposition 2.
Assume that ε (cid:28) and τ = τ c + ε ˆ τ . Then, the equations for p k and α k are approximatelygoverned by dp k dt = ∂p k ∂t + ε ∂p k ∂T + · · · = κα k + ε C ,k N (77a) dα k dt = ∂α k ∂t + ε ∂α k ∂T + · · · = C ,k N + ε (cid:18) ˆ τ κ α k − κα k (cid:19) . (77b) where C ,k is the k th element of C defined by Eq (75). The 2 N -dimensional ODEs (77) describe the motion of N spike solution observed in the PDEs (17) when thespikes are sufficiently close to the equilibrium and move slowly, with x k + εp k being the locations of the k th spike. Remark.
The terms C ,k (see Eq. 75), which are related to the Green’s function, serve as weak interactionsbetween the spikes even through they are far away from each other. In this section, we apply method of multiple time analysis for the reduced system (77). Numerical comparisonbetween ODEs and partial differential equations (PDEs) simulation is provided to validate our approximationfor the case of N = 1 , , In the case of one spike, after discard the high order term in the system (77), we obtain (cid:26) dpdt = κα − ε p Ddαdt = − p D + ε (cid:0) ˆ τ κ α − κα (cid:1) . (78)Numerical simulation of ODEs (78) and original PDE shows good agreement, see Fig. 2. The existence of smallorder nonlinear term in ODEs (78) make a further approximation possible. In fact, we can approximate (77) bythe following single equation of p up to O ( ε ).¨ p + κp D − ε (cid:20) (ˆ τ κ − D ) ˙ p − κ ˙ p (cid:21) = 0 (79)The system resembles a linear oscillator with weakly nonlinear damping. We proceed to use the multiple timescale analysis to construct uniformly valid approximations to the solution of Eq. (79). We define slow time scaleas T = εt, (80)and seek for a solution of the form: p = q ( t, T ) + εq ( t, T ) + · · · (81)Substituting expansion (81) into the Eq. (79) and separating at each order in ε , we get the sequence of problem, O (1) ∂ q dt + κq D = 0 (82) O ( ε ) ∂ q dt + κq D = − ∂ ∂t∂T q + (cid:20) (ˆ τ κ − D ) ˙ q − κ ˙ q (cid:21) (83)The O (1) solution is q = A ( T ) e iωt + A ∗ ( T ) e − iωt (84)where ω = κ D . Substituting Eq.(84) into the Eq.(83), we obtain ∂ q dt + κq D = (cid:20) − iω ∂ A ∂T e iωt + i ˜ τ κ ω A e iωt − κ (cid:16) − iω A e iωt + 3 A A ∗ iω e iωt (cid:17)(cid:21) + c.c (85) here c.c means the complex conjugate of the term inside the square bracket. To remove secular terms at order O ( ε ), the condition ∂ A ∂T = 12 (ˆ τ κ − D ) A − D A ∗ A (86)must be satisfied. Solving explicitly for A by setting A = B e iθ , we obtain separate equations for the amplitudeand phase of A , ∂ B ∂T = 12 (ˆ τ κ − D ) B − D B (87a) ∂θ∂T = 0 (87b)It is easy to see that the phase θ is a constant determined by the initial condition and the amplitude will convergeto (cid:0) ˆ τ κ − D (cid:1) as time approach infinity, that islim T →∞ B = max (cid:18) (cid:18) ˆ τ κ − D (cid:19) , (cid:19) . (88)A direct comparison between the PDE, ODE and our asymptotic prediction (88) of amplitude is shown in Fig. 3to validate our result. Figure 2: A direct comparison of the spike center location between PDE simulation of system (17),ODE simulation of system (78) and amplitude evolution from DE. (87a). Parameters are ˆ τ = 100 , ε =0 . , D = 0 . , κ = 0 .
2. The center of the spike oscillates as time goes on.
In this subsection, we use the same method to construct approximate solution to N-spike dynamics. In general,the dynamic system of N-spike motion can be written as d p dt ... d p N dt = κ M N p ... p N + ε (ˆ τ κ + M N ) ∂p ∂t ... ∂p N ∂t − κ ∂p ∂t ... ∂p N ∂t ◦ = 0 (89)where ∼ ◦ is the Hadamard power symbol and M is defined by Eq. (76).Define Q = ( q , · · · , q N ) (90) ε lim t →∞ B ) v.s ˆ τ between PDE simulation of (17), ODEsimulation of (78) and amplitude prediction from Eq. (88). Parameters are ε = 0 . , D = 0 . , κ = 0 . τ increases, the red circle is zoomed out in the rightfigure. The agreement are pretty well even when ˆ τ is big. where q = (cid:114) N (1 , − , , · · · , ( − N +1 ) (cid:48) (91) q i = ( q i, , · · · , q i,N ) (cid:48) , i = 2 , · · · , N (92) q i,j = (cid:114) N sin (cid:18) π ( j − N ( i −
12 ) (cid:19) (93)Then it is easy to verify that Q T Q = I, Q T M Q = Λ := λ . . . λ N (94)where λ i is defined by Eq. (15).Define ξ = Q T p , so that Q T M p = Λ Q T p = Λ ξ . Then ξ satisfies:¨ ξ = κ N Λ ξ + ε (cid:20) (ˆ τ κ + Λ3 N ) ˙ ξ − κ Q T ( Q ˙ ξ ) ◦ (cid:21) (95)Following the same procedure as in one spike case, the general form for the amplitude evolution of | ξ | is of theform: ∂ B i ∂t = B i (cid:34)
12 (ˆ τ κ + λ i N ) + 556 N (cid:88) j =1 b ij B j (cid:35) (96)where b ij = N (cid:16)(cid:80) Nk =1 Q kj (cid:17) λ i j = i N (cid:16) (cid:80) Nk =1 Q kj Q ki (cid:17) λ j j (cid:54) = i (97)Note that when D < D N , all the eigenvalues λ i are negative, which results in b ij <
0, thus the sign of term(ˆ τ κ + λ i N ) will determine whether the system (96) admit non-zero equilibrium points. We will further determinewhether non-zero equilibrium points are stable or not in the next subsection for N = 2 and 3. Remark.
When ˆ τ > − λ i Nκ , which corresponds to τ > τ N in the original variable, the zero equilibrium points ofthe system (96) becomes unstable and the system (96) admit at least i different non-zero equilibrium points. Inthis way, we recover the instability result in Proposition 1. In this subsection, we find and classify the fixed points of reduced ODE system in case of N = 2. When N = 2,the constants in Eqs. (96) are evaluated as: q = √
22 ( − , (cid:48) ; q = √
22 (1 , (cid:48) (98) = − D (99) λ = − D (cid:18) D (48 D − (cid:19) (100) b = (cid:18) λ λ λ λ (cid:19) . (101)Eq. (96) becomes: ∂ B ∂t = B (cid:34)
12 (ˆ τ κ + λ N ) + 556 (cid:88) j =1 b ij B j (cid:35) (102a) ∂ B ∂t = B (cid:34)
12 (ˆ τ κ + λ N ) + 556 (cid:88) j =1 b ij B j (cid:35) (102b)with ξ k ∼ B k cos( εω k t + θ k )and p = 1 √ ξ + ξ ) , p = 1 √ ξ − ξ ) . The equilibrium points satisfy: B (cid:34)
12 (ˆ τ κ + λ (cid:88) j =1 b ij B j (cid:35) = 0 (103a) B (cid:34)
12 (ˆ τ κ + λ (cid:88) j =1 b ij B j (cid:35) = 0 (103b)The system decouples when one of B or B is zero. When B (cid:54) = 0 and B = 0, the corresponding solutionexhibits “even” oscillations (out-of-phase oscillations), such as shown in Figure 1(e). When B (cid:54) = 0 and B = 0,the corresponding solution exhibits “odd” oscillations (in-phase oscillations), such as shown in Figure 1(d).Fig. 5 shows a numerical simulation of the ODEs (96) when N = 2 for various values of the bifurcation parameterˆ τ . The results are in very good agreement with the spike dynamics of the original PDE system w.r.t the amplitudeand frequency.By a conventional linear stability analysis, we obtain the following results w.r.t the equilibrium points and theirstability: • When ˆ τ < − D D (1 − D ) κ , Eqs. (103a) only admit one non-negative solution( B , B ) = (0 , , (104)and it is stable. • When − D D (1 − D ) κ < ˆ τ < Dκ , Eqs. (103a) admit two non-negative solutions( B , B ) = (0 ,
0) and , (cid:115) τ κ + λ )5 b . (105)It is easy to check that (0 ,
0) is unstable and (cid:32) , (cid:114) τκ + λ )5 b (cid:33) is stable. • When Dκ < ˆ τ < − D D (1 − D ) κ , Eqs. (103a) admit non-negative solutions( B , B ) = (0 , , (cid:115) τ κ + λ )5 b , , (cid:115) τ κ + λ )5 b , . (106)Again, (0 ,
0) is unstable, (cid:32) , (cid:114) τκ + λ )5 b (cid:33) is stable, and (cid:18) τκ + λ )5 b , (cid:19) is unstable . τ of Eqs. (102). The horizon axis is ˆ τ and the verticalaxis is ||B|| = B + B . Solid line is stable part and dash line is unstable part. BP indicate bifurcationbranch point. Parameters are D = , κ = 0 .
2. The three branch point are at ˆ τ = 82 . , . , . . < τ < .
28, there exist two types of oscillation, but only in-phase oscillations are stable. The coexistence of in-phase and out-of-phase oscillations happens when τ > . • when ˆ τ > − D D (1 − D ) κ , Eqs. (103a) admit four non-negative solutions( B , B ) = (0 , , (cid:115) τ κ + λ )5 b , , (cid:115) τ κ + λ )5 b , , (cid:32) (cid:114) b − (ˆ τ κ + λ , ˆ τ κ + λ (cid:33) . (107)Only (cid:32)(cid:114) τκ + λ )5 b , (cid:33) , (cid:32)(cid:114) τκ + λ )5 b , (cid:33) are stable.Fig. (4) is the bifurcation diagram for ˆ τ obtained via mathcont, [27]. The system has three branch pointsˆ τ = − D D (1 − D ) κ , Dκ , and Dκ . The first two branch points are the same as Hopf bifurcation pointsobtained from stability analysis of PDE in section §
2. The third branch point is the critical one to determinewhether both in-phase and out-of-phase oscillation can be observed or not at the same time. Below this point,only in phase oscillation is stable.
In this subsection, we numerically investigate ODEs and amplitude evolution DEs (96) for three spike dynamics.When N = 3, the constants in Eqs. (96) are evaluated as: q = 1 √ , − , (cid:48) ; q = (cid:114)
23 ( 12 , ,
12 ) (cid:48) ; q = (cid:114)
23 ( √ , , − √
32 ) (cid:48) (108) λ = − D (109) λ = − D − D (cid:18) − D (cid:19) − (110) λ = − D − D (cid:18) − D (cid:19) − (111)We choose D = , and obtain the bifurcation diagram, Fig. 6. Then we choose several special values of ˆ τ to validate our asymptotic result, see Fig. 7. The long term behaviour of amplitude and period from the PDEsimulation and reduced ODE simulation are in good agreement. Note that the amplitude from PDE and ODEoscillate around the amplitude from asymptotic equation in a small scale. A higher order approximation isrequired to further capture this small deviation. a) ˆ τ = 50 (b) ˆ τ = 500 (c) ˆ τ = 1000 (d) ˆ τ = 1000 Figure 5: (Color online) The first row shows the position of the spike center simulated by PDE (solidline) and ODE (dash line). The first column (a) is the simulation at ˆ τ = 50, when only ( B , B ) = (0 , |B | and |B | decays. The second column (b) is the simulationat ˆ τ = 500, although there exist three equilibrium point but only the one ( B , B ) = (0 , . B , B ) = (1 . , B , B ) = (0 , . τ = 1000,both ( B , B ) = (2 . ,
0) and (0 , . We have presented an extension of the Schnakenberg model, similar to the three-component gas discharge system[14, 15, 16]. For a solution consisting of N spikes, this extension exhibits for N distinct and nearly simultaneousHopf bifurcations in the spike positions, which leads to very complex oscillatory dynamics, many of which cannotbe observed in the usual two-component reaction-diffusion models. We have analysed the spike motion near theonset of oscillatory dynamics by first, deriving a system of 2 N ODEs for spike positions and their velocities, andsecond, using Multiple Scales techniques to further elucidate the dominant dynamics near the N-fold bifurcations.There are several differences between our model and gas discharge model, which make the dynamics of localizedpatterns different. First, unlike the gas discharge model, which models one activator and two inhibitors, oursystem belongs to a class of activator-substrate-inhibitor RD systems. Second, the interactions of pulses in gasdischarge system exponentially decay with respect to their distance. While the spikes in our system interactthrough the component ( v ) that do not localize, which is much stronger than the exponentially slow weakinteraction. Third, the motion of pulses in gas discharge system is driven through the drift bifurcation, while ourmoving spikes are triggered by the Hopf bifurcation.The reduction techniques (PDE → ODEs) are related to those used in a series of papers on the gas dischargemodel [14, 15, 16, 17]. Similarly, our analysis is only valid near the multi-Hopf bifurcation points. While thisis a rather limited parameter regime, it allows for a more complete description of the dynamics, including theuntangling of the complex interaction between the simultaneous oscillatory modes using multiple scales analysis.In a broader context, oscillatory localized patterns in RD system have been intensively studied. Spike oscillationshave been observed previously in the two-component RD systems such as Gray-Scott (of which Schnakenberg is alimiting case) [25, 28, 29, 23]. In these works, it was found that even oscillations is the dominant behaviour whentwo spikes oscillate (also the so-called breather oscillations). The oscillatory instability of single front solution for τ of Eqs (96) in three-spike case. the horizon axis is ˆ τ and the verticalaxis is ||B|| s = 2 B + B + 3 B . Solid line is stable part and dash line is unstable part. The red dot BP indicate bifurcation branch point. Parameters are D v = , κ = 0 .
2. The branch point from left toright are at ˆ τ = 72 . , . , . , . , . , . , .
90 respectively. A surprising resultis that three-spike pattern can at most have two stable oscillatory states. two component activator-inhibitor model was first studied analytically by [19]. For two-layer case, there occurs acompetition between in-phase and out-of-phase oscillations. The selection mechanism between them on a finiteinterval was discussed in [22]. Note that coexistence of two phases was not observed in [22]. In [21, 20] it wasshown that for a certain large class of RD systems, even oscillations dominate the dynamics. By contrast, forsystem (3) we found multiple coexisting periodic orbits, supporting both odd and even oscillations.Many open questions remain. Among them is to study the doubly-reduced ODE system (96) for N ≥ N = 2 spikes shows a rich bifurcation structure as well as coexistence and multiplefrequency oscillations. The question of multiple spikes ( N ≥
3) remains open. Further away from the Hopfbifurcation, there is a zoo of interesting dynamics. Some of these are shown in Fig. 8. In particular, for largerdomain sizes we observe “chaotic” dynamics and spike creation-destruction cycles. A similar creation-destructionprocess and chaotic motion for the two-component Gray-Scott model was analyzed numerically in [30]. See alsoa brief survey [31] in this direction. For moderate τ values, multiple oscillatory modes are seen to coexist leadingto rich dynamics. For larger values of τ, “zigzag” spike motion dominates. It would be very interesting to derivethe reduced equations of motion in this regime, far from the bifurcation points.Very rich dynamics are observed in two or higher dimensions, even for a single spike. Figure 9 shows complex“flower” orbits for a single spot in a square domain. Some of these are reminiscent of the trace of a meanderingtip of a spiral wave [32]. It is a completely open question to analyse these; however see [26] for analysis ofsimple (circular) orbit in the two-dimenisonal Schnakenberg model inside a disk. Also, It was shown in [33] thatrotational motion of spot can emerge as a result of a combination of drift and peanut instabilities.In conclusion, three-component RD systems exhibit very rich oscillatory spike motion. Combining PDE → ODEreduction and multiple scales techniques sheds light on the long-time behaviour of spikes in these systems.
References [1] Christopher A Klausmeier. Regular and irregular patterns in semiarid vegetation.
Science , 284(5421):1826–1828, 1999.[2] Jessica B McGillen, Eamonn A Gaffney, Natasha K Martin, and Philip K Maini. A general reaction–diffusionmodel of acidity in cancer invasion.
Journal of mathematical biology , 68(5):1199–1224, 2014.[3] J Schnakenberg. Simple chemical reaction systems with limit cycle behaviour.
Journal of theoretical biology ,81(3):389–400, 1979.[4] John E Pearson. Complex patterns in a simple system.
Science , 261(5118):189–192, 1993.[5] Shigeru Kondo, Rihito Asai, et al. A reaction-diffusion wave on the skin of the marine angelfish pomacanthus.
Nature , 376(6543):765–768, 1995.[6] RA Barrio, C Varea, JL Arag´on, and PK Maini. A two-dimensional numerical study of spatial patternformation in interacting turing systems.
Bulletin of mathematical biology , 61(3):483–505, 1999. a) ˆ τ = 800 (b) ˆ τ = 1200 (c) ˆ τ = 3000 (d) ˆ τ = 3000 Figure 7: The columns are the simulation at ˆ τ = 800 , , τ = 800 , τ = 3000. Note that column 1 and column 2 looks the same, but they locate at different branches, theyellow branch and sushi green branch respectively, which differ by the value of ξ . [7] Philip K Maini, Thomas E Woolley, Ruth E Baker, Eamonn A Gaffney, and S Seirin Lee. Turing’s modelfor biological pattern formation and the robustness problem. Interface focus , page rsfs20110113, 2012.[8] LJ Shaw and JD Murray. Analysis of a model for complex skin patterns.
SIAM Journal on AppliedMathematics , 50(2):628–648, 1990.[9] Debbie L Benson, Jonathan A Sherratt, and Philip K Maini. Diffusion driven instability in an inhomogeneousdomain.
Bulletin of mathematical biology , 55(2):365–384, 1993.[10] David Iron, Juncheng Wei, and Matthias Winter. Stability analysis of turing patterns generated by theschnakenberg model.
Journal of mathematical biology , 49(4):358–390, 2004.[11] Juncheng Wei and Matthias Winter. Stationary multiple spots for reaction–diffusion systems.
Journal ofmathematical biology , 57(1):53–89, 2008.[12] Theodore Kolokolnikov, Michael J Ward, and Juncheng Wei. Spot self-replication and dynamics for theschnakenburg model in a two-dimensional domain.
Journal of nonlinear science , 19(1):1–56, 2009.[13] Theodore Kolokolnikov and Juncheng Wei. Pattern formation in a reaction-diffusion system with space-dependent feed rate.
SIAM Review , 60(3):626–645, 2018.[14] CP Schenk, M Or-Guil, M Bode, and H-G Purwins. Interacting pulses in three-component reaction-diffusionsystems on two-dimensional domains.
Physical Review Letters , 78(19):3781, 1997.[15] M Bode, AW Liehr, CP Schenk, and H-G Purwins. Interaction of dissipative solitons: particle-like behaviourof localized structures in a three-component reaction-diffusion system.
Physica D: Nonlinear Phenomena ,161(1-2):45–66, 2002.[16] M Or-Guil, M Bode, CP Schenk, and H-G Purwins. Spot bifurcations in three-component reaction-diffusionsystems: The onset of propagation.
Physical Review E , 57(6):6432, 1998.[17] Takashi Teramoto, Kei-Ichi Ueda, and Yasumasa Nishiura. Phase-dependent output of scattering processfor traveling breathers.
Physical Review E , 69(5):056224, 2004. [18] Hitoshi Mahara, Nobuhiko J Suematsu, Tomohiko Yamaguchi, Kunishige Ohgane, Yasumasa Nishiura,and Masatsugu Shimomura. Three-variable reversible gray–scott model.
The Journal of chemical physics ,121(18):8968–8972, 2004.[19] Yasumasa Nishiura and Masayasu Mimura. Layer oscillations in reaction-diffusion systems.
SIAM Journalon Applied Mathematics , 49(2):481–514, 1989.[20] Shuangquan Xie and Theodore Kolokolnikov. Oscillations of many interfaces in the near-shadow regime oftwo-component reaction-diffusion systems.
Discrete & Continuous Dynamical Systems-Series B , 21(3):959–975, 2016.[21] Rebecca McKay, Theodore Kolokolnikov, and Paul Muir. Interface oscillations in reaction-diffusion systemsabove the hopf bifurcation.
Discrete & Continuous Dynamical Systems-Series B , 17(7), 2012.[22] Tsutomu Ikeda and Yasumasa Nishiura. Pattern selection for two breathers.
SIAM Journal on AppliedMathematics , 54(1):195–230, 1994.[23] Theodore Kolokolnikov, Michael J Ward, and Juncheng Wei. The existence and stability of spike equilibriain the one-dimensional gray–scott model on a finite domain.
Applied mathematics letters , 18(8):951–956,2005.[24] Michael J Ward and Juncheng Wei. The existence and stability of asymmetric spike patterns for theschnakenberg model.
Studies in Applied Mathematics , 109(3):229–264, 2002.[25] Arjen Doelman, Tasso J Kaper, and Wiktor Eckhaus. Slowly modulated two-pulse solutions in the gray–scottmodel i: Asymptotic construction and stability.
SIAM Journal on Applied Mathematics , 61(3):1080–1102,2000.[26] Shuangquan Xie and Theodore Kolokolnikov. Moving and jumping spot in a two-dimensional reaction–diffusion model.
Nonlinearity , 30(4):1536, 2017.[27] Annick Dhooge, Willy Govaerts, and Yu A Kuznetsov. Matcont: a matlab package for numerical bifurcationanalysis of odes.
ACM Transactions on Mathematical Software (TOMS) , 29(2):141–164, 2003.[28] CB Muratov and VV Osipov. Traveling spike autosolitons in the gray–scott model.
Physica D: NonlinearPhenomena , 155(1):112–131, 2001.[29] Wan Chen and Michael J Ward. Oscillatory instabilities and dynamics of multi-spike patterns for theone-dimensional gray-scott model.
European Journal of Applied Mathematics , 20(2):187–214, 2009. ∂ xx by ∂ xx + ∂ yy , on a square domain ( x, y ) ∈ ( − , . Here, ε = 0 . , D = 4 , κ = 1 and τ isas given in the caption. The spike is rotating clockwise. Higher values of τ lead to complex processionorbits. [30] Yasumasa Nishiura and Daishin Ueyama. Spatio-temporal chaos for the gray–scott model. Physica D:Nonlinear Phenomena , 150(3-4):137–162, 2001.[31] Yasumasa Nishiura. Dynamics of particle patterns in dissipative systems: Splitting · destruction · scattering. Sugaku expositions , 22(1):37–55, 2009.[32] Martin Golubitsky, Victor G LeBlanc, and Ian Melbourne. Meandering of the spiral tip: an alternativeapproach.
Journal of nonlinear science , 7(6):557–586, 1997.[33] Takashi Teramoto, Katsuya Suzuki, and Yasumasa Nishiura. Rotational motion of traveling spots in dissi-pative systems.
Physical Review E , 80(4):046208, 2009., 80(4):046208, 2009.