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

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.


Introduction
The behaviors of artificial and biological microswimmers such as oil droplets, bimetallic nanorods, catalytic Janus colloids, liposomes, flagellated bacteria, and Volvox 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][5][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][10][11][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].

Advances in Mathematical Physics
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][25][26][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.

Model Equation
We first consider the following three-component RD system with an activator , a chemotactic substance V, and an inhibitor : where () =  0 + ,  0 , , , , , , , and  are positive constants, and () is a step function satisfying  = 0 for  < 0 and  = 1 for  > 0. Throughout this study, we consider the system in a two-dimensional space, with ∇ = (  ,   ) 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 (, V, ) = (0, 0, −), and the excited state has spatially nonuniform values of , V, and .Between the rest and excited states, there appear boundary layers with thickness (), 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  is small and  and  are large, these bifurcations are suppressed by .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 ( − ()) − , 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 for  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 J d and the chemotaxis term J c .That is, J = J d + J c , where J d = − 2 ∇ and J c = ∇(V).It should be noted that the signs of these fluxes are different.The sign of J c suggests that the chemotaxis term provides a negative diffusion effect, which suppresses the expansion of .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) =    0 (V), where  0 (V) =  0 V 2 /(V 2 +  2 ) with  = 1.0.In order to satisfy the condition max | 0 /V| = 1, we choose  0 = 8 √ 3/9.We call   the intensity of chemotaxis [31].
In (3),  and  represent the relaxation time and diffusion constant of , respectively.Let us consider a situation in which  plays the role of feedback to suppress the static bifurcation and oscillatory bifurcation.For the rapid feedback mechanism,  and  must be small and large, respectively.In the limits  → 0 and  → ∞,  becomes a time-dependent but spatially independent variable, which is denoted by ⟨()⟩: where Ω 0 is the area of the entire system.Replacing  in () by ⟨()⟩, () becomes a global coupling term.In the case where  is very small and  is very large, we reduce the threecomponent RD system to the following two-component RD system with a global coupling term: where  {, V} = − +  ( −  {, V}) .
The functional {, V} represents a global coupling term given by where the integral is over the entire domain, and  and  are rescaled to absorb Ω 0 :   = /Ω 0 and   = Ω 0 , and the primes are dropped. corresponds to the intensity of the global coupling, and the value of  0 is chosen as  0 = 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  becomes an interface, as shown in Figure 1.The location of the interface is defined by the condition (r, ) = .In this limit,  and V satisfy the relation ( + 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 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.
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   for the bifurcation of stationary straight motion.Although the straight motion was stable for  >   , it was destabilized by rotation (straight-circular-motion) bifurcation for  ≤   .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 where  = 1 + .Setting V/ = 0 in (10), (10) becomes where  = |r| and  0 is the equilibrium radius of the symmetric spot.We impose the boundary conditions on V() as The radially symmetric equilibrium solution of (10) in two dimensions is given by where   and   are the modified Bessel functions of the first and second kinds of order , respectively.The corresponding Advances in Mathematical Physics solution of () is given by () = ( 0 − ) − V().The above equilibrium solution is employed to study the motion of the spot in the following sections.

Equation of Motion of the Interface
As  is small, the width of the boundary layer is small, and, therefore, the value of  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 where (ℎ) 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 ℎ is the value of V on the interface. is the normal component of the velocity, and  is the Lagrange multiplier for the constraint of domain area conservation (9).For the moving spot, this condition corresponds to where  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 (ℎ) is obtained as where  is a global coupling term (8).(ℎ) 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 along n is given by where   0 (V) =  0 (V)/V and |  represents the value of function  evaluated on the interface.n is a unit normal vector on the interface, and (V/n) 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/), 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  as where  is Green's function satisfying the equation In (18),  represents the integral 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  as in (18).

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 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   , which is given as where e  is a radial unit vector and (  ) is the distance between the center of gravity and a point on the interface, which is directed to an angle   with respect to the  axis.
Using the definition of R(),  is given by where   = (  )/  .For the deformed spot, the radius is given by with where i = √ −1,  0 is the equilibrium radius of the symmetric spot, and  corresponds to the deviation.In the expansion of ,  ± corresponds to a (2/)-periodic deformation.
In addition, we exclude the terms of  = ±1 because the translational motion of the spot is incorporated in .Thus, we set  ±1 = 0. Using (24), the curvature  and the normal component of the velocity  are given up to the first order of the deviations as respectively, where the overdot denotes the time derivative.For a general function (r, ), the Fourier transformation and its reverse transformation in two dimensions are defined as respectively, where ∫ q = ∫q/(2) 2 for brevity.For an isolated domain that forms a single loop without crossing, the step function ( − ) = ( − |r − |) is transformed by the Fourier transformation as Substituting ( 24) into (30), we obtain the expansion of  q up to the first order of  as where  = |q|,   is the angle between q and the -axis, and   is a Bessel function of the first kind of order .

