Reduction approach to the dynamics of interacting front solutions in a bistable reaction-diffusion system and its application to heterogeneous media

The dynamics of pulse solutions in a bistable reaction-diffusion system are studied analytically by reducing partial differential equations (PDEs) to finite-dimensional ordinary differential equations (ODEs). For the reduction, we apply the multiple-scales method to the mixed ODE-PDE system obtained by taking a singular limit of the PDEs. The reduced equations describe the interface motion of a pulse solution formed by two interacting front solutions. This motion is in qualitatively good agreement with that observed for the original PDE system. Furthermore, it is found that the reduction not only facilitates the analytical study of the pulse solution, especially the specification of the onset of local bifurcations, but also allows us to elucidate the global bifurcation structure behind the pulse behavior. As an application, the pulse dynamics in a heterogeneous bump-type medium are explored numerically and analytically. The reduced ODEs clarify the transition mechanisms between four pulse behaviors that occur at different parameter values.


Historical background to the present study
Reaction-diffusion systems have been widely utilized to model various spatiotemporal patterns observed in nature, such as chemical reactions, animal skin patterns, and morphogenesis. A class of spatially localized patterns is one of the most fundamental objects observed in many dissipative systems [58]. Nishiura and Mimura [42] studied a two-component bistable system in one-dimensional space of the form where 0 < ǫ ≪ 1 and τ, D > 0. This system, often referred to as a τ -ǫ system, includes the two key parameters τ and ǫ, which control the reaction creases [43]. These mathematical aspects of front solutions in bistable systems have been extensively studied [17][18] [19][26] [28][29] [35][39] [53]. Singularly perturbed three-component reaction-diffusion systems have also been investigated [4][8] [20][22] [24].
In this paper, we numerically and analytically study the pulse dynamics arising in the τ -ǫ system (1) with bistable-type nonlinearity. For the analysis, we do not deal with the original PDEs of (1), but with finite-dimensional ODEs that describe the pulse motion. Our study is inspired in part by the work of Ei et al. [11] [32]. In [11], they employed the following nonlinearities for f and g of (1): System (1) with (2) is of bistable type and has the three spatially uniform solutions P ± ≡ ( ±1/2, ±1/2 ), and P 0 ≡ ( 0, 0 ), where P ± are stable and P 0 is unstable. By the center manifold reduction, Ei where the variables l and r represent the location of the front interface and its velocity, respectively. The overdot ˙ denotes the time derivative d/dt, and interacting fronts in a bistable system in which the two front solutions are glued together to form a pulse shape. Note, however, that system (1) with (2) does not have a pulse solution, because the interaction between the front solutions is repulsive for the nonlinearities in (2). Hence, the interfaces monotonically repel each other because of the odd symmetry in f and g, and no pulse solution of finite width is formed. Such a repelling front motion has been extensively analyzed [3] [14][16] [31]. Thus, to construct a pulse solution with a finite width, Ei and Kusaka broke the odd symmetry by adding perturbation terms to Eq. (2): where 0 < δ 0 ≪ 1. The perturbation terms play the role of gluing the two fronts together. In fact, for particular choices of h 1 and h 2 , they numerically observed a stable pulse solution for system (1) with (4) The variables l 2 and l 1 represent the locations of the front interfaces, and r 2 and was destabilized by a pitchfork bifurcation at τ = τ p and a stable traveling pulse solution appeared instead. As τ decreased further, a Hopf bifurcation occurred as a secondary bifurcation at τ = τ H . As a result, a standing breather solution, which exhibited in-phase oscillation of the two interfaces, remained unstable after appearing via the Hopf bifurcation point. Regarding the order of the bifurcation points, the relation τ H < τ p was found analytically for ODE system (5). However, these analytic results disagreed with the pulse behavior observed by the direct numerical simulation of the original PDEs (1) with (4), which strongly suggested that the standing pulse solution underwent a Hopf bifurcation first and then a pitchfork bifurcation, implying that τ p < τ H . This was supported by a more elaborate numerical approach that clarified the global bifurcation structure of a pulse solution for a bistable reaction-diffusion system similar to ours [40]. Therefore, it seems most likely that the resulting ODEs (5) fail to properly reproduce the pulse dynamics for the original PDE system, and hence the underlying bifurcation structure.
One of the main goals of this paper is to derive equations for the interacting fronts that properly describe the pulse motion, thus resolving the discrepancy between the pulse dynamics for the original PDE system and those for the reduced ODE system. For the reduction, we start with the limit system in (7), which we introduced in our previous paper [41], and apply the multiple-scales method [17][18] [19] rather than the center manifold reduction. This approach yields the four-dimensional ODEs which are the same as in (16) in Section 2.3. The resulting equations are similar to (5) in appearance, but they successfully reproduce the pulse dynamics observed for the original reaction-diffusion system, including the aforementioned order of the Hopf and pitchfork bifurcation points. A key ingredient in remedying this discrepancy lies in the interaction terms +( G 0 − G 1 r 1 ) e − r 1 +φ(r 1 ) 2D h and  [59]. However, the mechanism for such heterogeneity-induced dynamics is still poorly understood considering the increasing number of experimental and numerical results [60]. In this paper, we deal with a bump-type spatial heterogeneity, for which a system parameter changes from one value to another in a certain interval of x. It is shown that our reduction method can be extended to the spatially heterogeneous case. First, we numerically obtain a complete phase diagram of the pulse dynamics, both for the original PDE system and the reduced ODE system. This diagram suggests four different kinds of behavior as the height and width of the bump are varied as bifurcation parameters. Based on the reduced ODEs, the transition mechanism for the pulse behavior is clarified from the perspective of a dynamical system. In particular, we characterize the difference between two kinds of behavior called pulse decomposition, whereby the pulse decomposes into two counter-propagating front solutions after encountering the heterogeneity. We find that the transition is caused by a change in the behavior of either pulse interface, which is unique to a pulse consisting of two interacting front solutions.
The remainder of this paper is organized as follows. In Section 2, we introduce the model system employed throughout the paper, namely, a system of two-component reaction-diffusion equations and its limit system of mixed ODE-PDE equations in the presence of a spatial heterogeneity. Next, we give a brief description of the derivation of the reduced four-dimensional ODEs for the motion of two interfaces of the pulse solution. The details of the derivation are given in Appendices A and B. In Section 3, we investigate the reaction-diffusion system and the reduced ODE system for the pulse dynamics in the bump-type heterogeneous medium, and study the transition mechanism for the pulse behavior from the perspective of a dynamical system. We conclude the paper with a summary and discussion in Section 4.

