In a reaction-diffusion-advection system, with a convectively unstable regime, a perturbation creates a wave train that is advected downstream and eventually leaves the system. We show that the convective instability coexists with a local absolute instability when a fixed boundary condition upstream is imposed. This boundary induced instability acts as a continuous wave source, creating a local periodic excitation near the boundary, which initiates waves travelling both up and downstream. To confirm this, we performed analytical analysis and numerical simulations of a modified Martiel-Goldbeter reaction-diffusion model with the addition of an advection term. We provide a quantitative description of the wave packet appearing in the convectively unstable regime, which we found to be in excellent agreement with the numerical simulations. We characterize this new instability and show that in the limit of high advection speed, it is suppressed. This type of instability can be expected for reaction-diffusion systems that present both a convective instability and an excitable regime. In particular, it can be relevant to understand the signaling mechanism of the social amoeba *Dictyostelium discoideum* that may experience fluid flows in its natural habitat.

In a reaction-diffusion-advection system, one or more species are carried away by a flowing medium with an externally imposed velocity. Therefore, the conditions of the system upstream become important to the phenomena observed downstream. In this work, we present the effects of adding an absorbing fixed boundary condition at the upstream end of the system. We focus on the convectively unstable regime, where a perturbation applied to the system dies out in the laboratory reference frame, while it grows in a moving one. By fixing the upstream boundary condition, the system becomes unstable, producing a trigger wave that travels upstream, and a wave train propagating downstream. The trigger wave is absorbed when it reaches the upstream boundary, then the system destabilizes again, and the phenomenon repeats. In 2-D simulations, the trigger wave propagating against the flow has a triangular shape, similar to the concentration profiles exhibiting a cusp in auto-catalytic advection reactions.^{1,2} The here reported mechanism can be expected to be applicable to other reaction-diffusion-advection systems in order to produce a continuous, periodic influx of wave trains.

## I. INTRODUCTION

Many out of equilibrium phenomena in nature can be described by reaction-diffusion systems. This includes the Belousov-Zhabotinsky reaction,^{3,4} electrical impulse dynamics in the heart,^{5} skin patterns in fish,^{6} calcium dynamics in oocytes,^{7} and slime mold aggregation,^{8} among others. In many cases, the active components of such reactions might be subjected to advective flows, which cause new kinds of instabilities.^{9} The most commonly studied types of these instabilities are of convective or absolute nature.^{10,11} Both types of instabilities have been observed in simulations,^{9,12} as well as in experiments such as the Belousov-Zhabotinsky reaction.^{13,14}

Due to the advective nature of the flow, the upstream boundary conditions have important consequences for the spatio-temporal dynamics downstream. Most studies have been performed with no-flux boundary conditions or periodic boundaries, which simplifies the analysis by going into a comoving reference frame. Under these boundaries, an initial perturbation creates a growing wave train^{15,16} whose wavelengths and velocities depend on the particular characteristics of the system. However, the comoving frame analysis is impossible with a Dirichlet (fixed) boundary condition. In particular, an absorbing (zero amplitude) boundary condition corresponds to a one dimensional defect and is the one dimensional equivalent of a spiral center in excitable systems.^{17–19} Up to now, the effects of this type of upstream condition on an advection-diffusion system have received little attention. Preliminary results on such a system were presented by Gholami *et al.*^{20,21} where a continuous influx of wave trains was observed.

Here, we show that in the reaction-diffusion-advection system under study (see below), a particular kind of boundary induced instability occurs when the advection velocities are below a threshold. This boundary condition creates waves periodically with a period dependent on the imposed flow velocity. Unlike the commonly emitted waves by a boundary, these waves do not travel in just one direction (either towards or away from the boundary as is usual in these systems^{22}), but instead two waves appear, one that travels towards and one that travels away from the boundary. In order for this to be possible, these waves do not grow directly at the boundary, but at a finite distance from it. The reaction of the system to this boundary driven instability is also different from the way it reacts to an external perturbation. In this system, a perturbation creates a growing wave train that is advected downstream, while in the absorbing boundary case, the growing instability produces not only a wave train downstream but also a wave travelling upstream.

