# Motion of a Spot in a Reaction Diffusion System under the Influence of Chemotaxis

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

---

## Abstract

We consider the motion of a spot under the influence of chemotaxis. We propose a two-component reaction diffusion system with a global coupling term and a Keller-Segel type chemotaxis term. For the system, we derive the equation of motion of the spot and the time evolution equation of the tensors. We show the existence of an upper limit for the velocity and a critical intensity for the chemotaxis, over which there is no circular motion. The chemotaxis suppresses the range of velocity for the circular motion. This braking effect on velocity originates from the refractory period behind the rear interface of the spot and the negative chemotactic velocity. The physical interpretation of the results and its plausibility are discussed.


---

## Body

## 1. Introduction

The behaviors of artificial and biological microswimmers such as oil droplets, bimetallic nanorods, catalytic Janus colloids, liposomes, flagellated bacteria, andVolvox have attracted widespread attention [1]. Under certain circumstances, some of these microswimmers are self-propelled particles, the mobility mechanism of which has been intensively studied [2, 3]. The motion of oil droplets, especially, has been studied in well-controlled experimental facilities with sufficient reproducibility. Although symmetric droplets cannot move in the absence of external force, the Marangoni effect can cause motion in the presence of an inhomogeneous chemical substance outside the droplet or a temperature gradient along the surface [4–6]. Numerical simulations and theoretical results support this mechanism and the existence of straight, circular, and complicated motions of droplets [7, 8], and experimental results qualitatively agree with the numerical results [9–12]. Droplet motion has also been the subject of a review article [13, 14].In a two-dimensional reaction diffusion (RD) system, the droplet is often referred to as a spot solution. In order to systematically describe the motion of spots in an RD system, the time evolution equation of the spot was derived and the mechanism of elastic collision of moving spots was clarified in a previous study [15]. This study was extended by studies on the drift and rotation bifurcations of spot solutions in RD systems [16, 17]. In order to describe the deformations of the spot, tensors were introduced. The bifurcation diagram of the spot suggested that, with increasing velocity of the spot, rotation bifurcation occurred causing the straight motion to become destabilized into circular motion.In addition to the Marangoni effect, which plays an important role in the motion of oil droplets, chemotaxis is an important property of cell migration; it is important in mass transfer and immunological response in biology. In inflammatory response, the neutrophils among blood cells have a remarkable migration potency (chemotaxis) and can change their form by generating pseudopods toward the antigen. In biophylaxis, several chemokines (chemoattractants) are released from the macrophages and mast cells. Then other immunocompetent cells (neutrophils) respond to the gradient of the chemoattractant. Consequently, the immunocompetent cells move unidirectionally to the source point of the antigen [18].The mathematical model for chemotaxis was first proposed by Keller and Segel [19], wherein the gradient of the chemoattractant was taken into consideration for the flow of amoeba. Neutrophil migration was considered with a Keller-Segel type chemotaxis term in [20, 21]. In these studies, the Cahn-Hilliard (CH) equation was employed, and the kinematic properties and morphological changes of the crawling cell distribution were shown. In addition to the chemotaxis of the neutrophil, cancer cell invasion under haptotaxis was modeled by the CH equation [22, 23]. The haptotactic response of cancer cells is described by the gradient of the haptoattractant. However, in the above studies, the gradients of chemoattractant and haptoattractant are assumed to be constant; there is no feedback between the cells and these chemical substances.As described above, with the recent increase in importance of chemotaxis in biology, medicine, and cytoengineering [24–27], many experimental and theoretical studies have been performed. Although there are model systems for the cell density and concentration of chemotactic substances, no mathematical analysis has been reported on the motion of the cell. Inspired by these points, we first propose an RD system including a naive Keller-Segel type chemotaxis term. The system is autonomous, the spot secretes a chemotactic substance, and the motion of the spot is influenced by it. For the proposed RD system, we apply the method reported in [16] to derive the equation of motion of the spot and time evolution equation of the tensors. Based on these equations, we study the bifurcation from straight motion to circular motion as well as the upper limit of the velocity of circular motion. In order to verify the theoretical result, we perform numerical simulations for the tensor model. The physical meaning and validity of the results are discussed.


## 2. Model Equation

We first consider the following three-component RD system with an activatoru, a chemotactic substance v, and an inhibitor w:(1)τϵ∂u∂t=ϵ2∇2u-ϵ∇·u∇χv+Hu-pw-u-v,(2)∂v∂t=∇2v+u-μv,(3)T∂w∂t=D∇2w+u+v-w-W,where p(w)=p0+αw,  p0, α, τ, ϵ, μ, T, D, and W are positive constants, and H(z) is a step function satisfying H=0 for z<0 and H=1 for z>0. Throughout this study, we consider the system in a two-dimensional space, with ∇=(∂x,∂y) and ϵ≪1. We choose μ such that the system is monostable. Here, we fix μ=0.3. In the above excitable system, there are two stationary states: a rest state and an excited state. The rest state is (u,v,w)=(0,0,-W), and the excited state has spatially nonuniform values of u, v, and w. Between the rest and excited states, there appear boundary layers with thickness O(ϵ), connecting the two different states.When the second term on the right hand side of (1) is absent, (1)–(3) describe an RD system with one activator and two inhibitors, which was studied in [28]. In that system, when τ is large, the localized domain (motionless spot solution) of an activator appears. With decreasing τ, the motionless spot is destabilized through static bifurcation or oscillatory bifurcation; however, when T is small and D and α are large, these bifurcations are suppressed by w. When τ is small, the motionless spot is primarily destabilized through translational bifurcation, causing the spot to move.In the presence of the second term on the right hand side of (1), the moving spot is influenced by the chemotaxis. A system similar to that described by (1)–(3), but with bistability, was studied in [29]. In that system, the nonlinear term in (1) was replaced by H(u-p(w))-u, and a front solution was obtained. Furthermore, maze patterns and branching from a front solution were observed. The stability analyses of the spot and front solutions were conducted by applying the singular perturbation method [30].The time evolution equation foru is obtained using the conservation equation. The diffusion term is derived from -∇·J, where the flux J is the sum of the normal diffusion (random motility) term Jd and the chemotaxis term Jc. That is, J=Jd+Jc, where Jd=-ϵ2∇u and Jc=ϵu∇χ(v). It should be noted that the signs of these fluxes are different. The sign of Jc suggests that the chemotaxis term provides a negative diffusion effect, which suppresses the expansion of u. The second term on the right hand side of (1) is the Keller-Segel type chemotaxis term; we express the chemotactic sensitivity function χ as χ(v)=kcχ0(v), where χ0(v)=f0θv2/(v2+θ2) with θ=1.0. In order to satisfy the condition max⁡dχ0/dv=1, we choose f0=83/9. We call kc the intensity of chemotaxis [31].In (3), T and D represent the relaxation time and diffusion constant of w, respectively. Let us consider a situation in which w plays the role of feedback to suppress the static bifurcation and oscillatory bifurcation. For the rapid feedback mechanism, T and D must be small and large, respectively. In the limits T→0 and D→∞, w becomes a time-dependent but spatially independent variable, which is denoted by wt:(4)wt=1Ω0∫drw=1Ω0∫dru+v-W,where Ω0 is the area of the entire system. Replacing w in p(w) by wt, p(w) becomes a global coupling term. In the case where T is very small and D is very large, we reduce the three-component RD system to the following two-component RD system with a global coupling term:(5)τϵ∂u∂t=ϵ2∇2u-ϵ∇·u∇χv+fu,v-v,(6)∂v∂t=∇2v+u-μv,where(7)fu,v=-u+Hu-pu,v.The functional p{u,v} represents a global coupling term given by(8)pu,v=p0+α∫dru+v-W,where the integral is over the entire domain, and α and W are rescaled to absorb Ω0:α′=α/Ω0 and W′=WΩ0, and the primes are dropped. α corresponds to the intensity of the global coupling, and the value of p0 is chosen as p0=0.275. Hereinafter, we consider the above two-component RD system to be described by (5) and (6).In the absence of chemotaxis, (5) and (6) describe the same system proposed by Krischer and Mikhailov [32]. This system had an activator and an inhibitor, and, for large τ, the motionless spot (localized particle-like structure) in two dimensions was stable. With decreasing τ under a large α, it was shown that the system had a stable moving spot. In order to understand intuitively the bifurcation from the motionless spot to the moving spot, we consider the limit ϵ→0. In the limit ϵ→0, the boundary layer of u becomes an interface, as shown in Figure 1. The location of the interface is defined by the condition u(r,t)=p. In this limit, u and v satisfy the relation (u+v)=1 (inside the domain) and 0 (outside the domain), which is obtained from (5). In the limits of α→∞ and ϵ→0, the area of the spot is conserved; using (8), we obtain(9)∫dru+v=W.This restriction prohibits the expansion and oscillation of the spot; the translational bifurcation firstly occurs with decreasing τ. Even for finite values of a large α and small ϵ, the area of the spot is approximately conserved by the feedback mechanism. This supports the existence of a moving spot for small τ with large α. In contrast, for small α, the static bifurcation or oscillatory bifurcation firstly occurs with decreasing τ, and the motionless spot is destabilized to form an expanding wave or is disintegrated by unstable oscillation.Figure 1
Boundary layer and interface. The solid and dotted curves correspond tou(r,t) and u0(r,t)=limϵ→0u(r,t) along a radial direction, respectively, that is, u(r,t) and u0(r,t). In the finite ϵ,  u(r,t) has a boundary layer around r~R0 with the thickness of ϵ, where R0 is determined by the condition u(r=R0)=p.  u0(r,t) does not have a boundary layer but has a sharp interface at r=R0.Although the bifurcation from a motionless spot to a straight moving spot and the collision of two moving spots were studied in [32], other types of moving spot were not studied. Recently, the bifurcation of the spot from straight motion to circular motion was theoretically analyzed for the system described by (5) and (6) without the chemotaxis term [16]. The result suggested that there was a critical value τc for the bifurcation of stationary straight motion. Although the straight motion was stable for τ>τc, it was destabilized by rotation (straight-circular-motion) bifurcation for τ≤τc. In other words, the straight motion was destabilized into circular motion with increasing velocity.In this study, we assume that the localized domain (spot) of an activator exists under global feedback in the system described by (5) and (6). In the system, the chemotactic substance is secreted from inside the spot, and the motion of the spot is influenced by the chemotaxis. In order to study the influence of chemotaxis on the motion of a spot, we apply the technique reported in [16] to the system described by (5) and (6). For this purpose, we first derive the radially symmetric equilibrium solution of (5) and (6) in the limit ϵ→0. In this limit, (6) becomes(10)∂v∂t=∇2v+Hu-p-βv,where β=1+μ. Setting ∂v/∂t=0 in (10), (10) becomes(11)d2dr2+1rddrv-βv+HR0-r=0,where r=r and R0 is the equilibrium radius of the symmetric spot. We impose the boundary conditions on v(r) as(12)dv0dr=0,vR0+=vR0-,dvR0+dr=dvR0-dr,v∞=0.The radially symmetric equilibrium solution of (10) in two dimensions is given by(13)v¯r=1β-R0βK1βR0I0βr,0<r<R0,R0βI1βR0K0βr,R0<r<∞,where In and Kn are the modified Bessel functions of the first and second kinds of order n, respectively. The corresponding solution of u(r) is given by u¯(r)=H(R0-r)-v¯(r). The above equilibrium solution is employed to study the motion of the spot in the following sections.


## 3. Equation of Motion of the Interface

Asϵ is small, the width of the boundary layer is small, and, therefore, the value of u changes sharply while that of v changes smoothly around the boundary layer. The velocity of the flat interface is a function of the value of v on the interface. In order to consider the dynamics of the spot in the system, we derive the equation of motion of the interface. The time evolution equation of the interface is described [16, 30, 33] by(14)τV=τCh+ϵκ+L,where C(h) is the velocity of a flat interface; κ is the curvature of the interface, the sign of which is chosen such that it is positive when the center of the curvature is outside of the excited domain; and h is the value of v on the interface. V is the normal component of the velocity, and L is the Lagrange multiplier for the constraint of domain area conservation (9). For the moving spot, this condition corresponds to(15)∫dωVω=0,where dω is the infinitesimal length of the interface and the integral is carried over the entire interface.We first consider the velocity of the interface in one dimension in (5) and (6). With the detailed derivations given in Appendix A, the velocity τC(h) is obtained as(16)τCh=21/2-h-ph+p1-h-p1/2+dχdvidvdxi,where p is a global coupling term (8). τC(h) consists of two components: the velocity of the traveling front solution in one dimension and the chemotactic velocity.In two dimensions, the velocity of a flat interface directed alongn is given by(17)τCh=21/2-h-ph+p1-h-p1/2+kcχ0′vdvdni≡Ph+Qhn·∇vi,where χ0′(v)=dχ0(v)/dv and Ai represents the value of function A evaluated on the interface. n is a unit normal vector on the interface, and (dv/dn) is a normal derivative.When the motion of the interface of the spot is slow compared with the relaxation rate of the chemotactic substance, the left hand side of (10), (∂v/∂t), is small. In this case, we deal with the time derivative of v in (10) as a perturbation. The asymptotic solution of (10) is written by perturbation expansion using the Green function G as(18)vr,t=GH-GG∂H∂t+GGG∂2H∂t2-GGGG∂3H∂t3+⋯,where G is Green’s function satisfying the equation(19)∇2-βGr-r′=-δr-r′.In (18), GH represents the integral(20)GH=∫dr1Gr-r1Hr1,t,GG∂H∂t=∫dr1dr2Gr-r1Gr1-r2∂Hr2,t∂t,⋯=⋯,for brevity. In this study, we consider the situation in which the motionless spot is destabilized supercritically and the spot moves with an arbitrarily small velocity. Then, we can safely apply the perturbative expansion of v in terms of G as in (18).


## 4. Deformed Spot Dynamics

In this section, we derive the equation of motion of the spot. In order to describe the deformations of the spot, tensors are introduced. The tensors depend on time, and the time evolution equation of the tensors are derived following [16].


### 4.1. Description of Deformed Spot