Equations of motion for two interfaces
2.1. Bistable reaction-diffusion system and its limit system We consider the two-component reaction-diffusion system in one-dimensional where u(x, t) and v(x, t) depend on space x ∈ R and t > 0. The function δ(x) represents a given spatial heterogeneity. For now, we assume that the system is homogeneous with δ(x) ≡ δ 0 , where δ 0 is some positive constant. In this case, we find that the system is bistable, with the two stable spatially uniform solutions (u ± (x, t), v ± (x, t)) ≡ (±1/2, ±1/2+δ 0 ) =: P ± and one unstable uniform solution P 0 . As shown in Figure 1(a), Eq. (6) not only has a front solution connecting the two spatially uniform solutions (u ± , v ± ), but also a pulse solution that exhibits four kinds of behavior as the parameter τ varies: (i) standing pulse (SP), (ii) standing breather (SB), (iii) traveling breather (TB), (iv) traveling pulse (TP) ( Fig. 1(b)).
In our previous paper [41], we took advantage of the fact that ǫ ≪ D and considered the limit system obtained as ǫ → 0, for which the equation for u is replaced by those for the two interfaces of the pulse solution located at x = l 2 , l 1 : where x = l 2 (t), l 1 (t) are the locations of the interfaces ( l 2 > l 1 ), the overdot denotes differentiation with respect to t, and the profile of the u component u(x; l 2 , l 1 ) is given by This type of mixed ODE-PDE system, which we call a hybrid system hereafter, typically arises in free boundary problems, such as the Stefan problem for phase transition [2], and has been used to study the dynamics of localized patterns in bistable reaction-diffusion systems [25][29] [50]. In our previous study [12][41], we examined the pulse dynamics for the hybrid system (7) both numerically and analytically. In particular, we followed the method described in [29] to investigate the eigenvalue problem for the stability of standing and traveling pulse solutions, confirming that the order of the Hopf and pitchfork bifurcations was consistent with that for the original PDEs, namely, τ p < τ H .
Numerical results indicated that the hybrid system (7) gave a good qualitative reproduction of the pulse dynamics observed in the original PDE system (6).
Comparing systems (6) and (7), we find that the equation for the u component in Eq. (6) is replaced by two ODEs for the interfaces, whereas the u term in the v component becomes piecewise constant. This allowed us to analytically show the existence and stability of the time-independent pulse solutions for both the heterogeneous and homogeneous cases, and hence clarify several mechanisms for the pulse dynamics observed in a jump-type heterogeneous medium [41]. However, the reduced hybrid system (7) still includes a PDE for the v component.
This has hindered the analytical study of more complicated behaviors, such as the standing breather and traveling breather shown in Figure 1(b-ii) and (b-iii).
We take a further step to reduce the hybrid system to finite-dimensional ODEs by perturbatively solving the equation for the v component. To this end, we employ the multiple-scales method, which was applied to a similar hybrid system by Hagberg et al. to study the front motion in a bistable reactiondiffusion system in both one-and two-dimensional space [17][18] [19]. In the present paper, we extend their approach to the dynamics of two interacting front solutions, and derive a system of ODEs that describes the interface motion observed in the original PDE system (6).

