# Pulse Dynamics in a Bistable Reaction-Diffusion System with Chemotaxis

Authors: Satoshi Kawaguchi
Journal: Advances in Mathematical Physics
 (2022)

---

## Abstract

We consider pulse dynamics in a bistable reaction-diffusion system with chemotaxis. We derive the ordinary differential equation of interfaces by applying the multiple scales method to the reaction-diffusion system for examining the effect of the chemotaxis on pulse dynamics in one dimension. The stability of the standing pulse is considered by two different methods, and the applicability of the methods is demonstrated. The chemotaxis influences 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 differential 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.


---

## Body

## 1. 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 solution 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–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–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–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–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–39]. Recently, pulse-shaped bacterial colonies, which traveled along one direction, have been observed [32, 40–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.


## 2. 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:
(1)τε∂u∂t=ε2∂2u∂x2+fu,v,∂v∂t=D∂2v∂x2+gu,v,with
(2)fu,v=u+1212−uu−12v,gu,v=u−v+θ0,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 fu,v and gu,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±vi/τ2, where vi 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 l1t and l2t, respectively. The equations of motion of l1 and l2 are
(3)l̇1=r1,l̇2=r2,m0ṙ1=2τc−τr1−g3r13−G0+G1r2e−−r2+ϕr2h/2D+θ0,m0ṙ2=2τc−τr2−g3r23+G0−G1r1e−r1+ϕr1h/2D−θ0,where overdot l̇ denotes the time derivative of variable lt, ϕr=r2+4D, and pulse width h=l2−l1 (l2>l1). The parameters in Equation (3) are m0=3/16D, τc=1/42D, g3=1/32D3, G0=1/2, and G1=1/4D. System (3) is called the reduced ODE system. The inhibitor v causes repulsive interaction between pulses. The reduced ODE system (3) is valid for τ−τc≪1, 0 <θ0<1/2, r1, r2≪1, and h≫1.The linear stability of the stationary state,r1=r2=0 with a pulse width h0=−Dlog2θ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+2θ0G1+G0h0/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−2θ0G1+G0h0/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 τ−τc≪1, the bifurcation points obtained by the two methods were different for larger values of θ0.


## 3. 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 byu 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:
(4)τε∂u∂t=ε2∂2u∂x2−ε∂∂xu∂χv∂x+fu,v,∂v∂t=D∂2v∂x2+gu,v,where f and g are the same as Equation (2); χ corresponds to the chemotactic sensitivity function given by χv=kcf0v2a/v2+a2. For a=1, we choose f0 such that the condition maxdf0v2a/v2+a2/dv=1 is satisfied; f0=83/9. kc (≥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 fu,v and gu,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.


### 3.1. Reduced ODE System of One Pulse

In the presence of chemotaxis, the equation of motion of interface is given ([31, 36, 47]) by
(5)dl1dt=1τvl12+kc2f0a3vl1vl12+a22dvdxx=l1,dl2dt=1τ−vl22+kc2f0a3vl2vl22+a22dvdxx=l2.Thus, it is determined by the values ofv 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
(6)l̇1=r1,l̇2=r2,m0ṙ1=2τc−τr1−g3r13−G0+G1r2e−−r2+ϕr2h/2D+Δ0,m0ṙ2=2τc−τr2−g3r23+G0−G1r1e−r1+ϕr1h/2D−Δ0,where the parameters are
(7)m0=316D+2f0kca3232D,τc=142D+2f0kca1−2θ08D,g3=132D3+2f0kca21−θ032D2,G0=12+2f0kca1+2θ022D,G1=14D+2f0kca28D,Δ0=1+2f0kca12Dθ0.The pulse width in the stationary stateh0 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).
(8)h0=−DlogΔ0G0,τH=τc+12Δ0G0G1+G0h02D,τd=τc−12Δ0G0G1+G0h02D.


### 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 h0,v¯x is determined by the roots of the equation F±z=0. Here, the width h0=−Dlog2θ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
(9)F±z=−z+1τ12+f0kc1−2θ0Da×12D1−2θ0−12Dz+11±e−h0z+1D.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,
(10)F−0=dF−0dz=0,the standing pulse is destabilized through the translational instability; condition (10) achieves τd. Here, τd is given by
(11)τd=14D12+f0kc1−2θ0Da1−2θ01+h0D.WhenF+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
(12)F+ik=0.


