Boundary-Induced Pattern Formation from Uniform Temporal Oscillation

Pattern dynamics triggered by fixing a boundary is investigated. By considering a reaction-diffusion equation that has a unique spatially-uniform and limit cycle attractor under a periodic or Neumann boundary condition, and then by choosing a fixed boundary condition, we found three novel phases depending on the ratio of diffusion constants of activator to inhibitor: transformation of temporally periodic oscillation into a spatially-periodic fixed pattern, travelling wave emitted from the boundary, and aperiodic spatiotemporal dynamics. The transformation into a fixed, periodic pattern is analyzed by crossing of local nullclines at each spatial point, shifted by diffusion terms. A spatial map, then, is introduced, whose temporal sequence can reproduce the spatially periodic pattern, by replacing the time with space. The generality of the boundary-induced pattern formation as well as its relevance to biological morphogenesis is discussed.


I. INTRODUCTION
Pattern formation in nonlinear-nonequilibrium systems has gathered much attention over decades both theoretically and experimentally. With the motivation to understand biological morphogeneis, Turing 1 , in his seminal paper, showed how spatial pattern is spontaneously organized from a spatially homogeneous state, without any external signal, as a result of symmetry breaking. One of the six possible pattern formation processes he classified is now known as the Turing pattern. Stationary periodic pattern with a finite wavelength is self-organized from a temporally stationary and spatially uniform state.
Since the pioneering study by Turing, spatiotemporal dynamics in chemical-reaction systems have been extensively investigated, including waves and spirals 2 , as well as spatiotemporal chaos 3 . On the other hand, understanding morphogenesis according to the Turing's thesis needed more time to be accepted in biology.
The Turing pattern observed so far in biology, however, is still restricted in few cases such as the patterns of shell 4 or fish skin 5 . In contrast, relevance of oscillatory dynamics to morphogenesis is now recognized by the discovery of oscillatory gene expressions in the somitogenesis of vertebrate development, and stem cell systems 6,7 . In chick development, temporal oscillation in the gene expression of c-hairy1 is uncovered in the presomitic mesoderm. Such temporal oscillation is later fixed into a spatial pattern 8,9 . For fixation into a spatial pattern, relevance of external gradient was discussed as clock-and-wave-front mechanism 4,10 .
In contrast, we recently found an alternative mechanism for the transformation from temporal oscillation to fixed spatial pattern, induced by diffusive interaction and boundary condition, even without any input of external gradient 11,12 . This finding can provide possible explanation for recent experimental results on somitogenesis, where cell-cell interaction is also suggested to be relevant 13,14 . As the proposed mechanism of ours requires just a fixed boundary and diffusion of chemicals, we expect that it can be applicable to a variety of chemical and biological systems.
Here, we first confirm the mechanism by investigating reaction-diffusion systems, focusing on the influence of boundary conditions to pattern formation. For the present study we consider a system in which spatially uniform and temporally oscillatory state exists as an attractor if a boundary condition allows, and then adopt a fixed boundary condition. In this case, temporal oscillation is ceased at the boundary point. Now, for some range of diffusion constants, homogeneous oscillation is recovered as the spatial point is departed from the boundary. In other cases, however, which we will mainly discuss in the present paper, the temporal oscillation is replaced by a spatially periodic pattern, induced by the fixed boundary condition. Note that the boundary sometime can influence on the global pattern, as was discussed in Ref. 15,16 , while selection of periodic pattern or other spatiotemporal pattern from a homogeneous periodic state is only recently found and will be analyzed here.
After introducing two-component reaction-diffusion equations in §2, we study the pattern formation induced by a fixed boundary condition in §3. Depending on the diffusion constants, we found fixation into a spatially periodic pattern, traveling wave, and spatiotemporal aperiodic dynamics. The wavelength of periodic pattern is not obtained by linear stability analysis as in the Turing pattern. Instead, in §4, we introduce local nullcline analysis that also includes the effect of diffusion. A four-dimensional spatial map is then introduced in §5, which successively determines the fixed pattern from the upper-to down-stream. The periodic attractor of the spatial map determines the organized spatial pattern. Relevance of the present boundary-induced pattern formation to biological problems is discussed in §6.