We firstly describe the position and velocity of a deformed spot. The motion of the spot consists of two components: the motion of the center of gravity and the motion of the interface relative to the center of gravity. The center of gravity is denoted byρ, and the velocity of the center of gravity is given by(21)dρdt=1Ω∫dωVωRω≡v,where Ω is the area of the spot and R(ω) is a position vector measured from the center of gravity. We consider the case that R(ω) is a single valued function of ϕr, which is given as(22)Rω=Rϕrer,where er is a radial unit vector and R(ϕr) is the distance between the center of gravity and a point on the interface, which is directed to an angle ϕr with respect to the x axis. Using the definition of Rω,  dω is given by(23)dω=dRdϕrdϕr=R2+R′2dϕr,where R′=dR(ϕr)/dϕr. For the deformed spot, the radius is given by(24)Rϕr=R0+δRϕr,twith(25)δRϕr,t=∑n=-∞∞cnteinϕr,where i=-1, R0 is the equilibrium radius of the symmetric spot, and δR corresponds to the deviation. In the expansion of δR, c±n corresponds to a (2π/n)-periodic deformation. According to the area conservation equation (9), c0=0. In addition, we exclude the terms of n=±1 because the translational motion of the spot is incorporated in ρ. Thus, we set c±1=0. Using (24), the curvature κ and the normal component of the velocity V are given up to the first order of the deviations as(26)κϕr,t=-1R0-1R02∑n=-∞∞n2-1cnteinϕr,(27)Vϕr,t=v·n+∑n=-∞∞c˙nteinϕr,respectively, where the overdot denotes the time derivative.For a general functionA(r,t), the Fourier transformation and its reverse transformation in two dimensions are defined as(28)Aqt=∫drAr,te-iq·r,(29)Ar,t=12π2∫dqAqteiq·r≡∫qAqteiq·r,respectively, where ∫q=∫dq/(2π)2 for brevity. For an isolated domain that forms a single loop without crossing, the step function H(u-p)=H(R-r-ρ) is transformed by the Fourier transformation as(30)Hq=∫drHu-pe-iq·r=∫r-ρ<Rdre-iq·r.Substituting (24) into (30), we obtain the expansion of Hq up to the first order of δR as(31)Hq=Hq0+Hq1,Hq0=2πR0qJ1qR0,Hq1=2πR0∑ncnti-neinϕqJnqR0,where q=q, ϕq is the angle between q and the x-axis, and Jn is a Bessel function of the first kind of order n.


### 4.2. Equation of Motion of a Spot

We derive the equation of motion of a spot from (14). τC is given as functions of h and h~ by (17), in which we used the notation h~=(dv/dn)i=(n·∇)vi. When the velocity of the spot is small, the deformation of the circle is small. For this situation, we derive the deviations of h and h~ from their stationary values as a power series of τC. Using (18), we first expand h and h~ up to the fourth-order time derivatives as(32)h=h0+h1+h2+h3+h4,(33)h~=h~0+h~1+h~2+h~3+h~4,respectively, where(34)h0=∫qGqHqeiq·Rω,h1=h1,1+h1,2=i∫qv·qGq2Hqeiq·Rω-∫qGq2∂Hq∂teiq·Rω,h2=h2,1+h2,2+h2,3=-i∫qv˙·qGq3Hqeiq·Rω-∫qv·q2Gq3Hqeiq·Rω+∫qGq3∂2Hq∂t2eiq·Rω,h3=-i∫qv·q3Gq4Hqeiq·Rω,h4=∫qv·q4Gq5Hqeiq·Rω,(35)h~0=i∫qn·qGqHqeiq·Rω,h~1=h~1,1+h~1,2=-∫qn·qv·qGq2Hqeiq·Rω-i∫qn·qGq2∂Hq∂teiq·Rω,h~2=h~2,1+h~2,2+h~2,3=∫qn·qv˙·qGq3Hqeiq·Rω-i∫qn·qv·q2Gq3Hqeiq·Rω+i∫qn·qGq3∂2Hq∂t2eiq·Rω,h~3=∫qn·qv·q3Gq4Hqeiq·Rω,h~4=i∫qn·qv·q4Gq5Hqeiq·Rωwith(36)Gq=1q2+β.In the above expansions, we omitted the term proportional to(q·v)(∂Hq/∂t) in the expressions of h2 and h~2. In the derivation of h3 and h~3, the terms proportional to (∂3Hq/∂t3), (q·v)(∂2Hq/∂t2), (q·v)2(∂Hq/∂t), (q·v˙)(∂Hq/∂t), (q·v¨)Hq, and (q·v)(q·v˙)Hq were omitted. In the derivation of h4 and h~4, all the terms of the derivatives of v and Hq with respect to time were omitted. These terms were omitted, because they include higher order derivatives of time and are smaller than the remaining terms, or they vanish in the later integral (1/Ω)∫dωR(ω)· due to the orthogonality relation of trigonometric functions. We remark that we include up to the third-order of v: h3 and h~3, for the derivation of the equation of motion of the spot. The justification for considering terms up to this order is discussed later. On the other hand, for the derivation of equations of motion of tensors, we include up to the fourth order of v: h4 and h~4. The justification for considering terms up to this order is discussed in Section 4.4.When the spot is in a stationary motionless state, we assume that the spot is a circle with radiusR0. The values of v and (n·∇)v at the interface are denoted as h0(0) and h~0(0), respectively. Using (34), h0(0) is given by h0 with the substitutions Hq=Hq(0) and R=R(0), where R(0)=R0er. Similarly, h~0(0) is given using (35). When ϵ is small but finite, substituting h=h0(0) and h~=h~0(0) into τC(h) with V=0, (14) becomes(37)Ph00+Qh00h~00-ϵR0+L=0,where κ=-1/R0 is used. We note that the Lagrange multiplier L in (37) can be absorbed into the constant term p0 in P(h0(0)). In the limit ϵ→0, we can calculate h0(0) and h~0(0) by using v¯(r), which is given by (13) (the validation is given in Appendix B). Using these expressions, (37) gives the dependence of R0 on kc and W. The numerical results obtained by using (13) and (17) are shown in Figure 2(a). When α=0, R0 monotonically decreases with kc. However, for large α, R0 depends weakly on kc and it is almost constant because of the global coupling; large global coupling suppresses variations in the area of the spot resulting in constant R0. The dependence of R0 on W when kc is large is shown in Figure 2(b). From (9), it can be seen that, for large α, R0 approximately satisfies πR02~W, and, therefore, R0 is proportional to the square root of W (see the dashed curve in Figure 2(b)).Figure 2
Dependence ofR0 on kc and W in the stationary state. (a) Dependence of R0 on kc.  W=1.0. The solid and dashed curves represent the cases of α=0 and 5, respectively. (b) Dependence of R0 on W.  kc=5.0. The solid and dashed curves represent the cases of α=0.1 and 5.0, respectively.
(a) (b)When the spot moves with an infinitesimal velocity,h and h~ deviate from h0(0) and h~0(0), respectively. By putting small deviations as δh=h-h0(0) and δh~=h~-h~0(0), we iteratively derive δh and δh~ in a power series of τC up to the third order. The procedure consists of two steps. We first expand τC given by (17) in a power series of δh and δh~. Then, we solve the equation in terms of δh and expand it by τC. When the stationary motionless spot is destabilized into a moving spot, τC is expanded up to the second-order of δh and δh~ as(38)τC=∂P∂h+∂Q∂hh~000δh+Q0δh~+12∂2P∂h2+∂2Q∂h2h~000δh2+2∂Q∂h0δhδh~,where, for a function A(h), (A)0 implies A(h0(0)). Up to the first order of δh and δh~, (38) results in τC=(∂P/∂h+∂Q/∂hh~0(0))0δh+(Q)0δh~. For higher order corrections, a supplementary relation between δh and δh~ is necessary. In order to relate δh with δh~, we consider the profile of v(r) in the radially symmetric function given by (13), where v(r) is a function of r  (=r) denoted by v¯(r). For this function, (n·∇)v¯=dv¯/dr. Then, v¯(R0+δr) and (n·∇)v¯(R0+δr) are expanded around v¯(R0) and dv¯(R0)/dr for small δr as(39)v¯R0+δr~v¯R0+dv¯drr=R0δr,(40)n·∇v¯R0+δr~dv¯R0dr+d2v¯dr2r=R0δr,respectively. In (39) and (40), (A(r))r=R0 implies A(R0). We make an ansatz that δh and δh~ are not independent but have a linear relation such that δh~=kδh with k=(v¯′′/v¯′)r=R0, where the prime corresponds to the derivative with respect to r. This assumption enables us to calculate δh as a function of τC. Substituting δh~=kδh into (38), τC becomes a quadratic equation of δh as(41)τC=a+kbδh+m1δh2,where(42)a=∂P∂h+∂Q∂hh~000,b=Q0,m1=12∂2P∂h2+∂2Q∂h2h~000+2k∂Q∂h0.As a result of numerical calculations, we remark that b>0 and a<0 for small h0(0), but a>0 for large h0(0). By solving (41) in terms of δh for small τC, the solution is expanded up to the second order of τC as(43)δh=τCa+kb-m1a+kbτCa+kb2.Here, it should be noted that although there are two solutions for (41), we chose one solution such that δh→τC/(a+kb) when the second-order term (τC)2 is neglected.In order to obtainδh up to the third order of τC, τC given by (17) is expanded up to the third order of δh and δh~. Using δh~=kδh, τC becomes a cubic equation of δh as(44)τC=a+kbδh+m1δh2+m2δh3,where m2 is given by(45)m2=16∂3P∂h3+∂3Q∂h3h~000+3k∂2Q∂h20.In order to obtain the third-order term of τC in δh, we add a correction Δ  (~O(τC)3) to δh given by (43) and substitute it into (44). On solving this equation in terms of Δ, δh up to the third order of τC is obtained as(46)δh=τCa+kb-m1a+kbτCa+kb2+1a+kb2m12a+kb-m2τCa+kb3.We remark that (46) is obtained by another method; by applying Cardano’s formula to (44) and expanding the solution up to the third order of τC, we obtain the same result as that of (46).In the above process, we iteratively derivedδh up to the third order of τC. We replace kδh~δh~ in (46), and finally the power series expansion of aδh+bδh~ in terms of τC is obtained as(47)aδh+bδh~=M1τC+M2τC2+M3τC3,where Mi  (i=1,2,3) is defined by(48)M1=1,M2=-m1a+bk2,M3=1a+bk32m12a+bk-m2.In the absence of chemotaxis (Q=0), (47) reproduces the result that was obtained in [16].Using (14) and (47), we derive the equation of motion of the spot by operating both sides of (47) with (1/Ω)∫dωR(ω). For the left hand side, we put(49)1Ω∫dωRωδh=δh=h1+h2+h3,1Ω∫dωRωδh~=δh~=h~1+h~2+h~3,where h4 and h~4 are neglected because these terms are higher than or equal to O(v4). For the calculation of the right hand side of (47), we use (14). The magnitudes of ϵκ and L were discussed in [16]; L is independent of ϕr, and, owing to the periodicity of the function, ∫dωR(ω)L=∫dωR(ω)L3=∫dωR(ω)V2L=0. In addition, L~O(v2), and, therefore, ∫dωR(ω)VL3 can be neglected up to the third order ofv. As the translational motion of the spot is incorporated in ρ, we chose c±1=0 in the expansion of δR(ϕr). This results in ϵ∫dωR(ω)κ=0 up to the first order of the deformation. The terms ϵ∫dωR(ω)κVL2 and ϵ∫dωR(ω)κL2, and the higher order terms of ϵ were small; therefore, they were neglected. Following the above discussion, we neglect the term (ϵκ+L) in (14).By carrying the integral overϕr, the following equation is obtained:(50)aδh+bδh~=τv+34M3τ3vv2,where(51)δh=h1,1+h2,1+h3,δh~=h~1,1+h~2,1+h~3.Here, h1,1=(1/Ω)∫dωR(ω)h1,1. The other terms on the right hand side of (51) are defined similarly. In the expansion (51), we neglected the terms h1,2 and h2,3, because these terms yielded only the terms c˙±1 and c¨±1, respectively; however, these terms are set to zero. h2,2 disappeared in the integral over ϕr owing to the orthogonality relation of trigonometric functions [15]. Owing to the same reasons, the terms h~1,2, h~2,3, and h~2,2 did not exist. In the integral over ϕr, there is no contribution from the second-order term of τC owing to the orthogonality relation of trigonometric function. This supports the justification to include terms up to the third order of τC in (47) for the minimal equation of motion of a spot. For the balance of order, it is sufficient to include terms up to the third order ofv in δh and δh~.After practical calculations ofδh and δh~, we finally obtain the equation of motion of the spot as(52)mv˙+τ-τcv+gvv2=aδh1,10+bδh~1,10,where δh1,1(0) is defined as δh1,1(0)=h1,1-h1,1(0). Here, h1,1(0) is given by h1,1 with the substitutions Hq=Hq(0) and R=R(0).  δh~1,1(0) is similarly defined with δh1,1(0). The coefficients m, τc, and g on the left hand side of (52) are given by(53)m=-aR0F1,3,1,1-bR0Fd2,3,1,1,(54)τc=-aR0F1,2,1,1-bR0Fd2,2,1,1,(55)g=34M3τ3-34aR0F3,4,1,1-34bR0Fd4,4,1,1,respectively, where we defined F(i,j,k,l) and Fd(i,j,k,l) as(56)Fi,j,k,l=∫0∞dqqiGqjJkJl,(57)Fdi,j,k,l=∫0∞dqqiGqjJk∂∂qR0Jl,respectively, where i,j,k, and l are integers and Jn(qR0) is written as Jn for brevity.Equation (52) has the same form as the one in the absence of chemotaxis [16]. It is a Newtonian equation of the spot; the effective mass m, damping coefficient τc, intensity of the cubic nonlinear term g, and the right hand side of (52) form the coupling term between deformation and velocity in the following subsection. When (τ-τc)>0, the motionless spot is stable; however, the motionless spot is destabilized into a moving spot for (τ-τc)<0. When g is positive, the third-order nonlinear coupling term gvv2 suppresses the divergence for large v. It is necessary that m, τc, and g are positive for the existence of a stable moving spot. The dependence of m, τc, and g on kc and W is shown in Figure 3. We see that these parameters are positive and monotonically increase with kc and W.Figure 3
Dependence ofm, τc, and g on kc and W. The solid, dashed, and dotted curves represent m, τc, and g, respectively. (a) Dependence of m, τc, and g on kc.  α=5.0 and W=1.0. In the calculation of g, τ is chosen as τ=τc. As M3 is positive, g for the case τ<(>)τc is smaller (larger) than this curve. (b) Dependence of m, τc, and g on W.  α=5.0 and kc=5.0. In the calculation of g, τ is chosen as τ=τc.
(a) (b)In the derivation of (47), we made an ansatz for the relation between δh and δh~ such that δh~=kδh with k=(v¯′′/v¯′)r=R0. Equation (55) suggests that g depends explicitly on k through M3. Through preliminary research, it is found that if k is chosen as a positive constant, the property that g is a positive and monotonically increasing function of kc and W holds. From the above, it is seen that m, τc, and g are positive for any value of kc and W, and, therefore, we fix W and examine the motion of the spot in the range kc≥0.