The downstream wave train is equivalent to the one observed with the no-flux boundary condition. We fully characterized this wave train using linear stability analysis in a moving reference frame and calculated the periodic travelling wave solutions. The upstream travelling wave is the novel feature of this process. This wave travels upstream until it reaches the fixed boundary where it is absorbed, and the process starts again. This process creates wave trains with a period dependent on the imposed flow velocity and thus provides a mechanism to continuously generate wave trains in the fixed reference frame.

To investigate this effect, we performed numerical simulations in one dimension of a model proposed by Martiel and Goldbeter^{23} which are reaction-diffusion equations, with the addition of an advection term due to an imposed external flow. To ensure accuracy in the simulations, we implemented a Runge-Kutta scheme with an adaptable time step based on the Merson error estimation.^{24} To complete the study of the convectively unstable regime, we also performed linear stability analysis of the system in a moving reference frame and periodic travelling wave calculations which we compared with the full nonlinear system solutions. Finally, we performed numerical simulations in 2-Dimensions to study the effect of the flow profile on the boundary induced oscillations. Similar to fronts in advected auto-catalytic reactions,^{1,2} we observed a strong triangular deformation of the trigger wave travelling upstream.

## II. THE REACTION-DIFFUSION-ADVECTION MODEL

Cellular slime moulds are unique organisms positioned between uni- and multi-cellular life in the evolutionary tree. The amoebae of the cellular slime mould *Dictyostelium discoideum* normally live as single cells in forest soil and feed on bacteria. They multiply by binary fission. Starvation induces a developmental program in which up to 10^{5} amoebae aggregate chemotactically to form a multicellular mass,^{25} the so-called slug, that behaves as a single organism and migrates to search for food and better environmental conditions. On failing to find nutrients, the slug culminates into a fruiting body consisting of a stalk and a mass of spores.^{26} Spores are dispersed by rain and small animals and under suitable conditions germinate to release amoebae and the whole cycle starts over again.

Cyclic adenosine monophosphate (cAMP) is the primary chemoattractant for the *D. discoideum* cells during early aggregation. cAMP is emitted from the aggregation centers in a pulsatile manner and surrounding cells detect it by highly specific cAMP receptors.^{27} When cAMP binds to the receptors, it triggers a series of intracellular reactions that activate an enzyme called Adenylate cyclase (ACA), which in turns consumes Adenosine triphosphate (ATP) to produce intracellular cAMP. The cAMP produced inside the cell is partially degraded by intracellular phosphodiesterase and partially transported to the extracellular medium. Phosphodiesterase secreted by the cells degrades extracellular cAMP and suppresses the accumulation of excessive cAMP in the aggregation field (Fig. 1). Since each cell responds to cAMP by moving towards the source of cAMP, by emitting a pulse of cAMP itself, and by a refractory period, the cAMP signal is relayed outward from the aggregation center as a wave.^{28,29} During the refractory period, the amoebas that have detected and produced cAMP do not react to it for a few minutes, and therefore, a new cAMP wave cannot pass during this period. This refractory phase is included 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 cAMP than the second one. Once the receptors bind with cAMP, they change their state to inactive and then slowly change back to their active state. This combination of relay and refractory phase is characteristic of excitable systems and produces target patterns or spirals in two dimensional systems. The geometry of propagating waves is analogous to the spatio-temporal pattern of chemical waves in the Belousov-Zhabotinski reaction.^{3,4}

The model we used for our study was initially proposed by Martiel and Goldbeter^{23} and extended by Tyson *et al.*^{30} In this reaction-diffusion system, the concentration of the signalling chemical cAMP is the activator, while the cAMP receptors on the cell membrane act as inhibitors. Since the inhibitor is cell bounded and we assume that the imposed flow is not strong enough to detach the cells from the substrate,^{31} we add the advection term only to the activator dynamics. A detailed model derivation and the biological correspondence of the model parameters can be found in Ref. 23.

