Stationary peaks in a multivariable reaction--diffusion system: Foliated snaking due to subcritical Turing instability
IIMA Journal of Applied Mathematics (2020) 1 − Stationary peaks in a multi-variable reaction–diffusion system:Foliated snaking due to subcritical Turing instability E DGAR K NOBLOCH AND A RIK Y OCHELIS , Department of Physics, University of California, Berkeley, California 94720, USA Department of Solar Energy and Environmental Physics, Swiss Institute for DrylandEnvironmental and Energy Research, Blaustein Institutes for Desert Research, Ben-GurionUniversity of the Negev, Sede Boqer Campus,Midreshet ben-gurion 8499000, Israel Department of Physics, Ben-Gurion University of the Negev,Beer Sheva 8410501, Israel [Received on August 7, 2020]
An activator-inhibitor model of side-branching used in the context of vascular and lung development isconsidered on the supposition that spatially localized concentrations of the activator trigger local branch-ing. The model consists of four coupled reaction-diffusion equations and its steady localized solutionstherefore obey an eight-dimensional spatial dynamical system in one dimension (1D). Stationary lo-calized structures within the model are found to be associated with a subcritical Turing instability andorganized within a distinct type of foliated snaking bifurcation structure. This behavior is in turn associ-ated with the presence of an exchange point in parameter space in which the complex leading eigenvaluesof the uniform concentration state are overtaken by a pair of real eigenvalues; this point plays the roleof a Belyakov-Devaney point in this system. The primary foliated snaking structure consists of periodicspike or peak trains with N identical equidistant peaks, N = , ,... , together with cross-links consistingof nonidentical, nonequidistant peaks. The structure is complicated by a multitude of multipulse states,some of which are also computed, and spans the parameter range from the primary Turing bifurcationall the way to the fold of the N = Keywords : localized states, homoclinic snaking, wavenumber selection, reaction–diffusion media
1. Introduction
Reaction–diffusion (RD) type models are widely studied as a laboratory for the study of pattern forma-tion in spatially extended systems that are driven far from equilibrium [18, 17], as exemplified by physic-ochemical autocatalytic systems [33, 22, 61, 34], biological and medical applications [48, 55, 36, 37]studies of vegetation cover in ecology [51] and even nonlinear optics [3]. Even though RD models can-not capture the complexity of many real-world applications, especially in the medical context, they sharemany generic similarities with them and so are often exploited to provide a qualitative understanding ofthe mechanistic basis of pattern formation and evolution in a wide variety of natural settings, includingperiodic pigmentation on animal skins, intracellular Ca + and actin waves, electrical impulses in neu-rons and in the heart, swarming phenomena in bacterial colonies or flocks of birds and fish, spatiallylocalized resonances in the inner ear, and the development of tissues and organs. In fact, this approachhas already been outlined in the seminal work of Turing on morphogenesis [70] and of Hodgkin andHuxley on action potentials in the giant squid axon [30]. Since then several prototypical RD models,such as the FitzHugh–Nagumo [24, 57], the Keller–Segel [38] and the Gierer–Meinhardt [27] equations, c (cid:13) Institute of Mathematics and its Applications 2020; all rights reserved. a r X i v : . [ n li n . PS ] A ug of 30 E. Knobloch and A. Yochelis have been subjected to detailed study in order to gain further insight into many of these systems.In the biomedical context, RD models are often referred to as activator–inhibitor (AI) systems, wherethe activator is locally expressed on a fast timescale but diffuses slowly while the inhibitor diffusesrapidly. These two kinetic ingredients are essential for the presence of the Turing instability [70] thatgives rise not only to periodic patterns [56, 14, 58, 41, 40, 25, 47, 42], such as striped or hexago-nal structures, but also to spatially localized structures that may range from stationary [85, 7, 64] orpropagating [21, 84] single peaks to groups of peaks. While a complete mathematical framework forstudying these structures has not yet been established, uncovering partial mechanisms behind biomed-ical pattern-formation phenomena remains of utmost importance for understanding functional aspectsand application design behind drug delivery systems or implants. On the other hand, mechanistic studiesof medical systems are a fertile source of new mathematical questions, especially in the context of modelsimplification and applicability as a useful description of highly detailed and multi-scale phenomena. Inthe present work our interest lies in spatially localized states that are symmetric under spatial reflections x → − x in one dimension (a symmetry broken by propagating excitable pulses) and their apparent rolein initiating side-branching in mammalian organ development [13, 66, 76, 53, 73, 29].1.1. Localized states and homoclinic snaking
Existence and stability of single steady-state localized states as well as groups of such states have beenintensively studied in the last two decades in the context of many application-driven model equationsincluding nonlinear optics [23, 60], Faraday-waves [65, 82, 12, 19, 1, 63], fluid convection [4, 44, 5],solidification [2, 69], electrically charged molecules [26], vegetation structures [87, 88, 89, 90, 86], hair-root formation in plants [8, 74] and vascular mesenchymal cells [85, 83]. The mathematical organizationof coexisting localized states in one space dimension (1D) is broadly attributed to the homoclinic snaking(HS) phenomenon [39] the commonest manifestation of which is the so–called snakes–and–laddersstructure of the snaking or pinning region [10]. Other types of organization, referred to as collapsedHS [9, 82, 83], defect mediated HS [46], slanted HS [19], and foliated HS [62, 59] also arise.These different cases are distinguished not only by the nonlinear properties of the localized solutionsin the parameter plane (e.g., types of homoclinic and hetereclinic connections) but also by their origin,i.e., the primary bifurcations through which the localized solutions form. The origin can be attributedto physical symmetry-breaking characteristics of the system, for example, due to a subcritical Turing(or finite wavenumber) instability from which the snakes–and–ladders HS emerges [10], while the or-ganization of the resulting localized structures may be slanted (instead of being vertical) in parameterspace due to the presence of a conserved quantity or other nonlocal effects [23, 19, 68, 63]. Althoughthe factors that affect the HS properties may be obscured in specific model equations, the use of spatialdynamics can reveal both its origin and organization [16, 39], an approach considerably advanced byPatrick Woods [77, 31].In the present work we study in detail one such model system from this point of view. The modelis of AI type and is motivated by increasing empirical evidence for the dominance of AI behavior asthe driving mechanism behind branching [32, 50] and specifically the development of the lungs [80]although the details of this mechanism are still unclear owing to the multi-scale nature of the processesinvolved that range from molecular to tissue levels [54, 73, 29].Beyond its potential applicability the model studied below also introduces new mathematical issues,arising from the fact that its spatial dynamics formulation is higher-dimensional than that of the pro-totypical models such as the Swift–Hohenberg and complex Ginzburg–Landau equation or applicationmotivated models (e.g., Lugiato–Lefever, Gray–Scott and Gierer–Meinhardt models) all of which lead tationary peaks in a multi-variable reaction–diffusion system
Localized states and side-branching
Meinhardt in 1976 [49] proposed an AI approach for studying branching phenomena, which has beenemployed in the context of pulmonary vascular and lung development [80]. The model is based onfour fields A , H , S and Y that represent, respectively, the concentrations of an activator (BMP), aninhibitor (MGP), and substrate (TGF- β /ALK1), as well as the concentration of an irreversible markerfor differentiated endothelial cells [80]: ∂ A ∂ t = c SA H − µ A + ρ A Y + D A ∇ A , (1a) ∂ H ∂ t = cSA − ν H + ρ H Y + D H ∇ H , (1b) ∂ S ∂ t = c − γ S − ε Y S + D S ∇ S , (1c) ∂ Y ∂ t = dA − eY + Y + fY + D Y ∇ Y . (1d)Although D Y = D Y as finite but much smaller than the smallest diffusion coefficient among all other fields, i.e, D Y (cid:28) D A . This increased spatial order is inconsequential: what matters for us is that the D Y = ρ H , asa control parameter while keeping all other parameters fixed and within the range of previous studies: c = . µ = . ρ A = . ν = . c = . γ = . ε = . d = . e = . f = D A = . D H = . D S = . D Y = − .Direct numerical simulations (DNS) of (1) in two space dimensions (2D) show that side-branchesdevelop via growing localized concentrations (top panel in Fig. 1) that are locked to the edge of apropagating front generated by the differentiation field (bottom panel in Fig. 1); the direction of prop-agation corresponds to the invasion of the light color region by the dark color region. Moreover, it hasbeen shown through analysis of (1) in 1D that the initiation of these side branches corresponds to theexistence of localized solutions that bifurcate subcritically via a Turing instability of an expanding dif-ferentiated state uniform in the variable x along the front [81]. In the following, we focus on the originand properties of the localized structures shown in Fig. 1 and so formulate the problem in one spacedimension, ignoring the motion of the front in the normal, y , direction. The latter will be the subject ofa future publication.1.3. Problem outline
We suppose that Eqs. (1) are posed on the real line with periodic boundary conditions on the domain x ∈ [ , L ] . The 1D domain and the choice of the boundary conditions mimic the planar interface thatconnects the P ∗ and P states prior to the formation of any peaks. Our aim is to reveal the mechanism of 30 E. Knobloch and A. YochelisF IG . 1. A representative snapshot of a solution obtained from a direct numerical simulations of (1) for ρ H = . · − withperiodic boundary conditions on x ∈ [ , ] and Neumann boundary conditions on y ∈ [ , ] . The top panel represents the A fieldwhile the bottom panel shows the Y field, the dark colors indicating higher values of the field. The contour lines in the bottompanel mark the locations of the localized A states shown in the top panel. responsible for the appearance of localized structures along this interface, prior to their expansion in thetransverse direction. This picture is predicated on the assumption that these structures form on a fastertimescale than that of the transverse growth, an assumption that will be tested in subsequent work onthe 2D problem.In the following we describe the nature of the stationary nonuniform solutions focusing on the prop-erties of spatially localized states when L =
60. Our results are supplemented with additional results forother values of L . The steady solutions of this system solve an eight-dimensional system in space. Be-cause of this higher dimension new behavior can arise that cannot occur in the four-dimensional systemshitherto studied. Thus our study of the system (1) represents a nontrivial extension of existing work onspatially localized states.Specifically, we show that the steady localized states of (1) are organized in a bifurcation structure,which we call foliated snaking [62, 59]. These states consist of trains of identical equispaced stationarypeaks we refer to as primary states, as well as secondary states consisting of unequal amplitude peakswith unequal separations. We use the numerical continuation software AUTO [20] to continue bothtypes of states as a function of the bifurcation parameter ρ H . We show that these states are generatedat small amplitude near the Turing bifurcation at ρ H = ρ T creating a subcritical Turing pattern andthen follow them through folds back to their termination near at finite amplitude, also near ρ T . Weexamine the role of spatial resonances in generating the secondary states and identify the existence ofan exchange point (EP), ρ H = ρ EP , which plays a role in the behavior of the system analogous to that ofthe Belyakov-Devaney (BD) point. We confirm our L =
60 results by performing similar computationsfor L =
20, 24 and 30, obtaining essentially identical results, and discuss the relation of finite domainnumerical continuation results to a theoretical analysis of the unfolding of a global bifurcation at BD inRef. [75]. We conjecture that the foliated snaking structure we identify here is both robust and universaland should therefore be present in other systems of partial differential equations supporting peak-like tationary peaks in a multi-variable reaction–diffusion system F IG . 2. (a) Bifurcation diagram showing the coexistence and linear stability (solid lines) of uniform solutions of (1) in terms theL norm (4). The red rectangle marks the region of interest; there are two distinct families of uniform solutions in this region, aprojection effect absent from panel (b). (b) Bifurcation diagram in terms of A (within the region of interest), showing the onset ofa Turing instability onset at ρ T (cid:39) . · − , with the inset depicting the dispersion relation, i.e., real part of the growth σ as afunction of the wavenumber k according to (2), where k T (cid:39) .
18 is the critical wavenumber. states, ranging from branching in physiology to formation of plant roots and to vegetation patches. Inparticular, we suggest that the structure we uncover arises generically in systems of two coupled spatiallyreversible reaction-diffusion equation on a line. The paper concludes with an outlook for future work.
2. Linear analysis of uniform solutions
We begin by solving (1) for steady spatially uniform states. In Fig. 2(a), we show that together with the‘trivial’ linearly stable solution P ≡ ( A , H , S , Y ) = ( , , c / γ , ) there are three or five additionalnontrivial solutions that coexist for different values of ρ H . Among these coexisting solutions, onlyone is linearly stable with respect to spatially uniform perturbations and we refer to this solution as P ∗ ≡ ( A ∗ , H ∗ , S ∗ , Y ∗ ) . Note that the front solutions that are shown in Fig. 1, connect the P and the P ∗ states. The P ∗ solution is thus the solution of interest in the present paper. of 30 E. Knobloch and A. YochelisF IG . 3. The six eigenvalues λ as a function of ρ H . Imaginary parts of λ , are indicated by the + symbol. The point ρ BD (cid:39) . · − corresponds to the Belyakov-Devaney (BD) point where λ , become real, while ρ EP (cid:39) . · − indicates theexchange point (EP) where the eigenvalue | λ | becomes smaller than | Re [ λ , ] | as ρ H increases. The eigenvalues ± λ ∼ ±
800 arenot shown. The ρ H axis here and in all subsequent figures is scaled by 10 , unless stated otherwise. The homogeneous states P ∗ are found on the bottom portion of a branch that resembles a backwards S , as shown in Fig. 2(b). The stability of the P ∗ states is described by linearized equations about P ∗ , AHSY − A ∗ H ∗ S ∗ Y ∗ ∝ e σ t + λ x . (2)Thus σ represents the growth rate in time of the perturbation while λ represents its growth rate in space.For our spatial analysis we take σ = λ , a consequence of the spatialreversibility of the system.In this case, because of the small value of D Y , there is always a pair of large real eigenvalues ± λ . Onthe lower branch to the right of the Turing bifurcation, ρ H > ρ T , there is a second pair of real eigenvalues, ± λ and a quartet of complex eigenvalues. At the Turing bifurcation ρ H = ρ T (cid:39) . · − , computedhere for an infinite system (for L =
60 the value is essentially identical), these complex eigenvaluescollide pairwise on the imaginary axis, forming a pair of imaginary eigenvalues λ , = ± ik T of doublemultiplicity (see Fig. 3), where k T (cid:39) .
18 is the critical wavenumber at the Turing onset (see inset inFig. 2(b)). For ρ H < ρ T these split into four purely imaginary eigenvalues, a situation that correspondsto temporal instability of P ∗ , i.e., in this region σ ( k ) = Re [ σ ( k )] > k T . With further decrease in ρ H the two imaginary eigenvalues closest to the origin collide at the origin,i.e., at λ =
0, indicating the presence of the left fold of the S-shaped branch and above this fold theeigenvalues are real (not shown here). On the middle part of the S-shaped branch there are three pairsof real eigenvalues and these become four pairs above the right fold.If we instead increase ρ H from ρ T we find that the imaginary parts of the eigenvalues ± λ , grad-ually decrease, and that these eigenvalues become real at the Belyakov-Devaney (BD) point ρ BD (cid:39) . · − implying the disapperance of an intrinsic wavelength from the behavior of the system.Such a picture is familiar from similar behavior found, for example, in the Lugiato-Lefever equationin nonlinear optics when the homogeneous state is modulationally unstable [59]. However, the steadystates of this equation represent a four-dimensional spatial system. In the following we shall see that the tationary peaks in a multi-variable reaction–diffusion system λ present in the system (1) plays a new and significant role in the observed behavior.Figure 3 indicates the essential new feature arising from the presence of λ : as ρ H increases | Re ( λ , ) | increases and at ρ EP | λ | falls below | Re ( λ , ) | , i.e., the complex eigenvalues ± λ , cease being leadingeigevalues and their role is taken by the real eigenvalues ± λ , a fact that has a profound influence onthe behavior of the system. This is because the leading eigenvalues generically determine the shape oflocalized structure near the background homogeneous state. In the following we refer to the transitionat ρ EP (cid:39) . · − as an exchange point (EP). Figure 3 shows that this point is associated with a spatial resonance .We mention that a Turing or modulational instability is not a prerequisite for the presence of a BDpoint, although in cases in which a Turing instability is absent the localized structures and the associatedfoliated snaking can be traced to other bifurcations, such as a fold [59] or a transcritical bifurcation [67].However, the BD point continues to play a similar role to that here whenever it is present. In higherdimensions the EP point is expected to do the same.
3. Nonlinear results
The Turing bifurcation at ρ H = ρ T (cid:39) . · − provides a key to much of the nonlinear behaviorobserved in this system since it gives rise to a subcritical branch of periodic states with wavenumber k T and corresponding wavelength (cid:96) T ≡ π / k T (cid:39) .
88. These periodic states, which extend towards ρ H > ρ T , are therefore temporally unstable, see [81] for details.Our interest is in the organization of the stationary localized states associated with this bifurcation.For this purpose, we employ numerical continuation [20] and compute solutions of the system 1 writtenas a set of eight first order ODEs in space: A (cid:48) = − a , (3a) a (cid:48) = D − (cid:18) c SA H − µ A + ρ A Y (cid:19) , (3b) H (cid:48) = − h , (3c) h (cid:48) = D − (cid:0) cSA − ν H + ρ H Y (cid:1) , (3d) S (cid:48) = − s , (3e) s (cid:48) = D − ( c − γ S − ε Y S ) , (3f) Y (cid:48) = − y , (3g) y (cid:48) = D − (cid:18) dA − eY + Y + fY (cid:19) . (3h)Here the prime indicates differentiation with respect to x . Our results are summarized in the next section.These resemble the phenomenon of foliated snaking originally described in [62] in the context of aspatially forced two-variable PDE (or four-variable spatial ODE), and subsequently identified in thehomogeneously forced Lugiato-Lefever equation in Ref. [59]. However, the present system is higher-dimensional and this fact is reflected in the observed behavior as discussed further below.3.1. Localized states on large domains: A new type of foliated snaking
In Fig. 4(a), we show the basic bifurcation diagram for a domain of length L =
60, showing the branchof homogeneous states near ρ T , as well as branches of N =
1, 2, 3 and 4 identical equispaced peaks in of 30
E. Knobloch and A. Yochelis the domain, as shown in the inset. The diagram shows the L norm of the solutions,L = (cid:115) L − (cid:90) L d x (cid:104) A + ( A (cid:48) ) + H + ( H (cid:48) ) + S + ( S (cid:48) ) + Y + ( Y (cid:48) ) (cid:105) , (4)as a function of the parameter ρ H . Observe that on sufficiently large domains ( L (cid:29) (cid:96) T ) the low N branches all bifurcate subcritically, and that all subsequently undergo a fold near essentially the samevalue ρ H = ρ SN (cid:39) . · − (see Fig. 5), before turning around towards smaller values of ρ H . In thefollowing we examine in detail the behavior of these primary branches near their point of origin nearthe Turing point at ρ T , their folds on the right, near ρ SN , and finally those on their left since all N -peakbranches appear to terminate back at the Turing point ρ T , although this time at finite amplitude.In fact Figs. 4(b,c) show that the primary branches do not terminate at ρ T (for ρ H < ρ T the leadingeigenvalues λ are imaginary thereby preventing approach to P ∗ ) but instead turn around in a narrowfold, and that in this region they start to develop small intervening peaks, i.e., in this region the solutionsno longer consist of equispaced identical peaks.To understand this behavior qualitatively recall that in the absence of locking between adjacentpeaks, i.e., beyond the so-called Belyakov-Devaney transition, ρ H > ρ BD , where all the spatial eigen-values λ are real. we expect all primary solutions to consist of equispaced peaks. In contrast when ρ T < ρ H < ρ BD the spatial eigenvalues become complex and such eigenvalues imply the presence ofoscillations in the tail of each peak and hence the presence of locking of adjacent peaks at distinct sep-arations. In other words, we expect a large variety of states for ρ T < ρ H < ρ BD owing to the presenceof spatial locking but a much simpler situation for ρ H > ρ BD where locking is absent. This is in fact thepicture that is consistent with the behavior observed in four-dimensional systems [59, 75].The present system is effectively six-dimensional, however, if we ignore the pair of large real eigen-values. Figure 3 shows the behavior of the six spatial eigenvalues along the branch of homogeneousstates as ρ H increases through ρ H = ρ T for our parameter values. We see that for ρ T < ρ H < ρ EP < ρ BD the complex eigenvalues have the smallest real part ( | λ | > | Re ( λ , ) | ) but that for ρ EP < ρ H < ρ BD thesituation is reversed ( | λ | < | Re ( λ , ) | ). Since generically the solution trajectory is expected to approachthe origin along the eigenvector corresponding the eigenvalue(s) closest to the origin, we see that for ρ EP < ρ H < ρ BD the trajectory approaches the origin along a direction associated with a real eigenvalue,implying absence of spatial locking in this parameter regime. Thus, in the present system the transi-tion from the presence of locking to its absence as ρ H increases occurs at ρ EP (cid:39) . · − insteadof ρ BD (cid:39) . · − and we expect, therefore, that states consisting of identical peaks with unequalspacing can only exist in the interval ρ T < ρ H < ρ EP . This result agrees with the continuation of theprimary branches of N -peak states with large interpeak separations (see Fig. 4(b,c)), such as N = , L =
60, for which the peak profiles develop oscillations once ρ H < ρ EP (cid:39) . · − leading, in general,to locking and hence non-equispaced peak states (if N >
1) as soon as ρ H < ρ EP . This is a significantnew feature inherent in the present problem.We focus first on the bifurcations creating the N =
1, 2, 3 and 4 branches near ρ H = ρ T , followedby the behavior near the right folds which are present in the region ρ H > ρ EP before returning to adiscussion of the behavior near the finite amplitude left folds which are also located in ρ T < ρ H < ρ EP .3.2. Origin of the primary states
The dispersion relation in Fig. 2(b) indicates that on the real line the homogeneous state loses stabilityas ρ H decreases at a Turing bifurcation that takes place at ρ H = ρ T . This bifurcation creates a subcritical tationary peaks in a multi-variable reaction–diffusion system T SN H L N=1N=2N=3N=4LLL LLL LLLLSSLSSL LLSLLSSSL LSLLSS LLSLSLS (b) T EP H L SSL
N=1
L LSS L (c) T EP H L N=2
LL LSLSLSSSLSSL F IG . 4. (a) Bifurcation diagram showing the primary N = , , , ρ T when L =
60 as well as additionalinterconnecting branches. Color-coded profiles at locations indicated by dots are provided alongside and represent periodic stateson the real line. (b,c) Detail of (a) showing that the finite amplitude termination points on the upper left of the (b) N = N = • (solid profiles) and (cid:7) (dashed profiles). TheLSLS branch in (a) is shown dashed since it has almost identical L norm to LLSS; branches of LLLS and SSSL are omitted (seeFig. 8(a)). E. Knobloch and A. Yochelis (a) T SN H L N=4N=3N=6N=1N=2N=8 (b) L S N N=1N=2N=3N=4N=6 F IG . 5. (a) Bifurcation diagram for L =
20 showing that the behavior of the system departs from its universal behavior (Fig. 4(a))when the number N of peaks within the domain becomes too large. (b) Location ρ SN of the right folds as a function of the domainsize L , obtained by continuation in L of the location of the right fold of the N = N -peak folds. The figure shows that ρ SN is close to its asymptotic value for L (cid:38) Turing pattern with wavelength (cid:96) T [81]. In a periodic domain of finite length L this bifurcation is slightlyshifted as is the Turing wavelength (cid:96) T in order that the domain accommodate an integer number n ofwavelengths ( n =
10 when L = L . These two states differ inwhether the maximum amplitude of the modulation coincides with a Turing peak or a Turing trough.Both of these branches generate states we refer to as N = N = N = L appear simultaneously via a secondary bifurcation from the Turing state, andso are strictly speaking secondary branches (purple curves in Fig. 6). However, since this secondarybifurcation occurs at very small amplitude (and is invisible in large scale bifurcation diagrams, such asFig. 4), we refer to these states in the following as primary states, as already done in Fig. 4.The single peak N = ρ H . In contrast, thetwo-peak N = tationary peaks in a multi-variable reaction–diffusion system
11 of 30 T H L TuringN=1 P * N=2 F IG . 6. Origin of the N -peak states for L =
30 showing the L norm of the subcritical Turing branch together with a pair ofsecondary branches bifurcating from it at small amplitude, all as functions of ρ H ; P ∗ denotes the homogeneous state. The firstsecondary states (purple curves) are the result of a modulational instability of the Turing state with period L and hence are referredto as N = N = N = ρ H until it reaches L / ρ EP where it connects to an N = L / N = N = N = ρ EP . The inset shows thecolor-coded solution profiles on each branch at ρ H = . · − (the solid profiles correspond to the lower branch of each pair,indicated by the • symbol, while the dashed profiles correspond to the upper branches, indicated by the (cid:7) symbol). gradually increases and approaches L / ρ H → ρ EP (Fig. 7). At this point the branch terminates on an N = L / N = N = N = ρ H and is also shown in Fig. 4(a).Note that these modulational instabilities correspond to spatial n : N resonances between the wavelength L / n Turing pattern ( n =
10) and the modulation wavelength L / N , viz. 10 : 1, 10 : 2 etc. We shall seein the next section that spatial resonances also play an important role at the right folds of the survivingprimary N -peak branches.We believe, but have not checked, that the other primary N -peak branches shown in Fig. 4(a) aregenerated via the same mechanism (modulation instability followed by spatial localization) as the N = , ρ H = ρ EP serves the “prune” the branchesemerging from successive modulational instabilities with wavelength L / N , N < n , leaving only theequispaced states beyond ρ EP .We mention that true N = , , . . . primary states are the result of primary bifurcations from thehomogeneous state that occur as ρ H decreases below ρ T , at parameter values where wavenumbers k = N π / L are themselves marginal ( σ = L =
20 the trueprimary N = ρ H (cid:39) . · − , a parameter value very close to the location of theleft fold of the homogeneous state, ρ T (cid:39) . · − . These states are completely distinct from the N -peak states created in the modulational mechanism just described and present in ρ H > ρ T . In particular,true primary states are present even when the Turing bifurcation is supercritical while the N -peak statesrequire the presence of a small amplitude modulation instability that in turn requires this bifurcationto be subcritical. Note that as the domain size increases all these secondary modulational instabilities2 of 30 E. Knobloch and A. Yochelis (a) T EP H L (b) T EP H L F IG . 7. Color-coded profiles along the N = N = ρ H → ρ EP and their separation approaches L / N . In the N = L / L =
30. The dashed line is the Turing branch. collapse on ρ T . Thus, when the Turing bifurcation is subcritical and the domain infinite, the Turingbifurcation is responsible for a large multiplicity of distinct states.3.3. Behavior near the right folds
Each fold on the right (see Fig. 4(a)) represents a fold of a branch of identical equispaced peaks createdas described in the preceding section. Near such folds there is always a Floquet multiplier that isclose to +
1, i.e., close to neutral stability with respect to amplitude perturbations. This fact allows,under appropriate circumstances, nearby bifurcations to secondary branches consisting of peaks with nonidentical heights. The number of new branches depends on the number N characterizing the primarybranch. Thus when N = N = N = N = x = ( n / N ) L , 0 (cid:54) n < N , modulo overall translation, even tationary peaks in a multi-variable reaction–diffusion system
13 of 30though they are no longer identical. However, when the number N of peaks within the domain becomestoo large the behavior of the system starts to depart from the universal behavior exhibited in Fig. 4(a),see Fig. 5.Figure 4(a) shows the details of this universal behavior. Each of the new states can be consideredto be the result of spatial modulation of a periodic train of identical peaks associated with the onsetof Eckhaus instability in the vicinity of the fold [6]. The bifurcation that takes place closest to thefold corresponds to modulation with the longest possible wavelength, and so leads to states with onemodulation wavelength in the domain; subsequent modulational instabilities generate states exhibitingtwo modulation wavelengths in the domain such as the state SLSL but these occur farther away fromthe fold. This is essentially the same as the mechanism responsible for the presence of the N -peak statesin the first place, except that here the relevant spatial resonances are of lower order since N is muchsmaller than n , the number of Turing wavelengths that fit in the domain. However, it is still true that thelong wave states are generated in N : 1 spatial resonances that take place nearest to the fold, while theother states mentioned are the consequence of additional spatial resonances that are present farther fromit [6].We mention that the local behavior near a N : 1 spatial resonance is described by amplitude equationsof the form ˙ B N = µ B N + c ¯ B N − N + c | B N | B N , (5)where B N is a complex modulation amplitude such that | B N | measures the height difference between Land S, µ represents the distance from the resonance and c , c are a pair of constants; ¯ B N is the complexconjugate of B N . Writing B N = r N exp i φ N we see that the new branches generated in the N : 1 resonancecorrespond to solutions of sin N φ N =
0, or φ N = π m / N , m = , , . . . , N −
1. Not all of these will havedistinct L norms, however. Note also that the strong resonances N = N = N = m =
1, while the N = m = m =
2, as seen inFig. 4(a). Note that this is (locally) a transcritical bifurcation only because the bifurcation does not infact coincide with the fold, although it is very close to it.The case N = µ =
0. The above theory predicts two distinct states with modulation wavelength L ,but there are in addition branches exhibiting modulation on smaller scales that also bifurcate very closeto the fold although slightly further from it. Given that there are four peaks in the domain, the possiblemodulated states are LLSS, LSLS, LLLS and LSSS. All of these have been found and are included inFig. 8(a), cf. [45]. Note that since the L norm of the LLSS and LSLS states is essentially the same thediagram appears to show 3 distinct branches bifurcating from the N = N − L = ρ H to the left of this fold the peak profiles present inthe observed unequal peak states do indeed accurately match the peak profiles on the upper and lowerbranches of the N -peak states present at this value of ρ H . Observe that in the LLSS state the SS peaksare significantly closer to one another than the LL peaks. We can understand this property in terms of aforce between adjacent peaks when ρ H > ρ EP . In this regime the preferred state is the N = N = repel one another [85, 74],and so inserting further peaks leads to N = , , ... primary states, all consisting of identical equispacedpeaks. The repulsive force depends on the amplitude of the peaks, however, and is stronger between4 of 30 E. Knobloch and A. Yochelis (a) T SN H L N=4N=1L LSLLSS LLLLLS N=3N=2LL LSLS LLLLLLLSLSSLSSS (b)(c) F IG . 8. (a) Details of the bifurcation diagram for L =
24 near the folds on the right confirming the universality of the behaviorshown in Fig. 4. (b) Solution profiles along the N = ρ H = . · − . Blue: LSLS. Purple: LLSS. Red: coexistingLLLL and SSSS states on the N = ρ H = . · − . L: solid line; SSL: dashed line. The profiles are shown on alogarithmic scale to evince the pair of small peaks S present in the latter case. taller peaks and weaker between shorter peaks. This fact accounts for the different separations withinthe LLSS state and the LSLS state; these separations depend on ρ H which in turn determines the height tationary peaks in a multi-variable reaction–diffusion system
15 of 30of the L and S peaks and hence the force between adjacent LL, LS and SS peaks. This is so for LLLS andLSSS as well (not shown here, but see Fig. 10(c) for L = N -peak states of the form SL ... L and LS ...
S to be present forany N . When N and the domain length L are both large, these states approximate homoclinic connec-tions from L ... L to itself and from S ...
S to itself, respectively, while states such as S ... SL ... L representheteroclinic cycles connecting S ...
S to L ...
L and back again, that is, connections between the small andlarge N -peak states that coexist at each ρ H to the left of the right fold of the N -peak state. Thus weexpect periodic trains of peaks to exhibit the same localization and subsequent snaking as that arisingfrom the subcritical Turing bifurcation itself.Note that from the point of view of the primary N = N = ρ EP for the small amplitude terminating branches) or due to adifference in the heights of adjacent peaks, as occurs at the right folds ρ SN ). A similar distinctionbetween phase and amplitude-generated period multiplication applies for other N : 1 resonances.3.4. Behavior near the upper left folds
As already mentioned, the apparent termination points of all multi-peak branches on the left (Fig. 4(b,c)for the case L =
60) are in fact sharp, cusp-like folds (see Fig. 9(a,b)). Similar behavior is present for L =
24 (Fig. 8(a)) as well as other domain sizes. We believe that these folds accumulate on ρ H = ρ T ,with N = N = ρ H < ρ T . Near these folds each primarybranch transitions into a state of nonidentical peaks. At the right folds this occurs when the periodicstate starts to develop differences in the peak heights, a transition that sets in via bifurcations generatedby spatial resonances. On the upper left the transition is instead continuous and occurs when new peaksstart to grow in the interpeak region between existing peaks. Figure 8(c) shows the transition from a1-peak state L to a 3-peak state SSL as one passes the left fold when L =
24: two identical new smallpeaks appear in the domain below the fold and grow in height as one passes the fold and beyond. Owingto the periodicity of the domain these small peaks may be thought of a growing sidebands of the mainpeak, and hence as an SLS state [75].Near the (upper) left fold the peaks in the resulting SSL (SLS) state are not equispaced, but asone follows the SSL branch through ρ EP and beyond, the amplitude of the small peaks continues togrow and their location changes in response. By the end of the branch the three peaks are identicaland equispaced, and the branch terminates in a 3:1 resonance near the right fold on the 3-peak primarybranch that emanates from the vicinity of the Turing bifurcation. The figure also shows a third branch,labeled SL. This branch bifurcates from the continuous branch near its fold and is characterized by theappearance of only a single new peak (Fig. 9(a)). Like the SSL branch, the SL branch can be continuedall the way to its termination near the right fold of the primary N = N = L the numerical continuation of the SSLbranch (and related branches) through the exchange point ρ EP as ρ H increases becomes exceedinglydifficult and one is tempted to conclude that the branch terminates at ρ EP . While such a termination is6 of 30 E. Knobloch and A. Yochelis (a)(b) F IG . 9. Detail of the L =
60 bifurcation diagram near the upper left for (a) N = N =
2, showing that the cusp-like featuresin Fig. 4 are in fact narrow folds with secondary bifurcations generating LS and LSLS branches. These folds are very close to ρ T (cid:39) . · − , the Turing bifucation that creates the N -peak states at the lower left. indeed possible on the real line where the peaks can move arbitrarily far apart as one approaches ρ EP this is not the case on finite domains. Instead, we find that while the peaks do indeed start to moveapart as ρ H increases the solution branch passes smoothly through ρ EP and in fact extends all the wayto the right fold of the N = , , . . . -peak states where all N peaks are of the same height and equidistant(Fig. 10). In fact numerical continuation decreasing ρ H through ρ EP is possible even when continuationin the opposite direction becomes impossible, an important fact that appears to have been missed inearlier work [74].We emphasize that the interval ρ T < ρ H < ρ EP is special in that it permits spatial locking owingto the presence of oscillations in the tails of the peaks while for ρ H > ρ EP the peak profiles becomemonotonic and spatial locking is absent. Our numerical continuation suggests that the new peaks thatappear along the N = , , . . . branches near the upper left folds arise from the growth of select peaks inthese tails, explaining why the N = , , . . . states start adding new peaks only once ρ H < ρ EP . Figure8(c) shows an example showing the creation of the SSL branch near the upper left fold when L = ρ EP the original and new peaks unlock, but do not become equispaced. Asalready explained, this is a consequence of unequal forces between small and large peaks. Once all thepeaks have acquired the same height the peaks become equispaced and connect to the corresponding tationary peaks in a multi-variable reaction–diffusion system
17 of 30(a) T SN H L SSLLSLSLSSS LLS (b) T SN H L SSLLSLSLSSS LLS (c) T SN H L LSSLLLSLSSS LS S F IG . 10. Evolution of the peak heights and spacing along (a) the SL branch, (b) the SSL branch and (c) the SSSL branch created inthe upper left N = L =
60. These branches terminate on the primary 2-peak, 3-peak and 4-peak branches, respectively.The peaks separate rapidly near their termination in order to reach the respective separations L / L / L / primary N -peak branch (Fig. 10), provided, of course, the domain length is sufficiently large ( L (cid:29) (cid:96) T ).8 of 30 E. Knobloch and A. Yochelis
Foliated snaking
The behavior described in the two preceding subsections resembles a phenomenon called foliated snakingthat was originally described in [62] in a spatially forced spatially two-dimensional system, and subse-quently identified in the homogeneously forced Lugiato-Lefever equation in Ref. [59]. However, thepresent system differs in two important ways. In our system the foliated snaking structure is createdfrom successive modulational instabilities of a subcritical Turing state. In contrast, the N -peak statespresent in the Lugiato-Lefever equation are created at a fold of the uniform states [59] while similarstates present in both vegetation [67] and Gierer-Meinhardt [85] models are the result of a transcriticalbifurcation of this state. Thus our system is closer to the systems studied in [43, 74, 75] than the Lugiato-Lefever equation. Moreover, it is higher-dimensional than these examples and this fact is reflected inthe behavior we observe.Foliated snaking is found in the present problem because the localized states created in the modula-tional instability of the Turing state extend into the region of no spatial locking where the leading spatialeigenvalues λ are real – in contrast to the Swift-Hohenberg equation in which the localized states remainconfined to the locking region where the eigenvalues λ form a quartet in the complex plane. We haveseen that outside the locking region the localized states take the form of N -peak states, i.e. a periodictrains of equispaced stationary peaks characterized by the number N of peaks in the domain.In general, we expect an absence of locking between adjacent peaks when ρ H > ρ BD , i.e., beyondthe so-called Belyakov-Devaney transition, where the four leading spatial eigenvalues λ are real. As aresult we expect only equispaced peaks when ρ H > ρ BD and non-equispaced peaks for ρ H < ρ BD . Inother words we expect a large variety of states when ρ H < ρ BD owing to the presence of spatial lockingbut a much simpler situation for ρ H > ρ BD where locking is absent. This is the picture that is consistentwith the behavior observed in 4-dimensional systems [59, 75] where the role played by the BD point isreliably established.The present system is 6-dimensional, however, if we ignore the pair of large real eigenvalues. Figure3 shows that for ρ T < ρ H < ρ EP < ρ BD the complex eigenvalues have the smallest real part ( | λ | > | Re ( λ , ) | ) but that for ρ EP < ρ H < ρ BD the situation is reversed ( | λ | < | Re ( λ , ) | ). Since generically thesolution trajectory approaches the origin along the eigenvector corresponding the eigenvalue(s) closest to the origin, we see that for ρ EP < ρ H < ρ BD the trajetory approaches the origin along a directionassociated with a real eigenvalue, implying absence of spatial locking in this parameter regime. Thus inthe present system the transition from the presence of locking to its absence as ρ H increases occurs at ρ EP (cid:39) . · − instead of ρ BD (cid:39) . · − and we therefore expect states with unequal peak spacingto exist in the interval ρ T < ρ H < ρ EP only. Our results are fully consistent with these conclusions.We now address a key issue, the relation between existing theory describing the unfolding of a globalbifurcation at BD point [75] and computations that are necessarily performed on finite (periodic) do-mains. On the real line identical peaks are expected to separate indefinitely as one approaches ρ EP fromabove, corresponding to the formation of an infinite number of homoclinic orbit(s) to the homogeneousstate [16]. Thus all primary N -peak ( N >
1) states would be expected to terminate at ρ EP . In a finiteperiodic domain, however, a global bifurcation cannot take place and all states remain periodic. Thepeak separation on the N -peak primary branches remains L / N until ρ EP as ρ H decreases from the rightfolds, but new intermediate peaks start to develop once ρ H < ρ EP and the resulting branches turn aroundat the upper left folds near ρ T beyond which the purely imaginary eigenvalues block approach to thehomogeneous state. Thus the peak states exist in the interval ρ T < ρ H (cid:54) ρ SN instead of the expectedinterval ρ EP < ρ H (cid:54) ρ SN (the single peak state N = ρ T ). Thusin finite periodic domains the BD point is no longer associated with global bifurcations but with period tationary peaks in a multi-variable reaction–diffusion system
19 of 30multiplication where the wavelength L / N jumps to the period L of the domain. Other states, composedof dissimilar peaks, such as those present after the upper left fold, or created in bifurcations from it, passthrough ρ EP without noticeable change since the peak separation is a function of the peak amplitudeswhich varies continuously along such branches until their termination in the vicinity of the right foldsof the primary N -peak states. These results are computationally difficult to establish unless one followsthese branches from their creation at the right folds down in ρ H . We conjecture that similar behavior ispresent in four-dimensional systems but may have been missed for this reason.We mention, finally, that branches such as the SL branch appear to be absent from the existingunfolding of the global bifurcation at the BD point [75]. We have seen that such states in fact play afundamental role in the overall structure of the foliated snaking bifurcation diagram identified here.3.6. Composite-peak isolas
Although the structure we have uncovered so far is already quite complex, it omits a large variety ofinteresting additional states. From studies of the quadratic-cubic Swift-Hohenberg equation it is knownthat the snaking or pinning region is populated by a variety of so-called multipulse states [11]. Thesimplest of these are equispaced 2-pulse states, consisting of two copies of a 1 , , , . . . -peak state,separated by L /
2. States of this type lie on a connected branch that snakes much like the 1-pulse states,except that the number of back-and-forth oscillations is half that of the 1-pulse branch. However, thetwo pulses need not be separated by L /
2: they can form bound groups, in which the pulses are separatedby an integer number p of half-wavelengths (cid:96) . For each p a 2-pulse state consisting of two identicalpulses, each with q peaks, lies on a closed branch called an isola, one for each pair of integers ( q , p ) .The smallest isola corresponds to p = p increases, accumulating onthe snakes-and-ladders structure of the 1-pulse bifurcation diagram in the pinning region, with higher q isolas corresponding to higher L norm. Thus the snaking or pinning region contains a stack ofnested isolas, one for each choice of ( p , q ) . This is not all, however, since bound states of two differentpulses are also possible. These also snake (see, e.g., Fig. 6, of [11] where the snaking of a bound stateconsisting of one 1-peak state and one 2-peak state is investigated). Similar 3-pulse states etc are alsopresent, altogether generating a structure of remarkable complexity.With this in mind we recall that for ρ T < ρ H < ρ EP the homogeneous state P ∗ of the present systemhas leading eigenvalues λ that form a quartet in the complex plane, and in this region we expect thepresent system to behave much like the Swift-Hohenberg equation. In particular, we expect that inaddition to the N -pulse equispaced states discussed thus far, the system also exhibits multipulse statesconsisting of groups of pulses with different intergroup separations, and that these are created in theinterval ρ T < ρ H < ρ EP in the same way as the corresponding states of the Swift-Hohenberg equation.Of course each of these states will exit the interval ρ T < ρ H < ρ EP and so participate in the foliatedsnaking structure of the primary 1-pulse states. In the following we present the results of detailedcomputations of one such sequence of states, fully cognizant of the limitations of numerical studies ofthis kind in revealing the full complexity of the system. We think of the branches of multipulse stateswe find as providing decoration of the basic foliated structure already described.Figure 11 provides a summary of the states computed on L =
24 and consisting of groups of 2-peakand 3-peak states. We first mention the branches labeled by N = , , . . . . In contrast to the primarybranches of equispaced 1-peak states forming the skeleton of the foliated snaking structure (not shownin this figure), these branches consist of N equispaced groups of two pulses. Profiles of N = N = x ∈ [ , L ] are shown in the inset using solid and dashed lines, respectively.The basic structure of the figure recapitulates that of Fig. 8(a), albeit with more branches. We now0 of 30 E. Knobloch and A. Yochelis T H L N=2N=4N=1L d F IG . 11. Bifurcation diagram for a subset of the multipulse states present in the system when L =
24. Here the integers N = , ,... refer to branches of N -pulse states consisting of N N = , discuss the behavior shown in this diagram in greater detail, focusing again on the folds.Figure 12 shows details of the N = N = d and L d , depending on whether the profiles lies below or abovethe N = d turn into L d at the fold. The N = d close to ρ T and does so in a fold (Fig. 12(b)). At the fold the branch turns smoothly into an S d SS state with theappearance of a pair of small single peaks; a third branch, of S d S states bifurcates from the fold, justas in the case of the 1-peak states. The figure also shows a second fold, where a branch labeled SS d Sturns into S d SSS through the addition of a single small peak. The S d SS and SS d S states are no longerthe same, as can be seen from the profiles on each branch, computed at the same value of ρ H , shown inthe inset.Each of these five branches can be traced through its fold on the right (Fig. 12(a)) and towardsthe associated upper left folds, a detail of which is shown in Fig. 12(c) where analogous behaviortakes place: L d turns smoothly into SL d S, while an L d S branch bifurcates directly from the fold. Twoadditional folds are also present, involving the transformation from L d S to L d SS and the transformationfrom a different L d SS into L d SSS. The boxed labels indicate branches that extend to large values of ρ H and so are created in the folds on the right. One of these is the L d S d branch that terminates near the N = ρ H the height of S d increases andthe separation of the two pulses adjusts in such a way that the state becomes an equispaced L d L d stateby the time the branch terminates.Near the upper left folds of Fig. 12(c) one also finds another set of folds, as shown in Fig. 13(a).These branches correspond to bound states of a 2-peak state L d and a 3-peak state S t below the right foldand L t above it. As shown in Fig. 13(b) these states are born in a manner that recapitulates the originof the N = N = ρ H , pass through folds (Fig. 13(a)) and then return to ρ H (cid:39) ρ T with S t having growninto L t (Fig. 13(c)). The main observation that merits a comment is that the right folds of these new2-pulse states do not align with those of the original 2-pulse state associated with N =
1. We believe thisto be a finite size effect: if one removes the L d portion of the state, the 3-peak component S t behavesas an isolated 3-peak state on a much smaller domain, something like half the domain length L = tationary peaks in a multi-variable reaction–diffusion system
21 of 30(a) T H L N=1 L d S d L d S d L d S d SSSS d SSS d L d SSS d SL d SL d SSSL d S d SSS (b) H L S d SSSS d SS d S S d SSS d S S d SS (c) H L L d S L d SSL d SL d SSL d L d S L d SSL d SS L d SSS F IG . 12. (a) Detail of the lower N = d SS and SS d S corresponding to the profiles shown in the inset. (c)Further detail of the region near the upper left folds. There are four folds and three distinct branches L d SS. The boxed labelsindicate branches that extend to large values of ρ H . (Fig. 13(c)). Since the wavelength of the oscillations and hence the inter-peak separation is somethinglike 2.8, the width of S t is of order 8.5 and hence quite comparable with the halflength of the domain.2 of 30 E. Knobloch and A. Yochelis (a)(b) H L L d S t SSL d S t L d S t S (c) H L L d L t S L d L t L d L t SS F IG . 13. (a) Detail of the upper N = d SS and SS d S corresponding to the profiles shown in the inset. (c)Further detail of the region near the upper left folds with two distinct branches L d L t S (see inset).
On such small effective domains the snaking of S t will depart from its asymptotic behavior on largedomains, pushing the right folds inwards.Note that the 2-pulse states based on L d and S t /L t form a completely distinct set of isolas from those tationary peaks in a multi-variable reaction–diffusion system
23 of 30based on L d and S d /L d . The latter do connect to the equispaced 2-pulse 2-peak state near the fold of the N = N = d SSS (Fig. 12(a,c)).It is clear that many more multipulse states of this type are present in the parameter range ρ T < ρ H < ρ SN but the above results suffice to illustrate the type of behavior expected. We mention that bothmultipulse and composite-peak states were computed in a four-dimensional non-Hamiltonian system byChampneys [15], lending support to the conjecture at the end of the previous section.
4. Discussion and conclusions
Isolated, spatially localized peaks or spikes (spots in two dimensions) play an increasingly importantrole in many systems of physical, chemical, and biological interest. We have studied here in some detailone such system that exhibits a large multiplicity of different multi-peak states and showed that thesestates are organized in a structure we call foliated snaking. We believe this structure to be universaland thus expect similar structures to be present in a variety of different systems. We have seen that inthe present system this structure is a associated with the presence of a subcritical Turing instability, andarises whenver the localized states familiar from other models exhibiting a subcritical pattern-forminginstability such as the Swift-Hohenberg equation extend into a parameter region with real leading eigen-values. In this case the localized structures cannot lock via oscillatory tails and instead repel one another,forming a sea of isolated structures within the domain whose dynamical properties in time remain to bestudied.A subcritical Turing instability is not the only one resulting in isolated localized structures. Earlierwork [85, 59, 67, 90] has shown that similar structures can be generated as a result of the presence ofa fold or transcritical bifurcation of a homogeneous state, provided only that a Belyakov-Devaney (BD)point is present. The important feature of all these mechanisms is that they generate such states over amuch larger range of parameter values than the better known homoclinic snaking mechanism. The lattergenerates related states but only within a typically narrow parameter range – the snaking or pinninginterval.In the present work we studied in detail the origin of the primary N -peak states – these form near ρ T as a consequence of a modulational instability of subcritical Turing states – and their subsequentevolution with increasing ρ H via folds on the right at ρ SN through to their transformation into trains ofunequal peaks near ρ T again. These states form the cross-links connecting the primary N -peak statesbetween the left and right folds. Spatial resonances were found to play an essential role both in the originof the primary N -peak states at small amplitude and in the generation of trains of unequal peaks near ρ SN and hence in the bifurcation structure associated with foliated snaking. Finally the exchange point(EP) point was found to play the role usually associated with the BD point, and this role was found to becrucial for the transformation from the primary N -peak states to secondary states consisting of unequalpeaks that form the cross-links in the bifurcation structure. We also made an extensive exploration ofthe so-called multipulse or subsidiary states that decorate the primary foliated snakind structures. Theseresults only scratch the surface of the complexity exhibited by the present system.We have employed several different domain sizes L to obtain these results. Normally the domain size L is largely immaterial when one studies well-localized states. Here, however, the peaks in a multipeakstate for ρ H > ρ EP want to be apart as far as possible and hence the domain size sets the length scalefor structures of this type, in contrast to the intrinsic wavelength (cid:96) = π / | Im ( λ , ) | that sets the lengthscale in the locking regime ρ T < ρ H < ρ EP . Despite this inevitable dependence of our results on L most aspects are essentially independent of it. This includes the locations ρ T , ρ EP and ρ BD which are4 of 30 E. Knobloch and A. YochelisF IG . 14. Representative snapshot of a solution obtained from a direct numerical simulations of (1) as in Fig. 1 but at a later time. determined to high accuracy by the infinite domain calculation. In the nonlinear regime, the right foldsaccumulate on ρ SN as N varies, a quantity that also almost independent of L when L is sufficiently large(Fig. 5(b)). Since the transitions in the foliated snaking structure take place near ρ T , ρ EP and ρ SN weconclude that this structure is almost independent of the domain size L after all.Temporal stability (determined via a standard eigenvalue formulation of the linearization of thesystem (1)) shows that all stationary peaks and groups of peaks are linearly unstable in 1D while in2D isolated spots in a background of P ∗ and away from any front such as that in Fig. 1 may gainlinear stability, facts that explain the oscillatory or decaying peak dynamics in 1D and persistence ofpeaks in 2D [81]. We are aware that in order to establish the further significance of our results it isalso necessary to investigate the temporal stability properties of the various states we have found, andthis, in addition to the results of elaborate time-integration, will form the topic of a future contribution.Nevertheless, nucleation of side-branches can already be regarded as a nontrivial wavenumber selectionproblem because of the multiplicity of solutions that may give rise to distinct separation distances atwhich the peaks form, as shown in Fig. 1. In practice, side-branches in lungs typically appear in themiddle of the domain (far from gradients imposed by the leading peak), in between the initial branchedpoint and the leading peak of the growing branch, and so the minimal length scale of the nucleationprocess depends also on the location of the right folds, i.e., for larger values of ρ H side-branches appearlater and at large distances from the growing tip [80, 72, 81]. Moreover, the observed self-repulsion ofthe peaks explains the so-called avoidance phenomenon, i.e., the absence of peak-peak collisions and themerging of counterpropagating peaks in 2D, as well as the potential inhibition of new peaks on nearbybranches [29]. Figure 14 shows that once the peaks feel a counterpropagating peak (perhaps reflectedfrom the boundary due to Neumann boundary conditions in x direction) their trajectory is deflected inthe y direction, something that cannot happen in 1D.Evidently, the conjectured universal organization of peaks in the foliated snaking structure generatedby a subcritical Turing instability provides a distinct mathematical description of homoclinic snaking EFERENCES
25 of 30in multi-variable reaction–diffusion media. The results suggest that the nonlinear mechanisms behindboth peak generation and wavenumber selection that appear to be involved in the side-branching modelstudied here lead to a robust phenomenon that is essential, among other examples, to studies of ecologi-cal systems, ranging from the initiation of root hairs [8] to the generation of vegetation patches [52] andthe design of agroforestry systems [71].
Acknowledgment : We have benefited from discussions with Nicol´as Verschueren. This work wassupported in part the National Science Foundation under grant DMS-1908891 (E.K.).
References [1] A. S. Alnahdi, J. Niesen, and A. M. Rucklidge. Localized patterns in periodically forced systems.SIAM Journal on Applied Dynamical Systems, 13(3):1311–1327, 2014.[2] A. J. Archer, M. C. Walters, U. Thiele, and E. Knobloch. Solidification in soft-core fluids: Disor-dered solids from fast solidification fronts. Physical Review E, 90(4):042404, 2014.[3] F. T. Arecchi, S. Boccaletti, and P. Ramazza. Pattern formation and competition in nonlinear optics.Physics Reports, 318(1-2):1–83, 1999.[4] O. Batiste, E. Knobloch, A. Alonso, and I. Mercader. Spatially localized binary-fluid convection.Journal of Fluid Mechanics, 560:149–158, 2006.[5] C. Beaume, A. Bergeon, and E. Knobloch. Convectons and secondary snaking in three-dimensional natural doubly diffusive convection. Physics of Fluids, 25(2):024105, 2013.[6] A. Bergeon, J. Burke, E. Knobloch, and I. Mercader. Eckhaus instability and homoclinic snaking.Physical Review E, 78(4):046201, 2008.[7] V. Bre˜na-Medina and A. Champneys. Subcritical Turing bifurcation and the morphogenesis oflocalized patterns. Physical Review E, 90(3):032923, 2014.[8] V. Bre˜na-Medina, A. R. Champneys, C. Grierson, and M. J. Ward. Mathematical modeling ofplant root hair initiation: Dynamics of localized patches. SIAM Journal on Applied DynamicalSystems, 13(1):210–248, 2014.[9] J. Burke and E. Knobloch. Localized states in the generalized Swift–Hohenberg equation. PhysicalReview E, 73(5):056211, 2006.[10] J. Burke and E. Knobloch. Snakes and ladders: Localized states in the Swift–Hohenberg equation.Physics Letters A, 360(6):681–688, 2007.[11] J. Burke and E. Knobloch. Multipulse states in the Swift–Hohenberg equation. Discrete andContinuous Dynamical Systems, Supplement 2009, 225:109–117, 2009.[12] J. Burke, A. Yochelis, and E. Knobloch. Classification of spatially localized oscillations in period-ically forced dissipative systems. SIAM Journal on Applied Dynamical Systems, 7(3):651–711,2008.[13] J. H. Caduff, L. C. Fischer, and P. H. Burri. Scanning electron microscope study of the developingmicrovasculature in the postnatal rat lung. The Anatomical Record, 216(2):154–164, 1986.6 of 30
REFERENCES [14] V. Castets, E. Dulos, J. Boissonade, and P. De Kepper. Experimental evidence of a sustainedstanding Turing-type nonequilibrium chemical pattern. Physical Review Letters, 64(24):2953–2956, 1990.[15] A. R. Champneys. Subsidiary homoclinic orbits to a saddle-focus for reversible systems.International Journal of Bifurcation and Chaos, 4(6):1447–1482, 1994.[16] A. R. Champneys. Homoclinic orbits in reversible systems and their applications in mechanics,fluids and optics. Physica D, 112(1-2):158–186, 1998.[17] M. Cross and H. Greenside. Pattern Formation and Dynamics in Nonequilibrium Systems. Cam-bridge University Press, 2009.[18] M. C. Cross and P. C. Hohenberg. Pattern formation outside of equilibrium. Reviews of ModernPhysics, 65:851–1112, 1993.[19] J. H. P. Dawes. Localized pattern formation with a large-scale mode: slanted snaking. SIAMJournal on Applied Dynamical Systems, 7(1):186–206, 2008.[20] E. J. Doedel, A. R. Champneys, T. Fairgrieve, Y. Kuznetsov, B. Oldeman, R. Paffenroth, B. Sand-stede, X. Wang, and C. Zhang. Auto07p: Continuation and bifurcation software for ordinarydifferential equations. Concordia University, http://indy.cs.concordia.ca/auto, 2012.[21] C. Elphick, E. Meron, and E. A. Spiegel. Patterns of propagating pulses. SIAM Journal on AppliedMathematics, 50(2):490–503, 1990.[22] I. R. Epstein and J. A. Pojman. An Introduction to Nonlinear Chemical Dynamics: Oscillations,Waves, Patterns, and Chaos. Oxford University Press, 1998.[23] W. J. Firth, L. Columbo, and A. J. Scroggie. Proposed resolution of theory-experiment discrepancyin homoclinic snaking. Physical Review Letters, 99(10):104503, 2007.[24] R. FitzHugh. Impulses and physiological states in theoretical models of nerve membrane.Biophysical Journal, 1(6):445–466, 1961.[25] A. Garfinkel, Y. Tintut, D. Petrasek, K. Bostr¨om, and L. L. Demer. Pattern formation by vascularmesenchymal cells. Proceedings of the National Academy of Sciences, 101(25):9247–9250, 2004.[26] N. Gavish, I. Versano, and A. Yochelis. Spatially localized self-assembly driven by electricallycharged phase separation. SIAM Journal on Applied Dynamical Systems, 16(4):1946–1968, 2017.[27] A. Gierer and H. Meinhardt. A theory of biological pattern formation. Kybernetik, 12:30–39,1972.[28] Y. A. Guo, M. Sun, A. Garfinkel, and X. Zhao. Mechanisms of side branching and tip splitting ina model of branching morphogenesis. PloS One, 9(7), 2014.[29] E. Hannezo and B. D. Simons. Multiscale dynamics of branching morphogenesis. Current Opinionin Cell Biology, 60:99–105, 2019.[30] A. L. Hodgkin and A. F. Huxley. A quantitative description of membrane current and its applicationto conduction and excitation in nerve. The Journal of Physiology, 117(4):500–544, 1952.
EFERENCES
27 of 30[31] G. W. Hunt, M. A. Peletier, A. R Champneys, P. D. Woods, M. A. Wadee, C. J. Budd, and G. J.Lord. Cellular buckling in long structures. Nonlinear Dynamics, 21(1):3–29, 2000.[32] D. Iber and D. Menshykau. The control of branching morphogenesis. Open Biology, 3(9):130088,2013.[33] R. Imbihl and G. Ertl. Oscillatory kinetics in heterogeneous catalysis. Chemical Reviews,95(3):697–733, 1995.[34] R. Kapral and K. Showalter. Chemical Waves and Patterns, volume 10. Springer Science &Business Media, 2012.[35] K. Kawasaki and T. Ohta. Kink dynamics in one-dimensional nonlinear systems. Physica A,116:573–593, 1982.[36] J. P. Keener and J. Sneyd. Mathematical Physiology. Part I: Cellular Physiology. Springer, NewYork, 1998.[37] J. P. Keener and J. Sneyd. Mathematical Physiology. Part II: Systems Physiology. Springer Sci-ence+Business Media, New York, 2008.[38] E. F. Keller and L. A. Segel. Model for chemotaxis. Journal of Theoretical Biology, 30(2):225–234, 1971.[39] E. Knobloch. Spatial localization in dissipative systems. Annual Review of Condensed MatterPhysics, 6(1):325–359, 2015.[40] S. Kondo. The reaction-diffusion system: a mechanism for autonomous pattern formation in theanimal skin. Genes to Cells, 7(6):535–541, 2002.[41] I. Lengyel and I. R. Epstein. A chemical approach to designing Turing patterns in reaction-diffusion systems. Proceedings of the National Academy of Sciences, 89(9):3977–3979, 1992.[42] C.-M. Lin, T. X. Jiang, R. E. Baker, P. K. Maini, R. B. Widelitz, and C.-M. Chuong. Spots andstripes: pleomorphic patterning of stem cells via p-ERK-dependent cell chemotaxis shown byfeather morphogenesis and mathematical simulation. Developmental Biology, 334(2):369–382,2009.[43] D. J. B. Lloyd and H. OFarrell. On localised hotspots of an urban crime model. Physica D,253:23–39, 2013.[44] D. Lo Jacono, A. Bergeon, and E. Knobloch. Magnetohydrodynamic convectons. Journal of FluidMechanics, 687:595–605, 2011.[45] D. Lo Jacono, A. Bergeon, and E. Knobloch. Spatially localized radiating diffusion flames.Combustion and Flame, 176:117–124, 2017.[46] Y.-P. Ma, J. Burke, and E. Knobloch. Defect-mediated snaking: A new growth mechanism forlocalized structures. Physica D, 239(19):1867–1883, 2010.[47] P. K. Maini, R. E. Baker, and C.-M. Chuong. The Turing model comes of molecular age. Science(New York, NY), 314(5804):1397–1398, 2006.8 of 30
REFERENCES [48] P. K. Maini, H. G. Othmer, et al. Mathematical Models for Biological Pattern Formation. Springer,2001.[49] H. Meinhardt. Morphogenesis of lines and nets. Differentiation, 6(2):117–123, 1976.[50] D. Menshykau, O. Michos, C. Lang, L. Conrad, A. P. McMahon, and D. Iber. Image-based mod-eling of kidney branching morphogenesis reveals GDNF-RET based Turing-type mechanism andpattern-modulating WNT11 feedback. Nature Communications, 10(1):1–13, 2019.[51] E. Meron. Nonlinear Physics of Ecosystems. CRC Press, 2015.[52] E. Meron. Vegetation pattern formation: The mechanisms behind the forms. Physics Today,72(11):30–36, 2019.[53] R. J. Metzger, O. D. Klein, G. R. Martin, and M. A. Krasnow. The branching programme of mouselung development. Nature, 453(7196):745–750, 2008.[54] E. E. Morrisey and B. L. M. Hogan. Preparing for the first breath: Genetic and cellular mechanismsin lung development. Developmental Cell, 18(1):8–23, 2010.[55] J. D. Murray. Mathematical Biology: I. An Introduction, volume 17. Springer Science & BusinessMedia, 2007.[56] B. N. Nagorcka. Evidence for a reaction-diffusion system as a mechanism controlling mammalianhair growth. Biosystems, 16(3-4):323–332, 1983.[57] J. Nagumo, S. Arimoto, and S. Yoshizawa. An active pulse transmission line simulating nerveaxon. Proceedings of the IRE, 50(10):2061–2070, 1962.[58] Q. Ouyang and H. L. Swinney. Transition from a uniform state to hexagonal and striped Turingpatterns. Nature, 352(6336):610–612, 1991.[59] P. Parra-Rivas, D. Gomila, L. Gelens, and E. Knobloch. Bifurcation structure of localized statesin the Lugiato–Lefever equation with anomalous dispersion. Physical Review E, 97(4):042204,2018.[60] P. Parra-Rivas, E. Knobloch, D. Gomila, and L. Gelens. Dark solitons in the Lugiato–Lefeverequation with normal dispersion. Physical Review A, 93(6):063839, 2016.[61] L. M. Pismen. Patterns and Interfaces in Dissipative Dynamics. Springer, 2006.[62] B. C. Ponedel and E. Knobloch. Forced snaking: Localized structures in the real Ginzburg–Landauequation with spatially periodic parametric forcing. The European Physical Journal Special Topics,225(13):2549–2561, 2016.[63] B. Pradenas, I. Araya, M. G. Clerc, C. Falc´on, P. Gandhi, and E. Knobloch. Slanted snaking oflocalized Faraday waves. Physical Review Fluids, 2(6):064401, 2017.[64] H.-G. Purwins, H. U. B¨odeker, and Sh. Amiranashvili. Dissipative solitons. Advances in Physics,59(5):485–701, 2010.
EFERENCES
29 of 30[65] R. Richter and I. V. Barashenkov. Two-dimensional solitons on the surface of magnetic fluids.Physical Review Letters, 94(18):184503, 2005.[66] M. Roth-Kleiner, T. M. Berger, M. R. Tarek, P. H. Burri, and J. C. Schittny. Neonataldexamethasone induces premature microvascular maturation of the alveolar capillary network.Developmental dynamics: An official publication of the American Association of Anatomists,233(4):1261–1271, 2005.[67] D. Ruiz-Reyn´es, L. Mart´ın, E. Hern´andez-Garc´ıa, E. Knobloch, and D. Gomila. Patterns, localizedstructures and fronts in a reduced model of clonal plant growth. arXiv:2001.00224, 2020.[68] U. Thiele, A. J. Archer, M. J. Robbins, H. Gomez, and E. Knobloch. Localized states in theconserved Swift–Hohenberg equation with cubic nonlinearity. Physical Review E, 87(4):042915,2013.[69] U. Thiele, T. Frohoff-H¨ulsmann, S. Engelnkemper, E. Knobloch, and A. J. Archer. First orderphase transitions and the thermodynamic limit. New Journal of Physics, 21(12):123021, 2019.[70] A. Turing. The chemical basis of morphogenesis. Philosophical Transactions of the Royal SocietyB, 237:37–72, 1952.[71] O. Tzuk, H. Uecker, and E. Meron. The role of spatial self-organization in the design of agro-forestry systems. PLoS One, 15(7):e0236325, 2020.[72] V. D. Varner and C. M. Nelson. Cellular and physical mechanisms of branching morphogenesis.Development, 141(14):2750–2759, 2014.[73] V. D. Varner and C. M. Nelson. Computational models of airway branching morphogenesis. InSeminars in Cell & Developmental Biology, volume 67, pages 170–176. Elsevier, 2017.[74] N. Verschueren and A. R. Champneys. A model for cell polarization without mass conservation.SIAM Journal on Applied Dynamical Systems, 16(4):1797–1830, 2017.[75] N. Verschueren and A. R. Champneys. Dissecting the snake: The transition from localized patternsto isolated spikes in pattern formation systems. arXiv:1809.07847, 2018.[76] D. Warburton. Order in the lung. Nature, 453(7196):733–734, 2008.[77] P. D. Woods and A. R. Champneys. Heteroclinic tangles and homoclinic snaking in the unfoldingof a degenerate reversible Hamiltonian–Hopf bifurcation. Physica D, 129(3-4):147–170, 1999.[78] H. Xu, M. Sun, and X. Zhao. Turing mechanism underlying a branching model for lung morpho-genesis. PloS One, 12(4), 2017.[79] Y. Yao, M. Jumabay, A. Wang, and K. I. Bostr¨om. Matrix GLA protein deficiency causes arteri-ovenous malformations in mice. The Journal of Clinical Investigation, 121(8):2993–3004, 2011.[80] Y. Yao, S. Nowak, A. Yochelis, A. Garfinkel, and K. I. Bostr¨om. Matrix GLA protein, an inhibitorymorphogen in pulmonary vascular development. Journal of Biological Chemistry, 282(41):30131–30142, 2007.0 of 30