### 4.3. Definitions of Tensors

In order to express the deformation of the spot, we introduce tensors. We rewrite the vector form (52) of each component, and, by denoting v=(v1,v2), the α component (α=1,2) of (52) becomes(58)mv˙α+τ-τcvα+gvαv2=aδh1,1α0+bδh~1,1α0,where δh1,1(α)(0) and δh~1,1(α)(0) correspond to the α component of δh1,1(0) and δh~1,1(0), respectively. Using the tensors, the right hand side of (58) is expressed in a simple form. The second rank tensor S and the third rank tensor U represent the elliptical deformation and the head-tail asymmetric deformation (2π/3-periodic deformation in the radial direction) of the spot, respectively. These tensors are traceless and symmetric. The detailed definitions are given in [16]. The deformation δR(ϕr) is expanded using the coefficients cn as given by (25); c±n corresponds to a (2π/n)-periodic deformation.We first introduce a second-rank tensorSαβ (α, β=1,2) using c±2 as follows:(59)S11=-S22=c2+c-2=δ22cos⁡2ϕ2,S12=S21=ic2-c-2=δ22sin⁡2ϕ2,where δ2 is a positive constant, which represents the radial deviation of the spot from R0, and ϕ2 is the angle between the long axis of the ellipse and the x-axis. The tensor elements (59) are represented by using a normal vector N=(cos⁡ϕ2,sin⁡ϕ2) along the long axis as(60)Sαβ=δ2NαNβ-δαβ2N2.The tensor S is the same as the nematic order parameter tensor in liquid crystals [34]. For an elliptical spot, R(ϕr) is represented as(61)Rϕr=R0+δ22cos⁡2ϕr-ϕ2.Next, in order to describe the head-tail asymmetric deformation, we first defineT1 and T2 using c±3 as(62)T1=c3+c-3=δ3cos⁡3ϕ3,(63)T2=ic3-c-3=δ3sin⁡3ϕ3,respectively, where δ3 is a positive constant and ϕ3 is the angle between one of the long axes of the deformed spot and the x-axis. In order to relate T1 and T2 with a tensor, we introduce the third-rank tensor Uαβγ (α, β, and γ=1 or 2). The tensor elements Uαβγ are represented by using vectors N(m) (m=1,2, and 3) as(64)Uαβγ=4δ33∑mNαmNβmNγm,where the normal vectors N(m) are defined by(65)N1=cos⁡ϕ3,sin⁡ϕ3,N2=cos⁡ϕ3+2π3,sin⁡ϕ3+2π3,N3=cos⁡ϕ3-2π3,sin⁡ϕ3-2π3.The tensor U is the same as the order parameter for banana (tetrahedral nematic) liquid crystals in two dimensions [35, 36]. We obtain the relations between the tensor elements Uαβγ and T1 and T2 as U111=T1, U222=-T2, and(66)U111=-U122=-U212=-U221,U222=-U112=-U121=-U211.The spot with head-tail asymmetric deformation is represented as(67)Rϕr=R0+δ3cos⁡3ϕr-ϕ3.The termsδh1,1(α)(0) and δh~1,1(α)(0) in (58) are expressed using tensors. The detailed calculations are given in Appendix C. The final form is as follows:(68)δh1,1α0=-a∗∑βSαβvβ,δh~1,1α0=-a~∗∑βSαβvβ,where(69)a∗=a1∗+a2∗,a1∗=-R02F2,2,1,2,(70)a2∗=-12R0F2,2,1,2-3F1,2,1,1,(71)a~∗=a~1∗+a~2∗,a~1∗=-R08-F3,2,1,2+3F3,2,1,1-34F2,2,0,1,(72)a~2∗=-2R04Fd2,2,2,1-Fd1,2,1,1.Thus, (58) is written using Sαβ in the form(73)mv˙α+τ-τcvα+gvαv2=-a¯∑βSαβvβ,where a¯=(aa∗+ba~∗).  δh1,1(α)(0) and δh~1,1(α)(0) yield the Sv term. This is due to the periodicity of the function in the integral over ϕr; in the expansion of δR(ϕr) and Hq(1)(ϕr), only c±2 terms contribute to the nonzero integral, resulting in the Sv term. Equation (73) suggests that the motionless spot is destabilized into a moving spot for τc>τ, and the velocity of the spot is approximately given by |v|2~(τc-τ)/g. The deformation and velocity coupling term Sv modifies the straight motion of the spot.


### 4.4. Time Evolution Equations of Tensors

In the previous subsection, we defined the tensors and described the equation of motion of the spot by using tensors, including the time-dependent tensorSαβ. In this subsection, we derive the time evolution equations of the tensors up to ~O(v4). We first discuss the order ofv, Sαβ, and Uαβγ, following [16]. From (73), the motionless spot is critical at τ=τc, and the moving spot occurs supercritically with increasing (τc-τ). We put δ=τc-τ for the smallness parameter. The time is scaled by t^=tδ, an all the terms in (73) are of the order O(δ3/2). Here v~O(δ1/2) and S~O(δ); therefore Sαβ~O(v2). In the later calculations, we can confirm that Sαβ~O(v2) and Uαβγ~O(v3). In the derivation of (73), we omitted the terms that include higher order derivatives of time. We estimate these terms, for example, d2v/dt2~O(δ5/2), SSv~O(δ5/2), (q·v)∂Hq/∂t~O(δ5/2), ∂3Hq/∂t3~O(δ4), (q·v)∂2Hq/∂t2~O(δ7/2), (q·v)2∂Hq/∂t~O(δ3), and (q·v˙)∂Hq/∂t~O(δ7/2) in the expansion of h and h~. From these estimates, the omission of those terms is justified.S and U are composed of c±2 and c±3, respectively, so that c±2~O(v2) and c±3~O(v3). The coupling terms such as cncm (n,|m|≥2) are much smaller. Therefore, we linearize (14) in terms of the deformation δR. Then, the time evolution equation of δR becomes(74)τ∂δR∂t=ϵR02∂2δR∂ϕr2+δR+aδh+bδh~+L.We consider terms up to the fourth order of v, where δh and δh~ are given by(75)δh=h-h00=δh0+δh1+δh2+δh3+δh4,(76)δh~=h~-h~00=δh~0+δh~1+δh~2+δh~3+δh~4,respectively. Here, δh0 and δhl (l=1,2,…,4) are defined by (34) as δh0=h0-h0(0) and δhl=hl, respectively. δh~l is similarly defined by (35). We expand δh1, δh2, δh~1, and δh~2 as δh1=δh1,1+δh1,2, δh2=δh2,2, δh~1=δh~1,1+δh~1,2, and δh~2=δh~2,2, respectively.In order to derive the equation of motion of tensors, we calculate the first-order time derivative ofcn. The linear combination of c˙n results in the equation of motion of tensors. In this process, v˙ terms in δh2,1 and δh~2,1 are neglected because these terms cause small order terms; ∑γUαβγv˙γ and Sαβv˙γ appear in the time evolution equation of Sαβ and Uαβγ, respectively. These terms are ∑γUαβγv˙γ~O(v6) and Sαβv˙γ~O(v5) so that they do not contribute to the equation of motion of tensors. In addition, (∂2Hq/∂t2) terms in δh2,3 and δh~2,3 are neglected because these terms result in the second-order time derivative of cn; these terms are c¨±2~O(v6) and c¨±3~O(v7) and they are small enough to be neglected.We first derive the equation of motion ofSαβ. Using (74), we can obtain the time evolution equation of c±2, and we derive the time evolution equation of Sαβ. The detailed derivation is given in Appendices D and E. The time evolution equation of Sαβ up to O(v4) is(77)Γ2dSαβdt=-K2Sαβ+b¯vαvβ-δαβ2v2+b¯svαvβ-δαβ2v2v2+b¯1∑γUαβγvγ+c¯v2Sαβ,where Γ2=τ+aE2+bE~2, K2=3ϵ/R02+aD2+bD~2, b¯=aG1+bG~1, b¯s=-aG1s-bG~1s, b¯1=-aB3/2+b(3A~∗-B~3)/2, and c¯=-a(2T2,2-V2)-b(2T~2,2-V~2). All the parameters in these expressions are given in Appendix E. The second term in (77) originates from δh2 and δh~2, the third term from δh4 and δh~4, the fourth term from δh1 and δh~1, and the fifth term from δh2 and δh~2. Thus, in order to consider the terms up to O(v4), δh4 and δh~4 are necessary in the expansion equations (75) and (76). From the first and second terms in (77), we note that Sαβ~O(v2).Next, we consider the time evolution equation ofUαβγ. Following procedures similar to those for Sαβ, the time evolution equation of Uαβγ up to O(v4) is(78)Γ3dUαβγdt=-K3Uαβγ+4d1vαvβvγ-v24δαβvγ+δβγvα+δγαvβ+2d23Sαβvγ+Sβγvα+Sγαvβ-∑ηvη2δαβSγη+δβγSαη+δγαSβη,where Γ3=τ+aE3+bE~3 and K3=8ϵ/R02+aD3+bD~3. d1 and d2 are defined as(79)d1=-R04aF3,4,1,3+bFd4,4,1,3,(80)d2=-aR02F2,2,2,3+Fd2,2,1,1-bR02Fd3,2,2,3+14F3,2,1,3-3F3,2,1,1-b2F2,2,0,1+F2,2,1,2,respectively. In the above equation, all the terms are ~O(v3). The terms of O(v4) disappeared because of the orthogonality relation of trigonometric functions. The next higher order term is O(v5), which is not included in expansions δh and δh~ in (75) and (76). The second term in (78) originates from δh3 and δh~3, and the third term from δh1 and δh~1. From the first and second terms in (78), we note that Uαβγ~O(v3).


## 4.1. Description of Deformed Spot

We firstly describe the position and velocity of a deformed spot. The motion of the spot consists of two components: the motion of the center of gravity and the motion of the interface relative to the center of gravity. The center of gravity is denoted byρ, and the velocity of the center of gravity is given by(21)dρdt=1Ω∫dωVωRω≡v,where Ω is the area of the spot and R(ω) is a position vector measured from the center of gravity. We consider the case that R(ω) is a single valued function of ϕr, which is given as(22)Rω=Rϕrer,where er is a radial unit vector and R(ϕr) is the distance between the center of gravity and a point on the interface, which is directed to an angle ϕr with respect to the x axis. Using the definition of Rω,  dω is given by(23)dω=dRdϕrdϕr=R2+R′2dϕr,where R′=dR(ϕr)/dϕr. For the deformed spot, the radius is given by(24)Rϕr=R0+δRϕr,twith(25)δRϕr,t=∑n=-∞∞cnteinϕr,where i=-1, R0 is the equilibrium radius of the symmetric spot, and δR corresponds to the deviation. In the expansion of δR, c±n corresponds to a (2π/n)-periodic deformation. According to the area conservation equation (9), c0=0. In addition, we exclude the terms of n=±1 because the translational motion of the spot is incorporated in ρ. Thus, we set c±1=0. Using (24), the curvature κ and the normal component of the velocity V are given up to the first order of the deviations as(26)κϕr,t=-1R0-1R02∑n=-∞∞n2-1cnteinϕr,(27)Vϕr,t=v·n+∑n=-∞∞c˙nteinϕr,respectively, where the overdot denotes the time derivative.For a general functionA(r,t), the Fourier transformation and its reverse transformation in two dimensions are defined as(28)Aqt=∫drAr,te-iq·r,(29)Ar,t=12π2∫dqAqteiq·r≡∫qAqteiq·r,respectively, where ∫q=∫dq/(2π)2 for brevity. For an isolated domain that forms a single loop without crossing, the step function H(u-p)=H(R-r-ρ) is transformed by the Fourier transformation as(30)Hq=∫drHu-pe-iq·r=∫r-ρ<Rdre-iq·r.Substituting (24) into (30), we obtain the expansion of Hq up to the first order of δR as(31)Hq=Hq0+Hq1,Hq0=2πR0qJ1qR0,Hq1=2πR0∑ncnti-neinϕqJnqR0,where q=q, ϕq is the angle between q and the x-axis, and Jn is a Bessel function of the first kind of order n.


## 4.2. Equation of Motion of a Spot