Equation of Motion of a Spot.
We derive the equation of motion of a spot from (14). is given as functions of ℎ and h by (17), in which we used the notation h = (V/n)|  = (n ⋅ ∇)V|  .When the velocity of the spot is small, the deformation of the circle is small.For this situation, we derive the deviations of ℎ and h from their stationary values as a power series of .Using (18), we first expand ℎ and h up to the fourth-order time derivatives as respectively, where with In the above expansions, we omitted the term proportional to (q ⋅ v)( q /) in the expressions of ℎ 2 and h2 .In the derivation of ℎ 3 and h3 , the terms proportional to ( 3  q / 3 ), (q⋅v)( 2  q / 2 ), (q⋅v) 2 ( q /), (q⋅ v )( q /), (q ⋅ v ) q , and (q ⋅ v)(q ⋅ v ) q were omitted.In the derivation of ℎ 4 and h4 , all the terms of the derivatives of v and  q 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/Ω) ∫ R()⋅ due to the orthogonality relation of trigonometric functions.We remark that we include up to the third-order of v: ℎ 3 and h3 , 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: ℎ 4 and h4 .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 radius  0 .The values of V and (n ⋅ ∇)V at the interface are denoted as ℎ (0) 0 and h(0) 0 , respectively.Using (34), ℎ (0) 0 is given by ℎ 0 with the substitutions  q =  (0) q and R = R (0) , where R (0) =  0 e  .Similarly, h(0) 0 is given using (35).When  is small but finite, substituting ℎ = ℎ (0) 0 and h = h(0) 0 into (ℎ) with  = 0, (14) becomes where  = −1/ 0 is used.We note that the Lagrange multiplier  in (37) can be absorbed into the constant term  0 in (ℎ (0) 0 ).In the limit  → 0, we can calculate ℎ (0) 0 and h(0) 0 by using V(), which is given by ( 13) (the validation is given in Appendix B).Using these expressions, (37) gives the dependence of  0 on   and .The numerical results obtained by using ( 13) and ( 17) are shown in Figure 2(a).When  = 0,  0 monotonically decreases with   .However, for large ,  0 depends weakly on   and it is almost constant because of the global coupling; large global coupling suppresses variations in the area of the spot resulting in constant  0 .The dependence of  0 on  when   is large is shown in Figure 2(b).From (9), it can be seen that, for large ,  0 approximately satisfies  2 0 ∼ , and, therefore,  0 is proportional to the square root of  (see the dashed curve in Figure 2(b)).
When the spot moves with an infinitesimal velocity, ℎ and h deviate from ℎ (0) 0 and h(0) 0 , respectively.By putting small deviations as ℎ = ℎ − ℎ (0) 0 and  h = h − h(0) 0 , we iteratively derive ℎ and  h in a power series of  up to the third order.The procedure consists of two steps.We first expand  given by ( 17) in a power series of ℎ and  h.Then, we solve the equation in terms of ℎ and expand it by .When the stationary motionless spot is destabilized into a moving spot,  is expanded up to the second-order of ℎ and  h as where, for a function (ℎ), () 0 implies (ℎ (0) 0 ).Up to the first order of ℎ and  h, (38) results in  = (/ℎ + /ℎ h(0) 0 ) 0 ℎ + () 0  h.For higher order corrections, a supplementary relation between ℎ and  h is necessary.In order to relate ℎ 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|) denoted by V().For this function, (n ⋅ ∇)V = V/.Then, V( 0 + ) and (n ⋅ ∇)V( 0 + ) are expanded around V( 0 ) and V( 0 )/ for small  as respectively.In ( 39) and ( 40), (()) = 0 implies ( 0 ).We make an ansatz that ℎ and  h are not independent but have a linear relation such that  h = ℎ with  = (V  /V  ) = 0 , where the prime corresponds to the derivative with respect to .This assumption enables us to calculate ℎ as a function of .Substituting  h = ℎ into (38),  becomes a quadratic equation of ℎ as where As a result of numerical calculations, we remark that  > 0 and  < 0 for small ℎ (0) 0 , but  > 0 for large ℎ (0) 0 .By solving (41) in terms of ℎ for small , the solution is expanded up to the second order of  as Here, it should be noted that although there are two solutions for (41), we chose one solution such that ℎ → /( + ) when the second-order term () 2 is neglected.
In order to obtain ℎ up to the third order of ,  given by ( 17) is expanded up to the third order of ℎ and  h.Using  h = ℎ,  becomes a cubic equation of ℎ as where  2 is given by In order to obtain the third-order term of  in ℎ, we add a correction Δ (∼ () 3 ) to ℎ given by (43) and substitute it into (44).On solving this equation in terms of Δ, ℎ up to the third order of  is obtained as 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 , we obtain the same result as that of (46).In the above process, we iteratively derived ℎ up to the third order of .We replace ℎ ∼  h in (46), and finally the power series expansion of ℎ +  h in terms of  is obtained as where   ( = 1, 2, 3) is defined by In the absence of chemotaxis ( = 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/Ω) ∫ R().For the left hand side, we put where h 4 and h4 are neglected because these terms are higher than or equal to (k 4 ).For the calculation of the right hand side of (47), we use (14).The magnitudes of  and  were discussed in [16];  is independent of   , and, owing to the periodicity of the function, ∫ R() = ∫ R() 3 = ∫ R() 2  = 0.In addition,  ∼ (k 2 ), and, therefore, ∫ R() 3 can be neglected up to the third order of v.
As the translational motion of the spot is incorporated in , we chose  ±1 = 0 in the expansion of (  ).This results in  ∫ R() = 0 up to the first order of the deformation.The terms  ∫ R() 2 and  ∫ R() 2 , and the higher order terms of  were small; therefore, they were neglected.Following the above discussion, we neglect the term ( + ) in (14).By carrying the integral over   , the following equation is obtained: where Here, h 1,1 = (1/Ω) ∫ R()ℎ 1,1 .The other terms on the right hand side of (51) are defined similarly.In the expansion Advances in Mathematical Physics (51), we neglected the terms h 1,2 and h 2,3 , because these terms yielded only the terms ċ±1 and c±1 , respectively; however, these terms are set to zero.h 2,2 disappeared in the integral over   owing to the orthogonality relation of trigonometric functions [15].Owing to the same reasons, the terms h1,2 , h2,3 , and h2,2 did not exist.In the integral over   , there is no contribution from the second-order term of  owing to the orthogonality relation of trigonometric function.This supports the justification to include terms up to the third order of  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 of v in h and  h.
After practical calculations of h and  h, we finally obtain the equation of motion of the spot as where h (0) 1,1 is defined as Here, h (0) 1,1 is given by h 1,1 with the substitutions  q =  (0) q and R = R (0) . h(0) is similarly defined with h (0) 1,1 .The coefficients ,   , and  on the left hand side of (52) are given by respectively, where we defined (, , , ) and   (, , , ) as respectively, where , , , and  are integers and   ( 0 ) is written as   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 , damping coefficient   , intensity of the cubic nonlinear term , and the right hand side of (52) form the coupling term between deformation and velocity in the following subsection.When ( −   ) > 0, the motionless spot is stable; however, the motionless spot is destabilized into a moving spot for ( −   ) < 0. When  is positive, the third-order nonlinear coupling term v|v| 2 suppresses the divergence for large |v|.It is necessary that ,   , and  are positive for the existence of a stable moving spot.The dependence of ,   , and  on   and  is shown in Figure 3.We see that these parameters are positive and monotonically increase with   and .
In the derivation of (47), we made an ansatz for the relation between ℎ and  h such that  h = ℎ with  = (V  /V  ) = 0 .Equation (55) suggests that  depends explicitly on  through  3 .Through preliminary research, it is found that if  is chosen as a positive constant, the property that  is a positive and monotonically increasing function of   and  holds.From the above, it is seen that ,   , and  are positive for any value of   and , and, therefore, we fix  and examine the motion of the spot in the range   ≥ 0.

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 where ℎ (0) 1,1() and  h(0) 1,1() correspond to the  component of h (0)  1,1 and  h(0) 1,1 , 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 (  ) is expanded using the coefficients   as given by (25);  ± corresponds to a (2/)-periodic deformation.
We first introduce a second-rank tensor   (,  = 1, 2) using  ±2 as follows: where  2 is a positive constant, which represents the radial deviation of the spot from  0 , and  2 is the angle between the long axis of the ellipse and the -axis.The tensor elements (59) are represented by using a normal vector N = (cos  2 , sin  2 ) along the long axis as The tensor S is the same as the nematic order parameter tensor in liquid crystals [34].For an elliptical spot, (  ) is represented as Next, in order to describe the head-tail asymmetric deformation, we first define  1 and  2 using  ±3 as respectively, where  3 is a positive constant and  3 is the angle between one of the long axes of the deformed spot and the axis.In order to relate  1 and  2 with a tensor, we introduce the third-rank tensor   (, , and  = 1 or 2).The tensor elements   are represented by using vectors N () ( = 1, 2, and 3) as where the normal vectors N () are defined by N (1) = (cos  3 , sin  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   and  1 and  2 as  111 =  1 ,  222 = − 2 , and The spot with head-tail asymmetric deformation is represented as The terms ℎ (0) 1,1() and  h(0) 1,1() in (58) are expressed using tensors.The detailed calculations are given in Appendix C. The final form is as follows: where Thus, (58) is written using   in the form where  = ( * + ã * ).ℎ (0) 1,1() and  h(0) 1,1() yield the Sv term.This is due to the periodicity of the function in the integral over   ; in the expansion of R(  ) and  (1)  q (  ), only  ±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   > , and the velocity of the spot is approximately given by |v| 2 ∼ (  −)/.The deformation and velocity coupling term Sv modifies the straight motion of the spot.

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 timedependent tensor   .In this subsection, we derive the time evolution equations of the tensors up to ∼ (k 4 ).We first discuss the order of v,   , and   , following [16].From (73), the motionless spot is critical at  =   , and the moving spot occurs supercritically with increasing (  − ).We put  =   −  for the smallness parameter.The time is scaled by t = , an all the terms in (73) are of the order ( 3/2 ).
In order to derive the equation of motion of tensors, we calculate the first-order time derivative of   .The linear combination of ċ results in the equation of motion of tensors.In this process, v terms in ℎ 2,1 and  h2,1 are neglected because these terms cause small order terms; ∑    V  and   V  appear in the time evolution equation of   and   , respectively.These terms are ∑    V  ∼ (k 6 ) and   V  ∼ (k 5 ) so that they do not contribute to the equation of motion of tensors.In addition, ( 2  q / 2 ) terms in ℎ 2,3 and  h2,3 are neglected because these terms result in the second-order time derivative of   ; these terms are c±2 ∼ (k 6 ) and c±3 ∼ (k 7 ) and they are small enough to be neglected.
We first derive the equation of motion of   .Using (74), we can obtain the time evolution equation of  ±2 , and we derive the time evolution equation of   .The detailed derivation is given in Appendices D and E. The time evolution equation of   up to (k 4 ) is where , and  = −(2 2,2 −  2 ) − (2 T2,2 − Ṽ2 ).All the parameters in these expressions are given in Appendix E. The second term in (77) originates from ℎ 2 and  h2 , the third term from ℎ 4 and  h4 , the fourth term from ℎ 1 and  h1 , and the fifth term from ℎ 2 and  h2 .Thus, in order to consider the terms up to (k 4 ), ℎ 4 and  h4 are necessary in the expansion equations ( 75) and (76).From the first and second terms in (77), we note that   ∼ (k 2 ).
Next, we consider the time evolution equation of   .Following procedures similar to those for   , the time evolution equation of   up to (k 4 ) is where Γ 3 =  +  3 +  Ẽ3 and  3 = 8/ 2 0 +  3 +  D3 . 1 and  2 are defined as respectively.In the above equation, all the terms are ∼ (k 3 ).The terms of (k 4 ) disappeared because of the orthogonality relation of trigonometric functions.The next higher order term is (k 5 ), which is not included in expansions ℎ and  h in (75) and (76).The second term in (78) originates from ℎ 3 and  h3 , and the third term from ℎ 1 and  h1 .From the first and second terms in (78), we note that   ∼ (k 3 ).
In the above full tensor model, we chose parameters such that / = 1 and put  = (  − )/ and   = / in (81).The other rescaled parameters in (82) and (83) were  =  2 /Γ 2 , .Henceforth, we drop the primes on the parameters in the full tensor model ( 81)-(83).In the full tensor model, when  is large and v is small, the effect of   on   is small in (82).In this case, we put  1 = 0, and the V  ,   , and   system described by ( 81)-( 83) is reduced to a V  and   system.We call this system the reduced tensor model.

Stationary Solution and Phase Diagrams
In the following subsections, we consider the stationary solution and phase diagrams.For this, we rewrite V  ,   , and   by introducing the following variables: where we choose V, ,  > 0. From (81), we obtain the time evolution equations of V and  as respectively.From (82), we obtain the time evolution equations of  and  as Advances in Mathematical Physics respectively.From (83), we obtain the time evolution equations of  and  as 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 Δ = 1.0 × 10 −4 and fix the parameters as  =  = −1,  = 0.3,  1 = 0.1, and  2 = 0.5.

Stationary Solution.
For the reduced and full tensor models, we consider the stationary solution of the moving spot, which moves straight in the  direction.This situation corresponds to  = 0.In the case of  1 =   =  = 0, the stationary solution of ( 85) and (86) becomes V 2 0 = /(1 + ) with  = /2.In order to avoid the divergence of V 2 0 , we assume that  > 0. For  ̸ = 0, we consider the dependence of the stationary solution of , , , and  on the signs of  and .In the stationary state of (86) with  = 0, there are two stationary solutions of  depending on the sign of  and : (I)  and  are positive with  0 =  and (II)  and  are negative with  0 = ( + 1/2), where  is an integer.Using (87)-( 90), the stationary solution of case (I) is where  0 satisfies sin 3 0 = 0 and V 0 is a stationary solution of V, which is determined later.On the other hand, for case (II), the stationary solution is From these results, it is seen that  0 has the same expression irrespective of the signs of  and .As  0 is positive,  0 satisfying sin 3 0 = 0 is chosen as For both cases ((I) and (II)), the stationary solution of  0 is given by In both cases ((I) and (II)), by substituting  0 into (85) in the stationary state, we obtain two stationary solutions V 0,± as In order to satisfy the condition that V 2 0,± → /(1 + ) in the limits  1 ,   , and  → 0, only V 2 0,− 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].