### 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 asl1t and l2t, respectively. Similar to pulse 1, we denote the positions of the left and right interfaces of pulse 2 as l3t and l4t, respectively. We consider the nearest neighbour interaction between the interfaces. The brief derivation is given in Appendix C, and the equations of motion of l1, l2, l3, and l4 are
(13)l̇1=r1,l̇2=r2,l̇3=r3,l̇4=r4,m0ṙ1=2τc−τr1−g3r13−G0+G1r2e−−r2+ϕr2h/2D+Δ0,m0ṙ2=2τc−τr2−g3r23+G0−G1r1e−r1+ϕr1h/2D−G0+G1r3e−−r3+ϕr3d/2D−Δ0,m0ṙ3=2τc−τr3−g3r33−G0+G1r4e−−r4+ϕr4h~/2D+G0−G1r2e−r2+ϕr2d/2D+Δ0m0ṙ4=2τc−τr4−g3r43+G0−G1r3e−r3+ϕr3h~/2D−Δ0,where ht=l2t−l1t, h~t=l4t−l3t, and dt=l3t−l2t represent the pulse width of pulses 1 and 2 and the distance between the right interface of pulse 1 and the left interface of pulse 2, respectively. Parameters m0, τc, g3, G0, G1, and Δ0 are the same as Equation (7). The reduced ODE system (13) is valid for τ−τc≪1, 0 <θ0<1/2, ri (i=1, 2, 3, 4) ≪1, and h, h~, d≫1.


## 3.1. Reduced ODE System of One Pulse

In the presence of chemotaxis, the equation of motion of interface is given ([31, 36, 47]) by
(5)dl1dt=1τvl12+kc2f0a3vl1vl12+a22dvdxx=l1,dl2dt=1τ−vl22+kc2f0a3vl2vl22+a22dvdxx=l2.Thus, it is determined by the values ofv 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
(6)l̇1=r1,l̇2=r2,m0ṙ1=2τc−τr1−g3r13−G0+G1r2e−−r2+ϕr2h/2D+Δ0,m0ṙ2=2τc−τr2−g3r23+G0−G1r1e−r1+ϕr1h/2D−Δ0,where the parameters are
(7)m0=316D+2f0kca3232D,τc=142D+2f0kca1−2θ08D,g3=132D3+2f0kca21−θ032D2,G0=12+2f0kca1+2θ022D,G1=14D+2f0kca28D,Δ0=1+2f0kca12Dθ0.The pulse width in the stationary stateh0 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).
(8)h0=−DlogΔ0G0,τH=τc+12Δ0G0G1+G0h02D,τd=τc−12Δ0G0G1+G0h02D.


## 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 h0,v¯x is determined by the roots of the equation F±z=0. Here, the width h0=−Dlog2θ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
(9)F±z=−z+1τ12+f0kc1−2θ0Da×12D1−2θ0−12Dz+11±e−h0z+1D.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,
(10)F−0=dF−0dz=0,the standing pulse is destabilized through the translational instability; condition (10) achieves τd. Here, τd is given by
(11)τd=14D12+f0kc1−2θ0Da1−2θ01+h0D.WhenF+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
(12)F+ik=0.


## 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 asl1t and l2t, respectively. Similar to pulse 1, we denote the positions of the left and right interfaces of pulse 2 as l3t and l4t, respectively. We consider the nearest neighbour interaction between the interfaces. The brief derivation is given in Appendix C, and the equations of motion of l1, l2, l3, and l4 are
(13)l̇1=r1,l̇2=r2,l̇3=r3,l̇4=r4,m0ṙ1=2τc−τr1−g3r13−G0+G1r2e−−r2+ϕr2h/2D+Δ0,m0ṙ2=2τc−τr2−g3r23+G0−G1r1e−r1+ϕr1h/2D−G0+G1r3e−−r3+ϕr3d/2D−Δ0,m0ṙ3=2τc−τr3−g3r33−G0+G1r4e−−r4+ϕr4h~/2D+G0−G1r2e−r2+ϕr2d/2D+Δ0m0ṙ4=2τc−τr4−g3r43+G0−G1r3e−r3+ϕr3h~/2D−Δ0,where ht=l2t−l1t, h~t=l4t−l3t, and dt=l3t−l2t represent the pulse width of pulses 1 and 2 and the distance between the right interface of pulse 1 and the left interface of pulse 2, respectively. Parameters m0, τc, g3, G0, G1, and Δ0 are the same as Equation (7). The reduced ODE system (13) is valid for τ−τc≪1, 0 <θ0<1/2, ri (i=1, 2, 3, 4) ≪1, and h, h~, d≫1.


