Drifting solitary waves in a reaction-diffusion medium with differential advection
aa r X i v : . [ n li n . PS ] M a y Drifting solitary waves in a reaction-diffusion medium with differential advection
Arik Yochelis and Moshe Sheintuch Department of Chemical Engineering, Technion – Israel Institute of Technology, Haifa 32000, Israel (Received October 31, 2018)Propagation of solitary waves in the presence of autocatalysis, diffusion, and symmetry breaking (differential)advection, is being studied. The focus is on drifting (propagating with advection) pulses that form via a convec-tive instability at lower reaction rates of the autocatalytic activator, i.e. the advective flow overcomes the fastexcitation and induces a drifting fluid type behavior. Using spatial dynamics analysis of a minimal case model,we present the properties and the organization of such pulses. The insights underly a general understanding oflocalized transport in simple reaction-diffusion-advection models and thus provide a background to potentialchemical and biological applications.
PACS numbers: 47.35.Fg, 82.40.Ck, 47.20.Ky, 47.54.-r
Solitary waves are prominent generic solutions to reaction-diffusion (RD) systems and basic to many applied science dis-ciplines [1]. In one spatial dimension, these spatially localizedpropagating pulses are qualitatively described by a fast exci-tation (leading front) from a rest state followed by a slow re-covery (rear front) to the same uniform state [1]. Thus, inisotropic RD media a single symmetric supra-threshold lo-calized perturbation results in counter-propagating pulses orwave trains [2].However, in chemical and biological media transport can befacilitated by both diffusion and advection and thus excitationproperties of solitary waves can be subjected to convectiveinstabilities [3]. Nevertheless, theoretical foundations of soli-tary waves in the presence of a symmetry breaking advectivetransport, i.e. dynamics in a differential reaction-diffusion-advection (RDA) media, are yet to be established. Among afew reported examples, it was only shown both experimen-tally and numerically (with no underlying theoretical basis),that excitable pulses can persist in RDA with a propagationdirection against the advective field [4], i.e. upstream.In this Letter we analyze an RDA case model and demon-strate that under low reaction rates, solitary waves may be-come convectively unstable and thus drift (see Fig. 1), i.e. theslow recovery becomes a leading front. We reveal the regionsand the properties of such drifting pulses and show that thephenomenon underlies a competition between a local kineticsof the activator and a differential advection. Our methods in-clude a bifurcation theory of coexisting spatial solutions (lin-ear analysis and numerical continuations) coupled to temporalstability; all the results well agree with direct numerical inte-grations. At the end, we discuss the potential applicability ofour findings to chemical and biological media.
Model setup.–
We start with a minimal RDA model thatincorporates local kinetics of activator v ( x, t ) and inhibitor u ( x, t ) type: u t + su x = f ( u, v ) − u, (1) Le v t + sv x = Bf ( u, v ) − αv + P e − v xx . These dimensionless equations describe a membrane (orcross-flow) reactor, with continuous feeding and cooling in which an exothermic reaction takes place f ( u, v ) ≡ Da (1 − u ) exp [ γv/ ( γ + v )] [6]. Eq. (1) admits a uniform rest state ( u, v ) = ( u , v ≡ Bu /α ) , where u obtained via Da = u (1 − u ) − exp [ − γu / ( γα/B + u )] . In what follows, weset P e = 15 , s = 1 , α = 4 , γ = 10000 , and use Le , Da and B as control parameters allowed to vary; parameter definitionsare given in [7].A standard linear stability analysis to periodic perturba-tions, shows that the uniform states ( u , v ) may loose sta-bility to two finite wavenumber Hopf instabilities, Da ± , thatemerge from ( B W , Da W ), as shown in Fig. 1; the instabili-ties are of a drifting type, i.e. in direction of advection. Thisanomaly arises due to the broken reflection symmetry of left-right traveling waves that is preserved in RD systems. Wenote that traveling waves T W − bifurcate (nonlinearly) sub-critically from Da − while traveling waves T W + bifurcate su-percritically from Da + [5]. While the region Da − < Da Propagation of solitary waves.– To reveal the propagationproperties and the regimes of solitary waves (see Fig. 1), welook at the steady state version of (1) in a comoving frame, ξ = x − ct [5]: u ξ = ( s − c ) − [ Daf ( u, v ) − u ] , v ξ = w, (2) w ξ = P e [( s − cLe ) w − BDaf ( u, v ) + αv ] . The advantage is that existence of nonuniform states can benow analyzed via spatial dynamics methods, i.e. where spaceis viewed as a time-like variable. Thus, solitary waves [in thecontext of (1)] become in (2) asymmetric homoclinic orbits B* 9 10 10.6 11 120 0.1Da*0.20.30.4 B D a (a)(b)(c)TW HopfDriftingPulses ExcitablePulsesSP Hopf Da W B W B Da + Da − FIG. 1: (color online) Top panel: Regions of excitable (thick dashedline) and drifting (thin dashed line) solitary waves (pulses) in a reac-tion rate parameter space ( B, Da ) at Le = 100 ; the thin dashed lineimplies zero velocity of a pulse. The solid line mark the onsets of fi-nite wavenumber instabilities of traveling waves, T W ± . The dottedline marks the criterion for stationary periodic ( SP ) solutions (seetext for details). The ( • ), marks the leftmost limits of homoclinicorbits ( B ∗ , Da ∗ ) ≃ (8 . , . and asymmetric finite wavenum-ber Hopf bifurcation ( B W , Da W ) ≃ (10 , . while ( (cid:4) ) marksthe leftmost limit of excitable pulses ( B , Da ∗ ) ≃ (10 . , . .Bottom panel: Space-time plots at B = 10 . and (a) Da = 0 . , (b) Da = 0 . , and (c) Da = 0 . . The plots show v ( x, t ) resultingfrom integration of Eq. (1) with no-flux boundary conditions, where x ∈ [0 , and t ∈ [0 , ; we used a top-hat initial conditionembedded in ( u , v ) at the respective Da values. ( HO ) and T W ± (which will be also discussed) correspondto periodic orbits undergoing Hopf bifurcations at Da ± (witha proper c ). In the following all these solutions will be com-puted numerically using a continuation package AUTO [8],where the speed c is obtained as a nonlinear eigenvalue prob-lem. Then temporal stability of such steady states will be cal-culated employing an eigenvalue problem via Eq. (1) in a co-moving frame.In Fig. 2(top panel), we present the branches of HO at Da = Da ∗ ≃ . (a horizontal cut in top panel in Fig. 1),resulting via a simultaneous variation of ( B, c ) . B = B ∗ ≃ . identifies a fold, where the stable branch corresponds to B* 9 10 11 12−4−202468 x 10 −3 B c B* 10 113456 B v m a x ξ v B b B b B 10 11 12−4−3.5−3−2.5−2−1.5−1−0.500.5 B c ξ v B FIG. 2: (color online) Bifurcation diagrams showing the branchesof homoclinic solutions as a function of B at Da ∗ ≃ . and Le = 100 (top panel) and Le = 1 (bottom panel). The branchesare plotted in terms of the propagation speed and the maximal valueof v ( ξ ) (top inset); solid lines indicate linear stability. Bottom insetsshow profiles of v ( ξ ) at the two folds, marked by ( • ). The brancheswere obtained via integration of Eq. (2) while the stable portionsof each branch coincide with solutions obtained by integration ofEq. (1); here the periodic domain is L = 24 but results are identicalalso on larger domains. large amplitude HO (see inset). The drifting pulses exist for B ∗ < B < B ≃ . since both branches have positivespeeds, and have similar profiles (see inset) as the standardexcitable pulses. Namely, drifting pulses propagate in the di-rection of the advection (downstream, c > ) where the lead-ing front is the oscillatory tail that was a trailing tail above B = B , for excitable pulses (upstream, c < ).Drifting pulses are expected at low reaction rate regimes ofthe activator, represented in (1) by dimensionless rate constant( Da ) and exothermocity ( B ). Under such conditions the exci-tation of nearest neighbors is suppressed due to the advectiveflow (a nonlinear convective instability) and thus the pulse af-ter speed reversal is no longer excitable since the leading frontnow develops from the rest state as a small amplitude pertur-bation. This scenario changes once the differential advectionis eliminated ( Le = 1 ), in this case a typical RD behavior is Re µ I m µ BB b −3 c L −3 cB=8.8 B=9.6 FIG. 3: (color online) Top panel: Schematic representation of typicaleigenvalue configurations about the uniform state ( u , v , corre-sponding to a saddle if B < B b and a saddle-focus if B > B b , where B b ≃ . is the Belyakov point. Bottom panel: Typical dispersionrelations that associated with the respective eigenvalues computedstarting from the stable homoclinic orbits (see top panel). Parame-ters as in the top panel of Fig 2. restored. While the c = 0 line for Le = 1 , in the ( B, Da )plane doesn’t change, we show in Fig. 2(bottom panel) thatnear the fold only a negative velocity region forms, i.e. a stan-dard excitable pulses are being restored (stability of the pulsesdoes not play a qualitative role).Nevertheless, drifting pulses in presence of a differentialadvection inherit the properties of excitable pulses, as demon-strated by monotonic and nonmonotonic dispersion relationsin Fig. 3. The latter are important characteristics of organi-zation and interaction of solitary waves [10] and are distin-guished here around B = B b ≃ . , a so called Belyakovpoint [11]. At this point and with an appropriate speed, thespatial eigenvalues [of Eq. (2)] correspond to one positive real(associated with ξ → −∞ ) and a degenerate pair of nega-tive reals (associated with ξ → ∞ ). Below B b , the degen-eracy is removed but the eigenvalues remain negative reals(a saddle) while above B b they become complex conjugatedcorresponding to a saddle-focus (a Shil’nikov type HO [12]),marked by ( × ) in top panel of Fig. 3. Importantly, such an in-terchange of eigenvalues implies a transition from monotonicto oscillatory dispersion relation [Fig. 3(bottom panel)] and amonotonic (in space) approach of the HO to the fixed point as ξ → ±∞ , which implies coexistence of bounded-pulse statesfor B > B b [10]. Organization of drifting states.– A standard theory of soli-tary waves qualitatively predicts an organization of HO to beaccompanied by periodic solutions [12]. Here the dispersionrelations obtained at B < B W [Fig. 3(bottom panel)], indeedimply existence of periodic orbits although the uniform state islinearly stable. These periodic solutions are in fact T W − that −3 B c B v m a x FIG. 4: (color online) Bifurcation diagram showing the branches oftraveling waves ( T W − ) as a function of B in terms of speed andthe maximal value of v ( ξ ) (in the inset), at B = 10 . (dark line)and B = 10 . , where Da − ≃ . , k c ≃ . , c ≃ . and Da − ≃ . , k c ≃ . , c ≃ . , respectively. Solid linesimply linear stability to long wave lengths perturbations [9], while( • ) marks the respective onsets of the linear finite wavenumber Hopfbifurcation to T W − . Integration details as in the top panel of Fig. 2but on distinct periodic domains. bifurcate subcritically from the locus of points Da = Da − for B > B W [with distinct critical wavenumbers and speedsobtained from the linear analysis of Eq. (1)], as shown by twoexamples in Fig. 4. Notably, there are infinite number of such T W − families. Unlike the HO , stability of T W − solutionsdo depend on domain size [9].The organization of all drifting nonuniform solutions canbe understood by varying Da at two representative B val-ues. Fig. 5(a), shows a bifurcation diagram of nonuniformsolutions at B ≃ . : while T W ± propagate downstream.The single pulse HO branch ends at the two rightmost ends(marked by dots), at which the profiles take the form of ho-moclinic tails (see bottom inset) [13]. Due to the proximityto the subcritical onset of T W − at Da − , the two rightmostends ever approach each other as domain ( L ) is increased, andconsequently, they inherit the propagation direction of the topand the bottom branches of T W − as discussed in [5]. As B is decreased below B W the HO and the T W − solutionsorganize in isolas and parts of their stability regions overlap[Fig. 5(b)], implying sensitivity to initial perturbations. Notethat the oscillations of the right tail in the profile had decreased(see bottom inset), which is consistent with the approach to-wards Belyakov point ( B = B b ). Conclusions and prospects.– We have showed that solitarywaves can propagate bidirectionally (without changing theirshape) due to a competition between activator autocatalysisand a symmetry breaking advection. Consequently, we distin-guish between excitable (upstream or against advection) anddrifting (downstream or with advection) propagations. Theformer is a characteristic behavior of RD systems and persistswhile the reaction rate of the activator is dominant (analogues Da v m a x −3 Da c ξ v TW − Da − Da + TW + TW + HO HO (a) (u ,v )TW − Da v m a x −3 Da c ξ v HO HO (b) (u ,v )TW − TW − FIG. 5: (color online) Bifurcation diagram showing the branches ofuniform states ( u , v ), homoclinic orbits ( HO ), and traveling waves( T W ± ) as a function of Da in terms of the maximal value of v ( ξ ) at (a) B = 10 . and (b) B = 9 . . Solid lines imply linear sta-bility, including stability of T W ± to long wave lengths perturba-tions [9], while Da ± mark the onsets of the linear finite wavenum-ber Hopf bifurcation to T W ± , respectively. The top inset representsthe nonuniform states in terms of speed while the large (small) isolacorresponds to the T W − family emerging from Da − at B = 10 . ( B = 10 . ). The bottom shows HO profiles at locations marked by( • ); in (a) the two dots mark also the two ends of the HO branch.Integration details as in the top panel of Fig. 2 but on distinct periodicdomains. to front dynamics [14]). While the latter is a consequenceof low excitation and thus subjected to a nonlinear convec-tive instability resulting in a fluid type behavior. Through abifurcation analysis of spatial extended steady states arisingin a minimal RDA model, we revealed the properties and theorganization of drifting pulses. Since the results center onhomoclinic orbits which known to act as organizing centersof spatial solutions, qualitative applicability to systems withother autocatalytic properties is naturally anticipated.Up-to-date only excitable (upstream) solitary waves havebeen observed experimentally in an autocatalytic RDA sys-tem [4], nevertheless chemical media operated in cross-flow(membrane) tubular reactors [6] or on a rotating disks [15], are the most natural setups to confirm our predictions and exploretechnological directions. Moreover, theoretical insights ex-plored here can be related to a profound puzzle of large intra-cellular particles (organelles) self-organization, in eucaryoticcells [16]. For example, localized aggregations of myosin-Xwithin the filopodia have been observed to propagate bidirec-tionally [17] and from the modeling point of view argued to bedriven by both diffusion and differential advection [18]. Con-sequently, a theoretical framework integrating autocatalytickinetics and distinct transport, is paramount to promoting amechanistic understanding of spatiotemporal trafficking of in-tracellular molecular aggregations.We thank to N. Gov and M. Naoz for the helpful discussionon molecular motors. This work was supported by the US-Israel Binational Science Foundation (BSF) and A.Y. was alsopartially supported by the Center for Absorption in Science,Israeli Ministry of Immigrant Absorbtion. M.S. is a memberof the Minerva Center of Nonlinear Dynamics and ComplexSystems. [1] E. Meron, Phys. Rep. , 1 (1992); M.C. Cross and P.C. Ho-henberg, Rev. Mod. Phys. , 851 (1993).[2] A. Yochelis, E. Knobloch, Y. Xie, Z. Qu, and A. Garfinkel, Eu-rophys. Lett. , 64005 (2008).[3] J.M. Chomaz, Phys. Rev. Lett. , 1931 (1992).[4] M. Kærn and M. Menzinger, Phys. Rev. E , 046202 (2002).[5] A. Yochelis and M. Sheintuch, e-print:arXiv:nlin.PS/0810.4690; e-print: arXiv:nlin.PS/0902.2688.[6] O. Nekhamkina, B.Y. Rubinstein, and M. Sheintuch, AIChE J. , 1632 (2000).[7] Le (Lewis number), is the ratio of solid- to fluid-phase heat ca-pacities, P e (P´eclet number), is the ratio of convective to con-ductive enthalpy fluxes, and Da (Damk¨ohler number), is thedimensionless rate constant [6].[8] E. Doedel et al. , AUTO2000: Continuation and bifurcationsoftware for ordinary differential equations (with HOMCONT) ,http://indy.cs.concordia.ca/auto/.[9] Temporal stability of T W ± was computed for large periodicdomains, L = nλ c , n > , λ c ≡ π/k c (until the onset didn’tchange with n ), via a standard numerical eigenvalue method us-ing Eq. (1) in a comoving frame. k c is the critical wavenumberat the Hopf onset.[10] C. Elphick, E. Meron, J. Rinzel, and E.A. Spiegel, J. Theor.Biol. , 249 (1990); M. Or-Guil, I.G. Kevrekidis, and M. B¨ar,Physica D , 154 (2000).[11] L.A. Belyakov, Mat. Zametki , 910 (1980).[12] P. Glendinning and C.T. Sparrow, J. Stat. Phys. , 645 (1984);N.J. Balmforth, Annu. Rev. Fluid Mech. 27, 335 (1995); A.R.Champneys et al. , SIAM J. Appl. Dyn. Syst. , 663 (2007).[13] J. Sneyd, A. LeBeaub, and D. Yule, Physica D , 158 (2000).[14] V. Yakhnin and M. Menzinger, Chem. Eng. Sci. , 4559(2002); M. Sheintuch, Y. Smagina, and O. Nekhamkina, Ind.Eng. Chem. Res. , 2136 (2002).[15] Y. Khazan and L.M. Pismen, Phys. Rev. Lett. , 4318 (1995).[16] M.A. Welte, Curr. Biol. , R525 (2004).[17] J.S. Berg and R.E. Cheney, Nature Cell. Biol. , 246 (2002).[18] D.A. Smith and R.M. Simmons, Biophys. J.80