We derive the equation of motion of a spot from (14). τC is given as functions of h and h~ by (17), in which we used the notation h~=(dv/dn)i=(n·∇)vi. When the velocity of the spot is small, the deformation of the circle is small. For this situation, we derive the deviations of h and h~ from their stationary values as a power series of τC. Using (18), we first expand h and h~ up to the fourth-order time derivatives as(32)h=h0+h1+h2+h3+h4,(33)h~=h~0+h~1+h~2+h~3+h~4,respectively, where(34)h0=∫qGqHqeiq·Rω,h1=h1,1+h1,2=i∫qv·qGq2Hqeiq·Rω-∫qGq2∂Hq∂teiq·Rω,h2=h2,1+h2,2+h2,3=-i∫qv˙·qGq3Hqeiq·Rω-∫qv·q2Gq3Hqeiq·Rω+∫qGq3∂2Hq∂t2eiq·Rω,h3=-i∫qv·q3Gq4Hqeiq·Rω,h4=∫qv·q4Gq5Hqeiq·Rω,(35)h~0=i∫qn·qGqHqeiq·Rω,h~1=h~1,1+h~1,2=-∫qn·qv·qGq2Hqeiq·Rω-i∫qn·qGq2∂Hq∂teiq·Rω,h~2=h~2,1+h~2,2+h~2,3=∫qn·qv˙·qGq3Hqeiq·Rω-i∫qn·qv·q2Gq3Hqeiq·Rω+i∫qn·qGq3∂2Hq∂t2eiq·Rω,h~3=∫qn·qv·q3Gq4Hqeiq·Rω,h~4=i∫qn·qv·q4Gq5Hqeiq·Rωwith(36)Gq=1q2+β.In the above expansions, we omitted the term proportional to(q·v)(∂Hq/∂t) in the expressions of h2 and h~2. In the derivation of h3 and h~3, the terms proportional to (∂3Hq/∂t3), (q·v)(∂2Hq/∂t2), (q·v)2(∂Hq/∂t), (q·v˙)(∂Hq/∂t), (q·v¨)Hq, and (q·v)(q·v˙)Hq were omitted. In the derivation of h4 and h~4, all the terms of the derivatives of v and Hq with respect to time were omitted. These terms were omitted, because they include higher order derivatives of time and are smaller than the remaining terms, or they vanish in the later integral (1/Ω)∫dωR(ω)· due to the orthogonality relation of trigonometric functions. We remark that we include up to the third-order of v: h3 and h~3, for the derivation of the equation of motion of the spot. The justification for considering terms up to this order is discussed later. On the other hand, for the derivation of equations of motion of tensors, we include up to the fourth order of v: h4 and h~4. The justification for considering terms up to this order is discussed in Section 4.4.When the spot is in a stationary motionless state, we assume that the spot is a circle with radiusR0. The values of v and (n·∇)v at the interface are denoted as h0(0) and h~0(0), respectively. Using (34), h0(0) is given by h0 with the substitutions Hq=Hq(0) and R=R(0), where R(0)=R0er. Similarly, h~0(0) is given using (35). When ϵ is small but finite, substituting h=h0(0) and h~=h~0(0) into τC(h) with V=0, (14) becomes(37)Ph00+Qh00h~00-ϵR0+L=0,where κ=-1/R0 is used. We note that the Lagrange multiplier L in (37) can be absorbed into the constant term p0 in P(h0(0)). In the limit ϵ→0, we can calculate h0(0) and h~0(0) by using v¯(r), which is given by (13) (the validation is given in Appendix B). Using these expressions, (37) gives the dependence of R0 on kc and W. The numerical results obtained by using (13) and (17) are shown in Figure 2(a). When α=0, R0 monotonically decreases with kc. However, for large α, R0 depends weakly on kc and it is almost constant because of the global coupling; large global coupling suppresses variations in the area of the spot resulting in constant R0. The dependence of R0 on W when kc is large is shown in Figure 2(b). From (9), it can be seen that, for large α, R0 approximately satisfies πR02~W, and, therefore, R0 is proportional to the square root of W (see the dashed curve in Figure 2(b)).Figure 2
Dependence ofR0 on kc and W in the stationary state. (a) Dependence of R0 on kc.  W=1.0. The solid and dashed curves represent the cases of α=0 and 5, respectively. (b) Dependence of R0 on W.  kc=5.0. The solid and dashed curves represent the cases of α=0.1 and 5.0, respectively.
(a) (b)When the spot moves with an infinitesimal velocity,h and h~ deviate from h0(0) and h~0(0), respectively. By putting small deviations as δh=h-h0(0) and δh~=h~-h~0(0), we iteratively derive δh and δh~ in a power series of τC up to the third order. The procedure consists of two steps. We first expand τC given by (17) in a power series of δh and δh~. Then, we solve the equation in terms of δh and expand it by τC. When the stationary motionless spot is destabilized into a moving spot, τC is expanded up to the second-order of δh and δh~ as(38)τC=∂P∂h+∂Q∂hh~000δh+Q0δh~+12∂2P∂h2+∂2Q∂h2h~000δh2+2∂Q∂h0δhδh~,where, for a function A(h), (A)0 implies A(h0(0)). Up to the first order of δh and δh~, (38) results in τC=(∂P/∂h+∂Q/∂hh~0(0))0δh+(Q)0δh~. For higher order corrections, a supplementary relation between δh and δh~ is necessary. In order to relate δh with δh~, we consider the profile of v(r) in the radially symmetric function given by (13), where v(r) is a function of r  (=r) denoted by v¯(r). For this function, (n·∇)v¯=dv¯/dr. Then, v¯(R0+δr) and (n·∇)v¯(R0+δr) are expanded around v¯(R0) and dv¯(R0)/dr for small δr as(39)v¯R0+δr~v¯R0+dv¯drr=R0δr,(40)n·∇v¯R0+δr~dv¯R0dr+d2v¯dr2r=R0δr,respectively. In (39) and (40), (A(r))r=R0 implies A(R0). We make an ansatz that δh and δh~ are not independent but have a linear relation such that δh~=kδh with k=(v¯′′/v¯′)r=R0, where the prime corresponds to the derivative with respect to r. This assumption enables us to calculate δh as a function of τC. Substituting δh~=kδh into (38), τC becomes a quadratic equation of δh as(41)τC=a+kbδh+m1δh2,where(42)a=∂P∂h+∂Q∂hh~000,b=Q0,m1=12∂2P∂h2+∂2Q∂h2h~000+2k∂Q∂h0.As a result of numerical calculations, we remark that b>0 and a<0 for small h0(0), but a>0 for large h0(0). By solving (41) in terms of δh for small τC, the solution is expanded up to the second order of τC as(43)δh=τCa+kb-m1a+kbτCa+kb2.Here, it should be noted that although there are two solutions for (41), we chose one solution such that δh→τC/(a+kb) when the second-order term (τC)2 is neglected.In order to obtainδh up to the third order of τC, τC given by (17) is expanded up to the third order of δh and δh~. Using δh~=kδh, τC becomes a cubic equation of δh as(44)τC=a+kbδh+m1δh2+m2δh3,where m2 is given by(45)m2=16∂3P∂h3+∂3Q∂h3h~000+3k∂2Q∂h20.In order to obtain the third-order term of τC in δh, we add a correction Δ  (~O(τC)3) to δh given by (43) and substitute it into (44). On solving this equation in terms of Δ, δh up to the third order of τC is obtained as(46)δh=τCa+kb-m1a+kbτCa+kb2+1a+kb2m12a+kb-m2τCa+kb3.We remark that (46) is obtained by another method; by applying Cardano’s formula to (44) and expanding the solution up to the third order of τC, we obtain the same result as that of (46).In the above process, we iteratively derivedδh up to the third order of τC. We replace kδh~δh~ in (46), and finally the power series expansion of aδh+bδh~ in terms of τC is obtained as(47)aδh+bδh~=M1τC+M2τC2+M3τC3,where Mi  (i=1,2,3) is defined by(48)M1=1,M2=-m1a+bk2,M3=1a+bk32m12a+bk-m2.In the absence of chemotaxis (Q=0), (47) reproduces the result that was obtained in [16].Using (14) and (47), we derive the equation of motion of the spot by operating both sides of (47) with (1/Ω)∫dωR(ω). For the left hand side, we put(49)1Ω∫dωRωδh=δh=h1+h2+h3,1Ω∫dωRωδh~=δh~=h~1+h~2+h~3,where h4 and h~4 are neglected because these terms are higher than or equal to O(v4). For the calculation of the right hand side of (47), we use (14). The magnitudes of ϵκ and L were discussed in [16]; L is independent of ϕr, and, owing to the periodicity of the function, ∫dωR(ω)L=∫dωR(ω)L3=∫dωR(ω)V2L=0. In addition, L~O(v2), and, therefore, ∫dωR(ω)VL3 can be neglected up to the third order ofv. As the translational motion of the spot is incorporated in ρ, we chose c±1=0 in the expansion of δR(ϕr). This results in ϵ∫dωR(ω)κ=0 up to the first order of the deformation. The terms ϵ∫dωR(ω)κVL2 and ϵ∫dωR(ω)κL2, and the higher order terms of ϵ were small; therefore, they were neglected. Following the above discussion, we neglect the term (ϵκ+L) in (14).By carrying the integral overϕr, the following equation is obtained:(50)aδh+bδh~=τv+34M3τ3vv2,where(51)δh=h1,1+h2,1+h3,δh~=h~1,1+h~2,1+h~3.Here, h1,1=(1/Ω)∫dωR(ω)h1,1. The other terms on the right hand side of (51) are defined similarly. In the expansion (51), we neglected the terms h1,2 and h2,3, because these terms yielded only the terms c˙±1 and c¨±1, respectively; however, these terms are set to zero. h2,2 disappeared in the integral over ϕr owing to the orthogonality relation of trigonometric functions [15]. Owing to the same reasons, the terms h~1,2, h~2,3, and h~2,2 did not exist. In the integral over ϕr, there is no contribution from the second-order term of τC owing to the orthogonality relation of trigonometric function. This supports the justification to include terms up to the third order of τC in (47) for the minimal equation of motion of a spot. For the balance of order, it is sufficient to include terms up to the third order ofv in δh and δh~.After practical calculations ofδh and δh~, we finally obtain the equation of motion of the spot as(52)mv˙+τ-τcv+gvv2=aδh1,10+bδh~1,10,where δh1,1(0) is defined as δh1,1(0)=h1,1-h1,1(0). Here, h1,1(0) is given by h1,1 with the substitutions Hq=Hq(0) and R=R(0).  δh~1,1(0) is similarly defined with δh1,1(0). The coefficients m, τc, and g on the left hand side of (52) are given by(53)m=-aR0F1,3,1,1-bR0Fd2,3,1,1,(54)τc=-aR0F1,2,1,1-bR0Fd2,2,1,1,(55)g=34M3τ3-34aR0F3,4,1,1-34bR0Fd4,4,1,1,respectively, where we defined F(i,j,k,l) and Fd(i,j,k,l) as(56)Fi,j,k,l=∫0∞dqqiGqjJkJl,(57)Fdi,j,k,l=∫0∞dqqiGqjJk∂∂qR0Jl,respectively, where i,j,k, and l are integers and Jn(qR0) is written as Jn for brevity.Equation (52) has the same form as the one in the absence of chemotaxis [16]. It is a Newtonian equation of the spot; the effective mass m, damping coefficient τc, intensity of the cubic nonlinear term g, and the right hand side of (52) form the coupling term between deformation and velocity in the following subsection. When (τ-τc)>0, the motionless spot is stable; however, the motionless spot is destabilized into a moving spot for (τ-τc)<0. When g is positive, the third-order nonlinear coupling term gvv2 suppresses the divergence for large v. It is necessary that m, τc, and g are positive for the existence of a stable moving spot. The dependence of m, τc, and g on kc and W is shown in Figure 3. We see that these parameters are positive and monotonically increase with kc and W.Figure 3
Dependence ofm, τc, and g on kc and W. The solid, dashed, and dotted curves represent m, τc, and g, respectively. (a) Dependence of m, τc, and g on kc.  α=5.0 and W=1.0. In the calculation of g, τ is chosen as τ=τc. As M3 is positive, g for the case τ<(>)τc is smaller (larger) than this curve. (b) Dependence of m, τc, and g on W.  α=5.0 and kc=5.0. In the calculation of g, τ is chosen as τ=τc.
(a) (b)In the derivation of (47), we made an ansatz for the relation between δh and δh~ such that δh~=kδh with k=(v¯′′/v¯′)r=R0. Equation (55) suggests that g depends explicitly on k through M3. Through preliminary research, it is found that if k is chosen as a positive constant, the property that g is a positive and monotonically increasing function of kc and W holds. From the above, it is seen that m, τc, and g are positive for any value of kc and W, and, therefore, we fix W and examine the motion of the spot in the range kc≥0.


## 4.3. Definitions of Tensors

In order to express the deformation of the spot, we introduce tensors. We rewrite the vector form (52) of each component, and, by denoting v=(v1,v2), the α component (α=1,2) of (52) becomes(58)mv˙α+τ-τcvα+gvαv2=aδh1,1α0+bδh~1,1α0,where δh1,1(α)(0) and δh~1,1(α)(0) correspond to the α component of δh1,1(0) and δh~1,1(0), respectively. Using the tensors, the right hand side of (58) is expressed in a simple form. The second rank tensor S and the third rank tensor U represent the elliptical deformation and the head-tail asymmetric deformation (2π/3-periodic deformation in the radial direction) of the spot, respectively. These tensors are traceless and symmetric. The detailed definitions are given in [16]. The deformation δR(ϕr) is expanded using the coefficients cn as given by (25); c±n corresponds to a (2π/n)-periodic deformation.We first introduce a second-rank tensorSαβ (α, β=1,2) using c±2 as follows:(59)S11=-S22=c2+c-2=δ22cos⁡2ϕ2,S12=S21=ic2-c-2=δ22sin⁡2ϕ2,where δ2 is a positive constant, which represents the radial deviation of the spot from R0, and ϕ2 is the angle between the long axis of the ellipse and the x-axis. The tensor elements (59) are represented by using a normal vector N=(cos⁡ϕ2,sin⁡ϕ2) along the long axis as(60)Sαβ=δ2NαNβ-δαβ2N2.The tensor S is the same as the nematic order parameter tensor in liquid crystals [34]. For an elliptical spot, R(ϕr) is represented as(61)Rϕr=R0+δ22cos⁡2ϕr-ϕ2.Next, in order to describe the head-tail asymmetric deformation, we first defineT1 and T2 using c±3 as(62)T1=c3+c-3=δ3cos⁡3ϕ3,(63)T2=ic3-c-3=δ3sin⁡3ϕ3,respectively, where δ3 is a positive constant and ϕ3 is the angle between one of the long axes of the deformed spot and the x-axis. In order to relate T1 and T2 with a tensor, we introduce the third-rank tensor Uαβγ (α, β, and γ=1 or 2). The tensor elements Uαβγ are represented by using vectors N(m) (m=1,2, and 3) as(64)Uαβγ=4δ33∑mNαmNβmNγm,where the normal vectors N(m) are defined by(65)N1=cos⁡ϕ3,sin⁡ϕ3,N2=cos⁡ϕ3+2π3,sin⁡ϕ3+2π3,N3=cos⁡ϕ3-2π3,sin⁡ϕ3-2π3.The tensor U is the same as the order parameter for banana (tetrahedral nematic) liquid crystals in two dimensions [35, 36]. We obtain the relations between the tensor elements Uαβγ and T1 and T2 as U111=T1, U222=-T2, and(66)U111=-U122=-U212=-U221,U222=-U112=-U121=-U211.The spot with head-tail asymmetric deformation is represented as(67)Rϕr=R0+δ3cos⁡3ϕr-ϕ3.The termsδh1,1(α)(0) and δh~1,1(α)(0) in (58) are expressed using tensors. The detailed calculations are given in Appendix C. The final form is as follows:(68)δh1,1α0=-a∗∑βSαβvβ,δh~1,1α0=-a~∗∑βSαβvβ,where(69)a∗=a1∗+a2∗,a1∗=-R02F2,2,1,2,(70)a2∗=-12R0F2,2,1,2-3F1,2,1,1,(71)a~∗=a~1∗+a~2∗,a~1∗=-R08-F3,2,1,2+3F3,2,1,1-34F2,2,0,1,(72)a~2∗=-2R04Fd2,2,2,1-Fd1,2,1,1.Thus, (58) is written using Sαβ in the form(73)mv˙α+τ-τcvα+gvαv2=-a¯∑βSαβvβ,where a¯=(aa∗+ba~∗).  δh1,1(α)(0) and δh~1,1(α)(0) yield the Sv term. This is due to the periodicity of the function in the integral over ϕr; in the expansion of δR(ϕr) and Hq(1)(ϕr), only c±2 terms contribute to the nonzero integral, resulting in the Sv term. Equation (73) suggests that the motionless spot is destabilized into a moving spot for τc>τ, and the velocity of the spot is approximately given by |v|2~(τc-τ)/g. The deformation and velocity coupling term Sv modifies the straight motion of the spot.