II. REACTION-DIFFUSION MODEL
We consider a spatially reaction-diffusion system of two components X(r, t) and Y (r, t) on a one dimensional space r.
where f (X, Y ) and g(X, Y ) are the reaction functions for X and Y , D X (D Y ) is the diffusion constant of X (Y ), respectively. In order to study the fixation into spatial pattern from a homogeneous oscillatory state, we assume that the attractor of the dynamical system without the diffusion term, i.e., the dynamical system dX/dt = f (X, Y ), dY /dt = g(X, Y ) is a unique limit cycle. It is generated by a Hopf bifurcation from the fixed point (X * , Y * ). Hence the Jacobi matrix at the fixed point satisfies a + d > 0 and −2 < d−a √ −bc < 2 1 . In all the examples to be studied here, the spatiallyuniform limit cycle is an attractor as long as the boundary condition allows for its existence.
For example, by taking Neumann or periodic boundary conditions, the uniform limit cycle state is the only attractor in the reaction-diffusion equation (1). Throughout the present paper, we discuss such case and how the boundary condition alters an attractor.
To be specific we study the following two systems. The first one (model A) is a simplified from of the gene regulation network, while the second one is Brusselator (model B). The first one is given by where X and Y are the protein concentrations. X inhibits the expression of the protein Y while Y activates the expressions of both X and Y 17,18 . When β > 8, this system has one unstable fixed point (X * , Y * ) = (1/2, 1/2), and the limit cycle is an attractor. Unless otherwise mentioned, we choose β = 40 throughout the paper.
The second one, the so-called Brusselator is given by Ref. 19 : As is well known, the systems has a limit cycle attractor for A 2 + 1 < B < (A + 1) 2 . Here we choose A = 1 and B = 3.

III. PATTERN FORMATION BY BOUNDARY CONDITION
Now we adopt a fixed boundary condition X(0, t) = X 0 and Y (0, t) = Y 0 for the above reaction-diffusion models. The fixed boundary is set at the leftmost point r = 0, and we study how this choice alters an attractor by taking the models (A) and (B). In order to exclude the effect of initial condition to pattern formation dynamics, the uniform initial condition with the boundary values, i.e., X(r, 0) = X 0 and Y (r, 0) = Y 0 is adopted. Unless otherwise mentioned the right boundary is given by Neumann condition, so that the uniform state is not destabilized from that boundary. Now by the fixed boundary at the left side r = 0, we found that the spatially uniform oscillation is not always reached, even if the system goes sufficiently far to the right side (which we call downflow).
Depending on the values of diffusion constants D X and D Y , the behaviors downflow are characterized into the following phases.
I. The oscillation downflow is replaced by a temporally-fixed and spatially-periodic pattern.
II. As the spatial location goes downflow, X and Y oscillate in time. The amplitudes of oscillation increase with r. As r goes larger, the original spatially-uniform, limit cycle attractor is reached.
III. Neither spatially-uniform nor temporally-fixed pattern is reached. Instead, travelling wave is emitted from the boundary, which replaces the original uniform oscillatory state.
IV. Neither temporally periodic nor spatially periodic state is stable. Spatiotemporal aperiodic dynamics replaces the original uniform oscillatory state.
By using the model (A), the typical spatiotemporal patterns in the four phases are shown in Fig.1, which appear depending on the diffusion constants D X and D Y . The phase diagram against the diffusion constants is given in Fig.2, where the phase is determined basically by Now we describe the dynamics in each phase in detail.
, the diffusion of inhibitor is much larger than that of the activator, spatially-uniform limit cycle is replaced by a spatially periodic pattern. When D X is sufficiently larger than D Y , one stripe is shaped from one periodic cycle (see Fig.1a). Indeed this transformation from temporal oscillation to spatial pattern is recently discussed for the extreme case with D Y = 0, i.e., the case in which only the inhibitor diffuses 12 . Here, even if D Y is not zero but very small, this transformation from temporal oscillation to spatial pattern works (see the next section for the mechanism of this transformation). As D Y is increased, the formation of stripe is slowed When one stripe is formed by one periodic cycle, the proportionality between the period of the limit cycle T and the wavelength of fixed pattern λ is expected against the change in the diffusion constants (at least approximately). In this case, by using the dimension analysis 12 , the relationship between the wavelength λ and the oscillation period T is λ ∝ √ D X T when D Y is fixed at a small value. This relationship is shown in Fig.3, against the change in diffusion constants D X . Note that when D X /D Y is decreased and 1-to-1 time-space transformation is replaced by n-to-m locking (n > 1, m ≥ 1) between cycle period and wavelength, the above proportionality is deviated to prolong the wavelength.
II. As D Y is further increased, then the diffusion of activator is sufficiently fast, so that the fixation to a spatial pattern no longer works, as the attraction to a homogeneous state is fast (Fig1d). Near the boundary the oscillation is suppressed, but as the space goes downflow, the (X, Y ) values are shifted and start to oscillate. After such "spatial transient" the state reaches the original limit cycle with a spatially uniform state (see  .1e). The period in the wave is shorter than that of the limit cycle. Near the fixed boundary, the oscillation amplitude is small, and in the present model the oscillation period, then, is shorter than the original limit cycle.
With the diffusion of activator, this fast oscillation propagates to downflow, and hence the travelling wave is generated.
IV. With the further increase in D Y , the global travelling wave turns to be unstable, and aperiodic spatiotemporal pattern replaces the uniform periodic pattern (Fig.1f). The  Then spatiotemporal homogeneous limit cycle is reached after long transient (Fig.5b).
This transient increases rapidly with the system size, suggesting the super-transient behavior 24,25 .
The transition between homogeneous limit cycle state and the fixed pattern is generally observed in reaction-diffusion systems. In the model A, the uniform periodic attractor remains to exist near the Hopf bifurcation point β = 8, but for β > ∼ 12, the fixation to spatially periodic state is observed. In Fig.6, we have plotted the case with model B. Again one stripe is generated from one cycle when D X ≫ D Y . In this case, however, we have observed only the phases I and II.

