Convective Instability and Boundary Driven Oscillations in a Reaction-Diffusion-Advection Model
Estefania Vidal-Henriquez, Vladimir Zykov, Eberhard Bodenschatz, Azam Gholami
CConvective instability and boundary driven oscillations in a reaction-diffusion-advection model
Estefania Vidal-Henriquez, Vladimir Zykov, Eberhard Bodenschatz, and Azam GholamiCitation: Chaos , 103110 (2017); doi: 10.1063/1.4986153View online: http://dx.doi.org/10.1063/1.4986153View Table of Contents: http://aip.scitation.org/toc/cha/27/10Published by the American Institute of Physics onvective instability and boundary driven oscillationsin a reaction-diffusion-advection model Estefania Vidal-Henriquez,
Vladimir Zykov, Eberhard Bodenschatz, and Azam Gholami Max Planck Institute for Dynamics and Self-Organization, Am Fassberg 17, D-37077 G € ottingen, Germany Institute for Nonlinear Dynamics, University of G € ottingen, D-37073 G € ottingen, Germany Laboratory of Atomic and Solid-State Physics and Sibley School of Mechanical and Aerospace Engineering,Cornell University, Ithaca, New York 14853, USA (Received 2 June 2017; accepted 25 September 2017; published online 9 October 2017)In a reaction-diffusion-advection system, with a convectively unstable regime, a perturbationcreates a wave train that is advected downstream and eventually leaves the system. We show thatthe convective instability coexists with a local absolute instability when a fixed boundary conditionupstream is imposed. This boundary induced instability acts as a continuous wave source, creatinga local periodic excitation near the boundary, which initiates waves travelling both up anddownstream. To confirm this, we performed analytical analysis and numerical simulations of amodified Martiel-Goldbeter reaction-diffusion model with the addition of an advection term. Weprovide a quantitative description of the wave packet appearing in the convectively unstableregime, which we found to be in excellent agreement with the numerical simulations. We charac-terize this new instability and show that in the limit of high advection speed, it is suppressed. Thistype of instability can be expected for reaction-diffusion systems that present both a convectiveinstability and an excitable regime. In particular, it can be relevant to understand the signalingmechanism of the social amoeba
Dictyostelium discoideum that may experience fluid flows in its natu-ral habitat. V C .https://doi.org/10.1063/1.4986153 In a reaction-diffusion-advection system, one or morespecies are carried away by a flowing medium with anexternally imposed velocity. Therefore, the conditions ofthe system upstream become important to the phenom-ena observed downstream. In this work, we present theeffects of adding an absorbing fixed boundary conditionat the upstream end of the system. We focus on the con-vectively unstable regime, where a perturbation appliedto the system dies out in the laboratory reference frame,while it grows in a moving one. By fixing the upstreamboundary condition, the system becomes unstable, pro-ducing a trigger wave that travels upstream, and a wavetrain propagating downstream. The trigger wave isabsorbed when it reaches the upstream boundary, thenthe system destabilizes again, and the phenomenonrepeats. In 2-D simulations, the trigger wave propagatingagainst the flow has a triangular shape, similar to theconcentration profiles exhibiting a cusp in auto-catalyticadvection reactions.
The here reported mechanism canbe expected to be applicable to other reaction-diffusion-advection systems in order to produce a continuous, peri-odic influx of wave trains.
I. INTRODUCTION
Many out of equilibrium phenomena in nature can bedescribed by reaction-diffusion systems. This includes theBelousov-Zhabotinsky reaction, electrical impulse dynam-ics in the heart, skin patterns in fish, calcium dynamics inoocytes, and slime mold aggregation, among others. Inmany cases, the active components of such reactions mightbe subjected to advective flows, which cause new kinds ofinstabilities. The most commonly studied types of theseinstabilities are of convective or absolute nature.
Bothtypes of instabilities have been observed in simulations, as well as in experiments such as the Belousov-Zhabotinskyreaction.
Due to the advective nature of the flow, the upstreamboundary conditions have important consequences for thespatio-temporal dynamics downstream. Most studies havebeen performed with no-flux boundary conditions or periodicboundaries, which simplifies the analysis by going into acomoving reference frame. Under these boundaries, an initialperturbation creates a growing wave train whose wave-lengths and velocities depend on the particular characteris-tics of the system. However, the comoving frame analysis isimpossible with a Dirichlet (fixed) boundary condition. Inparticular, an absorbing (zero amplitude) boundary conditioncorresponds to a one dimensional defect and is the onedimensional equivalent of a spiral center in excitable sys-tems.
Up to now, the effects of this type of upstream a) Electronic mail: [email protected] b) Electronic mail: [email protected] V C Author(s) 2017. , 103110-1 CHAOS , 103110 (2017) ondition on an advection-diffusion system have received lit-tle attention. Preliminary results on such a system were pre-sented by Gholami et al. where a continuous influx ofwave trains was observed.Here, we show that in the reaction-diffusion-advectionsystem under study (see below), a particular kind of bound-ary induced instability occurs when the advection velocitiesare below a threshold. This boundary condition createswaves periodically with a period dependent on the imposedflow velocity. Unlike the commonly emitted waves by aboundary, these waves do not travel in just one direction(either towards or away from the boundary as is usual inthese systems ), but instead two waves appear, one thattravels towards and one that travels away from the boundary.In order for this to be possible, these waves do not growdirectly at the boundary, but at a finite distance from it. Thereaction of the system to this boundary driven instability isalso different from the way it reacts to an external perturba-tion. In this system, a perturbation creates a growing wavetrain that is advected downstream, while in the absorbingboundary case, the growing instability produces not only awave train downstream but also a wave travelling upstream.The downstream wave train is equivalent to the oneobserved with the no-flux boundary condition. We fullycharacterized this wave train using linear stability analysis ina moving reference frame and calculated the periodic travel-ling wave solutions. The upstream travelling wave is thenovel feature of this process. This wave travels upstreamuntil it reaches the fixed boundary where it is absorbed, andthe process starts again. This process creates wave trainswith a period dependent on the imposed flow velocity andthus provides a mechanism to continuously generate wavetrains in the fixed reference frame.To investigate this effect, we performed numerical sim-ulations in one dimension of a model proposed by Martieland Goldbeter which are reaction-diffusion equations, withthe addition of an advection term due to an imposed externalflow. To ensure accuracy in the simulations, we implementeda Runge-Kutta scheme with an adaptable time step based onthe Merson error estimation. To complete the study of the convectively unstable regime, we also performed linear sta-bility analysis of the system in a moving reference frame andperiodic travelling wave calculations which we comparedwith the full nonlinear system solutions. Finally, we per-formed numerical simulations in 2-Dimensions to study theeffect of the flow profile on the boundary induced oscilla-tions. Similar to fronts in advected auto-catalytic reac-tions, we observed a strong triangular deformation of thetrigger wave travelling upstream.
II. THE REACTION-DIFFUSION-ADVECTION MODEL
Cellular slime moulds are unique organisms positionedbetween uni- and multi-cellular life in the evolutionary tree.The amoebae of the cellular slime mould
Dictyostelium dis-coideum normally live as single cells in forest soil and feedon bacteria. They multiply by binary fission. Starvation indu-ces a developmental program in which up to 10 amoebaeaggregate chemotactically to form a multicellular mass, the so-called slug, that behaves as a single organism andmigrates to search for food and better environmental condi-tions. On failing to find nutrients, the slug culminates into afruiting body consisting of a stalk and a mass of spores. Spores are dispersed by rain and small animals and undersuitable conditions germinate to release amoebae and thewhole cycle starts over again.Cyclic adenosine monophosphate (cAMP) is the primarychemoattractant for the
D. discoideum cells during earlyaggregation. cAMP is emitted from the aggregation centersin a pulsatile manner and surrounding cells detect it byhighly specific cAMP receptors. When cAMP binds to thereceptors, it triggers a series of intracellular reactions thatactivate an enzyme called Adenylate cyclase (ACA), whichin turns consumes Adenosine triphosphate (ATP) to produceintracellular cAMP. The cAMP produced inside the cell ispartially degraded by intracellular phosphodiesterase andpartially transported to the extracellular medium.Phosphodiesterase secreted by the cells degrades extracellu-lar cAMP and suppresses the accumulation of excessivecAMP in the aggregation field (Fig. 1). Since each cell
FIG. 1. Schematic representation ofthe model proposed by Martiel andGoldbeter for production and relay ofcAMP in
D. discoideum . The extracel-lular cAMP binds with the membranereceptors at a rate K R for those in theactive state, activating ACA which inturn produces intracellular cAMP thatis transported to the extracellularmedia at a rate k t . After binding, thereceptors change to their inactive format a rate k which has a lower probabil-ity of binding with cAMP ( K D ). Theextracellular cAMP is then degradedvia phosphodiestrase at a rate k e .Reproduced from Ref. 20. et al. Chaos , 103110 (2017) esponds to cAMP by moving towards the source of cAMP,by emitting a pulse of cAMP itself, and by a refractoryperiod, the cAMP signal is relayed outward from the aggre-gation center as a wave. During the refractory period,the amoebas that have detected and produced cAMP do notreact to it for a few minutes, and therefore, a new cAMPwave cannot pass during this period. This refractory phase isincluded in the model in terms of the membrane receptors.These receptors are present in two states, active and inactive.The first one has a higher probability to bind with cAMPthan the second one. Once the receptors bind with cAMP,they change their state to inactive and then slowly changeback to their active state. This combination of relay andrefractory phase is characteristic of excitable systems andproduces target patterns or spirals in two dimensional sys-tems. The geometry of propagating waves is analogous tothe spatio-temporal pattern of chemical waves in theBelousov-Zhabotinski reaction.
The model we used for our study was initially pro-posed by Martiel and Goldbeter and extended by Tyson et al. In this reaction-diffusion system, the concentrationof the signalling chemical cAMP is the activator, while thecAMP receptors on the cell membrane act as inhibitors.Since the inhibitor is cell bounded and we assume that theimposed flow is not strong enough to detach the cells fromthe substrate, we add the advection term only to the acti-vator dynamics. A detailed model derivation and the bio-logical correspondence of the model parameters can befound in Ref. 23.The main equations of the model in its three compo-nent version are as follows, where q stands for the per-centage of active receptors on the cell membrane, c , theextracellular concentration of cAMP, and b , the intracel-lular amount of cAMP. The receptor dynamics are givenby @ t q ¼ (cid:2) k f ð c Þ q þ k f ð c Þð (cid:2) q Þ with f ð c Þ ¼ þ jc þ c ; f ð c Þ ¼ L þ j L c c þ c c ; where f controls the receptor desensitization (change fromactive to inactive state) and f , the resensitization. The intra-cellular cAMP is increased by the cAMP production, whichin turn depends on the extracellular cAMP and the activereceptors. This production is tuned through the rate r atwhich the activated ACA produces cAMP. The intracellularcAMP is diminished through degradation by intracellularphosphodiesterase at a rate k i and passive transport outsideof the cell at a rate k t @ t b ¼ q ra U ð q ; c Þ = ð þ a Þ (cid:2) ð k i þ k t Þ b with U ð q ; c Þ ¼ k þ Y k þ Y ; Y ¼ qc þ c ; k ¼ ð þ ah Þ = ð (cid:2) ð þ a ÞÞ , and k ¼ kh =(cid:2) . The extracellularconcentration of cAMP c is degraded at a rate k e by theextracellular phosphodiesterase and is increased by the trans-port of cAMP from the intracellular medium @ t c ¼ D r c (cid:2) v (cid:3) r c þ k t b = h (cid:2) k e c : We nondimensionalize the system by introducing dimen-sionless time and space as t ¼ t (cid:3) k and x ¼ x (cid:3) k = ffiffiffiffiffiffiffiffi k e D p .Dropping primes and setting (cid:2) ¼ k = k e ; (cid:2) ¼ k = ð k i þ k t Þ ,we arrive at @ t q ¼ (cid:2) f ð c Þ q þ f ð c Þð (cid:2) q Þ ; (1a) (cid:2) @ t b ¼ q ra U ð q ; c Þ = ðð þ a Þð k i þ k t ÞÞ (cid:2) b ; (1b) @ t c ¼ (cid:2) r c (cid:2) v (cid:3) r c þ ð k t b = ð hk e Þ (cid:2) c Þ =(cid:2) : (1c)Finally, we reduce this system to a two component modelwhich simplifies its theoretical treatment. For this, we assume (cid:2) small, which means that the intracellular cAMP is instanta-neously transmitted to the outside media (for a discussion onthe validity of this approximation, refer to Refs. 23 and 30).We then arrive at the two component Martiel-Goldbeter,which we will use during the rest of this paper @ t c ¼ (cid:2) r c (cid:2) v (cid:3) r c þ ð s U ð q ; c Þ (cid:2) c Þ =(cid:2) ; (2a) @ t q ¼ (cid:2) f ð c Þ q þ f ð c Þð (cid:2) q Þ ; (2b)where s ¼ qk t ar = ð k e ð k t þ k i Þ h ð þ a ÞÞ . All used parametersare listed in Table I and were selected as suggested byLauzeral et al. because of their good agreement withexperimental measurements. We selected r and k e as con-trol parameters since they account for the productionand degradation of extracellular cAMP, respectively.Depending on these two parameters, this system can haveone, two, or three steady state solutions, as is shown inthe phase diagram in Fig. 2. We focused on the rangewhere only one steady state exists (green, yellow, and bluein Fig. 2). We performed linear stability analysis aroundthis steady state solution ð c ; q Þ by setting c ¼ c þ c ; q ¼ q þ q ; linearizing, dropping primes, and performingFourier transform ð c ; q Þ ¼ ð ð c k ; q k Þ e x ð k Þ t þ ikx dk ; we arrive at the dispersion relation0 ¼ x þ x ð(cid:2) T þ (cid:2) k þ i v k Þ þ D (cid:2) a ð (cid:2) k þ i v k Þ ; (3)where D ¼ a a (cid:2) a a ; T ¼ a þ a , TABLE I. Parameters used for simulations of Eq. (2). c ¼ h ¼ k ¼ :
09 min (cid:2) k ¼ :
665 min (cid:2) K R ¼ (cid:2) M k i ¼ : (cid:2) k t ¼ : (cid:2) L ¼ L ¼ : q ¼ (cid:2) ¼ k ¼ : h ¼ : a ¼ D ¼ (cid:3) min (cid:2) et al. Chaos , 103110 (2017) ¼ s (cid:2) q c ð k (cid:2) k Þð þ c Þ ð k þ Y Þ (cid:2) (cid:2) ; a ¼ s (cid:2) c q ð k (cid:2) k Þð þ c Þ ð k þ Y Þ ; a ¼ ð (cid:2) q Þ j L c (cid:2) c L ð þ c c Þ (cid:2) q ð j (cid:2) Þð þ c Þ ; and a ¼ (cid:2) f ð c Þ (cid:2) f ð c Þ : From here, the different regimes can be distinguished.Starting with the left green area of the phase diagram ofFig. 2 and with no imposed flow velocity v ¼
0, the systemhas T < D >
0; therefore, Re ð x Þ < k andthe system is stable. Increasing k e , the system has a Hopfbifurcation (at the boundary between yellow and blue area inFig. 2) and a limit cycle appears (Oscillatory regime). When v
0, part of the stable regime becomes convectively unsta-ble (yellow in Fig. 2). In this area, a is positive and we cancalculate the minimum imposed velocity at which the systembecomes unstable, by calculating when the real part of x becomes positive. This gives the following relation: v ð k Þ ¼ ð(cid:2) D = a þ (cid:2) k Þð (cid:2) k (cid:2) T Þ k ð a (cid:2) (cid:2) k Þ : (4)This is a convex curve dependent on k with asymptotes at k ¼ k ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffi a =(cid:2) p . Its global minimum corresponds tothe critical velocity v c at which the system destabilizes. Thistype of instability is of the convective type, which meansthat although a perturbation applied to the system will dieout in the laboratory reference frame, it will grow in a refer-ence frame moving with a speed v , when the system isadvected with a flow higher than v c . All our simulations wereperformed in this regime.Before proceeding to the characterization of the bound-ary driven instability, we perform a general description ofthe wave trains present in this system. III. NO-FLUX BOUNDARY CONDITION
In the convectively unstable regime, when the advectionvelocity v is above the critical value v c [calculated as the min-imum of Eq. (4)], a perturbation creates a peak that isadvected downstream. This peak creates further peaks behindit, producing a wave train, as can be observed in Fig. 3. Thefront of this wave train travels with a speed v f higher than theimposed flow v , while the rear of the wave train travels with avelocity v b < v . This difference between v b and v f translatesinto the wave train growing in size and having more peaks astime passes. These velocities are indicated by colored lines inFig. 3. The characteristics of these wave trains can be esti-mated by taking the Fourier transform in a moving referenceframe y ¼ x (cid:2) v t , where v is a free parameter ð c ; q Þ ¼ ð ð c k ; q k Þ e t ð x þ ik v iky dk ; with k C and x ð k Þ given by the dispersion relation, Eq.(3). According to the method of steepest descents, the longterm behavior of this integral is given by the saddle point ofthe term accompanying t , i.e., ddk x ð k Þ þ ik v (cid:3) (cid:4) ¼ : Since x is also complex, we can use the Cauchy-RiemannEquations @ x r @ k r ¼ @ x i @ k i ¼ @ x r @ k i (cid:2) v ¼ @ x i @ k r þ v ¼ ; (5)where k ¼ k r þ ik i and x ¼ x r þ i x i . This gives pairs ofsolutions ð k ; v Þ , each with its growing rate k r ¼ x r (cid:2) k i v . Atypical curve k r vs v is shown in Fig. 4. The maximum ofthis curve corresponds to the group velocity of the wavetrain, it is the fastest growing mode and has k i ¼ ; k r FIG. 2. Phase diagram of the system described by Eq. (2). Stable regime ingreen, in this regime only one solution exists and is stable. In the yellowarea labeled CU exists one steady state that is convectively unstable. In theAU labeled blue area exists one unstable steady state surrounded by a limitcycle. The orange regime marked as excitable presents three steady states,one of which is excitable, while the light blue bistable regime has threesteady states, two of which are stable. FIG. 3. Space-time plot of a simulation performed in the convectively unsta-ble regime using no-flux (Neumann) boundary condition. The wave train isgenerated by an initial perturbation and is advected downstream (to theright) by the imposed flow. The relevant velocities present in the wave trainare highlighted, these are the velocity of wave train rear v b in black,front velocity v f in white, and individual peak velocity v p in yellow. Allnumerical simulations were performed using k e ¼ : (cid:2) and r ¼ :
45 min (cid:2) . et al. Chaos , 103110 (2017) re the pairs with zero growing rate, because these will cor-respond to the first and last points at which the system desta-bilizes and therefore mark the boundaries of the velocityrange at which the wave train can be observed. There aretwo velocities v with zero growing rate, the lower corre-sponds to v b and the higher to v f . This linear calculation has avery good agreement with the velocities calculated from thenumerical simulations of the full nonlinear system, Eq. (2).This is shown in Fig. 5 where these two data sets arecompared.It has been shown that for some systems, the previouslyused method may not catch the fastest growing mode in amoving reference frame. For this, the more reliableBriggs collision criterion is recommended. However, inour system the function x ð k Þ has only two local maxima andone unstable branch for real k , and under these conditions,the saddle point approach is enough to find all the unstablepoints. To connect to other results in literature, it is worth men-tioning that in our calculation, v f is equivalent to the spread-ing speed to the right of the system, which means that it is the supreme of the velocities v , such that the system isunstable in the comoving frame moving at v .To characterize each individual peak velocity v p , westudied the periodic travelling wave solutions of this system.These waves are characteristic of oscillatory systems andhave the property c ð z þ T Þ ¼ c ð z Þ with z ¼ x (cid:2) ct for a cer-tain combination of propagation velocity c and period T . Thewave calculation and stability analysis were performed usingthe software Wavetrain. We found a range of velocities c at which the periodictravelling wave solutions exist. Inside this range, there is aband of velocities c where they are stable. The velocities ofeach individual peak fall into this band as shown in Fig. 6.The selection of a particular wave solution depends on theinitial conditions.The velocity of each particular peak v p is higher than thefront velocity; therefore, each peak moves forward in thetrain until it approaches the front, where it has to slow downuntil it matches v f , the velocity of the front of the wave train.Since wavelength and velocity are uniquely linked, the peakscloser to the front of the wave train have a smaller wave-length than the rest of the train. This creates a traffic jamwhere more peaks start to accumulate in this shorter wave-length area at the front of the train. A similar process hasbeen observed in other reaction-diffusion systems. IV. FIXED UPSTREAM BOUNDARY CONDITION
We performed numerical simulations with a Dirichlet(fixed) boundary condition upstream c ð x ¼ Þ ¼ q ð x ¼ Þ ¼ d ¼ (cid:2) = v with the timeindependent version of Eq. (2) as d @ x x c ¼ @ x c (cid:2) ð s U ð q ; c Þ (cid:2) c Þ =(cid:2) ; q ¼ f ð c Þ = ð f ð c Þ þ f ð c ÞÞ ; FIG. 4. Growth rate in the different reference systems v for v ¼ x -axis mark the back and front velocities of thewave train. For these parameters v b ¼ :
00 mm/min and v f ¼ :
13 mm = min.FIG. 5. Dependence of the wave train velocities on the imposed advectionflow. The lower value corresponds to the back of the train, that is, the firstpoint that destabilises, while the higher value corresponds to the front of thewave train, that is, the last point that destabilises. The continuous line corre-sponds to the prediction obtained by the linear analysis, and the dots are thevalues obtained from the simulations of the full nonlinear equations. FIG. 6. Region of existence of periodic travelling wave solutions. v standsfor the imposed advection velocity, and c , the velocity of the periodic travel-ling wave. The solutions exist above the black line and are stable in the bandbetween the red lines. The dots correspond to the solution selected by thesystem in the middle of the wave train in our numerical simulations. et al. Chaos , 103110 (2017) here x ¼ x = v and c ð x ¼ Þ ¼
0. The first two terms of theexpansion were calculated taking c ¼ u þ du ¼ @ x u (cid:2) ð s U ð u Þ (cid:2) u Þ =(cid:2) ; (6a) @ x x u ¼ @ x u (cid:2) u s d U d c (cid:5)(cid:5)(cid:5) c ¼ u (cid:2) (cid:6) (cid:7). (cid:2) : (6b)This solution connects smoothly the zero boundary conditionwith the steady state of the system. This approximationmatches quite well with the full solution as it is shown inFig. 7.We performed numerical linear analysis of this mono-tone solution and found that it becomes unstable at smallervelocities (when d gets larger). The eigenvalues cross thereal axis with non-zero imaginary part when the imposedflow velocity v is lowered below a threshold. This bifurcationis shown in Fig. 8. The fastest growing eigenvector has theshape of a peak centered close to the fixed border, the dis-tance between the peak and the border increasing withincreasing imposed flow velocity.To study this instability, we performed numerical simu-lations with Dirichlet boundary condition upstream andsmall imposed flow velocities. We observed that the systeminitially reaches a state similar to the one showed in Fig. 7,that is, a smooth connection between the boundary and thesteady state. However, this solution becomes unstable producing a peak which, as it grows, divides into two peaks.One of the peaks travels downstream and produces a wavetrain as was previously described in Sec. III. The secondpeak travels upstream until it reaches the boundary. Once theupstream travelling peak has been absorbed by the boundary,the system goes back to the smooth solution, which thenagain becomes unstable and repeats the cycle. This wholeprocess generates periodically wave trains propagatingdownstream, as shown in Fig. 9.The period of these perturbations is hard to measuredownstream due to the wave train that it generates, whoseperiod is given by the periodic travelling wave solution. Tosolve this, we measured the period of the initial destabiliza-tion peak at its point of creation, as shown in white in Fig. 9.This nucleation location moves farther away from the bound-ary as the imposed flow velocity increases. This relation isshown in Fig. 10.This period T does not appear to have a relation to any ofthe periods in the train wave previously studied. This, com-bined with the difference in the back and front velocities v b and v f , produces phase slips. The phase slips occur when thefront of the newly generated wave train catches up with theback of the previous wave train, thus forming downstream FIG. 7. (a) High speed solution withDirichlet boundary condition. Advectionvelocity v ¼ v ¼ u solution of Eq. (6a)in black dotted lined. Scaled space x ¼ x = v . (b) Comparison for approxi-mation at smaller speed. Full solutionwith Dirichlet boundary condition and v ¼ u in green and firstorder approximation c ¼ u þ du indashed black line.FIG. 8. Frequencies of the linear analysis of the monotone profile c ¼ u þ du showing the oscillatory bifurcation when the imposed flow v islowered. v ¼ v ¼ v ¼ v ¼ np , where the destabilization occurs, T ,the oscillations period, and v u , the velocity of the upstream travelling peak.Inset with a zoom of the wave generation area with previously defined quan-tities of the wave train highlighted, v b , velocity of the back of the wave train, v f , velocity of the front of the wave train, and v p , velocity of each individualpeak. The imposed flow velocity is 1.2 mm/min. et al. Chaos , 103110 (2017) ne larger wave train with phase slips. This process ishighlighted on the inset of Fig. 9.As expected, the velocity of the upstream travellingpeak v u decreases with imposed flow velocity. Since the newwave does not appear until the previous one has travelled upto the boundary, the instability period is directly affected bythe velocity of the upstream peak. Therefore, the periodincreases with increasing imposed flow velocities. All thesedependencies are shown in Fig. 10. The periodical travelling wave solutions selected by theboundary condition could not be measured for every velocityvalue, because of the interaction between the new wave trainand the old one. This interaction produces numerous phaseslips that change the wavelength along the wave train. Forthose values of v where it was measured, the selected travel-ling wave falls into the stable range shown in Fig. 6.We understand the upstream travelling peak as a triggerwave, analogous to the ones present in the excitable regimein this system. Trigger waves are non-linear excitation wavesthat propagate in excitable media when a perturbation abovea threshold is applied. In these systems, small perturbationsdamp out but supra-threshold ones are amplified and excitethe neighboring area allowing for wave propagation. Atrigger wave has a velocity which is nonlinearly selected bythe system. Another important characteristic is that a newtrigger wave cannot enter the system until some recoverytime has elapsed. In the case of the upstream travellingwave, the system cannot sustain another upstream travellingpeak until the old one has reached the boundary and thecAMP close to the boundary has been washed away.Schematically, the wave works as follows. The cells closerto the boundary have been exposed to very small amounts ofcAMP because it is initially washed away due to the bound-ary. As a result, they have a very high percentage of activereceptors on the cell membrane. Therefore, they quicklyreact to the small perturbation of cAMP produced by thegrowing peak, emitting cAMP themselves and producing atrigger wave. It has been shown that trigger waves can travelagainst imposed flows when the advection is not too strong,experimentally in the Belousov-Zhabotinsky reaction andnumerically in the excitable regime of the Martiel-Goldbetermodel and in the FitzHugh-Nagumo model. V. 2-DIMENSIONAL RESULTS
To study the instability already investigated in onedimension, we performed numerical simulations in a 2-Dimensional system. The dimensions were chosen followingthe
D. discoideum experiments of Gholami et al. In thismicrofluidic setup, the amoebas were placed in a 30 mm (cid:4) (cid:4) l m channel, where a constant flow wasapplied along the longest axis. Because of the small heightand velocities of this system, the flow can be assumed to belaminar and constant in the long channel axis ( x -axis), thusmaking a Poiseuille flow. We solved the Navier-Stokes equa-tion under these assumptions and used this flow as ourimposed advection for the simulations. The resulting flow isparabolic in the short axis ( z -axis). This is the direction overwhich we averaged to have a 2-Dimensional system. In the xy -plane, the flow is almost planar in the center with a sharpboundary layer of the order of 50 l m on the top ( y ¼ y ¼ q ; c ð x ¼ Þ ¼
0] boundary condition upstream.The simulations confirmed our previous observations in one
FIG. 10. Different properties of the boundary driven oscillations as functionof the imposed flow velocity v . (a) Period of the oscillations, as measured atthe point of nucleation of the instability. (b) Velocity of the upstream travel-ling peak. (c) Destabilization point position. et al. Chaos , 103110 (2017) imension: When a small advection flow is applied, an insta-bility appears, which creates a wave train downstream and atravelling peak upstream. This process can be observed inFig. 11, the destabilization peak begins to appear in Fig.11(b), creating a train wave. The back travelling wave isalready visible in Fig. 11(d) and more clear in Fig. 11(f).Remarkable in comparison with the 1-Dimensional sim-ulations are the range of existence of the instability and theshape of the upstream travelling peak. In the 2-Dimensionalsimulations, we observed that the system becomes stable ata higher speed v ¼ v ¼ r ¼ :
45 and k e ¼ : where the advection flow velocity is much smaller in awider region, thus making the instability range of existencemuch larger.Of particular interest is the shape that the upstream trav-elling peak acquires while it travels towards the boundary.Since this peak travels against the flow, its shape getsdeformed due to the different speeds along the perpendicularaxis. When the peak originally appears, it has a much flattershape, similar to the imposed flow, as can be observed at thefar right of Fig. 12. As the peak travels upstream (towardsthe left), it gets increasingly deformed until it acquires a triangular shape. Contours of the peak taken every 0.5 minare displayed in Fig. 12 showing this process.The triangular deformation of a front due to an adverseflow was theoretically predicted by Edwards and experi-mentally confirmed by Leconte et al. for an auto-catalyticreaction. The main difference with our system is that in ourreaction-diffusion-advection system, only the activator c isadvected, while the inhibitor q remains static. Like in thosesystems, the deformation of the wave is larger at largerimposed flows. This is shown in Fig. 13 for three differentadvection velocities. This wave deformation makes the char-acterization of the system difficult, because it produces dif-ferent arrival times at the boundary. More work is needed inthis direction to fully characterize this system in 2-D. VI. CONCLUSIONSA. No-flux boundary
We have analyzed and characterized the convectivelyunstable regime in the model proposed by Martiel andGoldbeter for cAMP production in
D. discoideum . In thisregime, an initial perturbation generates a wave train ofgrowing size (i.e., it contains more peaks as time passes) thattravels downstream. In particular, the speed of the peakslocated near the wave front (back) is higher (lower) than theadvection flow, thus causing the growing size. These twovelocities were numerically characterized through linear sta-bility analysis and have an excellent agreement with thevelocities measured in the nonlinear simulations of themodel. The growing mode on the center of the wave traincorresponds to one of the periodic travelling wave solutionsof the system and moves faster than the front of the train.Therefore, a peak will move towards the front of the train,where then it will decrease its speed to match the front veloc-ity. As a result of this smaller speed, the wavelength near thefront of the train is smaller than in the center of the wavetrain, thus producing a traffic jam. These wave trains are sim-ilar to the differential flow induced convective instability(DIFICI) waves which were first predicted in Ref. 9, experi-mentally observed in Ref. 13 and further investigated inRefs. 14–16.
FIG. 11. Colormap of the concentration c every 0.9 min starting at t ¼ v ¼ c ¼ : v ¼ t ¼
12 min on the right until t ¼
16 min on the left. Contours not in theiroriginal positions but spatially separated for better visualization. FIG. 13. Shape of the upstream travelling peak for different velocities takenas a contour at c ¼ :
0. Top: v ¼ v ¼ v ¼ x- and y- axis forbetter visualisation. et al. Chaos , 103110 (2017) . Fixed boundary When slow advection speeds are applied along with aDirichlet (absorbing) boundary condition, an instabilityappears that periodically produces wave trains. This instabil-ity initially generates a peak that divides into two, with onepeak travelling upstream towards the boundary and the otherone producing a wave train downstream. Once the peaktravelling upstream has reached the absorbing boundary, theprocess starts again, thus acting as a continuous source ofwaves. The velocity of the wave travelling upstream isaffected by the imposed flow velocity. As expected, it travelsslower at higher advection, and since the instability does notappear until the peak reaches the boundary, this affects theperiod of the oscillation. The faster the imposed flow, thelonger the period. The location of appearance of this instabil-ity also increases with the advected flow velocity. Weemphasize that these upstream travelling waves are nonli-nearly selected and depend solely on the system parameters.This instability was also observed in 2-Dimensionalsimulations, where the upstream travelling peak acquires thetriangular shape of fronts propagating against adverse flows. This triangular shape increases its height with increasingadvection flow. The instability persists up until higher veloc-ities than in one dimension and similarly increases periodwith increased imposed flow.The observed phenomena is different from other wavetrains emitted by Dirichlet boundary conditions in that thewaves are not directly emitted or absorbed by the boundary,but instead appear as a pair of waves from a nucleation pointwhich exists downstream from the boundary. From this pairof waves, one travels upstream and is absorbed by the bound-ary while the other creates a wave train downstream. Weexpect that a similar mechanism may exist in systems wherethe convective or absolute unstable regime exists close to anexcitable regime, thus facilitating the creation of an upstreamtravelling peak. This mechanism can then be used to producea constant wave influx. ACKNOWLEDGMENTS
The authors acknowledge A. Bae for fruitful discussionsand unknown referees for insightful comments. E.V.H.thanks the Deutsche Akademische Austauschdienst (DAAD),Research Grants—Doctoral Programs in Germany. A.G.acknowledges MaxSynBio Consortium, which is jointlyfunded by the Federal Ministry of Education and Research ofGermany and the Max Planck Society.
APPENDIX: POSEUILLE FLOW CALCULATION
To estimate the flow profile inside the channel, we usedthe Navier-Stokes equations and assumed incompressibleflow in a 3D rectangular geometry ð x ; L (cid:5) ; y c ; c (cid:5) ; z b ; b (cid:5)Þ , with zero velocity as the boundary conditionalong the two shortest directions, u ð y ¼ c Þ ¼ u ð z ¼ b Þ ¼
0, thus obtaining q D u Dt ¼ q g (cid:2) r p þ l r u ; where bold text denotes vectors, l is the system viscosity, q its density, p the pressure, and u the fluid velocity. We fur-ther simplified by assuming that the flow is constant overtime, only exists in the ^ x -direction, and is constant over thisdirection, that is, u ¼ u ^ x and @ x u ¼ @ t u ¼
0; therefore, theprevious equation reduces to l @ u @ y þ @ u @ z ! ¼ @ p @ x (cid:6) (cid:2) G ; where G is an externally applied pressure difference. Thiscan be solved by setting the auxiliary function F ¼ u (cid:2) G ð b (cid:2) z Þ l ; which reduces the system to solve r F ¼
F ð z ¼ b Þ ¼
F ð y ¼ c Þ ¼ (cid:2) G ð b (cid:2) z Þ = l . Using variable separation F ð y ; z Þ ¼ F y ð y Þ F z ð z Þ and con-sidering the symmetry of the system, we obtain F z ¼ cos ð k z z Þ ; F y ¼ cosh ð k y y Þ with k z ¼ k y . The boundary condition F z ð z ¼ b Þ ¼ k z ¼ k y ¼ m p b with m an odd integer. The other boundary condition is ful-filled by using Fourier series, finally obtaining u ¼ G ð b (cid:2) z Þ l þ X odd A n cosh n p b y (cid:6) (cid:7) cos n p b z (cid:6) (cid:7) ; where A n ¼ (cid:2) Gb lp n sin n p (cid:6) (cid:7) cosh n p c b (cid:6) (cid:7) with c ¼ b ¼ l m. Since in our case the z direc-tion is much shorter than the others, this is the length thatsets the boundary layer in the system. Therefore, the flowlooks parabolic in the z -axis, but almost planar in the y -axis,with a very sharp drop to zero close to the boundaries. B. F. Edwards, “Poiseuille advection of chemical reaction fronts,” Phys.Rev. Lett. , 104501 (2002). M. Leconte, J. Martin, N. Rakotomalala, and D. Salin, “Pattern of reactiondiffusion fronts in laminar flows,” Phys. Rev. Lett. , 128302 (2003). A. Zaikin and A. Zhabotinsky, “Concentration wave propagation in two-dimensional liquid-phase self-oscillating system,” Nature , 535–537(1970). A. T. Winfree, “Spiral waves of chemical activity,” Science , 634–636(1972). M. A. Allessie, F. I. Bonke, and F. J. Schopman, “Circus movement in rab-bit atrial muscle as a mechanism of tachycardia,” Circ. Res. , 54–62(1973). S. Kondo and R. Asai, “A reaction-diffusion wave on the skin of themarine angelfish pomacanthus,” Nature , 765 (1995). J. Lechleiter, S. Girard, E. Peralta, and D. Clapham, “Spiral calcium wavepropagation and annihilation in xenopus laevis oocytes,” Science ,123–126 (1991). P. Devreotes, “Dictyostelium discoideum: A model system for cell-cellinteractions in development,” Science , 1054–1058 (1989). et al.
Chaos , 103110 (2017) A. B. Rovinsky and M. Menzinger, “Chemical instability induced by a dif-ferential flow,” Phys. Rev. Lett. , 1193 (1992). R. J. Deissler, “Noise-sustained structure, intermittency, and theGinzburg-Landau equation,” J. Stat. Phys. , 371–395 (1985). B. Sandstede and A. Scheel, “Absolute and convective instabilities ofwaves on unbounded and large bounded domains,” Phys. D: NonlinearPhenom. , 233–277 (2000). R. Satnoianu, J. Merkin, and S. Scott, “Interaction between hopf and con-vective instabilities in a flow reactor with cubic autocatalator kinetics,”Phys. Rev. E , 3246 (1998). A. B. Rovinsky and M. Menzinger, “Self-organization induced by the dif-ferential flow of activator and inhibitor,” Phys. Rev. Lett. , 778 (1993). R. Toth, A. Papp, V. Gaspar, J. Merkin, S. Scott, and A. Taylor, “Flow-driven instabilities in the Belousov-Zhabotinsky reaction: Modelling andexperiments,” Phys. Chem. Chem. Phys. , 957–964 (2001). R. A. Satnoianu, J. H. Merkin, and S. K. Scott, “Spatio-temporal structuresin a differential flow reactor with cubic autocatalator kinetics,” Phys. D:Nonlinear Phenom. , 345–367 (1998). R. A. Satnoianu, J. H. Merkin, and S. K. Scott, “Pattern formation in a dif-ferential–flow reactor model,” Chem. Eng. Sci. , 461–469 (2000). P. S. Hagan, “Spiral waves in reaction-diffusion equations,” SIAM J.Appl. Math. , 762–786 (1982). I. S. Aranson and L. Kramer, “The world of the complex Ginzburg-Landau equation,” Rev. Mod. Phys. , 99 (2002). L. M. Pismen,
Vortices in Nonlinear Fields: From Liquid Crystals toSuperfluids, from Non-Equilibrium Patterns to Cosmic Strings (OxfordUniversity Press, 1999), Vol. 100. A. Gholami, O. Steinbock, V. Zykov, and E. Bodenschatz, “Flow-driveninstabilities during pattern formation of Dictyostelium discoideum,” NewJ. Phys. , 063007 (2015). A. Gholami, V. Zykov, O. Steinbock, and E. Bodenschatz, “Flow-driventwo-dimensional waves in colonies of Dictyostelium discoideum,” New J.Phys. , 093040 (2015). J. A. Sherratt, “Periodic travelling wave selection by dirichlet boundaryconditions in oscillatory reaction-diffusion systems,” SIAM J. Appl. Math. , 1520–1538 (2003). J.-L. Martiel and A. Goldbeter, “A model based on receptor desensitizationfor cyclic AMP signaling in Dictyostelium cells,” Biophys. J. , 807(1987). R. Merson, “An operational method for the study of integration proc-esses,” in
Proceedings of the Symposium on Data Processing (1957), pp.1–25. J. Bonner, “Evidence for the formation of aggregates by chemotaxis in thedevelopment of the slime mold Dictyostelium discoideum,” J. Exp. Zool. , 1 (1947). W. B. Loomis,
The Development of Dictyostelium Discoideum (Elsevier,1982). C. Parent and P. Devreotes, “Molecular genetics of signal transduction inDictyostelium,” Annu. Rev. Biochem. , 411–440 (1996). M. Dinauer, T. Steck, and P. Devreotes, “Cyclic, 3 ,5 -AMP relay inDictyostelium discoideum. V. Adaptation of the cAMP signaling responseduring cAMP stimulation,” J. Cell Biol. , 554 (1980). B. M. Shaffer, “Secretion of cyclic AMP induced by cyclic AMP in thecellular slime mould Dictyostelium discoideum,” Nature , 549–552(1975). J. J. Tyson, K. A. Alexander, V. Manoranjan, and J. Murray, “Spiral wavesof cyclic AMP in a model of slime mold aggregation,” Phys. D: NonlinearPhenom. , 193–207 (1989). A. Gholami, O. Steinbock, V. Zykov, and E. Bodenschatz, “Flow-drivenwaves and phase-locked self-organization in quasi-one-dimensional colo-nies of Dictyostelium discoideum,” Phys. Rev. Lett. , 018103 (2015). J. Lauzeral, J. Halloy, and A. Goldbeter, “Desynchronization of cells onthe developmental path triggers the formation of spiral waves of cAMPduring Dictyostelium aggregation,” Proc. Natl. Acad. Sci. , 9153–9158(1997). L. Brevdo, “A study of absolute and convective instabilities with an applica-tion to the eady model,” Geophys. Astrophys. Fluid Dyn. , 1–92 (1988). L. Brevdo, P. Laure, F. Dias, and T. J. Bridges, “Linear pulse structure andsignalling in a film flow on an inclined plane,” J. Fluid Mech. , 37–71(1999). R. J. Briggs,
Electron-Stream Interaction with Plasmas (MIT, 1964). L. Brevdo, “Convectively unstable wave packets in the blasius boundarylayer,” Z. Angew. Math. Mech. , 423–436 (1995). W. Van Saarloos, “Front propagation into unstable states,” Phys. Rep. , 29–222 (2003). M. Holzer and A. Scheel, “Criteria for pointwise growth and their role ininvasion processes,” J. Nonlinear Sci. , 661–709 (2014). N. Kopell and L. Howard, “Plane wave solutions to reaction-diffusionequations,” Stud. Appl. Math. , 291–328 (1973). J. A. Sherratt, “Numerical continuation methods for studying periodictravelling wave (wavetrain) solutions of partial differential equations,”Appl. Math. Comput. , 4684–4694 (2012). J. D. Rademacher, B. Sandstede, and A. Scheel, “Computing absolute andessential spectra using continuation,” Phys. D: Nonlinear Phenom. ,166–183 (2007). J. A. Sherratt, “Numerical continuation of boundaries in parameter spacebetween stable and unstable periodic travelling wave (wavetrain) solutionsof partial differential equations,” Adv. Comput. Math. , 175–192 (2013). C. T. Hamik and O. Steinbock, “Shock structures and bunching fronts inexcitable reaction-diffusion systems,” Phys. Rev. E , 046224 (2002). N. Manz and O. Steinbock, “Dynamics of excitation pulses with attractiveinteraction: Kinematic analysis and chemical wave experiments,” Phys.Rev. E , 066213 (2004). V. S. Zykov and A. T. Winfree,
Simulation of Wave Processes inExcitable Media (John Wiley & Sons, Inc., 1992). K. Agladze, M. Braune, H. Engel, H. Linde, and V. Krinsky, “Autowavepropagation in a Belousov-Zhabotinsky medium with immobilized catalystand stationary flow of reagents,” Z. Phys. Chem. , 79–85 (1991). J. Lindner, H. (cid:2)
Sevcˇ (cid:3) ıkov (cid:3) a, and M. Marek, “Influence of an external electricfield on cAMP wave patterns in aggregating Dictyostelium discoideum,”Phys. Rev. E , 041904 (2001). E. A. Ermakova, E. E. Shnol, M. A. Panteleev, A. A. Butylin, V. Volpert,and F. I. Ataullakhanov, “On propagation of excitation waves in movingmedia: The FitzHugh-Nagumo model,” PloS One , e4454 (2009). et al. Chaos27