## 4.4. Time Evolution Equations of Tensors

In the previous subsection, we defined the tensors and described the equation of motion of the spot by using tensors, including the time-dependent tensorSαβ. In this subsection, we derive the time evolution equations of the tensors up to ~O(v4). We first discuss the order ofv, Sαβ, and Uαβγ, following [16]. From (73), the motionless spot is critical at τ=τc, and the moving spot occurs supercritically with increasing (τc-τ). We put δ=τc-τ for the smallness parameter. The time is scaled by t^=tδ, an all the terms in (73) are of the order O(δ3/2). Here v~O(δ1/2) and S~O(δ); therefore Sαβ~O(v2). In the later calculations, we can confirm that Sαβ~O(v2) and Uαβγ~O(v3). In the derivation of (73), we omitted the terms that include higher order derivatives of time. We estimate these terms, for example, d2v/dt2~O(δ5/2), SSv~O(δ5/2), (q·v)∂Hq/∂t~O(δ5/2), ∂3Hq/∂t3~O(δ4), (q·v)∂2Hq/∂t2~O(δ7/2), (q·v)2∂Hq/∂t~O(δ3), and (q·v˙)∂Hq/∂t~O(δ7/2) in the expansion of h and h~. From these estimates, the omission of those terms is justified.S and U are composed of c±2 and c±3, respectively, so that c±2~O(v2) and c±3~O(v3). The coupling terms such as cncm (n,|m|≥2) are much smaller. Therefore, we linearize (14) in terms of the deformation δR. Then, the time evolution equation of δR becomes(74)τ∂δR∂t=ϵR02∂2δR∂ϕr2+δR+aδh+bδh~+L.We consider terms up to the fourth order of v, where δh and δh~ are given by(75)δh=h-h00=δh0+δh1+δh2+δh3+δh4,(76)δh~=h~-h~00=δh~0+δh~1+δh~2+δh~3+δh~4,respectively. Here, δh0 and δhl (l=1,2,…,4) are defined by (34) as δh0=h0-h0(0) and δhl=hl, respectively. δh~l is similarly defined by (35). We expand δh1, δh2, δh~1, and δh~2 as δh1=δh1,1+δh1,2, δh2=δh2,2, δh~1=δh~1,1+δh~1,2, and δh~2=δh~2,2, respectively.In order to derive the equation of motion of tensors, we calculate the first-order time derivative ofcn. The linear combination of c˙n results in the equation of motion of tensors. In this process, v˙ terms in δh2,1 and δh~2,1 are neglected because these terms cause small order terms; ∑γUαβγv˙γ and Sαβv˙γ appear in the time evolution equation of Sαβ and Uαβγ, respectively. These terms are ∑γUαβγv˙γ~O(v6) and Sαβv˙γ~O(v5) so that they do not contribute to the equation of motion of tensors. In addition, (∂2Hq/∂t2) terms in δh2,3 and δh~2,3 are neglected because these terms result in the second-order time derivative of cn; these terms are c¨±2~O(v6) and c¨±3~O(v7) and they are small enough to be neglected.We first derive the equation of motion ofSαβ. Using (74), we can obtain the time evolution equation of c±2, and we derive the time evolution equation of Sαβ. The detailed derivation is given in Appendices D and E. The time evolution equation of Sαβ up to O(v4) is(77)Γ2dSαβdt=-K2Sαβ+b¯vαvβ-δαβ2v2+b¯svαvβ-δαβ2v2v2+b¯1∑γUαβγvγ+c¯v2Sαβ,where Γ2=τ+aE2+bE~2, K2=3ϵ/R02+aD2+bD~2, b¯=aG1+bG~1, b¯s=-aG1s-bG~1s, b¯1=-aB3/2+b(3A~∗-B~3)/2, and c¯=-a(2T2,2-V2)-b(2T~2,2-V~2). All the parameters in these expressions are given in Appendix E. The second term in (77) originates from δh2 and δh~2, the third term from δh4 and δh~4, the fourth term from δh1 and δh~1, and the fifth term from δh2 and δh~2. Thus, in order to consider the terms up to O(v4), δh4 and δh~4 are necessary in the expansion equations (75) and (76). From the first and second terms in (77), we note that Sαβ~O(v2).Next, we consider the time evolution equation ofUαβγ. Following procedures similar to those for Sαβ, the time evolution equation of Uαβγ up to O(v4) is(78)Γ3dUαβγdt=-K3Uαβγ+4d1vαvβvγ-v24δαβvγ+δβγvα+δγαvβ+2d23Sαβvγ+Sβγvα+Sγαvβ-∑ηvη2δαβSγη+δβγSαη+δγαSβη,where Γ3=τ+aE3+bE~3 and K3=8ϵ/R02+aD3+bD~3. d1 and d2 are defined as(79)d1=-R04aF3,4,1,3+bFd4,4,1,3,(80)d2=-aR02F2,2,2,3+Fd2,2,1,1-bR02Fd3,2,2,3+14F3,2,1,3-3F3,2,1,1-b2F2,2,0,1+F2,2,1,2,respectively. In the above equation, all the terms are ~O(v3). The terms of O(v4) disappeared because of the orthogonality relation of trigonometric functions. The next higher order term is O(v5), which is not included in expansions δh and δh~ in (75) and (76). The second term in (78) originates from δh3 and δh~3, and the third term from δh1 and δh~1. From the first and second terms in (78), we note that Uαβγ~O(v3).


## 5. Tensor Model

Using (73), (77), and (78), we discuss the stationary solution, the stability conditions of straight motion, and the critical velocity of circular motion. For the sake of convenience, we rescale (73), (77), and (78) into(81)dvαdt=γvα-vαv2-a¯′∑βSαβvβ,(82)dSαβdt=-κSαβ+b¯′vαvβ-12v2δαβ+b¯s′vαvβ-12v2δαβv2+b¯1′∑γUαβγvγ+c¯′v2Sαβ,(83)dUαβγdt=-KUαβγ+d1′vαvβvγ-v24δαβvγ+δβγvα+δγαvβ+d2′3Sαβvγ+Sβγvα+Sγαvβ-∑ηvη2δαβSγη+δβγSαη+δγαSβη,respectively. We name (81)–(83) as the full tensor model [37]. In the above full tensor model, we chose parameters such that g/m=1 and put γ=(τc-τ)/m and a¯′=a¯/m in (81). The other rescaled parameters in (82) and (83) were κ=K2/Γ2, b¯′=b¯/Γ2, b¯s′=b¯s/Γ2, b¯1′=b¯1/Γ2, c¯′=c¯/Γ2, K=K3/Γ3, d1′=2d1/Γ3, and d2′=d2/Γ3.We summarize the coefficients in (81)–(83) as follows. γ is a damping coefficient in the time evolution equation of vα; although the motionless spot is stable for γ<0, the motionless spot is destabilized for γ>0 and there appears a moving spot. κ and K are damping coefficients in the time evolution equations of Sαβ and Uαβγ, respectively. In this study, we fix K as a positive constant and consider κ in the range κ>0. b¯′ and b¯s′ are coefficients of the quadratic and fourth order of v, respectively, in the time evolution equation of Sαβ. The term b¯s′ can be absorbed into b¯′ with the replacement b¯′→(b¯′+b¯s′v2). From the results, it can be seen that the term b¯′ is enhanced (suppressed) when the sign of b¯s′ is the same as (different from) that of b¯′. d1′ is the coefficient of the cubic term of v in the time evolution equation of Uαβγ. These terms cause deformations with increasing velocity. The other parameters are coupling coefficients of deformation and velocity. a¯′ and d2′ are the coupling coefficients of Sαβvβ and Sαβvγ, respectively, and these terms influence the time evolution of vα and Uαβγ, respectively. b¯1′ and c¯′ are the coupling coefficients of Uαβγvγ and v2Sαβ, respectively. When b¯1′=0, Uαβγ is decoupled from vα and Sαβ. c¯′v2Sαβ can be absorbed into -κSαβ with the replacement κ→(κ-c¯′v2). Then, the term κ is suppressed (enhanced) when c¯′ is positive (negative), with v2≠0.The dependence ofa¯′, b¯′, b¯s′, c¯′, b¯1′, d1′, and d2′ on kc is shown in Figure 4. In each figure, these parameters (shown on the vertical axis) are scaled by R02. We see that both a¯′ and b¯′ are negative, and b¯1′, c¯′≠0 even at kc=0. Although b¯1′ is a monotonically increasing function of kc, c¯′ is a monotonically decreasing function of kc. d1′ and d2′ are positive but d1′ is much smaller than d2′.Figure 4
Dependence ofa-′, b-′, b-s′, c-′, b-1′, d1′, and d2′ on kc. α=5.0 and W=1.0. All the parameters in the vertical axis are scaled by R02. (a) a-′ and b-′. The solid and dashed curves represent a-′ and b-′, respectively. (b) b-s′. (c) c-′. (d) b-1′. (e) d1′ and d2′. The solid and dashed curves represent d1′ and d2′, respectively.
(a) (b) (c) (d) (e)Henceforth, we drop the primes on the parameters in the full tensor model (81)–(83). In the full tensor model, when K is large and v is small, the effect of Uαβγ on Sαβ is small in (82). In this case, we put b¯1=0, and the vα, Sαβ, and Uαβγ system described by (81)–(83) is reduced to a vα and Sαβ system. We call this system the reduced tensor model.


## 6. Stationary Solution and Phase Diagrams

In the following subsections, we consider the stationary solution and phase diagrams. For this, we rewritevα, Sαβ, and Uαβγ by introducing the following variables:(84)v1=vcos⁡ϕ,v2=vsin⁡ϕ,S11=s2cos⁡2θ,S12=s2sin⁡2θ,U111=c3+c-3,U222=-ic3-c-3,c±3=δ32e∓i3ϕ3≡z4e∓i3φ,where we choose v, s, z>0. From (81), we obtain the time evolution equations of v and ϕ as(85)ddtv=vγ-v2-12a¯svcos⁡2θ-ϕ,(86)ddtϕ=-12a¯ssin⁡2θ-ϕ,respectively. From (82), we obtain the time evolution equations of s and θ as(87)ddts=-κs+b¯v2cos⁡2θ-ϕ+b¯sv4cos⁡2θ-ϕ+b¯1zvcos⁡2θ+ϕ-3φ+c¯v2s,(88)ddtθ=-b¯v22ssin⁡2θ-ϕ-b¯sv42ssin⁡2θ-ϕ+b¯1zv2ssin⁡3φ-ϕ-2θ,respectively. From (83), we obtain the time evolution equations of z and φ as(89)ddtz=-Kz+d1v3cos⁡3φ-ϕ+d2svcos⁡3φ-2θ-ϕ,(90)ddtφ=-d13zv3sin⁡3φ-ϕ-d23zsvsin⁡3φ-2θ-ϕ.Using (85)–(90), we theoretically derive the stationary solution and phase diagrams. In order to verify the theoretical result, we perform simulations of the reduced tensor model. For the numerical calculations, we employ the fourth-order Runge-Kutta algorithm with the time increment Δt=1.0×10-4 and fix the parameters as a¯=b¯=-1, K=0.3, d1=0.1, and d2=0.5.


### 6.1. Stationary Solution

For the reduced and full tensor models, we consider the stationary solution of the moving spot, which moves straight in thex direction. This situation corresponds to ϕ=0. In the case of b¯1=b¯s=c¯=0, the stationary solution of (85) and (86) becomes v02=γ/(1+B) with B=a¯b¯/2κ. In order to avoid the divergence of v02, we assume that a¯b¯>0. For c¯≠0, we consider the dependence of the stationary solution of s, z, θ, and φ on the signs of a¯ and b¯. In the stationary state of (86) with ϕ=0, there are two stationary solutions of θ depending on the sign of a¯ and b¯: (I) a¯ and b¯ are positive with θ0=lπ and (II) a¯ and b¯ are negative with θ0=(l+1/2)π, where l is an integer. Using (87)–(90), the stationary solution of case (I) is(91)s0=b¯+b¯sv02K+b¯1d1v02Kκ-c¯v02-d2b¯1v02v02,z0=cos⁡3φ0d1κ-c¯v02+d2b¯+b¯sv02Kκ-c¯v02-d2b¯1v02v03,where φ0 satisfies sin⁡3φ0= 0 and v0 is a stationary solution of v, which is determined later. On the other hand, for case (II), the stationary solution is(92)s0=-b¯+b¯sv02K+b¯1d1v02Kκ-c¯v02-d2b¯1v02v02,z0=cos⁡3φ0d1κ-c¯v02+d2b¯+b¯sv02Kκ-c¯v02-d2b¯1v02v03.From these results, it is seen that z0 has the same expression irrespective of the signs of a¯ and b¯. As z0 is positive, φ0 satisfying sin⁡3φ0=0 is chosen as(93)iφ0=2m3π,ifd1κ-c¯v02+d2b¯+b¯sv02Kκ-c¯v02-v02d2b¯1>0,iiφ0=2m+13π,ifd1κ-c¯v02+d2b¯+b¯sv02Kκ-c¯v02-v02d2b¯1<0.For both cases ((I) and (II)), the stationary solution of z0 is given by(94)z0=d1κ-c¯v02+d2b¯+b¯sv02Kκ-c¯v02-d2b¯1v02v03.In both cases ((I) and (II)), by substituting s0 into (85) in the stationary state, we obtain two stationary solutions v0,± as(95)v0,±2=1f1f2±f3withf1=2Kc¯+d2b¯1-12a¯b¯1d1+Kb¯s,f2=12a¯b¯+κ+γc¯K+γd2b¯1,f3=12a¯b¯-κ+γc¯K+γd2b¯12+2a¯κb¯+b¯sγK2+γb¯1d1K.In order to satisfy the condition that v0,±2→γ/(1+B) in the limits b¯1, b¯s, and c¯→0, only v0,-2 exists. In the physical images of the moving spot in cases (I) and (II) with φ0=0, the long axes of the elliptical domain are parallel and perpendicular to the moving direction, respectively [16].