IV. TRANSFORMATION FROM TEMPORAL OSCILLATION TO SPATIAL PATTERN: EXPLANATION BY NULLCLINE SHIFT
Now we analyze how oscillation is fixed into a periodic spatial pattern. Note that a uniform oscillatory state exists for r > 0 initially. Thus the pattern formation progresses far from a fixed point state. Accordingly, the standard analysis for Turing instability cannot be applied, since it treats the instability of a spatially-uniform and temporally fixed point state against perturbation.
Intuitively, the fixation mechanism we found is understood as follows. In a homogeneous state, the temporal evolution is given by the local dynamics dX/dt = f (X, Y ), dY /dt = g(X, Y ). There are no stable fixed points given by the crossing point of X− and Y − nullclines (see Fig.7). Due to the fixed boundary, however, oscillation stops at r ∼ 0, and there, the values X, Y cannot remain identical with those at the limit cycle attractor. Then the diffusion gives an additional term to each of the local dynamics for a downflow point r > ∼ 0. Accordingly, the nullclines for the local dynamics including such diffusion terms could be shifted from the original, so that they can cross to generate a novel stable fixed point.
Thus a fixed point replaces the oscillatory state at r ∼ 0. This gives a spatial gradient again to further downflow, so that diffusion term gives shift of nullclines for the local dynamics at the downflow. If this shift is sufficient to make the nullclines cross with each other, a fixed point is again generated at the downflow. With the repetition of this process, the oscillatory state is replaced by a fixed spatial pattern. Now we discuss if this scenario really works in the reaction-diffusion systems we studied here. The fixed point pattern should satisfy Eq.(5) describes local nullclines of X and Y for the fixed pattern, at each location r. Now, consider the fixed boundary at r = 0. In the model (A), the term D X ∂ 2 X ∂r 2 at r > 0 gives a horizontal shift in X-nullcline, and the term D Y ∂ 2 Y ∂r 2 gives a diagonal shift in Y-nullcline. In Fig.8, we plotted how these nullclines are shifted. With this shift, a pair of cross-points is generated, one of which give a stable fixed point. The diffusion term successively allows for the emergence of a stable fixed point that varies with r.
To sum up, a stable fixed point is generated by the shift of nullclines, while the chain of such novel fixed points along space generate the diffusion term which self-consistently generates the shift. Now, we study how the transition to the homogeneous limit cycle state occurs in terms of nullcline shift. In the above case, as D X ≫ D Y , the horizontal shift of X-nullcline is The abscissa is X, and the ordinate is Y . The two dashed lines are the original nullclines, i.e., f (X, Y ) = 0 and g(X, Y ) = 0 (see also Fig.7). Distance between each local nullcline and the original nullcline is provided by the local diffusion terms, i.e., D X ∂ 2 X ∂r 2 The nullcline of X is shifted horizontally and the nullcline of Y is shifted diagonally. The stationary state from the initial condition X(r, 0) = Y (r, 0) = 0 is also depicted as the red circle, which lies at the cross point of two nullclines.
dominant, which keeps on generating a pair of novel crossing points of the nullclines, one of which corresponds to a stable fixed point. As the diffusion of activator is increased, however, such cross-point is not sustained with time, as shown in Fig.9. Hence the temporal oscillation remains even at upper flow. Further downflow, the spatial gradient in X(r) or Y (r) decreases, and such small shift can no longer generate a crossing point. At further down flow, the spatial uniformity is reached, and no more shifts of nullclines exist, so that the original flow in Fig.7 is recovered, giving rise to the original uniform limit cycle attractor (see also Fig.5).
With the further increase in D Y , the diagonal shift of nullclines is larger. Again, crossingpoints of nullclines appear at a certain moment at the upper flow. The actual value of (X(r, t), Y (r, t)), however, is deviated from the crossing point, so that the self-consistent condition for the crossing point and the actual diffusion term to give the shift is not satisfied.
In the travelling wave phase (III), however, even though (X(r, t), Y (r, t)) is no longer fixed in time for given r, the solution written as (X(R), Y (R)) with R = r − vt and v as a speed of travelling wave. Then, instead of the condition ∂X/∂t = f + D X ∂ 2 X/∂r 2 = 0 to be satisfied for the fixed pattern, −vdX/dR = f + D X ∂ 2 X/dR 2 has to be satisfied (Y has to satisfy the corresponding equation). In other words, even though the shifted nullclines do not keep on generating a requested crossing point, instead, the fixed point condition for nullcline crossing is expected to be satisfied in an appropriate inertial frame that moves with the speed of the traveling wave.
As D Y is further increased, X-and Y -nullclines can cross at some spatial location but do not remain crossing globally over the space. No moving frame can keep the crossing points of nullclines, either. This leads to the aperiodic spatiotemporal dynamics.
For the model (B) again, the shifts of X-and Y -nullclines by the diffusion term generate a crossing point and accordingly a fixed pattern is formed.