The main equations of the model in its three component version are as follows, where *ρ* stands for the percentage of active receptors on the cell membrane, *γ*, the extracellular concentration of cAMP, and *β*, the intracellular amount of cAMP. The receptor dynamics are given by

with

where *f*_{1} controls the receptor desensitization (change from active to inactive state) and *f*_{2}, the resensitization. The intracellular cAMP is increased by the cAMP production, which in turn depends on the extracellular cAMP and the active receptors. This production is tuned through the rate *σ* at which the activated ACA produces cAMP. The intracellular cAMP is diminished through degradation by intracellular phosphodiesterase at a rate *k _{i}* and passive transport outside of the cell at a rate

*k*

_{t}with

$\lambda 2=(1+\alpha \theta )/(\u03f5(1+\alpha ))$, and $\lambda 1=\lambda \theta /\u03f5$. The extracellular concentration of cAMP *γ* is degraded at a rate *k _{e}* by the extracellular phosphodiesterase and is increased by the transport of cAMP from the intracellular medium

We nondimensionalize the system by introducing dimensionless time and space as $t\u2032=t\xb7k1$ and $x\u2032=x\xb7k1/keD$. Dropping primes and setting $\u03f51=k1/ke,\u2009\u03f5\u2032=k1/(ki+kt)$, we arrive at