### 6.2. Stability of Straight Motion in the Reduced Tensor Model

WhenK is large in (89), z is rapidly damped to zero. In addition, when v is small, the influences of z on s and θ in (87) and (88) are small. For this case, we discuss the stability of the straight motion of the spot by setting b¯1=0 in (87) and (88). Using (86) and (88), we can derive the time evolution equation of ψ=θ-ϕ as(96)dψdt=-12-a¯s+b¯v2s+b¯sv4ssin⁡2ψ.For the stationary statesv=v0,-, s=s0, we discuss the stability condition of the straight motion. From (96), the stationary solution ψ0 satisfies 2ψ0=2lπ or (2l+1)π, and their corresponding stability conditions are (-a¯s0+b¯v0,-2/s0+b¯sv0,-4/s0)>0 or (-a¯s0+b¯v0,-2/s0+b¯sv0,-4/s0)<0, respectively. Since a¯ and b¯ are negative in our RD system, as shown in Figure 4(a), we consider the case of ψ0=π/2 as follows. In this case, the stability condition is calculated using the s0 given by (92). The bifurcation of straight-circular motion, when b¯s=c¯=0, was already studied in [16], and the critical value for the straight motion is expressed as(97)γc=κ2a¯b¯+κ2.The phase diagram in the κ-γ plane is shown by the solid curve in Figure 5(a). When γ≤0, the motionless spot is stable. When γ is in the range 0<γ≤γc, the spot moves straight. However, when γ>γc, the straight motion become unstable and changes into circular motion. On the other hand, when b¯s≠0 and c¯≠0, the stability condition for straight motion is(98)a¯b¯v0,-2+b¯sv0,-4κ-c¯v0,-2<κ-c¯v0,-2.The phase diagram is a curved surface in κ-γ-c¯ space; the boundary of different motions of spot is obtained by using (98). The phase diagram in the κ-γ plane is shown in Figure 5(a) with dashed and dotted curves. Comparing the curves with that for the case when b¯s=c¯=0 (the solid curve), it can be seen that when c¯>0, the parameter region for the circular motion is larger. However, when c¯<0, the parameter region for the circular motion is smaller. We can explain this result using (87); the term c¯v2s essentially reduces the parameter κ of the damping term -κs such that the damping term becomes -(κ-c¯v2)s. When c¯ is positive (negative), the term c¯v2s effectively reduces (increases) κ so that the deformation s becomes large (small). This leads to positive (negative) c¯ results in the larger parameter region of the circular motion. In our RD system, c¯ is negative and c¯ monotonically increases with kc (see Figure 4(c)), and the parameter region for circular motion becomes smaller as the intensity of chemotaxis increases. The phase diagrams in the c¯-γ plane and c¯-κ plane are shown in Figures 5(b) and 5(c), respectively. Figure 5(b) suggests that the parameter region for the circular motion is smaller for a larger value of κ. Equation (87) suggests that a larger value of κ damps s strongly such that the deformation is small, resulting in a smaller region of the circular motion. On the other hand, Figure 5(c) suggests that the parameter region for the circular motion is larger for a larger value of γ. Equation (85) suggests that a large positive value of γ enhances v such that the deformation becomes large, resulting in a larger region of circular motion. In each figure, the simulation results for the reduced tensor model are shown; they agree well with the theoretical results.Figure 5
Phase diagram of spot motion.a-=b-=-1 and b-s=0.05. Data are obtained by the analysis of (96). M, S, and C in the figure represent regions of no motion, straight motion, and circular motion, respectively. The mark ∗ represents the data obtained by the simulation of the reduced tensor model. (a) Phase diagram in the κ-γ plane. The solid, dashed, and dotted curves correspond to the cases of c-=0, 0.2, and -0.2, respectively. The upper (lower) region of each curve corresponds to the region of circular (straight) motion, and the region γ≤0 corresponds to the region of no motion. (b) Phase diagram in the c--γ plane. The solid, dashed, and dotted curves correspond to cases of κ=0.2, 0.5, and 0.8, respectively. The region across each curve is the same as in (a). (c) Phase diagram in the c--κ plane. The solid, dashed, and dotted curves correspond to the cases of γ=0.25, 0.5, and 1.0, respectively. The upper (lower) region of each curve corresponds to the region of straight (circular) motion.
(a) (b) (c)


### 6.3. Dependence of Critical Velocity on Parameters

In the previous subsection, we showed that a spot appears in the circular motion whenγ>γc. In this parameter region, we examine the stationary circular motion of the spot with a constant angular frequency ω and velocity vr. For this, we set b¯1=0 and substitute v=vr, s=sr, θ=ωt+ζ/2, and ϕ=ωt into (85)–(88), and obtain the relations among ζ, sr, vr, and ω. The calculation is straightforward; the final expressions are shown.We first consider the case whenb¯s=c¯=0, where up to O(v3) terms are considered, and obtain the relations cos⁡ζ=κ/a¯sr, sr2=b¯vr2/a¯, vr2=(γ-κ/2), ω2=(a¯b¯/4)(vr2-vc2), and vc2=κ2/a¯b¯. When the spot rotates with an angular frequency ω, ω satisfies ω2≥ 0, and this condition leads to vr≥vc [16]. The dependence of the critical velocity vc for the stable angular frequency on κ  (≥0) is shown in Figure 6(a). In this case, there is a lower limit of velocity vc for the circular motion but no upper limit of velocity for the same.Figure 6
Dependence of the critical velocityvc, vc,±, Δvc2, and c-∗ on κ and c-. a-=b-=-1. The mark ∗ represents data obtained by the simulation of the reduced tensor model. (a) Dependence of vc2 on κ.  b-s=c-=0. (b) Dependence of vc,+2 on κ.  b-s=-0.2 and c-=-0.25. (c) Dependence of vc,±2 on κ.  b-s=0.05. The curves with filled (•) and open (∘) circles correspond to the cases of c-=-0.25 and c-=-0.5, respectively. The solid and dashed curves represent vc,-2 and vc,+2, respectively. (d) Dependence of Δvc2 on c-.  b-s=0.05. The solid, dashed, and dotted curves correspond to κ=0.1, 0.5, and 0.8, respectively. (e) Dependence of c-∗ on κ.  b-s=0.05.
(a) (b) (c) (d) (e)Next, we consider the case whenb¯s≠0 and c¯≠0, where up to O(v4) terms are considered. After similar calculations to the case of b¯s=c¯=0, we obtain the relations cos⁡ζ=(κ-c¯vr2)/a¯sr, sr2=(b¯vr2+b¯svr4)/a¯, vr2=(γ-κ/2)/(1-c¯/2), and ω2=(a¯b¯s-c¯2)vr4/4+(a¯b¯+2κc¯)vr2/4-κ2/4. In this case, ω2 is a quartic function of vr. For the condition ω2≥0, there are two cases depending on the sign of (a¯b¯s-c¯2). When (a¯b¯s-c¯2)>0, ω2 exists for vr2≥vc,+2, where vc,±2 are the solutions of ω2=0, which are given by vc,±2=[-(a¯b¯+2κc¯)±Dr]/[2(a¯b¯s-c¯2)] with Dr=(a¯b¯+2κc¯)2+4κ2(a¯b¯s-c¯2). The dependence of vc,+2 on κ is shown in Figure 6(b); there is only a lower limit of velocity vc,+. This property is the same with the case when b¯c=c¯=0. On the other hand, when (a¯b¯s-c¯2)<0, ω2 exists in a certain range of vr2; vc,+2≤vr2≤vc,-2 under the condition Dr≥0. In this case, there is a lower limit and an upper limit of velocity, vc,+ and vc,-, respectively. The gap between vc,-2 and vc,+2 is defined by Δvc2=vc,-2-vc,+2=Dr/(c¯2-a¯b¯s). Let us consider two limiting cases. In the limit c¯→∞, vc,±2 and Δvc2→0. That is, the circular motion cannot exist. On the other hand, for small c¯ and b¯s≪1, we expand vc,±2 in terms of c¯ and b¯s and obtain vc,+2~κ2/(a¯b¯+2κc¯) and vc,-2~(a¯b¯+2κc¯)/(c¯2-a¯b¯s). Thus, in the limits b¯s and c¯→0, vc,+2→κ2/a¯b¯, vc,-2→∞, and Δvc2→∞. This suggests that vc,+→vc and there is no upper limit of velocity, which are the same with the case when b¯s=c¯=0. From the numerical results shown in Figures 4(a), 4(b), and 4(c), (a¯b¯s-c¯2)<0 in our RD system; therefore, it turns out that vr2 exists only in a finite range, vc,+2≤vr2≤vc,-2. For larger c¯, the gap Δvc,±2 is smaller, as shown in Figure 6(c). The dependence of Δvc2 on c¯ is shown in Figure 6(d), which suggests that Δvc2 monotonically decreases with increasing c¯. Since |c¯| increases monotonically as kc increases (Figure 4(c)), Δvc2 decreases as the chemotaxis increases. When the condition Dr=0 is satisfied, vc,-=vc,+ and Δvc2=0, which yields the critical value c¯∗ as(99)c¯∗=-a¯b¯2+4κ2a¯b¯s4κa¯b¯.When c¯>c¯∗, Dr<0, and, therefore, vc,±2 does not exist. There is no parameter region of γ for a positive value of κ; the circular motion does not occur for any values of γ. The dependence of c¯∗ on κ is shown in Figure 6(e). Equation (99) can be positive for large κ. However, when c¯ is chosen as a positive value, the velocities v and s increase monotonically with increasing γ, and, therefore, the system diverges in the simulation. For this reason, c¯∗ is shown in the range c¯∗≤0.In order to verify the above theoretical results, the simulation results for the reduced tensor model are shown in each figure. The results obtained via simulation agree well with the theoretical results. In Figures6(a) and 6(b), for given κ and c¯, there is only a lower limit of γ for the circular motion of the spot, resulting in vc and vc,+, respectively. On the other hand, in Figure 6(c), for given κ and c¯, there is not only a lower limit of γ but also an upper limit of γ for the circular motion of the spot, resulting in vc,+ and vc,-, respectively. However, in Figure 6(e), for given κ and c¯∗, there is no parameter region of γ for the circular motion of the spot.The dependence of the stationary velocity and radius of the circular motion onc¯ are shown in Figure 7. The velocity v increases monotonically with increasing c¯, as shown in Figure 7(a). The corresponding radius of the spot motion Rc is shown in Figure 7(b). For fixed values of γ and κ, v decreases and Rc increases as c¯ increases in the range c¯≤0. The effect of chemotaxis is incorporated in c¯, and c¯ is negative; c¯ increases as the intensity of chemotaxis kc increases in our RD system (Figure 4(c)). Thus, it is seen that the chemotaxis reduces v and increases Rc.Figure 7
Dependence of velocity and radius onc- in the stationary circular motion. Data are obtained by the simulation of the reduced tensor model. b-s=0.1. The solid, dashed, and dotted curves correspond to (γ,κ)=(0.25,0.1), (0.5,0.2), and (1.0,0.3), respectively. (a) Dependence of v on c-. (b) Dependence of Rc on c-.
(a) (b)From the above analyses on (85)–(88) with b¯1=0, we conclude four points on the properties of stationary circular motion. (i) There is a lower limit of velocity. (ii) Although there is no upper limit for the velocity when up to O(v3) terms are considered, there will be an upper limit for the velocity if up to O(v4) terms are considered. (iii) The range of velocity decreases as the chemotaxis increases. (iv) There is a critical value c¯∗ (corresponding critical intensity of the chemotaxis is denoted by kc∗) such that when c¯>c¯∗ (kc>kc∗), the circular motion does not occur for any values of γ.


### 6.4. Simulation of the Full Tensor Model