V. SPATIAL MAP
When the fixed pattern exists, ∂X/∂t = 0 and ∂Y /∂t = 0 are satisfied. This condition corresponds to the self-consistent shift of nullclines, and with this condition, the pattern is successively determined from the boundary point. To formulate this fixed-pattern condition, we adopt spatial discretization of the reaction-diffusion equation. Then the condition for a fixed spatial pattern is given by where ℓ is a discretized space index with a as the discretized unit length, i.e., r = aℓ and X ℓ = X(r/a), Y ℓ = Y (r/a). Here, D ′ X,Y is the rescaled diffusion constant which satisfies D X,Y ′ = D X,Y /a 2 . From these equations, we can derive the four-dimensional map as follows: Previously we considered the system with D Y = 0, in which case the mapping is reduced to a two-dimensional map (X ℓ−1 , X ℓ ) → X ℓ+1 , as Y ℓ is determined from g(X ℓ .Y ℓ ) = 0 12 .
This two-dimensional map is more tractable than the above four-dimensional mapping, but instead, additional selection criterion was needed to choose an appropriate branch in the mapping in obtaining Y ℓ from g(X ℓ .Y ℓ ) = 0 (The use of a "spatial map" was introduced earlier in the analysis of a spatial pattern, see also 26,27 ).
From this four-dimensional mapping, the spatial pattern is provided as a time-series generated by the mapping (7), and then by replacing ℓ to spatial sequence r = aℓ, the spatial pattern is obtained.
Examples of spatial sequences derived from Eq.7, for given three different initial conditions are plotted in Fig.10, which reproduces the fixed patterns obtained from direct numerical simulation of the reaction diffusion equation. Now note that the spatially periodic pattern is obtained as an attractor of the spatial map (7), hence, independent of the value of the fixed boundary, the periodic pattern of the same wavelength (only with the phase difference) is obtained, as is observed in the direct simulation.(see Fig. 10).
As D Y is increased, the spatial period increases from 1 to 2, and 3 (see Fig.1c,d). With further increase in D Y , the periodic attractor in the spatial map disappears, and the fixed pattern is no longer generated in the original equation (1), either.
Also, when β in the model (A) is decreased, the pattern predicted by the spatial map can be quasiperiodic or chaotic, and again, a fixed spatial pattern is not generated, and is replaced by aperiodic spatiotemporal dynamics.
Finally, the spatial map works well in the model (B) also.