Finally, we reduce this system to a two component model which simplifies its theoretical treatment. For this, we assume $\u03f5\u2032$ small, which means that the intracellular cAMP is instantaneously transmitted to the outside media (for a discussion on the 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

*et al.*

^{32}because of their good agreement with experimental measurements. We selected

*σ*and

*k*as control parameters since they account for the production and degradation of extracellular cAMP, respectively. Depending on these two parameters, this system can have one, two, or three steady state solutions, as is shown in the phase diagram in Fig. 2. We focused on the range where only one steady state exists (green, yellow, and blue in Fig. 2). We performed linear stability analysis around this steady state solution $(\gamma 0,\rho 0)$ by setting $\gamma =\gamma 0+\gamma \u2032,\rho =\rho 0+\rho \u2032,$ linearizing, dropping primes, and performing Fourier transform

_{e}we arrive at the dispersion relation

where $\Delta =a11a22\u2212a12a21,\u2009T=a11+a22$,

From here, the different regimes can be distinguished. Starting with the left green area of the phase diagram of Fig. 2 and with no imposed flow velocity *v* = 0, the system has *T* < 0 and $\Delta >0$; therefore, $Re(\omega )<0$ for every *k* and the system is stable. Increasing *k _{e}*, the system has a Hopf bifurcation (at the boundary between yellow and blue area in Fig. 2) and a limit cycle appears (Oscillatory regime). When $v\u22600$, part of the stable regime becomes convectively unstable (yellow in Fig. 2). In this area,

*a*

_{11}is positive and we can calculate the minimum imposed velocity at which the system becomes unstable, by calculating when the real part of

*ω*becomes positive. This gives the following relation:

This is a convex curve dependent on *k* with asymptotes at *k* = 0 and $k=a11/\u03f51$. Its global minimum corresponds to the critical velocity *v _{c}* at which the system destabilizes. This type of instability is of the convective type, which means that although a perturbation applied to the system will die out in the laboratory reference frame, it will grow in a reference frame moving with a speed $v\u2032$, when the system is advected with a flow higher than

*v*. All our simulations were performed in this regime.

_{c}c = 10 | h = 5 | $k1=0.09\u2009min\u22121$ |

$k2=1.665\u2009min\u22121$ | $KR=10\u22127M$ | $ki=1.7\u2009min\u22121$ |

$kt=0.9\u2009min\u22121$ | $\u21131=10$ | $\u21132=0.005$ |

q = 4000 | ϵ = 1 | $\lambda =0.01$ |

$\theta =0.01$ | α = 3 | D = 0.024 $mm2\xb7min\u22121$ |

c = 10 | h = 5 | $k1=0.09\u2009min\u22121$ |

$k2=1.665\u2009min\u22121$ | $KR=10\u22127M$ | $ki=1.7\u2009min\u22121$ |

$kt=0.9\u2009min\u22121$ | $\u21131=10$ | $\u21132=0.005$ |

q = 4000 | ϵ = 1 | $\lambda =0.01$ |

$\theta =0.01$ | α = 3 | D = 0.024 $mm2\xb7min\u22121$ |

Before proceeding to the characterization of the boundary driven instability, we perform a general description of the wave trains present in this system.

## III. NO-FLUX BOUNDARY CONDITION

In the convectively unstable regime, when the advection velocity *v* is above the critical value *v _{c}* [calculated as the minimum of Eq. (4)], a perturbation creates a peak that is advected downstream. This peak creates further peaks behind it, producing a wave train, as can be observed in Fig. 3. The front of this wave train travels with a speed

*v*higher than the imposed flow

_{f}*v*, while the rear of the wave train travels with a velocity $vb<v$. This difference between

*v*and

_{b}*v*translates into the wave train growing in size and having more peaks as time passes. These velocities are indicated by colored lines in Fig. 3. The characteristics of these wave trains can be estimated by taking the Fourier transform in a moving reference frame $y=x\u2212v\u2032t$, where $v\u2032$ is a free parameter

_{f}with $k\u2208\u2102$ and $\omega (k)$ given by the dispersion relation, Eq. (3). According to the method of steepest descents,^{10} the long term behavior of this integral is given by the saddle point of the term accompanying *t*, i.e.,

Since *ω* is also complex, we can use the Cauchy-Riemann Equations

where $k=kr+iki$ and $\omega =\omega r+i\omega i$. This gives pairs of solutions $(k,v\u2032)$, each with its growing rate $\lambda r=\omega r\u2212kiv\u2032$. A typical curve *λ _{r}* vs $v\u2032$ is shown in Fig. 4. The maximum of this curve corresponds to the group velocity of the wave train, it is the fastest growing mode and has $ki=0,kr\u22600$. To calculate the edges of the wave train, the relevant values are the pairs with zero growing rate, because these will correspond to the first and last points at which the system destabilizes and therefore mark the boundaries of the velocity range at which the wave train can be observed. There are two velocities $v\u2032$ with zero growing rate, the lower corresponds to

*v*and the higher to

_{b}*v*. This linear calculation has a very good agreement with the velocities calculated from the numerical simulations of the full nonlinear system, Eq. (2). This is shown in Fig. 5 where these two data sets are compared.

_{f}It has been shown that for some systems, the previously used method may not catch the fastest growing mode in a moving reference frame.^{33,34} For this, the more reliable Briggs collision criterion^{35} is recommended. However, in our system the function $\omega (k)$ has only two local maxima and one unstable branch for real *k*, and under these conditions, the saddle point approach is enough to find all the unstable points.^{36}

To connect to other results in literature, it is worth mentioning that in our calculation, *v _{f}* is equivalent to the spreading speed to the right of the system,

^{37,38}which means that it is the supreme of the velocities $v\u2032$, such that the system is unstable in the comoving frame moving at $v\u2032$.

To characterize each individual peak velocity *v _{p}*, we studied the periodic travelling wave solutions of this system. These waves are characteristic of oscillatory systems

^{39}and have the property $\gamma (z+T)=\gamma (z)$ with $z=x\u2212ct$ for a certain combination of propagation velocity

*c*and period

*T*. The wave calculation and stability analysis were performed using the software Wavetrain.

^{40–42}

We found a range of velocities *c* at which the periodic travelling wave solutions exist. Inside this range, there is a band of velocities *c* where they are stable. The velocities of each individual peak fall into this band as shown in Fig. 6. The selection of a particular wave solution depends on the initial conditions.

The velocity of each particular peak *v _{p}* is higher than the front velocity; therefore, each peak moves forward in the train until it approaches the front, where it has to slow down until it matches

*v*, the velocity of the front of the wave train. Since wavelength and velocity are uniquely linked, the peaks closer to the front of the wave train have a smaller wavelength than the rest of the train. This creates a traffic jam where more peaks start to accumulate in this shorter wavelength area at the front of the train. A similar process has been observed in other reaction-diffusion systems.

_{f}^{43,44}

## IV. FIXED UPSTREAM BOUNDARY CONDITION

We performed numerical simulations with a Dirichlet (fixed) boundary condition upstream $\gamma (x=0)=\rho (x=0)=0$ in the convectively unstable regime. We found that for very high flow speeds, the advection dominates over the diffusion and the system reaches a stable extended steady state. This state can be approximated in powers of $\delta =\u03f51/v2$ with the time independent version of Eq. (2) as

where $x\u2032=x/v$ and $\gamma (x\u2032=0)=0$. The first two terms of the expansion were calculated taking $\gamma =\phi 0+\delta \phi 1$

This solution connects smoothly the zero boundary condition with the steady state of the system. This approximation matches quite well with the full solution as it is shown in Fig. 7.

We performed numerical linear analysis of this monotone solution and found that it becomes unstable at smaller velocities (when *δ* gets larger). The eigenvalues cross the real axis with non-zero imaginary part when the imposed flow velocity *v* is lowered below a threshold. This bifurcation is shown in Fig. 8. The fastest growing eigenvector has the shape of a peak centered close to the fixed border, the distance between the peak and the border increasing with increasing imposed flow velocity.

To study this instability, we performed numerical simulations with Dirichlet boundary condition upstream and small imposed flow velocities. We observed that the system initially reaches a state similar to the one showed in Fig. 7, that is, a smooth connection between the boundary and the steady 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 wave train as was previously described in Sec. III. The second peak travels upstream until it reaches the boundary. Once the upstream travelling peak has been absorbed by the boundary, the system goes back to the smooth solution, which then again becomes unstable and repeats the cycle. This whole process generates periodically wave trains propagating downstream, as shown in Fig. 9.

The period of these perturbations is hard to measure downstream due to the wave train that it generates, whose period is given by the periodic travelling wave solution. To solve this, we measured the period of the initial destabilization peak at its point of creation, as shown in white in Fig. 9. This nucleation location moves farther away from the boundary as the imposed flow velocity increases. This relation is shown in Fig. 10.

This period *T* does not appear to have a relation to any of the periods in the train wave previously studied. This, combined with the difference in the back and front velocities *v _{b}* and

*v*, produces phase slips. The phase slips occur when the front of the newly generated wave train catches up with the back of the previous wave train, thus forming downstream one larger wave train with phase slips. This process is highlighted on the inset of Fig. 9.

_{f}As expected, the velocity of the upstream travelling peak *v _{u}* decreases with imposed flow velocity. Since the new wave does not appear until the previous one has travelled up to the boundary, the instability period is directly affected by the velocity of the upstream peak. Therefore, the period increases with increasing imposed flow velocities. All these dependencies are shown in Fig. 10.

The periodical travelling wave solutions selected by the boundary condition could not be measured for every velocity value, because of the interaction between the new wave train and the old one. This interaction produces numerous phase slips that change the wavelength along the wave train. For those values of *v* where it was measured, the selected travelling wave falls into the stable range shown in Fig. 6.

We understand the upstream travelling peak as a trigger wave, analogous to the ones present in the excitable regime in this system. Trigger waves are non-linear excitation waves that propagate in excitable media when a perturbation above a threshold is applied. In these systems, small perturbations damp out but supra-threshold ones are amplified and excite the neighboring area allowing for wave propagation.^{45} A trigger wave has a velocity which is nonlinearly selected by the system. Another important characteristic is that a new trigger wave cannot enter the system until some recovery time has elapsed. In the case of the upstream travelling wave, the system cannot sustain another upstream travelling peak until the old one has reached the boundary and the cAMP close to the boundary has been washed away. Schematically, the wave works as follows. The cells closer to the boundary have been exposed to very small amounts of cAMP because it is initially washed away due to the boundary. As a result, they have a very high percentage of active receptors on the cell membrane. Therefore, they quickly react to the small perturbation of cAMP produced by the growing peak, emitting cAMP themselves and producing a trigger wave. It has been shown that trigger waves can travel against imposed flows when the advection is not too strong, experimentally in the Belousov-Zhabotinsky reaction^{46} and numerically in the excitable regime of the Martiel-Goldbeter model^{47} and in the FitzHugh-Nagumo model.^{48}

## V. 2-DIMENSIONAL RESULTS

To study the instability already investigated in one dimension, we performed numerical simulations in a 2-Dimensional system. The dimensions were chosen following the *D. discoideum* experiments of Gholami *et al.*^{31} In this microfluidic setup, the amoebas were placed in a 30 mm× 2 mm × 100 *μ*m channel, where a constant flow was applied along the longest axis. Because of the small height and velocities of this system, the flow can be assumed to be laminar and constant in the long channel axis (*x*-axis), thus making a Poiseuille flow. We solved the Navier-Stokes equation under these assumptions and used this flow as our imposed advection for the simulations. The resulting flow is parabolic in the short axis (*z*-axis). This is the direction over which we averaged to have a 2-Dimensional system. In the *xy*-plane, the flow is almost planar in the center with a sharp boundary layer of the order of 50 *μ*m on the top (*y* = 2 mm) and on the bottom (*y* = 0 mm) boundaries, where the velocity quickly drops to zero. The detailed flow profile calculation is presented in the Appendix.

We performed numerical simulations with no-flux boundary conditions on the top and bottom boundaries and Dirichlet [$\rho ,\gamma (x=0)=0$] boundary condition upstream. The simulations confirmed our previous observations in one dimension: When a small advection flow is applied, an instability appears, which creates a wave train downstream and a travelling peak upstream. This process can be observed in Fig. 11, the destabilization peak begins to appear in Fig. 11(b), creating a train wave. The back travelling wave is already visible in Fig. 11(d) and more clear in Fig. 11(f).

Remarkable in comparison with the 1-Dimensional simulations are the range of existence of the instability and the shape of the upstream travelling peak. In the 2-Dimensional simulations, we observed that the system becomes stable at a higher speed *v* = 1.75 mm/min compared to the 1-Dimensional ones *v* = 1.33 mm/min, when measured at the center of the channel ($\sigma =0.45$ and $ke=3.0$). We attribute this difference to the smaller advection speeds at the boundary layer which are enough to destabilize the whole system. This phenomenon was also observed in some preliminary simulations using a parabolic advection flow,^{21} where the advection flow velocity is much smaller in a wider region, thus making the instability range of existence much larger.

Of particular interest is the shape that the upstream travelling peak acquires while it travels towards the boundary. Since this peak travels against the flow, its shape gets deformed due to the different speeds along the perpendicular axis. When the peak originally appears, it has a much flatter shape, similar to the imposed flow, as can be observed at the far right of Fig. 12. As the peak travels upstream (towards the left), it gets increasingly deformed until it acquires a triangular shape. Contours of the peak taken every 0.5 min are displayed in Fig. 12 showing this process.

The triangular deformation of a front due to an adverse flow was theoretically predicted by Edwards^{1} and experimentally confirmed by Leconte *et al.*^{2} for an auto-catalytic reaction. The main difference with our system is that in our reaction-diffusion-advection system, only the activator *γ* is advected, while the inhibitor *ρ* remains static. Like in those systems, the deformation of the wave is larger at larger imposed flows. This is shown in Fig. 13 for three different advection velocities. This wave deformation makes the characterization of the system difficult, because it produces different arrival times at the boundary. More work is needed in this direction to fully characterize this system in 2-D.

## VI. CONCLUSIONS

### A. No-flux boundary

We have analyzed and characterized the convectively unstable regime in the model proposed by Martiel and Goldbeter for cAMP production in *D. discoideum*. In this regime, an initial perturbation generates a wave train of growing size (i.e., it contains more peaks as time passes) that travels downstream. In particular, the speed of the peaks located near the wave front (back) is higher (lower) than the advection flow, thus causing the growing size. These two velocities were numerically characterized through linear stability analysis and have an excellent agreement with the velocities measured in the nonlinear simulations of the model. The growing mode on the center of the wave train corresponds to one of the periodic travelling wave solutions of 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 velocity. As a result of this smaller speed, the wavelength near the front of the train is smaller than in the center of the wave train, thus producing a traffic jam. These wave trains are similar to the differential flow induced convective instability (DIFICI) waves which were first predicted in Ref. 9, experimentally observed in Ref. 13 and further investigated in Refs. 14–16.

### B. Fixed boundary

When slow advection speeds are applied along with a Dirichlet (absorbing) boundary condition, an instability appears that periodically produces wave trains. This instability initially generates a peak that divides into two, with one peak travelling upstream towards the boundary and the other one producing a wave train downstream. Once the peak travelling upstream has reached the absorbing boundary, the process starts again, thus acting as a continuous source of waves. The velocity of the wave travelling upstream is affected by the imposed flow velocity. As expected, it travels slower at higher advection, and since the instability does not appear until the peak reaches the boundary, this affects the period of the oscillation. The faster the imposed flow, the longer the period. The location of appearance of this instability also increases with the advected flow velocity. We emphasize that these upstream travelling waves are nonlinearly selected and depend solely on the system parameters.

This instability was also observed in 2-Dimensional simulations, where the upstream travelling peak acquires the triangular shape of fronts propagating against adverse flows.^{1} This triangular shape increases its height with increasing advection flow. The instability persists up until higher velocities than in one dimension and similarly increases period with increased imposed flow.

The observed phenomena is different from other wave trains emitted by Dirichlet boundary conditions^{22} in that the waves are not directly emitted or absorbed by the boundary, but instead appear as a pair of waves from a nucleation point which exists downstream from the boundary. From this pair of waves, one travels upstream and is absorbed by the boundary while the other creates a wave train downstream. We expect that a similar mechanism may exist in systems where the convective or absolute unstable regime exists close to an excitable regime, thus facilitating the creation of an upstream travelling peak. This mechanism can then be used to produce a constant wave influx.

## ACKNOWLEDGMENTS

The authors acknowledge A. Bae for fruitful discussions and 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 jointly funded by the Federal Ministry of Education and Research of Germany and the Max Planck Society.

### APPENDIX: POSEUILLE FLOW CALCULATION

To estimate the flow profile inside the channel, we used the Navier-Stokes equations and assumed incompressible flow in a 3D rectangular geometry $(x\u2208[0,L],y\u2208[\u2212c,c],z\u2208[\u2212b,b])$, with zero velocity as the boundary condition along the two shortest directions, $u(y=\xb1c)=0$ and $u(z=\xb1b)=0$, thus obtaining

where bold text denotes vectors, *μ* is the system viscosity, *ρ* its density, *p* the pressure, and $u$ the fluid velocity. We further simplified by assuming that the flow is constant over time, only exists in the $x\u0302$-direction, and is constant over this direction, that is, $u=ux\u0302$ and $\u2202xu=\u2202tu=0$; therefore, the previous equation reduces to

where *G* is an externally applied pressure difference. This can be solved by setting the auxiliary function

which reduces the system to solve $\u22072F=0$ with boundary conditions $F(z=\xb1b)=0$ and $F(y=\xb1c)=\u2212G(b2\u2212z2)/2\mu $. Using variable separation $F(y,z)=Fy(y)Fz(z)$ and considering the symmetry of the system, we obtain $Fz=cos\u2009(kzz),Fy=cosh(kyy)$ with *k _{z}* =

*k*. The boundary condition $Fz(z=\xb1b)=0$ sets

_{y}with *m* an odd integer. The other boundary condition is fulfilled by using Fourier series, finally obtaining

where

with *c* = 1 mm and *b* = 50 *μ*m. Since in our case the *z* direction is much shorter than the others, this is the length that sets the boundary layer in the system. Therefore, the flow looks parabolic in the *z*-axis, but almost planar in the *y*-axis, with a very sharp drop to zero close to the boundaries.