In this subsection, we examine the effect ofUαβγ on vα and Sαβ. For this, we consider the case of b¯1≠0. The data are obtained by the simulation of the full tensor model.The phase diagrams of the spot motion are shown in Figure8. The phase diagrams in the c¯-γ plane and c¯-κ plane are shown in Figures 8(a) and 8(b), respectively. We note that the behaviors in the full tensor model are similar to that in the reduced model (Figures 5(b) and 5(c)). For positive (negative) b¯1, the region of the circular motion is larger (smaller). The phase diagrams in the b¯1-γ plane and b¯1-κ plane are shown in Figures 8(c) and 8(d), respectively. The behaviors are similar to the phase diagrams in the c¯-γ plane and c¯-κ plane. This result suggests that the effect of ∑γUαβγvγ on the time evolution equation of Sαβ is similar to that of v2Sαβ.Figure 8
Phase diagram of spot motion. Data are obtained by the simulation of the full tensor model.b-s=0.1. M, S, and C in the figure represent regions of no motion, straight motion, and circular motion, respectively. (a) Phase diagram in the c--γ plane. κ=0.8. The solid, dashed, and dotted curves correspond to the cases of b-1=0, 0.25, and −0.25, respectively. The upper (lower) region of each curve corresponds to the region of circular (straight) motion, and the region γ≤0 corresponds to the region of no motion. (b) Phase diagram in the c--κ plane. γ=0.5. The solid, dashed, and dotted curves correspond to the cases of b-1=0, 0.25, and −0.25, respectively. The upper (lower) region of each curve corresponds to the region of straight (circular) motion. (c) Phase diagram in the b-1-γ plane. κ=0.8. The solid, dashed, and dotted curves correspond to the cases of c-=0, 0.25, and −0.25, respectively. The region across each curve is the same as in (a). (d) Phase diagram in the b-1-κ plane. γ=0.5. The solid, dashed, and dotted curves correspond to the cases of c-=0, 0.25, and −0.25, respectively. The region across each curve is the same as in (b).
(a) (b) (c) (d)The phase diagram of spot motion in thec¯-b¯1 plane is shown in Figure 9. It can be seen that the circular motion occurs in the parameter region of positive and large values of c¯ and b¯1. Sαβ and Uαβγ measure the elliptical deformation and head-tail asymmetric deformation of the spot, respectively; these deformations lead to the circular motion.Figure 9
Phase diagram of spot motion in thec--b-1 plane. Data are obtained by the simulation of the full tensor model. b-s=0.1. The marks S and C in the figure represent regions of straight motion and circular motion, respectively. The solid, dashed, and dotted curves correspond to the cases of (γ,κ)=(0.25,0.3), (0.5,0.5), and (1.0,0.8), respectively.The effects ofc¯ and b¯1 on the stationary velocity v and the corresponding radius of circular motion Rc are shown in Figure 10. On comparing these figures, it can be seen that although v and Rc depend on c¯ and b¯1, the influences of b¯1 on v and Rc are much smaller than those of c¯.Figure 10
Dependence of velocity and radius onc- and b-1 in the stationary circular motion. Data are obtained by the simulation of the full tensor model. b-s= 0.1, γ= 0.5, and κ=0.2. (a) Dependence of v on c-. (b) Dependence of Rc on c-, where, in (a) and (b), the solid, dashed, and dotted curves correspond to b-1=0, 0.5, and −0.5, respectively. (c) Dependence of v on b-1. (d) Dependence of Rc on b-1, where, in (c) and (d), the solid, dashed, and dotted curves correspond to c-=0, −0.25, and −0.5, respectively.
(a) (b) (c) (d)The dependence of the critical velocityvc,± and c¯∗ on κ is shown in Figure 11. There exists a critical velocity vc± for b¯1≠0 as shown in Figure 11(a). For larger b¯1, the gap Δvc2 is larger. This suggests that b¯1∑γUαβγvγ enhances the circular motion, which is consistent with the results shown in Figure 8. For positive b¯1, c¯∗ becomes larger than that in the reduced tensor model (solid curve) as shown in Figure 11(b). On the other hand, for negative b¯1, there is no c¯∗. In our system, b¯1 is positive and increases as the intensity of chemotaxis kc increases as shown in Figure 4(d). Thus, c¯∗ still exists in the full tensor model, and c¯∗ is larger than that in the reduced model.Figure 11
Dependence of the critical velocityvc,± and c-∗ on κ. Data are obtained by the simulation of the full tensor model. b-s=0.1. (a) Dependence of vc,±2 on κ.  c-=-0.5. The curves with filled (•) and open (∘) circles correspond to the cases of b-1=0.25 and b-1=0.5, respectively. The solid and dashed curves represent vc,-2 and vc,+2, respectively. (b) Dependence of c-∗ on κ.  b-s=0.1. The solid, dashed, and dotted curves correspond to the cases of b-1=0, 0.25, and 0.5, respectively.
(a) (b)


## 6.1. Stationary Solution

For the reduced and full tensor models, we consider the stationary solution of the moving spot, which moves straight in thex direction. This situation corresponds to ϕ=0. In the case of b¯1=b¯s=c¯=0, the stationary solution of (85) and (86) becomes v02=γ/(1+B) with B=a¯b¯/2κ. In order to avoid the divergence of v02, we assume that a¯b¯>0. For c¯≠0, we consider the dependence of the stationary solution of s, z, θ, and φ on the signs of a¯ and b¯. In the stationary state of (86) with ϕ=0, there are two stationary solutions of θ depending on the sign of a¯ and b¯: (I) a¯ and b¯ are positive with θ0=lπ and (II) a¯ and b¯ are negative with θ0=(l+1/2)π, where l is an integer. Using (87)–(90), the stationary solution of case (I) is(91)s0=b¯+b¯sv02K+b¯1d1v02Kκ-c¯v02-d2b¯1v02v02,z0=cos⁡3φ0d1κ-c¯v02+d2b¯+b¯sv02Kκ-c¯v02-d2b¯1v02v03,where φ0 satisfies sin⁡3φ0= 0 and v0 is a stationary solution of v, which is determined later. On the other hand, for case (II), the stationary solution is(92)s0=-b¯+b¯sv02K+b¯1d1v02Kκ-c¯v02-d2b¯1v02v02,z0=cos⁡3φ0d1κ-c¯v02+d2b¯+b¯sv02Kκ-c¯v02-d2b¯1v02v03.From these results, it is seen that z0 has the same expression irrespective of the signs of a¯ and b¯. As z0 is positive, φ0 satisfying sin⁡3φ0=0 is chosen as(93)iφ0=2m3π,ifd1κ-c¯v02+d2b¯+b¯sv02Kκ-c¯v02-v02d2b¯1>0,iiφ0=2m+13π,ifd1κ-c¯v02+d2b¯+b¯sv02Kκ-c¯v02-v02d2b¯1<0.For both cases ((I) and (II)), the stationary solution of z0 is given by(94)z0=d1κ-c¯v02+d2b¯+b¯sv02Kκ-c¯v02-d2b¯1v02v03.In both cases ((I) and (II)), by substituting s0 into (85) in the stationary state, we obtain two stationary solutions v0,± as(95)v0,±2=1f1f2±f3withf1=2Kc¯+d2b¯1-12a¯b¯1d1+Kb¯s,f2=12a¯b¯+κ+γc¯K+γd2b¯1,f3=12a¯b¯-κ+γc¯K+γd2b¯12+2a¯κb¯+b¯sγK2+γb¯1d1K.In order to satisfy the condition that v0,±2→γ/(1+B) in the limits b¯1, b¯s, and c¯→0, only v0,-2 exists. In the physical images of the moving spot in cases (I) and (II) with φ0=0, the long axes of the elliptical domain are parallel and perpendicular to the moving direction, respectively [16].


## 6.2. Stability of Straight Motion in the Reduced Tensor Model

WhenK is large in (89), z is rapidly damped to zero. In addition, when v is small, the influences of z on s and θ in (87) and (88) are small. For this case, we discuss the stability of the straight motion of the spot by setting b¯1=0 in (87) and (88). Using (86) and (88), we can derive the time evolution equation of ψ=θ-ϕ as(96)dψdt=-12-a¯s+b¯v2s+b¯sv4ssin⁡2ψ.For the stationary statesv=v0,-, s=s0, we discuss the stability condition of the straight motion. From (96), the stationary solution ψ0 satisfies 2ψ0=2lπ or (2l+1)π, and their corresponding stability conditions are (-a¯s0+b¯v0,-2/s0+b¯sv0,-4/s0)>0 or (-a¯s0+b¯v0,-2/s0+b¯sv0,-4/s0)<0, respectively. Since a¯ and b¯ are negative in our RD system, as shown in Figure 4(a), we consider the case of ψ0=π/2 as follows. In this case, the stability condition is calculated using the s0 given by (92). The bifurcation of straight-circular motion, when b¯s=c¯=0, was already studied in [16], and the critical value for the straight motion is expressed as(97)γc=κ2a¯b¯+κ2.The phase diagram in the κ-γ plane is shown by the solid curve in Figure 5(a). When γ≤0, the motionless spot is stable. When γ is in the range 0<γ≤γc, the spot moves straight. However, when γ>γc, the straight motion become unstable and changes into circular motion. On the other hand, when b¯s≠0 and c¯≠0, the stability condition for straight motion is(98)a¯b¯v0,-2+b¯sv0,-4κ-c¯v0,-2<κ-c¯v0,-2.The phase diagram is a curved surface in κ-γ-c¯ space; the boundary of different motions of spot is obtained by using (98). The phase diagram in the κ-γ plane is shown in Figure 5(a) with dashed and dotted curves. Comparing the curves with that for the case when b¯s=c¯=0 (the solid curve), it can be seen that when c¯>0, the parameter region for the circular motion is larger. However, when c¯<0, the parameter region for the circular motion is smaller. We can explain this result using (87); the term c¯v2s essentially reduces the parameter κ of the damping term -κs such that the damping term becomes -(κ-c¯v2)s. When c¯ is positive (negative), the term c¯v2s effectively reduces (increases) κ so that the deformation s becomes large (small). This leads to positive (negative) c¯ results in the larger parameter region of the circular motion. In our RD system, c¯ is negative and c¯ monotonically increases with kc (see Figure 4(c)), and the parameter region for circular motion becomes smaller as the intensity of chemotaxis increases. The phase diagrams in the c¯-γ plane and c¯-κ plane are shown in Figures 5(b) and 5(c), respectively. Figure 5(b) suggests that the parameter region for the circular motion is smaller for a larger value of κ. Equation (87) suggests that a larger value of κ damps s strongly such that the deformation is small, resulting in a smaller region of the circular motion. On the other hand, Figure 5(c) suggests that the parameter region for the circular motion is larger for a larger value of γ. Equation (85) suggests that a large positive value of γ enhances v such that the deformation becomes large, resulting in a larger region of circular motion. In each figure, the simulation results for the reduced tensor model are shown; they agree well with the theoretical results.Figure 5
Phase diagram of spot motion.a-=b-=-1 and b-s=0.05. Data are obtained by the analysis of (96). M, S, and C in the figure represent regions of no motion, straight motion, and circular motion, respectively. The mark ∗ represents the data obtained by the simulation of the reduced tensor model. (a) Phase diagram in the κ-γ plane. The solid, dashed, and dotted curves correspond to the cases of c-=0, 0.2, and -0.2, respectively. The upper (lower) region of each curve corresponds to the region of circular (straight) motion, and the region γ≤0 corresponds to the region of no motion. (b) Phase diagram in the c--γ plane. The solid, dashed, and dotted curves correspond to cases of κ=0.2, 0.5, and 0.8, respectively. The region across each curve is the same as in (a). (c) Phase diagram in the c--κ plane. The solid, dashed, and dotted curves correspond to the cases of γ=0.25, 0.5, and 1.0, respectively. The upper (lower) region of each curve corresponds to the region of straight (circular) motion.
(a) (b) (c)


## 6.3. Dependence of Critical Velocity on Parameters

In the previous subsection, we showed that a spot appears in the circular motion whenγ>γc. In this parameter region, we examine the stationary circular motion of the spot with a constant angular frequency ω and velocity vr. For this, we set b¯1=0 and substitute v=vr, s=sr, θ=ωt+ζ/2, and ϕ=ωt into (85)–(88), and obtain the relations among ζ, sr, vr, and ω. The calculation is straightforward; the final expressions are shown.We first consider the case whenb¯s=c¯=0, where up to O(v3) terms are considered, and obtain the relations cos⁡ζ=κ/a¯sr, sr2=b¯vr2/a¯, vr2=(γ-κ/2), ω2=(a¯b¯/4)(vr2-vc2), and vc2=κ2/a¯b¯. When the spot rotates with an angular frequency ω, ω satisfies ω2≥ 0, and this condition leads to vr≥vc [16]. The dependence of the critical velocity vc for the stable angular frequency on κ  (≥0) is shown in Figure 6(a). In this case, there is a lower limit of velocity vc for the circular motion but no upper limit of velocity for the same.Figure 6
Dependence of the critical velocityvc, vc,±, Δvc2, and c-∗ on κ and c-. a-=b-=-1. The mark ∗ represents data obtained by the simulation of the reduced tensor model. (a) Dependence of vc2 on κ.  b-s=c-=0. (b) Dependence of vc,+2 on κ.  b-s=-0.2 and c-=-0.25. (c) Dependence of vc,±2 on κ.  b-s=0.05. The curves with filled (•) and open (∘) circles correspond to the cases of c-=-0.25 and c-=-0.5, respectively. The solid and dashed curves represent vc,-2 and vc,+2, respectively. (d) Dependence of Δvc2 on c-.  b-s=0.05. The solid, dashed, and dotted curves correspond to κ=0.1, 0.5, and 0.8, respectively. (e) Dependence of c-∗ on κ.  b-s=0.05.
(a) (b) (c) (d) (e)Next, we consider the case whenb¯s≠0 and c¯≠0, where up to O(v4) terms are considered. After similar calculations to the case of b¯s=c¯=0, we obtain the relations cos⁡ζ=(κ-c¯vr2)/a¯sr, sr2=(b¯vr2+b¯svr4)/a¯, vr2=(γ-κ/2)/(1-c¯/2), and ω2=(a¯b¯s-c¯2)vr4/4+(a¯b¯+2κc¯)vr2/4-κ2/4. In this case, ω2 is a quartic function of vr. For the condition ω2≥0, there are two cases depending on the sign of (a¯b¯s-c¯2). When (a¯b¯s-c¯2)>0, ω2 exists for vr2≥vc,+2, where vc,±2 are the solutions of ω2=0, which are given by vc,±2=[-(a¯b¯+2κc¯)±Dr]/[2(a¯b¯s-c¯2)] with Dr=(a¯b¯+2κc¯)2+4κ2(a¯b¯s-c¯2). The dependence of vc,+2 on κ is shown in Figure 6(b); there is only a lower limit of velocity vc,+. This property is the same with the case when b¯c=c¯=0. On the other hand, when (a¯b¯s-c¯2)<0, ω2 exists in a certain range of vr2; vc,+2≤vr2≤vc,-2 under the condition Dr≥0. In this case, there is a lower limit and an upper limit of velocity, vc,+ and vc,-, respectively. The gap between vc,-2 and vc,+2 is defined by Δvc2=vc,-2-vc,+2=Dr/(c¯2-a¯b¯s). Let us consider two limiting cases. In the limit c¯→∞, vc,±2 and Δvc2→0. That is, the circular motion cannot exist. On the other hand, for small c¯ and b¯s≪1, we expand vc,±2 in terms of c¯ and b¯s and obtain vc,+2~κ2/(a¯b¯+2κc¯) and vc,-2~(a¯b¯+2κc¯)/(c¯2-a¯b¯s). Thus, in the limits b¯s and c¯→0, vc,+2→κ2/a¯b¯, vc,-2→∞, and Δvc2→∞. This suggests that vc,+→vc and there is no upper limit of velocity, which are the same with the case when b¯s=c¯=0. From the numerical results shown in Figures 4(a), 4(b), and 4(c), (a¯b¯s-c¯2)<0 in our RD system; therefore, it turns out that vr2 exists only in a finite range, vc,+2≤vr2≤vc,-2. For larger c¯, the gap Δvc,±2 is smaller, as shown in Figure 6(c). The dependence of Δvc2 on c¯ is shown in Figure 6(d), which suggests that Δvc2 monotonically decreases with increasing c¯. Since |c¯| increases monotonically as kc increases (Figure 4(c)), Δvc2 decreases as the chemotaxis increases. When the condition Dr=0 is satisfied, vc,-=vc,+ and Δvc2=0, which yields the critical value c¯∗ as(99)c¯∗=-a¯b¯2+4κ2a¯b¯s4κa¯b¯.When c¯>c¯∗, Dr<0, and, therefore, vc,±2 does not exist. There is no parameter region of γ for a positive value of κ; the circular motion does not occur for any values of γ. The dependence of c¯∗ on κ is shown in Figure 6(e). Equation (99) can be positive for large κ. However, when c¯ is chosen as a positive value, the velocities v and s increase monotonically with increasing γ, and, therefore, the system diverges in the simulation. For this reason, c¯∗ is shown in the range c¯∗≤0.In order to verify the above theoretical results, the simulation results for the reduced tensor model are shown in each figure. The results obtained via simulation agree well with the theoretical results. In Figures6(a) and 6(b), for given κ and c¯, there is only a lower limit of γ for the circular motion of the spot, resulting in vc and vc,+, respectively. On the other hand, in Figure 6(c), for given κ and c¯, there is not only a lower limit of γ but also an upper limit of γ for the circular motion of the spot, resulting in vc,+ and vc,-, respectively. However, in Figure 6(e), for given κ and c¯∗, there is no parameter region of γ for the circular motion of the spot.The dependence of the stationary velocity and radius of the circular motion onc¯ are shown in Figure 7. The velocity v increases monotonically with increasing c¯, as shown in Figure 7(a). The corresponding radius of the spot motion Rc is shown in Figure 7(b). For fixed values of γ and κ, v decreases and Rc increases as c¯ increases in the range c¯≤0. The effect of chemotaxis is incorporated in c¯, and c¯ is negative; c¯ increases as the intensity of chemotaxis kc increases in our RD system (Figure 4(c)). Thus, it is seen that the chemotaxis reduces v and increases Rc.Figure 7
Dependence of velocity and radius onc- in the stationary circular motion. Data are obtained by the simulation of the reduced tensor model. b-s=0.1. The solid, dashed, and dotted curves correspond to (γ,κ)=(0.25,0.1), (0.5,0.2), and (1.0,0.3), respectively. (a) Dependence of v on c-. (b) Dependence of Rc on c-.
(a) (b)From the above analyses on (85)–(88) with b¯1=0, we conclude four points on the properties of stationary circular motion. (i) There is a lower limit of velocity. (ii) Although there is no upper limit for the velocity when up to O(v3) terms are considered, there will be an upper limit for the velocity if up to O(v4) terms are considered. (iii) The range of velocity decreases as the chemotaxis increases. (iv) There is a critical value c¯∗ (corresponding critical intensity of the chemotaxis is denoted by kc∗) such that when c¯>c¯∗ (kc>kc∗), the circular motion does not occur for any values of γ.


