Pulse Dynamics in a Bistable Reaction-Diffusion System with Chemotaxis

We consider pulse dynamics in a bistable reaction-di ﬀ usion system with chemotaxis. We derive the ordinary di ﬀ erential equation of interfaces by applying the multiple scales method to the reaction-di ﬀ usion system for examining the e ﬀ ect of the chemotaxis on pulse dynamics in one dimension. The stability of the standing pulse is considered by two di ﬀ erent methods, and the applicability of the methods is demonstrated. The chemotaxis in ﬂ uences the Hopf and drift bifurcations and the collision of two traveling pulses. It also enlarges the bifurcation point and enhances the repulsive force between pulses so that the parameter region of the elastic collision becomes large. Although the ordinary di ﬀ erential equation of interfaces can describe the elastic collision, it cannot describe the pair annihilation of pulses caused by the collision. The conditions for the reliable calculation of pulse collision are discussed.


Introduction
Reaction-diffusion (RD) systems are employed in wide fields. The normal diffusion equation, which underlies the Brownian motion of particles, is derived by using the equation of continuity (law of conservation of mass) and Fick's law. Although the mean squared displacement (MSD) of particles is linear in time for the normal diffusion, there are cases where the MSD is described by a power law [1]. For these cases, the anomalous diffusion is modeled using fractional or fractal derivatives, where the usual second derivative in space is replaced by a fractional derivative of order 0 < α < 2 [2,3]. The fractional derivative is employed in many fields such as control engineering, condensed matter physics, and electronics. For example, the vibration of viscoelastic bodies, thermal conduction in amorphous semiconductors, and discharge of Lithium-ion batteries are studied. On the other hand, the fractal derivative is introduced on the basis of the fractal concept to model power law scaling phenomena such as the turbulence, seepage in porous media, and groundwater flow in aquifer [4,5]. Although it is almost impossible to obtain an analytic solu-tion of anomalous diffusion systems, methods to solve them have been proposed such as the variational approach [6], homotopy perturbation method [7], and reproducing kernel method [8]. By contrast, with the use of normal diffusion, many RD systems have been proposed to describe the spatiotemporal patterns in physics, chemistry, biology, and neuroscience [9][10][11][12]. A homogeneous state is destabilized through the Turing instability, resulting in inhomogeneous states such as front and pulse [13].
The pulse dynamics and motion of interfaces have been considered in many systems [14][15][16][17]. For example, center manifold reduction with a weak interaction theory is employed in RD systems when front or pulse solutions of the system interact weakly [18][19][20][21]. Additionally, Kusaka derived the equation of motion of interfaces, that is, the ordinary differential equation (ODE) system by applying the method to a bistable RD system [22]. However, significant differences were observed when comparing the numerical results obtained using the RD system with that obtained using the ODE system. In the RD system, the standing pulse was destabilized through the Hopf bifurcation by changing a bifurcation parameter, and this was followed by drift bifurcation. By contrast, in the ODE system, the order of the bifurcation was reversed. Nishi et al. applied the multiple scales method to the RD system to overcome this discrepancy [23]. The multiple scales method has been employed to examine the bifurcatoin of a stationary solution in many systems [24][25][26]. By applying the multiple scales method to the bistable RD system, an extra term was added in addition to the equations obtained by Kusaka in ref. [22]; this term corrected the above discrepancy. However, the stability analysis of the ODE system did not agree well with the numerical simulation results. Therefore, the ODE system was reliable only within a limited parameter region.
Chemotaxis is an important mechanism for moving and transporting biological substances in biology [27,28]. Living organisms efficiently utilize chemotaxis to survive in an external environment [29,30]. The process of chemotaxis involves the detection of the direction of the gradient of the chemical substances in two or three dimensions in a biological organ; after the chemotactic substance is detected, the cell moves to or away from the source depending on whether the substances are chemoattractant or chemorepellent, respectively [31,32]. Chemotaxis is proportional to the gradient of the chemical substance in the medium. The mathematical model of chemotaxis was proposed by Keller and Segel, wherein the flow of amoeba was modeled by the gradient of the chemoattractant [33,34]. The collective movement of bacteria caused by chemotaxis has been observed and modeled based on the Keller-Segel model [35][36][37][38][39]. Recently, pulse-shaped bacterial colonies, which traveled along one direction, have been observed [32,[40][41][42]. Additionally, two types of collision of the pulse-shaped bacterial colonies were considered in a model under experimental conditions [43]. Although the mathematical models of traveling waves were considered, the stability and bifurcation of the solution were not sufficiently explained [44,45]. Even if the stability and bifurcation of the spatially uniform state were considered in ref. [46], the traveling pulse did not exist in their model.
To clarify the mechanism of emergence of bacterial traveling pulse, it is necessary to examine the stability of the stationary pulse in a mathematically tractable model. Thus, in this study, we consider the pulse dynamics in a bistable RD system with chemotaxis. Following ref. [23], we employ the multiple scales method for the model and derive an ODE system. First, we consider the stability of the standing pulse to examine the effect of chemotaxis on the system. Then, we consider the collision of two traveling pulses and compare the mechanisms of the repulsion of pulses observed in other models. Finally, reliable conditions of traveling pulse and pulse collision in the ODE system are discussed.

RD System and Reduced ODE System without Chemotaxis
In this section, we introduce the preceding study reported in ref. [23] and briefly explain the properties of the bistable RD and ODE systems. A two-component RD system without chemotaxis in one dimension was proposed [21,22]. Activator u and inhibitor v satisfy the time evolution equation: where we choose ε and θ 0 such that 0 < ε ≪ 1 and 0 ≤ θ 0 ≤ 1/2. Here, τ and D represent the positive constants. Let us briefly summarize the solution of system (1). When θ 0 = 0, system (1) has stationary front solutions connecting u = ±1/2 for a large τ. The stationary front solution bifurcates to the traveling front solutions when τ is decreased to less than the critical value τ c . However, system (1) does not have a pulse solution because of the odd symmetry of the reaction terms f ðu, vÞ and gðu, vÞ. For θ 0 > 0, the odd symmetry is broken, and system (1) has pulse solutions. The pulse solution is called the front-back pulse because the system is bistable; the pulse is caught by two interfaces that connect two stable points. There are four types of solutions depending on the value of τ: standing pulse, standing breather, traveling breather, and traveling pulse.
The velocity of the interface is obtained as ±v i /ðτ ffiffi ffi 2 p Þ, where v i corresponds to the value of v at the interface [47]. The equations of motion of the left and right interfaces were derived by applying the multiple scales method to Equation (1) [23]. We consider one pulse in system (1) and denote the positions of the left and right interfaces of the pulse as l 1 ðtÞ and l 2 ðtÞ, respectively. The equations of motion of l 1 and l 2 are where overdot _ l denotes the time derivative of variable lðtÞ, ϕðrÞ= ffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi r 2 + 4D p , and pulse width h=l 2 − l 1 (l 2 >l 1 ). The parameters in Equation (3)  Advances in Mathematical Physics reduced ODE system (3) is valid for jτ − τ c j≪1, 0 <θ 0 <1/2, jr 1 j, jr 2 j≪1, and h≫1: The linear stability of the stationary state, r 1 =r 2 = 0 with a pulse width h 0 =− ffiffiffi ffi D p log ð2θ 0 Þ, was studied [23]. Although the standing pulse was stable for a large τ, it was destabilized through the Hopf bifurcation as τ decreased at τ H =τ c + ffiffi ffi 2 p θ 0 ðG 1 + G 0 h 0 /ð2DÞÞ, and a standing breather appeared. The standing breather was destabilized through pitchfork bifurcation with a further decrease in τ; then, the standing breather started to move while retaining interface oscillation (i.e., a traveling breather). With a further decrease in τ, the interface oscillation ceased, and the traveling breather changed to a traveling pulse at τ d =τ c − ffiffi ffi 2 p θ 0 ðG 1 + G 0 h 0 /ð2DÞÞ. There exists a small parameter region of coexistence between the standing and traveling breathers. Although the bifurcation points are obtained by the RD, and reduced ODE systems agreed well for small values of θ 0 with jτ − τ c j≪1, the bifurcation points obtained by the two methods were different for larger values of θ 0 .

RD System and Reduced ODE System with Chemotaxis
We now consider a two-component RD system with chemotaxis in one dimension. We denote the population density of biological individuals and concentration of chemotactic substance by u and v, respectively. We assume the case wherein the chemotactic substance is released from the biological individual, and it diffuses outward and influences the movement of each individual. We consider the following time evolution equations of u and v: where f and g are the same as Equation (2); χ corresponds to the chemotactic sensitivity function given by χðvÞ= . k c (≥0) represents the intensity of the chemotaxis [48]. We choose ε and θ 0 such that 0<ε≪1 and 0≤θ 0 ≤1/2. Here, τ and D represent positive constants. We note three important points about system (4). First, because we use f ðu, vÞ and gðu, vÞ as presented in Equation (2), the ranges of u and v are −1/2≤u, v≤1/2. From the practical perspective, the population density of biological individuals u and concentration of chemotactic substance v should satisfy u, v≥0. To overcome this problem, although we adopt system (4) in the calculation, the values u and v shall be replaced with ðu + 1/2Þ and ðv + 1/2Þ, respectively, in practical situations (experiments). Second, v represents not only the concentration of chemotactic substance but also the concentration of inhibitor in system (4). This dual role of v is a new aspect of the system. Third, we note the difference between our model (4) and other chemotaxis models [40,41,43]. In these models, three variables are introduced: the density of bacteria, concentration of cheomotactic substance, and concentration of nutrition. In our model, the variable of the concentration of nutrition is not considered; it corresponds to the scenario that the nutrition is sufficiently supplied, and therefore, we do not need to consider its distribution. In this study, we consider the effect of the chemotaxis on the pulse solution. The dual role of v as the chemotactic substance and inhibitor changes the stability of the standing pulse, and it reduces the velocity of pulses during the collision. We consider the stabilty of the standing pulse solution by two different methods: the singular perturbation method of the interface and the linear stability of the reduced ODE system. Further, we consider the collision of the two traveling pulses. Then, we compare the numerical results obtained by the direct calculation of the RD and reduced ODE systems to obtain reliable conditions for the reduced ODE system.

Reduced ODE System of One Pulse.
In the presence of chemotaxis, the equation of motion of interface is given ( [31,36,47]) by Thus, it is determined by the values of v and ðdv/dxÞ at the interface. Following the process reported in ref. [23], we derive the reduced ODE system of the RD system with chemotaxis. The brief derivation is given in Appendix A. The equations of motion of interface have the same form as in Equation (3); however, the parameters are different. The final expressions are where the parameters are ,   Advances in Mathematical Physics

Advances in Mathematical Physics
, The pulse width in the stationary state h 0 and the Hopf and drift bifurcation points, τ H and τ d , respectively, are obtained in a similar manner as that without chemotaxis by using the parameters given in Equation (7).
3.2. Stability Analysis Using a Singular Perturbation Method.
The stability of the standing pulse is examined by the singular perturbation method for the RD system (4) [49]. The exponential growth of the perturbation to the standing pulse solution ðh 0 , vðxÞÞ is determined by the roots of the equation F ± ðzÞ= 0. Here, the width h 0 =− ffiffiffi ffi D p log ð2θ 0 Þ and vðxÞ represents the stationary solution of the RD system (4) in the limit ε↓0, which is symmetric with respect to x= 0. The brief derivation of the stability formula F ± ðzÞ is given in Appendix B; the final expression is Here, F + and F − correspond to the cases of symmetric and antisymmetric perturbation with respect to the center of the pulse, respectively.
F − ð0Þ= 0 corresponds to the translational invariance. When F − ðzÞ is degenerated at z= 0, that is, the standing pulse is destabilized through the translational instability; condition (10) achieves τ d . Here, τ d is given by When F + ðzÞ= 0 has a pair of imaginary solutions ±ik with real number k, the standing pulse is destabilized through the Hopf bifurcation; the interfaces oscillate symmetrically with respect to the center of the pulse. This condition achieves τ H ; here, τ H is given by the solution 3.3. Reduced ODE System of Two Pulses. We derive the reduced ODE system of two pulses to consider the collision of two traveling pulses. We consider the case with two pulses: pulse 1 and pulse 2. Let us denote the positions of the left and right interfaces of pulse 1 as l 1 ðtÞ and l 2 ðtÞ, respectively. Similar to pulse 1, we denote the positions of the left and right interfaces of pulse 2 as l 3 ðtÞ and l 4 ðtÞ, respectively. We consider the nearest neighbour interaction between the interfaces. The brief derivation is given in Appendix C, and the equations of motion of l 1 , l 2 , l 3 , and l 4 are  Advances in Mathematical Physics

Numerical Results
In this section, we consider the numerical results. In the calculation of the RD system (4), we fix ε= 5 × 10 −2 and D= 1 and choose differences Δx= 2:5 × 10 −2 and Δt= 5 × 10 −4 . The spatial direction is calculated using the implicit method under Neumann boundary conditions unless otherwise noted. Additionally, in the calculation of the reduced ODE systems (6) and (13), we employ the Runge-Kutta scheme with Δt= 5 × 10 −4 . The profiles of u and v for the standing and traveling pulses in system (4) are shown in Figure 1. These profiles are almost the same with that in system (1). In the limit ε↓0, the profile of u becomes rectangular, connecting u= 1/2 with −1/2; the sharp edge is called the interface.
For the finite value of ε (≪1), the interface becomes a rather dull layer (inner layer), and its width is approximately~ε. By decreasing τ in system (4), four different types of solutions appear, as shown in Figure 2. Similar four types of solution are observed in system (1). Thus, solutions of the RD system Equations (4) and (1) have similar properties; the influence of chemotaxis on the solution is not evident. The area of the pulse is shown in Figure 3, where the area is defined by S= Ð l 2 l 1 dxðuðx, tÞ + 1/2Þ; this clearly indicates the oscillation of the interface of the traveling breather (τ= 0.4743). Four pulse solutions observed in the RD system (4) occur in the reduced ODE system (6), as shown in Figure 4.

Advances in Mathematical Physics
These pulse solutions are also observed in the reduced ODE system (3). In our ODE system (6), the coexistence of the standing and traveling breathers are observed in a small parameter region. Thus, the reduced ODE system (6) describes similar properties observed in the RD system (4).
We consider the dependence of the Hopf and drift bifurcations on θ 0 using the singular perturbation method. Comparing Figure 5(a) with Figure 5(b), it suggests that the chemotaxis enlarges the bifurcation points; τ d and τ H become large under the influence of chemotaxis. To confirm the validity of the reduced ODE system (6), we calculate the bifurcation diagram using the linear stability analysis, as shown in Figure 6. When k c = 0, τ d and τ H are symmetric with respect to τ~0:1823. By contrast, when k c = 1, τ d and τ H are not symmetric. In the bifurcation diagram obtained by the singular perturbation method ( Figure 5), τ d and τ H approach 0 with an increase in θ 0 . However, when using the linear stability analysis of the reduced ODE system, τ H does not approach 0 with an increase in θ 0 . Additionally, for larger θ 0 , although the Hopf bifurcation exists, the drift bifurcation does not occur. This is attributed to the fact that the amplitude of the oscillation grows larger and h becomes negative with decreasing τ; it becomes impossible to continue the calculation. Thus, τ H obtained by the reduced ODE system is unreliable for the larger values of θ 0 . For the small values of θ 0 , τ d~τH . In this case, the dependence of τ H on k c is calculated by two different methods: the linear stability analysis and singular perturbation method. τ H obtained by the different methods agree well with each other, as shown in Figure 7. Thus, the reduced ODE system is reliable for small θ 0 and jτ − τ c j≪1. The dependence of h 0 , τ c , and τ d on k c are calculated using Equations (7) and (8), as shown in Figure 8. This suggests that h 0 slightly depends on k c , and it is almost constant. By contrast, τ c and τ d increase linearly with k c . These properties hold for the small and large values of θ 0 .
Let us consider the collision of two traveling pulses using the RD system (4). When the traveling velocity is low, they collide elastically; the two traveling pulses approach slowly and stop before collision with finite acceleration. Then, the pulses travel backward. By contrast, when the traveling velocity is high, the two traveling pulses collide inelastically; the two traveling pulses do not stop before collision, resulting in pair annihilation (Figure 9(a), i-ii). The parameter region of the elastic collision is shown in Figure 9(b). The parameter region of elastic collision with chemotaxis (vertically striped region) is wider than that without chemotaxis (horizontally striped region). Thus, the chemotactic substance emitted inside the pulse reduces the speed during the collision, resulting in an elastic collision. The time evolution of h and d during collision is shown in Figure 9(c), where we define the position of the interface as the position satisfying u= 0. Here, h and d are apparently positive; when h is locally the smallest value (t~54), d becomes the locally maximum value, and the velocity of the two pulses becomes zero around this period.
The collision of two traveling pulses is considered using the reduced ODE system of the two pulses (13). The left and right interfaces of each pulse are traced, as shown in Figure 10(a). The interaction between the interfaces is considered among the nearest neighbours, as in Equation (13). The interface is repelled to the backward direction when an interface approaches the nearest neighbour interfaces. The pulse width h and the distance between two interfaces d can be negative under a weak repulsive force when the traveling velocity is high. In fact, when the braking effect on pulse speed is small, h and d are close to zero, as shown in Figure 10 9 Advances in Mathematical Physics includes a repulsive force resulted from the inhibitor (and chemotactic substance) v; therefore, the elastic collision of pulses is reasonable when the pulse velocity is low. By contrast, for a smaller τ, the pulse velocity is high, and h and d can have negative values by collision. It is impossible to calculate the time evolution for the negative values of h and d because of the divergence of exponential terms in Equation (13). Comparing Figure 10(b) with Figure 9(c), although the qualitative properties are similar, h and d become smaller in the reduced ODE system (13). For the RD system (4), with a smaller τ, the velocity of the pulse is high, and we can observe the pair annihilation of two pulses, as shown in Figure 9(a), ii. However, the reduced ODE system (13) cannot describe pair annihilation by collision.

Conclusions
We considered the effect of chemotaxis on pulse stability and dynamical properties. In the RD system (4), v played the dual role as a chemotactic substance and inhibitor, and influenced the motion of biological individuals. We derived the equation of motion of two interfaces, called the reduced ODE system, by applying the multiple scales method to this RD system. The equation included the linear, cubic nonlinear, and weak interacting terms and had the same form as that without chemotaxis derived in ref. [23]. The standing pulse was destabilized through the Hopf and drift bifurcations by decreasing the bifurcation parameter τ. These properties were the same as those without chemotaxis. The parameter region of the Hopf and drift bifurcations became wider under the presence of chemotaxis. To verify the validity of the reduced ODE system, we examined the stability of the standing pulse by two different methods: the singular perturbation method and linear stability of reduced ODE system. Definite differences were observed: when θ 0 was small, the reduced ODE system described the RD system well. However, for a large θ 0 , the Hopf bifurcation point increased, which did not agree with the results obtained using the singular perturbation method. Thus, the reduced ODE system (6) is reliable only for jτ − τ c j≪1, 0<θ 0 ≪1/2, jr i j≪1 (i= 1 and 2), and h≫1: We derived the reduced ODE system for two pulses to consider the collision of two traveling pulses. The nearest neighbour interaction was considered, and the term of repulsive force appeared in the equation that defined the motion of interfaces. Although this term decayed exponentially with the distance between the two interfaces, it caused a strong force that changed the direction of motion when the distance between the two interfaces was small. Chemotaxis enhanced the repulsive force, resulting in a wider parameter region of the elastic collision. The elastic collision of pulses with low speed was theoretically considered in an RD system without chemotaxis [19,[50][51][52][53]. The equation of motion of a pulse was derived for the pulse with low speed; this had a similar form as that of Equation (3). When the traveling velocity became zero with finite acceleration (corresponding to the third term in the left-hand side of Equation (3)), the pulse could travel in the backward direction. Although the traces of the interfaces of pulses were similar between the RD system and the reduced ODE system, the pair annihilation of pulses with a smaller τ could not be observed in the reduced ODE system. In the reduced ODE system, although the repulsive term existed, there was no restriction on the time evolution of the interfaces; this resulted in negative h, h, and d. In the RD system (4), the inhibitor (and chemotactic substance) v caused the pair annihilation of traveling pulses with high velocity. By contrast, in the reduced ODE system (13), h,h, and d could be negative, and the reduced ODE system could not describe the pair annihilation of traveling pulses with high velocity. More generally, the reduced ODE system (13) is valid under the conditions jτ − τ c j≪1, 0<θ 0 ≪1/2, jr i j≪1 (i= 1, 2, 3, 4), and h,h, d≫1: To overcome the limitation of the multiple scales method, other perturbation methods are necessary. In this study, we employed the multiple scales method to derive the reduced ODE system, and we expanded h and r using small parameters under the assumption jτ − τ d j≪1. There are perturbation methods which do not require any small parameters, for example, the homotopy perturbation and variational iteration method [54][55][56]. Recently, He's multiple scales method was proposed, where the homotopy perturbation method was incorporated into the multiple scales method [57,58]. Although this method is promising and may be applicable to our RD system, the discussion on detailed calculations is beyond the scope of the present study.
The collision of two traveling pulses was considered in [43]. In their model, two types of collisions were considered: an elastic collision and a fusion of two traveling pulses (inelastic collision) followed by a single traveling pulse. The concentrations of chemoattractant and nutrition were considered under the constraint of constant total bacteria. The relative strength of the attractive force caused by these two substances determined the elastic or inelastic collisions; when the attractive force caused by the chemoattractant was larger than that of the nutrition, elastic collision occurred. By contrast, when the attractive force caused by the nutrition was larger than that of chemoattractant, the fusion of two traveling pulses occurred under the constraint of constant total bacteria, followed by a single traveling pulse with the help of noise or asymmetries of the system. We note that, in their model, the pair annihilation of pulses did not appear due to the constraint of constant total bacteria. Thus, the elastic and inelastic collisions resulted from the distributions of chemoattractant and nutrition. Conversely, in our system, the chemotactic substance played the role of an inhibitor; the inhibitor in system (1) was replaced by the chemotactic substance in our system (4), and the chemotactic substance had a repulsive force on the other pulse. When the velocities of the traveling pulses were low, the repulsive force led to elastic collision. By contrast, when the velocities of the traveling pulses were high, one traveling pulse moved into the refractory zone of the other traveling pulse and resulted in the pair annihilation of the two pulses. Thus, in our system, the elastic collision and pair annihilation of traveling pulses resulted from the property of the inhibitor. In addition, the property of the chemotactic substance enlarged the parameter region of the elastic collision. A. Reduced ODE System of One Pulse In this appendix, we briefly derive the reduced ODE system in the presence of chemotaxis. The equations of motion of interface are ðA:1Þ Following ref. [23], we employ the multiple scales method to derive the equations of motion of interface. The motion of the interface is slow at the onset of the bifurcation from the standing pulse. To describe the slow motion of the interfaces, we use the scaled time T=μt, where μ is a small parameter in the range 0 <μ≪1. Thus, l i (i= 1 and 2) and v are l i =l i ðTÞ as v=vðx, TÞ. Then, vðx, TÞ is expanded with μ as where a prime represents the derivative of a function with respect to x, and U ð1Þ i and U i (i= 1 and 2) are U  where ϕðrÞ= ffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 4D + r 2 p . The first term of (A.1) is already estimated in ref. [23]. We show the final result of the second term, which is the chemotaxis term. Let us consider the equation of motion of l 2 . When _ l 2 is small, the denominator ðvðl 2 Þ 2 + a 2 Þ 2 can be approximated as a 4 ðA:6Þ Substituting the above equations into the second term of (A.1) and putting μ= 1, we can derive the equation of motion of l 2 . Similarly, we can derive the equation of motion of l 1 ; the final expressions of the equation of motion of the interfaces are given in (6) with the parameters given in (7).

B. Stability Formula Using a Singular Perturbation Method
In this appendix, following the singular perturbation method [49], we briefly derive the stability formula F ± given in Section 3. Let us consider the symmetric solution with respect to x= 0, which is denoted by ðh 0 , vðxÞÞ, where the