Stability of Straight Motion in the Reduced Tensor Model.
When  is large in (89),  is rapidly damped to zero.In addition, when V is small, the influences of  on  and  in (87) and (88) are small.For this case, we discuss the stability of the straight motion of the spot by setting  1 = 0 in (87) and (88).Using ( 86) and (88), we can derive the time evolution equation of  =  −  as For the stationary states V = V 0,− ,  =  0 , we discuss the stability condition of the straight motion.From (96), the stationary solution  0 satisfies 2 0 = 2 or (2 + 1), and their corresponding stability conditions are (− 0 + V 2 0,− / 0 +   V 4 0,− / 0 ) > 0 or (− 0 + V 2 0,− / 0 +   V 4 0,− / 0 ) < 0, respectively.Since  and  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  0 given by (92).The bifurcation of straight-circular motion, when   =  = 0, was already studied in [16], and the critical value for the straight motion is expressed as 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 <  ≤   , the spot moves straight.However, when  >   , the straight motion become unstable and changes into circular motion.On the other hand, when   ̸ = 0 and  ̸ = 0, the stability condition for straight motion is The phase diagram is a curved surface in -- 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   =  = 0 (the solid curve), it can be seen that when  > 0, the parameter region for the circular motion is larger.However, when  < 0, the parameter region for the circular motion is smaller.We can explain this result using (87); the term V 2  essentially reduces the parameter  of the damping term − such that the damping term becomes −( − V 2 ).When  is positive (negative), the term V 2  effectively reduces (increases)  so that the deformation  becomes large (small).This leads to positive (negative)  results in the larger parameter region of the circular motion.In our RD system,  is negative and || monotonically increases with   (see Figure 4(c)), and the parameter region for circular motion becomes smaller as the intensity of chemotaxis increases.The phase diagrams in the - plane and - 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  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.