VI. DISCUSSION
In the present paper, we have studied the boundary-condition-induced pattern dynamics from a uniform oscillatory state in a reaction-diffusion system. By imposing a fixed boundary condition, the spatially-uniform limit cycle state is replaced by the following three phases, depending on the diffusion constants, : (I) fixed stripe pattern (III) travelling wave from the boundary (IV) aperiodic spatiotemporal pattern with the mixture of local travelling-wave and locally uniform oscillatory states, instead of the (II) original spatially-uniform limit cycle state. Note that in all cases, uniform oscillatory states are stable as long as the adopted boundary condition allows for such state globally, as given by, say, a periodic or Neumann boundary condition. In contrast, just by fixing a value at the boundary, the uniform state is replaced by the spatiotemporal patterns (I), (III) and (IV).
To analyze how the limit cycle is transformed to spatial pattern in the case (I), we introduced a diffusion term into the (local) nullcline analysis at each spatial point. Each The present mechanism of boundary-induced pattern formation is expected to work in a variety of systems that have a spatially uniform limit cycle attractor when the values of diffusion constants are appropriate. As such limit cycle attractor is common in a twovariable system with a negative-feedback interaction, say with an activator and an inhibitor, there will be a variety of systems that suit the present mechanism. Also, extension to a two or three-dimensional space is straightforward, where patterns are formed depending on the conditions on boundaries.
Further, the boundary-condition dependence, which is not fully explored as compared with the initial-condition dependence, will be important in a broad area of nonlinear systems.
By using the analysis by spatial map, the spatial pattern is obtained by the evolution of the mapping, and the initial condition in the mapping corresponds to the boundary condition.
Hence the standard analysis of dynamical systems is applicable. The spatial map will provide a useful tool to study the boundary-induced formation of fixed patterns.
Considering that the real system in nature is always finite in space, the boundarycondition can influence the system's behavior. The present result implies that the boundary can influence not only on its vicinity but globally throughout the space. Besides the stripe formation, spatiotemporal aperiodic behavior (chaos) is generated from a stable period oscillation just by fixing the boundary. It will be interesting to examine experimentally the observed behaviors in reaction-diffusion and other spatially extended systems.
Note that we first found the present fixation of oscillation in the evolution simulation to fit the protein expression pattern generated by gene-regulation-network with a requested spatial pattern 11 . Instead of the standard Turing pattern, we observed the present oscillationfixation mechanism frequently in these evolution simulations. Furthermore, in biological morphogenesis, oscillatory gene expressions are often observed, while the boundary condition is essential to morphogenesis.
In this context, it is interesting to recall that oscillation in protein expression is observed in somitogenesis 8,9 and in neural and other stem cells 6,28 . In the somitogenesis, such oscillation is transformed into spatial pattern, where relevance of cell-cell interaction has recently been discussed 14 . Also under in vitro condition, the differentiation from such oscillation is recently observed experimentally, where both the gradient due to cellular heterogeneity and cell-cell interaction seem to be important to induce the differentiation from the oscillatory dynamics 29,30 . In the present study, the fixed boundary generates a gradient, which causes the fixed pattern with differentiation of protein expressions X and Y through the diffusive interaction, as is consistent with the experiment. Further comparison with biological experiments on morophogenesis will be important both in biology and in nonlinear dynamics.
The present manuscript is dedicated to the memory of John Hudson. I (K.K.) first read his paper on chaos in BZ reaction 31 , about 36 years ago, and since then I have always been impressed by his marvelous studies unveiling beautiful structure from chemical-reaction experiments that could be potentially complicated. It was my fortune that our studies come closer, through his beautiful experiment on clustering 32 in electrochemical reaction which exhibits amazing correspondence with the behaviors in globally coupled maps 33 . I regret that I could no longer hear his invaluable opinion on the present study, based on his deep insight as experimentalist.