## 4. 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 ofu 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=∫l1l2dxux,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. 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).Figure 1
Profiles of standing and traveling pulses in the RD system (4). Solid and dashed curves represent u and v, respectively. System size L=100, kc=1, and θ0=1×10−2. (a) Standing pulse. τ=0.6. (b) Traveling pulse. τ=0.4.
(a)(b)Figure 2
Four types of interface motion in the RD system (4). System size L=100, kc=1, and θ0=1×10−2. (a) Standing pulse. τ=0.6. (b) Standing breather. τ=0.48. (c) Traveling breather. τ=0.4743. (d) Traveling pulse. τ=0.4. Periodic boundary condition is imposed.
(a)(b)(c)(d)Figure 3
Oscillation of the area of traveling breather. Data are obtained by the RD system (4). kc=1, θ0=1×10−2, and τ=0.4743.Figure 4
Four types of the interface motion of the reduced ODE system (6). θ0=1×10−2 and kc=1. (a) Standing pulse. τ=0.7. (b) Standing beather. τ=0.55. (c) Traveling breather. τ=0.487. (d) Traveling pulse. τ=0.35.
(a)(b)(c)(d)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 kc=0, τd and τH are symmetric with respect to τ~0.1823. By contrast, when kc=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 kc 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 τ−τc≪1. The dependence of h0, τc, and τd on kc are calculated using Equations (7) and (8), as shown in Figure 8. This suggests that h0 slightly depends on kc, and it is almost constant. By contrast, τc and τd increase linearly with kc. These properties hold for the small and large values of θ0.Figure 5
Bifurcation diagram of the standing pulse using the singular perturbation method. Solid and dashed curves represent Hopf and drift bifurcations, respectively. (a)kc=0. (b) kc=1.
(a)(b)Figure 6
Bifurcation diagram of the standing pulse by the linear stability analysis for the reduced ODE system. Solid and dashed curves represent Hopf and drift bifurcations, respectively. (a)kc=0. (b) kc=1.
(a)(b)Figure 7
Dependence of the Hopf bifurcation pointτH on kc in the case of θ0=1×10−3. Solid and dashed curves represent data obtained using the singular perturbation method and linear analysis of the reduced ODE system, respectively.Figure 8
Dependence ofh0, τc, and τd on kc. (a) Dependence of h0 on kc. Solid and dashed curves represent the cases θ0=1×10−3 and θ0=1×10−1, respectively. (b) Dependences of τc and τd on kc. Solid and dashed curves represent τc and τd in the case of θ0=1×10−3, respectively. Dash-dotted and dotted curves represent τc and τd in the case of θ0=1×10−1, respectively.
(a)(b)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.Figure 9
Collision of two traveling pulses obtained by the RD system (4). (a) Elastic collision and pair annihilation by collision. System size L=400, θ0=1×10−3, and kc=1. (a) (i) Elastic collision. τ=0.51575. (ii) Pair annihilation. τ=0.25157. (b) Parameter region of elastic collision in τ−θ0 space. Horizontally and vertically striped regions represent the parameter region, in which elastic collisions are observed in the case of kc=0 and 1, respectively. (c) Time evolution of h and d. τ=0.51575. Solid and dashed curves represent h and d, respectively. h~=h.
(a)(b)(c)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(b). The equation of motion of the interface 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.Figure 10
Collision of two traveling pulses obtained by the reduced ODE system (13). System size L=100, τ=0.45, θ0=1×10−3, and kc=1. (a) Collision of two traveling pulses. Solid, dashed, dash-dotted, and dotted curves represent l1, l2, l3, and l4, respectively. (b) Time evolution of h and d. Solid and dashed curves represent h and d, respectively. h~=h.
(a)(b)


## 5. 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 τ−τc≪1, 0<θ0≪1/2, ri≪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–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 τ−τc≪1, 0<θ0≪1/2, ri≪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 τ−τd≪1. There are perturbation methods which do not require any small parameters, for example, the homotopy perturbation and variational iteration method [54–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.



---

*Source: 1637071-2022-09-12.xml*