Dependence of Critical Velocity on Parameters.
In the previous subsection, we showed that a spot appears in the circular motion when  >   .In this parameter region, we examine the stationary circular motion of the spot with a constant angular frequency  and velocity V  .For this, we set  1 = 0 and substitute V = V  ,  =   ,  =  + /2, and  =  into (85)-(88), and obtain the relations among ,   , V  , and .The calculation is straightforward; the final expressions are shown.
We first consider the case when   =  = 0, where up to (k 3 ) terms are considered, and obtain the relations cos ), and V 2  =  2 /.When the spot rotates with an angular frequency ,  satisfies  2 ≥ 0, and this condition leads to V  ≥ V  [16].The dependence of the critical velocity V  for the stable angular frequency on  (≥ 0) is shown in Figure 6(a).In this case, there is a lower limit of velocity V  for the circular motion but no upper limit of velocity for the same.
Next, we consider the case when   ̸ = 0 and  ̸ = 0, where up to (k 4 ) terms are considered.After similar calculations to the case of   =  = 0, we obtain the relations cos In this case,  2 is a quartic function of V  .For the condition  2 ≥ 0, there are two cases depending on the sign of (  − 2 ).When (  − 2 ) > 0,  2 exists for V 2  ≥ V 2 ,+ , where V 2 ,± are the solutions of  2 = 0, which are given by V ,+ on  is shown in Figure 6(b); there is only a lower limit of velocity V ,+ .This property is the same with the case when   =  = 0. On the other hand, when (  −  2 ) < 0,  2 exists in a certain range of V 2  ; V 2 ,+ ≤ V 2  ≤ V 2 ,− under the condition   ≥ 0. In this case, there is a lower limit and an upper limit of velocity, V ,+ and V ,− , respectively.The gap between V 2 ,− and V 2 ,+ is defined by Let us consider two limiting cases.In the limit || → ∞, V 2 ,± and ΔV 2  → 0. That is, the circular motion cannot exist.On the other hand, for small || and in terms of  and   and obtain V 2 ,+ ∼  2 /( + 2) and V 2 ,− ∼ ( + 2)/( 2 −   ).Thus, in the limits   and  → 0, V 2 ,+ →  2 /, V 2 ,− → ∞, and ΔV 2  → ∞.This suggests that V ,+ → V  and there is no upper limit of velocity, which are the same with the case when   =  = 0. From the numerical results shown in Figures 4(a), 4(b), and 4(c), (  −  2 ) < 0 in our RD system; therefore, it turns out that V 2  exists only in a finite range, decreases as the chemotaxis increases.When the condition   = 0 is satisfied, V ,− = V ,+ and ΔV 2  = 0, which yields the critical value  * as When || > | * |,   < 0, and, therefore, V 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  * on  is shown in Figure 6(e).Equation (99) can be positive for large .However, when  is chosen as a positive value, the velocities V and  increase monotonically with increasing , and, therefore, the system diverges in the simulation.For this reason,  * is shown in the range  * ≤ 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 Figures 6(a) and 6(b), for given  and , there is only a lower limit of  for the circular motion of the spot, resulting in V  and V ,+ , respectively.On the other hand, in Figure 6(c), for given  and , there is not only a lower limit of  but also an upper limit of  for the circular motion of the spot, resulting in V ,+ and V ,− , respectively.However, in Figure 6(e), for given  and  * , 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 on  are shown in Figure 7.The velocity V increases monotonically with increasing , as shown in Figure 7(a).The corresponding radius of the spot motion   is shown in Figure 7(b).For fixed values of  and , V decreases and   increases as || increases in the range  ≤ 0. The effect of chemotaxis is incorporated in , and  is negative; || increases as the intensity of chemotaxis   increases in our RD system (Figure 4(c)).Thus, it is seen that the chemotaxis reduces V and increases   .
From the above analyses on (85)-( 88) with  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 (k 3 ) terms are considered, there will be an upper limit for the velocity if up to (k 4 ) terms are considered.(iii) The range of velocity decreases as the chemotaxis increases.(iv) There is a critical value  * (corresponding critical intensity of the chemotaxis is denoted by  *  ) such that when || > | * | (  >  *  ), the circular motion does not occur for any values of .