## 6.4. Simulation of the Full Tensor Model

In this subsection, we examine the effect ofUαβγ on vα and Sαβ. For this, we consider the case of b¯1≠0. The data are obtained by the simulation of the full tensor model.The phase diagrams of the spot motion are shown in Figure8. The phase diagrams in the c¯-γ plane and c¯-κ plane are shown in Figures 8(a) and 8(b), respectively. We note that the behaviors in the full tensor model are similar to that in the reduced model (Figures 5(b) and 5(c)). For positive (negative) b¯1, the region of the circular motion is larger (smaller). The phase diagrams in the b¯1-γ plane and b¯1-κ plane are shown in Figures 8(c) and 8(d), respectively. The behaviors are similar to the phase diagrams in the c¯-γ plane and c¯-κ plane. This result suggests that the effect of ∑γUαβγvγ on the time evolution equation of Sαβ is similar to that of v2Sαβ.Figure 8
Phase diagram of spot motion. Data are obtained by the simulation of the full tensor model.b-s=0.1. M, S, and C in the figure represent regions of no motion, straight motion, and circular motion, respectively. (a) Phase diagram in the c--γ plane. κ=0.8. The solid, dashed, and dotted curves correspond to the cases of b-1=0, 0.25, and −0.25, respectively. The upper (lower) region of each curve corresponds to the region of circular (straight) motion, and the region γ≤0 corresponds to the region of no motion. (b) Phase diagram in the c--κ plane. γ=0.5. The solid, dashed, and dotted curves correspond to the cases of b-1=0, 0.25, and −0.25, respectively. The upper (lower) region of each curve corresponds to the region of straight (circular) motion. (c) Phase diagram in the b-1-γ plane. κ=0.8. The solid, dashed, and dotted curves correspond to the cases of c-=0, 0.25, and −0.25, respectively. The region across each curve is the same as in (a). (d) Phase diagram in the b-1-κ plane. γ=0.5. The solid, dashed, and dotted curves correspond to the cases of c-=0, 0.25, and −0.25, respectively. The region across each curve is the same as in (b).
(a) (b) (c) (d)The phase diagram of spot motion in thec¯-b¯1 plane is shown in Figure 9. It can be seen that the circular motion occurs in the parameter region of positive and large values of c¯ and b¯1. Sαβ and Uαβγ measure the elliptical deformation and head-tail asymmetric deformation of the spot, respectively; these deformations lead to the circular motion.Figure 9
Phase diagram of spot motion in thec--b-1 plane. Data are obtained by the simulation of the full tensor model. b-s=0.1. The marks S and C in the figure represent regions of straight motion and circular motion, respectively. The solid, dashed, and dotted curves correspond to the cases of (γ,κ)=(0.25,0.3), (0.5,0.5), and (1.0,0.8), respectively.The effects ofc¯ and b¯1 on the stationary velocity v and the corresponding radius of circular motion Rc are shown in Figure 10. On comparing these figures, it can be seen that although v and Rc depend on c¯ and b¯1, the influences of b¯1 on v and Rc are much smaller than those of c¯.Figure 10
Dependence of velocity and radius onc- and b-1 in the stationary circular motion. Data are obtained by the simulation of the full tensor model. b-s= 0.1, γ= 0.5, and κ=0.2. (a) Dependence of v on c-. (b) Dependence of Rc on c-, where, in (a) and (b), the solid, dashed, and dotted curves correspond to b-1=0, 0.5, and −0.5, respectively. (c) Dependence of v on b-1. (d) Dependence of Rc on b-1, where, in (c) and (d), the solid, dashed, and dotted curves correspond to c-=0, −0.25, and −0.5, respectively.
(a) (b) (c) (d)The dependence of the critical velocityvc,± and c¯∗ on κ is shown in Figure 11. There exists a critical velocity vc± for b¯1≠0 as shown in Figure 11(a). For larger b¯1, the gap Δvc2 is larger. This suggests that b¯1∑γUαβγvγ enhances the circular motion, which is consistent with the results shown in Figure 8. For positive b¯1, c¯∗ becomes larger than that in the reduced tensor model (solid curve) as shown in Figure 11(b). On the other hand, for negative b¯1, there is no c¯∗. In our system, b¯1 is positive and increases as the intensity of chemotaxis kc increases as shown in Figure 4(d). Thus, c¯∗ still exists in the full tensor model, and c¯∗ is larger than that in the reduced model.Figure 11
Dependence of the critical velocityvc,± and c-∗ on κ. Data are obtained by the simulation of the full tensor model. b-s=0.1. (a) Dependence of vc,±2 on κ.  c-=-0.5. The curves with filled (•) and open (∘) circles correspond to the cases of b-1=0.25 and b-1=0.5, respectively. The solid and dashed curves represent vc,-2 and vc,+2, respectively. (b) Dependence of c-∗ on κ.  b-s=0.1. The solid, dashed, and dotted curves correspond to the cases of b-1=0, 0.25, and 0.5, respectively.
(a) (b)


## 7. Discussions

In this section, we discuss the physical origins of the braking effect observed in the previous section. For our RD system, we derived the time evolution equation ofSαβ up to O(v4) (see (77)). There were three terms of O(v4): b¯s(vαvβ-δαβv2/2)v2, b¯1∑γUαβγvγ, and c¯v2Sαβ. The b¯s term can be absorbed into the b¯(vαvβ-δαβv2/2) term by replacing b¯ with (b¯+b¯sv2). Then, with increasing velocity, the b¯s term changes the value of only b¯. In our theoretical analysis, considering that the influence of Uαβγ is small, we set b¯1=0 and examined the effect of c¯v2Sαβ on the stationary solution in Sections 6.2 and 6.3. When up to O(v3) terms were considered, there was no upper limit of velocity for the circular motion of the spot, as shown in Figure 6(a). However, on including the O(v4) terms, even in the absence of chemotaxis (kc=0), c¯≠0. The term c¯v2Sαβ influenced the upper limit of velocity of the circular motion, as shown in Figure 6(c). One of the physical origins of the upper limit of velocity is the refractory period behind the rear interface. When kc=0, the RD system described by (5) and (6) is considered to have an activator and an inhibitor, as reported by Krischer and Mikhailov [32]. When τ is large, the motionless pulse (localized domain) in one dimension is stable. With decreasing τ, the motionless pulse is destabilized into a traveling pulse in one dimension. For a traveling pulse in one dimension, there is a refractory period behind the rear interface. A repulsive interaction occurs between pulses in the wave train; the repulsion results from the refractory period imposed on the leading interface of a traveling pulse by the tail of the preceding pulse [38]. In two dimensions, when a spot moves along a circle, there is an upper limit of velocity due to the refractory period behind the rear interface; the refractory period causes a braking effect on the velocity of the spot traveling along a circle.In the absence of chemotaxis, that is,kc=0, we estimate the upper limit of velocity. Let us consider the traveling pulse with velocity c in one dimension of the system described by (5) and (6) in the limit ϵ→0. The spatial dependence of the inhibitor before the leading interface and that behind the rear interface are given by vlead(z)=Cleadexp⁡(λ-z) (z>0) and vrear(z)=Crearexp⁡(λ+z) (z<0), respectively, where z=x-ct, λ±=[-c±c2+4β]/2, and Clead(rear) is a coefficient. Using these expressions, we estimate the relaxation length (refractory period) of the inhibitor behind the rear interface by the linear approximation of 1/λ+~c/β for large c. For the spot moving along a circle with a given radius r0, the upper limit of velocity cmax will approximately satisfy the relation 2πr0≫1/λ+~cmax/β: cmax≪2πr0β. This suggests that, for large c, the refractory period is proportional to c. As the velocity of the spot increases, the inhibitor at the leading interfaces becomes large due to the overlap, and, therefore, the velocity becomes small (because of the first term in (16)). These effects on velocity are manifested when the spot moves rapidly along a circle with a small radius.In the presence of chemotaxis, that is,kc>0, we discuss another physical origin of the braking effect. The chemotaxis is caused by the gradient of the chemotactic substance at the interface. In our RD system, the chemotactic substance is autosecretion; it is produced inside the spot and diffuses outward. The chemotactic substance is distributed around the center of the spot and monotonically decreases away from the center. For the traveling pulse, the chemotactic substance at the rear interface is larger than that at the leading interface. The gradient of the chemotactic substance at the leading (rear) interface is negative (positive). The absolute value of the gradient of the chemotactic substance at the leading interface is larger than that at the rear interface. In addition, (dχ/dv) at the leading and rear interfaces are positive. Overall, the chemotactic velocity is negative (because of the second term in (16)). Thus, the chemotaxis in our system essentially reduces velocity.


## 8. Conclusions

In this study, we considered the motion of a spot solution in two dimensions under the influence of chemotaxis. Starting from a three-component RD system, we proposed a two-component (an activator and a chemotactic substance) RD system with a global coupling term. We remark that, in our model system, the spot secretes the chemotactic substance from the inside and the motion of the spot is influenced by the chemotaxis. Thus, the model is an autonomous system. The chemotaxis term is of the Keller-Segel type, and the chemotactic velocity is opposite to the traveling direction. The reason for the opposite direction is that the system involves autosecretion, and the gradient of the chemotactic substance at the leading interface is negative. Although there have been several studies on the motion of spots under the influence of chemotaxis, the chemoattractant was supplied from the outside [20–23]. The spot in these models is driven toward the source (higher concentration) point of the chemoattractant. These models are nonautonomous systems, and, therefore, different from our model.For the RD system, by employing the method proposed in [16], we derived the equation of motion of the spot and the time evolution equations of the tensors. Terms up to the fourth order of v were considered, and we found that the terms v2Sαβ and ∑γUαβγvγ played an important role in the motion of the spot. Our numerical results suggested the existence of an upper limit of velocity for the circular motion of the spot due to the braking effect. There are two physical origins for the braking effect: the refractory period behind the rear interface and the chemotactic velocity opposite to the moving direction. The former is unique to the circular motion and is not applicable for a spot in straight motion. On the other hand, the latter is general in our autonomous chemotactic system and is applicable for all moving directions. Owing to the former reason, an upper limit of velocity for the spot moving in a circle exists even in the absence of chemotaxis. Based on the results of our chemotaxis system, we conclude three points. (i) A spot moving on a straight line is destabilized to circular motion with increasing velocity. (ii) The velocity of a spot moving along a circle is restricted to a certain range, which is narrower than that in the absence of the chemotaxis. (iii) There is a critical intensity of chemotaxis kc∗, and the spot does not undergo circular motion for kc>kc∗.In practical experiments, the candidates of autonomous systems including chemotaxis areE. coli and macrophage [39, 40]. However, there is an essential difference between these biological systems and our model system;E. coli is a rod-shaped cell with a flagellum and macrophage has an amoeba-like shape with pseudopods, and the shape of the motionless state of these cells is not described by a circle.E. coli swims by rotating the flagellum, and macrophage crawls by changing its size and shape. In contrast, in our model system, the spot was a solution of an RD system and the shape of the motionless state was assumed to be a circle. The motionless spot was primarily destabilized through drift bifurcation, and it changed into a moving spot with straight motion. The straight moving spot was secondly destabilized through rotation bifurcation, and it changed into a moving spot with circular motion. Due to these differences, the bifurcation points of the straight and circular motions and the critical velocity for the circular motion in these biological systems must be different from those of our model system. However, the collective motions of bacteria and amoebas are modeled by RD systems, and many spatiotemporal patterns are shown [41, 42]: the diffusion-limited aggregation, Eden-like, uniform disk, branched, target, spiral patterns, and so on. This suggests that the behavior of individual cell in these biological systems can be described by RD systems. From these considerations, the above biological systems are promising candidates of autonomous systems including chemotaxis, and experimental observations of moving cell with circular motion are expected.



---

*Source: 6152961-2018-05-21.xml*