Reduction to ODE system
In this section, a reduction from the hybrid system (7) to ODEs is performed by perturbatively solving the PDE for the v component in Eq. (7). For the reduction, a multiple-scales method is utilized, which proceeds as follows.
We introduce a slow time scale T = µ t, where µ is some infinitesimally small constant 0 < µ ≪ 1, and assume that the variables l 2 , l 1 , and v in Eq. (7) have two timescales, t and T . However, we also assume that, after a sufficiently long time, these variables are described by the slow timescale T alone and become independent of the fast time scale t as t → ∞. Now, we expand the solution v in power series of µ as Substituting the above ansatz into the third equation in (7) and collecting terms with equal powers of µ, we can write equations for the functions v i (x, t) ( i = 0, 1, 2, · · · ) as where the order of the heterogeneity function δ(x) is assumed to be O(µ 0 ), which does not affect the final result of the reduction. In general, the following relation holds for n ≥ 1: Each equation in (10) is a linear ODE with respect to x and can be solved iteratively to give explicit solutions v n (x, T ). Substituting these v(x, t) into the first two equations in (7) yields closed-form ODEs for l 1 and l 2 ( see Appendix A for the details ). Neglecting the higher-order terms, we obtain the following four-dimensional ODEs for ( l 2 (t), l 1 (t), r 2 (t), r 1 (t) ): where h(t) := l 2 (t)−l 1 (t) and φ(r) := √ r 2 + 4D. The function ∆ 0 (x) originates from the heterogeneity term δ(x) in Eq. (6), which solves As l 2 and l 1 represent the locations of the two interfaces, the first two equations in (12) indicate that the variables r 2 and r 1 correspond to the velocities of the interfaces. The prefactor functions are defined as Note that Ohta et al. [49] and Mimura et al. [37][38] also studied a hybrid system similar to Eq. (7) with two ODEs for the interface motion and one PDE for the diffusion field, and reduced this to a four-dimensional system. The resulting ODEs were similar to Eq. (12) in appearance, but their reduction method was different from the aforementioned multiple-timescale technique. They first applied a Fourier transform to the PDE and solved for v(x, t) in an integral form.
Next, to evaluate the integral, they expanded the interface locations l 1 (t) and l 2 (t) in series of t up to the third order, thus obtaining closed-form second-order ODEs for l 1 (t) and l 2 (t). They also utilized AUTO [5] to study the bifurcation structure of the resulting second-order ODEs numerically, which will be mentioned in Section 2.3.
In this paper, we do not deal with the reduced ODE system (12), but make further simplifications by assuming the interfaces move very slowly and are sufficiently far apart.
Proof. Expand the functions m(r), g(r), and G(r) in (14) as power series of r: where the coefficients are given by m 0 = 3/( 16 . Substituting these series into (12) yields Truncating the higher-order terms including O( , we have the equations in (15). ✷ The four-dimensional ODEs (15) describe the motion of the two interfaces of the pulse solution observed in the PDEs (6), with the terms +( h serving as exponentially weak interactions between the interfaces. In the limit h → ∞, the interaction terms vanish and the system reduces to equations for the separate motion of two front solutions, for which, in the absence of the heterogeneity terms, l = l 2 = l 1 and r = r 2 = r 1 are exactly the same as Eq. (3), which was derived by center manifold reduction around a drift bifurcation point of a front solution [11].
Furthermore, the interaction terms simplify to ± G 0 e −h/ √ D when both r 2 and r 1 are set to zero. In weak interaction theory [10], reduced ODEs usually have this simple form of interaction term. In fact, the simple interaction term also appears in the ODEs of (5), which are quite similar to those in (15) in appearance, except for the additional terms ±M 0 e −h/ √ D and ± δ 0 β 1 , and the form of the interaction terms. However, the form of the interaction term qualitatively changes the behavior of the solution to the reduced ODE system, and hence its bifurcation structure, as we shall see in the next section.

Analysis of ODE system
Section 2.1 presented numerical results for the homogeneous case in which ODEs (15) can be found by solving Eq. (13), yielding ∆ 0 (l) ≡ δ 0 . In this section, we investigate the resulting ODEs (15) for the homogeneous case ∆ 0 (l) ≡ δ 0 : Using the relation h = l 2 − l 1 , this reduces to the three-dimensional system Figure 2(a) shows a numerical simulation of the ODEs (17) for various values of the bifurcation parameter τ . Compared with the numerical results in Figure   1(b), the four kinds of pulse behavior (SP, SB, TB, and TP) are successfully reproduced, and are in qualitatively very good agreement with the pulse dynamics of the original PDE system. Now that the system has been reduced to ODEs, the continuation and bifurcation software AUTO [5] can be utilized to investigate the global bifurcation structure for the four kinds of solutions ( Fig. 2(b)). The SP solution is stable for large values of τ , but loses its stability via a Hopf bifurcation as τ decreases.
From the Hopf bifurcation point, a stable SB solution branch appears. As τ decreases further, the SP solution undergoes a pitchfork bifurcation, at which point an unstable TP solution appears. The TP solution recovers its stability via a Hopf bifurcation, and a stable TB solution appears instead. The TB solution branch turns around at a fold bifurcation. Note that both the SB and TB solutions appear via Hopf bifurcations and their maximal h values seem to diverge as the parameter τ approaches a particular value. However, the TB solution is only stable before the fold bifurcation, whereas the SB solution remains stable after the Hopf bifurcation, yielding a small coexistence regime of SB and TB solutions, as seen in Figure 2(a).
As mentioned at the end of Section 2.2, Mimura et al. [37] studied secondorder ODEs similar to Eq. (12) using AUTO. They also obtained a qualitatively identical bifurcation structure for the four solution (SP, SB, TP, and TB) by varying a parameter that corresponds to τ in our ODEs. They observed that the pitchfork bifurcation on the SP solution branch became subcritical as another parameter varied, producing a saddle-node point on the TP solution branch.
This was confirmed for Eq. (17) with a larger value of δ 0 . Nagayama et al. [40] numerically studied a two-component τ -ǫ system with a bistable nonlinearity to clarify the bifurcation structure of the pulse solution. They obtained a bifurcation structure that is qualitatively the same as in Figure 2(b), except that the SB solution branch, unlike in the case of the ODEs (17), did not remain stable after the Hopf bifurcation point, but was destabilized via a saddle-node bifurcation. Considering that the PDE system investigated by Nagayama et al.
was quite similar to Eq. (6), it seems plausible to assume that the SB solution also undergoes a saddle-node bifurcation in our PDE system (6), implying that the reduced ODEs (17) do not completely reproduce the pulse dynamics of the original PDE system. Remarkably, similar pulse dynamics have been found in an excitatory-inhibitory neural field model [15]. The author explored the behavior of traveling pulse solutions both numerically and analytically for a wide range of parameters, and revealed the bifurcation properties of stationary pulse solutions. The neural field model has a non-local interaction across a wide spatial region, in contrast to our system of diffusion-induced local interactions, and the similarities in the pulse behavior and associated bifurcation properties seem to indicate that the bistable reaction-diffusion system is a good approximation of the neural model system, at least for some parameter regimes.
The stationary solution ( h, r 2 , r 1 ) = ( h * , r * , r * ) to Eq. (17) corresponds to either the TP solution (r * = 0) or the SP solution (r * = 0). The values of τ at which the SP solution undergoes the Hopf and pitchfork bifurcations can also be found analytically.
is a stationary solution to Eq. (17). Furthermore, the stationary solution undergoes a pitchfork bifurcation at τ = τ d and a Hopf bifurcation at τ = τ H as τ varies. These points are explicitly given by Substituting is a stationary solution to (17). Linearizing Eq. (17) around this stationary solution, we obtain the Jacobi matrix The constants φ 0 and φ 1 are the coefficients in the series expansion of where φ(r) = √ r 2 + 4D, and are given by The cubic equation for λ in (22) has one real root λ R and a pair of imaginary roots λ I ,λ I for some appropriate values of P , Q, and R, which depend on the parameters in the ODEs (17). In particular, the eigenvalues λ R and λ I change depending on the parameter τ through Q. Suppose that λ R = 0 for τ = τ d , The existence of the TP solution near the pitchfork bifurcation point τ d can also be shown by the regular perturbation method. where The constants φ 0 , φ 1 , and φ 2 are the coefficients in the series expansion of The constantg 3 is defined as Proof. The stationary solution ( h, r 2 , r 1 ) = ( h * , r * , r * ) to (17) satisfies We now introduce a small parameter 0 < η ≪ 1, and assume that the solution has the form where τ 0 , h (i) , and r (i) ( i = 0, 1, 2, 3 ) are unknown constants to be determined. Note that rescaling η allows us to set the coefficient of η 2 in the first equation to 1. Substituting the ansatz (26) into the first equation in (25) and collecting powers of η, we obtain where the constants Q i ( i = 0, · · · , 3 ) are defined as with The first equation in (27) yields e −φ0h (0) = δ 0 /G 0 or, equivalently, Substituting (28) and (29) into the second equation for O(η 1 ) in (27) yields By substituting the ansatz (26) into the second equation in (25), we have a similar equation: which is also obtained by the change in sign Equations (31) and (32) for h (1) and r (1) lead to Similarly, from the equation for O(η 2 ) in (27) together with the sign inversion (33), we obtain Finally, from the equation for O(η 3 ) in (27), we obtain whereg 3 is a constant defined bỹ From the first equation in (36), the value of r (1) = 0 is determined as Now that we have obtained r (1) and h (i) ( i = 0, 1, 2 ), we construct the approximate solution in (26). The first equation in (26) The second and third equations in (26) yield where h (0) and h (2) are defined in Eqs. (30) and (35), respectively. ✷ For later use, let us define the SP and TP solutions more rigorously based on the above two propositions.
In particular, as the constants on the right-hand side of (19) are all positive, we find that τ H > τ d , which agrees with the numerical results shown in Figure 2(b).
As mentioned at the end of the previous section, an interaction term of the form also reduce to ± G 0 e − h √ D when r 1 = r 2 = 0. In this case, both G 1 and φ 1 vanish in (19), and we are led to the following corollary. Thus, the above analysis suggests that higher-order interaction terms must be taken into account for our particular case of the front-back pulse dynamics observed in PDE system (6).

Numerical simulation of PDE system
So far, we have assumed that the function δ(x) in Eq. (6) is constant with δ(x) ≡ δ 0 , and hence ∆(l) ≡ δ 0 in Eq. (15).  [23]. The geometric singular perturbation theory has been effectively used to clarify the existence and stability of stationary pulse solutions induced by heterogeneity, and in our previous paper [41], we considered a mixed ODE-PDE system (7) and applied a center manifold reduction to derive ODEs describing the pulse dynamics in the presence of a jump-type heterogeneity.
In the present paper, we deal with a bump-type heterogeneity of the form and consider the situation in which a traveling pulse coming from the left infinity encounters this bump-type heterogeneity ( Fig. 3(a)). In Eq. (40), the parameter ǫ 0 represents the bump height, d 0 is the bump width, and γ is the steepness of the bump interface. Through the numerical simulation of the PDEs in (6) with Eq. (40), we find four kinds of pulse behavior as d 0 and ǫ 0 vary as bifurcation parameters: penetration (PEN), decomposition 1 (DEC1), rebound (REB), and decomposition 2 (DEC2). These are illustrated in Figure 3(b). For PEN, the traveling pulse passes the bump region and travels on ( Fig. 3(b-i)). As the bump height is increased, the pulse exhibits either DEC1 or REB behavior, depending on the bump width d 0 . When d 0 is relatively small, DEC1 is observed, whereby the interfaces of the pulse separate after encountering the bump (Fig. 3(b-ii)).
When d 0 is sufficiently large, the DEC1 behavior is no longer observed. Instead, the REB regime appears, whereby the pulse traveling to the right encounters the bump and turns back to the left (Fig. 3(b-iii)). For negative values of ǫ 0 , DEC2 is observed, whereby, as for DEC1, the two interfaces move away from each other and the pulse width diverges (Fig. 3(b-iv)).
Note the difference between the two types of pulse decomposition behavior,

Mechanism for PEN-DEC and DEC-REB transitions
The reduced four-dimensional ODEs shown in Section 2.2 are and these can also be applied to analytically investigate the pulse dynamics of the PDEs, as we saw in the previous section. For the sake of analysis, we assume that the bump interfaces are sufficiently steep and take the limit γ → ∞, so that the function δ 0 (x) simplifies to For the piecewise constant function (42), (13) is solved analytically as follows.
Lemma 1. Suppose that the function δ(x) is given by the piecewise constant function in (42). Then, a bounded C 1 solution ∆ 0 (x) to (13) is given by Proof. As (13) with (42) is an ODE with constant coefficients, it is readily solved as where the unknown constants C i ( i=1, 2, 3, 4 ) are determined by the condition that both ∆ 0 (x) and ∆ ′ 0 (x) are continuous at x = ± d 0 /2. Thus, we have Eq. (43). ✷ Now, we define penetration, rebound, and decomposition in terms of the solution orbit of the ODE system. Definition 2. When the orbit of (41) with (43) starting from TP + converges to TP + as t → +∞, we refer to the dynamics as penetration.
Definition 3. When the orbit of (41) with (43) starting from TP + converges to TP − as t → +∞, we refer to the dynamics as rebound.
Definition 4. When the orbit of (41) with (43) starting from TP + converges to the solution consisting of TF + 2 and TF − 1 as t → +∞, we refer to the dynamics as decomposition.
Remark For the definitions of TF + 2 and TF − 1 corresponding to the right-and left-moving traveling front solutions, see Definitions 5 and 6, respectively. Note that, for the decomposition behavior, the distance between the two interfaces diverges as t → +∞, and hence h → ∞. As mentioned at the end of Section 2.2, the limit h → ∞ decouples the four-dimensional ODEs (41) into two ODE systems, each describing the separate motions of the left and right interfaces. ODEs provide a good qualitative reproduction of the pulse dynamics in the original PDE system (Fig. 3(c)). In the following, we utilize the ODEs to investigate the mechanism for the PEN-DEC1, DEC1-REB, and PEN-DEC2 transitions from the viewpoint of bifurcation theory. In particular, we characterize the difference between DEC1 and DEC2, which exhibit similar asymptotic behavior

as shown in Figures 3(c-ii) and 3(c-iv).
To this end, the parameters were set very close to the phase boundaries; d 0 was fixed and ǫ 0 was varied ( Fig. 4(a)). The results reveal that the pulse interfaces exhibit different kinds of transient behavior right after the pulse encounters the heterogeneity. For parameters close to the PEN-DEC1 boundary, the left interface bounces off the bump heterogeneity and transiently moves with almost a constant velocity to the left before eventually traveling either to the left or to the right, resulting in the DEC1 or PEN behavior, respectively ( Fig. 4(a-i)). For those close to the DEC1-REB boundary, however, the right interface transiently moves to the right ( Fig. 4(a-ii)). In contrast, for the PEN-DEC2 transition, the left interface stays near the right edge of the bump heterogeneity for a certain period of time ( Fig. 4(a-iii)).
These numerical findings lead us to infer that one of the interfaces transiently approaches some unstable traveling front solution for the PEN-DEC1 and DEC1-REB transitions, whereas it approaches some unstable stationary front solution for the PEN-DEC2 transition. This can be confirmed analytically using the ODEs in (41). First, for the PEN-DEC1 transition, the spatiotemporal plot in Figure 4(a-i) suggests that the left interface moves away from both the right interface and the heterogeneity immediately after it bounces off, and that the influence of these two is negligible. Therefore, we consider the following ODEs for the motion of the left interface by setting h → ∞ and ∆ 0 (l 1 ) = δ 0 in Eq. (41): The solution with r 1 = const in (44) corresponds to the traveling front solution in PDE system (6).
Lemma 2. Assume that 0 < δ 0 ≪ 1, g 3 > 0, and τ < τ c . Then, the cubic equation for r 1 , has three roots, r 1 = r ± 1 and r 0 1 . These are approximately given by Proof. By assuming a solution of the form r 1 = a 0 + a 1 δ 0 + a 2 δ 2 0 + · · · and substituting this into (45), we have lecting the a i terms up to O(δ 2 0 ), we have the three solutions in (46). ✷ The solutions r 1 = r + 1 , r − 1 , and r 0 1 approximate the velocities of TF + 1 , TF − 1 , and TF 0 1 , respectively. Similarly, the solution that corresponds to the traveling front for the right interface, whose motion is governed by is defined as follows.
Definition 6. The solution ( l 2 , r 2 ) to Eq. (47) withṙ 2 = 0 is defined as The approximate value of r 2 in Definition 6 is obtained by inverting the sign of δ 0 to −δ 0 in Eq. (44).
Corollary 2. Assume that 0 < δ 0 ≪ 1, g 3 > 0, and τ < τ c . Then, the cubic equation for r 2 , has three roots, r 2 = r ± 2 and r 0 2 . These are approximately given by  For the PEN-DEC2 transition for the bump-down case ǫ 0 < 0, Figure 4(aiii) shows that the left interface transiently stays close to the right edge of the bump, meaning that the influence from the heterogeneity may not be neglected in this case. Therefore, we consider the following ODEs for the motion of the left interface: For the explicit form of ∆ 0 (x) in (43), the stationary solution to Eq. (50) is solved as follows. (50) is given by r * 1 = 0 and Furthermore, the stationary solution is unstable for τ < τ c .
Proof. Settingl 1 =ṙ 1 = 0 and ( l 1 , r 1 ) = ( l * 1 , r * 1 ) yield r * 1 = 0 and ∆ 0 ( l * 1 ) = 0. As ∆ 0 ( x ) in (43) is symmetric about x = 0, we assume l * 1 ≥ 0. Then, ∆ 0 ( l * 1 ) = 0 ( l * 1 ≥ 0 ) is solved as where z * 1 := e holds in the first equation in (52). Furthermore, the interval z 0 ≤ z * 1 ≤ 1 with in the first equation can be rewritten in terms of ǫ 0 as ǫ 0 ≥ −2δ 0 / ( 1 − z 2 0 ). Similarly, the interval z * 1 ≤ z 0 in the second equation is rewritten as ǫ 0 ≤ −2δ 0 / ( 1 − z 2 0 ). Finally, by the relation l * 1 = − √ D log z * 1 , we obtain the equations for l * 1 . The stability of the stationary solution is determined by the eigenvalues of the Jacobi matrix As tr J f = √ 2 ( τ c − τ )/m 0 > 0 and det J f = −p 0 > 0, we find that the eigenvalues are both positive for τ < τ c , and hence the stationary solution is unstable. ✷ Note that l 1 = l * 1 corresponds to the location of the stationary front solution. This is unstable for τ < τ c , the parameter range we now consider.        (5) and (16) In this study, the dynamics of pulse solutions have been examined both numerically and analytically for a bistable reaction-diffusion system (6). Applying the multiscale method to the hybrid system (7), we formally derived four-dimensional ODEs. These were found to successfully reproduce the pulse dynamics observed in the original PDE system, both for the homogeneous and heterogeneous cases. In particular, they correctly preserved the order of the pitchfork and Hopf bifurcations of the SP solution, which the previously derived ODEs had failed to do within the framework of weak interaction theory  ODEs. In the present study, we made a detour via the limit system (7) to the four-dimensional ODE system (16).

Comparison between four dimensional ODEs in
The pulse dynamics considered in this paper fall into the category of a semistrong front-front interaction [6][7] [9]. Away from the interfaces, the pulse solution comes very close to the equilibrium value for the u component, but not for the v component, which is largely deformed in between the two interfaces, as shown in Figure 1(a) for a typical profile of the pulse solution. Hence, the front interaction is essentially determined by the slowly varying component of v. In contrast, Ei and Kusaka [32] constructed the pulse solution by combining the front solutions in the framework of weak interaction theory. In this case, the fronts were assumed to be so far apart that they could only interact through their exponentially decaying tails, both for the u and v components. Although the two ODEs (5) and (16)  center manifold reduction may work well if the higher-order interaction terms are appropriately taken into account, which is left as future work. One may also be inclined to apply a renormalization group method for the reduction of pulse behavior, as in [24], in which the authors not only derived ODEs for the motion of N -fronts by the geometric singular perturbation technique, but also ensured their validity based on the idea of the renormalization group method.
As the resulting ODEs only describe monotonic attracting/repelling motion of the fronts, it would be a challenge to extend their rigorous approach to admit time-periodic and more complicated motions, as shown in Figures 2 and 8.

Sliding motion of standing breather in heterogeneous media
As an application of the reduced ODE system (15), we also considered the propagation manner of the traveling pulse that encountered the bump-type heterogeneity, and revealed the influence of parameter variations on the pulse behavior. In particular, we focused on the behavioral changes that occurred near the phase boundaries to identify the unstable solutions, called scattors, which played a central role in determining the pulse behavior. The related unstable traveling and stationary solutions were explicitly given based on the reduced system.
In this study, the bump height ǫ 0 and the bump width d 0 were varied as bifurcation parameters. Note that the bump width becomes infinitely large as d 0 → ∞. In this limit, the bump-type heterogeneity can be regarded as a jump, which was studied in our previous paper [41]. As in the bump case, the PEN and REB behavior were also observed for the jump-type heterogeneity, for which the traveling pulse coming from the left eventually converges to another stable traveling pulse that exists for the homogeneous system either in the right infinity x → ∞ (PEN) or the left infinity x → −∞ (REB). However, the DEC1 regime appeared between the PEN and REB regimes (Figs. 3(b),(c)); this was not observed for the jump-type case. The DEC1 regime occupies a relatively large parameter space when the bump width d 0 is comparable to the pulse size, and becomes narrower as d 0 increases. Eventually, the DEC1 regime seems to disappear for large values of d 0 , which indicates that the PEN behavior directly changes to REB behavior as ǫ 0 increases. We numerically confirmed this for the reduced ODE system (Figs. 6(a)(b)). After the pulse entered the bump region, it started to oscillate, and its amplitude gradually grew to that of a stable oscillatory solution for the homogeneous system δ(x) = δ 0 + ǫ 0 .
When the bump width was not sufficiently large, the pulse could not sustain the oscillation inside the bump region. After crossing the bump edges, the pulse interfaces decomposed into two counter-propagating front solutions. For sufficiently large values of d 0 , however, the oscillatory pulse did not immediately decompose into front solutions, but rather exhibited slow sliding motion that lasted exponentially long as the bump height ǫ 0 remained close to the PEN-REB boundary at ǫ PEN-REB 0 ( Fig. 6(c)). The mechanism for this sliding motion was analyzed for the jump case in [41]. As a result of the sliding motion, the oscillatory pulse slowly approached one of the bump edges, and its interface eventually crossed the bump edge, leading to either the PEN or REB behavior.
These numerical observations strongly suggest that there exists an unstable oscillatory pulse solution located at the center of the bump, playing the role of a scattor for the PEN and REB behavior. As we showed in Section 3, the reduced ODEs (15) can also be used to study unstable oscillatory pulse solutions. As these oscillatory pulse solutions often appear following the Hopf bifurcation of a stationary pulse solution, we first seek the corresponding stationary solution seemed to diverge at τ ∞ , indicating that the oscillatory pulse solution exists for τ ∞ < τ < τ H . By varying ǫ 0 as well, we obtained the τ -ǫ 0 diagram for the existence of the oscillatory pulse solution when d 0 = 40 ( Fig. 7(b)). In This implies that the unstable oscillatory pulse solution described above is not responsible for the PEN-REB transition. It remains to be elucidated whether there is another unstable oscillatory solution that plays the role of the scattor for the PEN-REB transition.

Designing trapped motion by square-well-type heterogeneity
Finally, we remark that the reduction method presented here is not limited to the dynamics of the two interacting fronts along the whole line, but can readily be extended to the pulse dynamics in finite domains [54] [55], as well as to dynamics that involve more than two interfaces [48] [51]. Another application of our study is the design of spatial heterogeneity based on the reduced system.
In the present paper, we have dealt with the bump-type heterogeneity, partly because this is one of the most fundamental and ubiquitous heterogeneities seen in nature. In principle, however, the function δ(x) in Eq. (6) can take any type of spatial heterogeneity, and may be designed to give the desired dynamics of moving localized patterns.
For instance, the pinning and trapping of moving localized patterns, for which pulses or spots become trapped in a finite region of space, are one of the most dramatic effects of heterogeneity. They play a central role in controlling the motion of localized patterns, together with the reversal and change of propagation direction [34] [36]. Such trapped motion of traveling pulses and fronts has been observed for reaction-diffusion systems with jump-or bump-type heterogeneities [23][21] [27] [62], where, unlike in our case, the heterogeneity was introduced to the activator (u component) rather than the inhibitor (v component). In fact, trapped motion was not observed in our system, where the bump-type heterogeneity was introduced to the v component.
Regardless, we can realize trapped motion by manipulating the form of δ(x) in Eq. (6), based on the information obtained for the bump-type heterogeneity.
To this end, we note that the traveling pulse exhibited the PEN and REB behavior depending on the bump height (Fig. 3), which readily leads to the idea that the traveling pulse may by trapped in between the two bumps by appropriately adjusting the bump height. Thus, we applied the square-welltype heterogeneity for δ(x) in (13): and observed trapped oscillatory motions of the traveling pulse for the ODE system ( Fig. 8) by changing the values of ǫ 1 and ǫ 2 . These motions were also found for the original PDE system. Note that the swaying motion in Figure 8(biv) does not occur when ǫ 2 = 0, where the two identical bumps are simply placed side by side.
Thus, the pulse dynamics for the original PDE system may be explored through the reduced ODEs, which facilitates numerical and analytical study of more complex behavior for other types of heterogeneity, as well as the design of external perturbations for controlling the traveling pulse.  Appendix A. Derivation of reduced ODE system (12) We introduce a slow time T = µ t and consider a perturbation series for v in where 0 < µ ≪ 1 is some infinitesimally small parameter. Then, the total time derivative can be rewritten as ∂/∂t → ∂/∂ t + µ ∂/∂ T . Substituting this into the third equation in (7), we obtain Collecting terms with equal powers of µ yields equations for v i ( i = 0, 1, 2, · · · ): Here, we focus on the asymptotic pulse behavior governed by the slow timescale T , and assume that v, l 2 , and l 1 are independent of t as t → ∞. Hence, where u(x ; l 2 , l 1 ) = F (x − l 1 ) − F (x − l 2 ) − 1/2, as shown in (8). This is readily where ∆ 0 (x) and v 0 (x) satisfy which is solved by where v 1 satisfies Similarly, the second-order equation in Eq. (A.3) can be rewritten as where v 20 and v 11 satisfy Thus, calculating up to O( µ 4 ), we obtain Here, to avoid redundancy, the notation [ terms of l 1 ] is introduced to denote the preceding v -terms with l 2 replaced by l 1 . Note that each equation for v i is a nonhomogeneous linear ODE with constant coefficients. In particular, the nonhomogeneous term F (x) for O(µ 0 ) is piecewise linear in our case, which allows us to derive explicit formulas for not only v 0 , but all v i ( i ≥ 1 ). Substituting the v i in (A.13) into the expansion (A.1) and rewriting the differentiation with (A.14) Among the many terms in Eq. (A.14), we choose the following and truncate the others as higher-order terms: v(x, T ) where we have defined Note that both U (1) i (x, T ) and U (2) i (x, T ) are expanded in power series of µ l i , which can be written as follows in renormalized form. First, applying the operator L to U (1) i and using the relations between v given by the lemma in Appendix B, we obtain (in the following, the argument x for v i is sometimes omitted, unless this may cause confusion): 17) or, equivalently,L This is again a nonhomogeneous linear differential equation with constant coef-ficients, which can be solved explicitly as where we have defined φ(l i ) := (µl i ) 2 + 4D.