Simulation of the Full Tensor Model.
In this subsection, we examine the effect of   on V  and   .For this, we consider the case of  1 ̸ = 0.The data are obtained by the simulation of the full tensor model.
The phase diagrams of the spot motion are shown in Figure 8.The phase diagrams in the - plane and - 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)  1 , the region of the circular motion is larger (smaller).The phase diagrams in the  1 - plane and  1 - plane are shown in Figures 8(c) and 8(d), respectively.The behaviors are similar to the phase diagrams in the - plane and - plane.This result suggests that the effect of ∑    V  on the time evolution equation of   is similar to that of |k| 2   .
The phase diagram of spot motion in the - 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  and  1 .  and   measure the elliptical deformation and head-tail asymmetric deformation of the spot, respectively; these deformations lead to the circular motion.
The effects of  and  1 on the stationary velocity V and the corresponding radius of circular motion   are shown in Figure 10.On comparing these figures, it can be seen that although V and   depend on  and  1 , the influences of  1 on V and   are much smaller than those of .
The dependence of the critical velocity V ,± and  * on  is shown in Figure 11.There exists a critical velocity V ± for  1 ̸ = 0 as shown in Figure 11(a).For larger  1 , the gap ΔV 2  is larger.This suggests that  1 ∑    V  enhances the circular motion, which is consistent with the results shown in Figure 8.For positive  1 , | * | becomes larger than that in the reduced tensor model (solid curve) as shown in Figure 11(b).On the other hand, for negative  1 , there is no  * .In our system,  1 is positive and increases as the intensity of chemotaxis   increases as shown in Figure 4(d).Thus,  * still exists in the full tensor model, and | * | is larger than that in the reduced model.

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 of   up to (k 4 ) (see (77)).There were three terms of (k 4 ):   (V  V  −   |k| 2 /2)|k| 2 ,  1 ∑    V  , and |k| 2   .The   term can be absorbed into the (V  V  −   |k| 2 /2) term by replacing  with ( +   |k| 2 ).Then, with increasing velocity, the   term changes the value of only .In our theoretical analysis, considering that the influence of   is small, we set  1 = 0 and examined the effect of |k| 2   on the stationary solution in Sections 6.2 and 6.3.When up to (k 3 ) terms were considered, there was no upper limit of velocity for Advances in Mathematical Physics the circular motion of the spot, as shown in Figure 6(a).However, on including the (k 4 ) terms, even in the absence of chemotaxis (  = 0),  ̸ = 0.The term |k| 2   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   = 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,   = 0, we estimate the upper limit of velocity.Let us consider the traveling pulse with velocity  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 V lead () =  lead exp( − ) ( > 0) and V rear () =  rear exp( + ) ( < 0), respectively, where  =  − ,  ± = [− ± √ 2 + 4]/2, and  lead(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/ + ∼ / for large .For the spot moving along a circle with a given radius  0 , the upper limit of velocity  max will approximately satisfy the relation 2 0 ≫ 1/ + ∼  max /:  max ≪ 2 0 .This suggests that, for large , the refractory period is proportional to .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,   > 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, (/V) 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.

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 twocomponent (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][21][22][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 |k| 2   and ∑    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  respectively.We first derive ℎ (0) 0 from the definition (34).Substituting  q =  (0) q and R = R (0) into (34), we obtain )  0 ( 0 )  1 ( 0 ) .Next, we derive h(0) 0 from the definition (35).Substituting  q =  (0) q and R = R (0) into (35), we obtain h(0)

D. Expansions of 𝛿ℎ and 𝛿 h in (75) and (76)
The expansions of ℎ and  h in (75) and ( 76 where V ± = V 1 ± iV 2 and all the coefficients in (D.1) and (D.2) are given in Appendix E. In order to calculate the time

Figure 2 :
Figure 2: Dependence of  0 on   and  in the stationary state.(a) Dependence of  0 on   . = 1.0.The solid and dashed curves represent the cases of  = 0 and 5, respectively.(b) Dependence of  0 on .  = 5.0.The solid and dashed curves represent the cases of  = 0.1 and 5.0, respectively.

Figure 3 :
Figure 3: Dependence of ,   , and  on   and .The solid, dashed, and dotted curves represent ,   , and , respectively.(a) Dependence of ,   , and  on   . = 5.0 and  = 1.0.In the calculation of ,  is chosen as  =   .As  3 is positive,  for the case  < (>)   is smaller (larger) than this curve.(b) Dependence of ,   , and  on . = 5.0 and   = 5.0.In the calculation of ,  is chosen as  =   .

Figure 5 :
Figure 5: Phase diagram of spot motion. =  = −1 and   = 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  = 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 - 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 - 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.

Figure 8 :
Figure 8: Phase diagram of spot motion.Data are obtained by the simulation of the full tensor model.  = 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 - plane. = 0.8.The solid, dashed, and dotted curves correspond to the cases of  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 - plane. = 0.5.The solid, dashed, and dotted curves correspond to the cases of  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  1 - plane. = 0.8.The solid, dashed, and dotted curves correspond to the cases of  = 0, 0.25, and −0.25, respectively.The region across each curve is the same as in (a).(d) Phase diagram in the  1 - plane. = 0.5.The solid, dashed, and dotted curves correspond to the cases of  = 0, 0.25, and −0.25, respectively.The region across each curve is the same as in (b).

Figure 9 :
Figure 9: Phase diagram of spot motion in the - 1 plane.Data are obtained by the simulation of the full tensor model.  = 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.

Figure 10 :Figure 11 :
Figure 10: Dependence of velocity and radius on  and  1 in the stationary circular motion.Data are obtained by the simulation of the full tensor model.  = 0.1,  = 0.5, and  = 0.2.(a) Dependence of V on .(b) Dependence of   on , where, in (a) and (b), the solid, dashed, and dotted curves correspond to  1 = 0, 0.5, and −0.5, respectively.(c) Dependence of V on  1 .(d) Dependence of   on  1 , where, in (c) and (d), the solid, dashed, and dotted curves correspond to  = 0, −0.25, and −0.5, respectively.