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.
1. 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–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 activator , a chemotactic substance , and an inhibitor :where , , , , , , , , and are positive constants, and is a step function satisfying for and for . Throughout this study, we consider the system in a two-dimensional space, with and . We choose such that the system is monostable. Here, we fix . In the above excitable system, there are two stationary states: a rest state and an excited state. The rest state is , and the excited state has spatially nonuniform values of , , 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 , where the flux is the sum of the normal diffusion (random motility) term and the chemotaxis term . That is, , where and . It should be noted that the signs of these fluxes are different. The sign of 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 , where with . In order to satisfy the condition , we choose . 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 and , becomes a time-dependent but spatially independent variable, which is denoted by :where 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 three-component RD system to the following two-component RD system with a global coupling term:whereThe functional represents a global coupling term given bywhere the integral is over the entire domain, and and are rescaled to absorb and , and the primes are dropped. corresponds to the intensity of the global coupling, and the value of is chosen as . 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 . In the limit , the boundary layer of becomes an interface, as shown in Figure 1. The location of the interface is defined by the condition . In this limit, and satisfy the relation (inside the domain) and 0 (outside the domain), which is obtained from (5). In the limits of and , the area of the spot is conserved; using (8), we obtainThis 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 . In this limit, (6) becomeswhere . Setting in (10), (10) becomeswhere and is the equilibrium radius of the symmetric spot. We impose the boundary conditions on asThe radially symmetric equilibrium solution of (10) in two dimensions is given bywhere and are the modified Bessel functions of the first and second kinds of order , respectively. The corresponding solution of is given by . 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 changes sharply while that of changes smoothly around the boundary layer. The velocity of the flat interface is a function of the value of 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] bywhere 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 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 towhere 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 aswhere 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 is given bywhere and represents the value of function evaluated on the interface. is a unit normal vector on the interface, and 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), , is small. In this case, we deal with the time derivative of in (10) as a perturbation. The asymptotic solution of (10) is written by perturbation expansion using the Green function aswhere is Green’s function satisfying the equationIn (18), represents the integralfor 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 in terms of 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 bywhere is the area of the spot and is a position vector measured from the center of gravity. We consider the case that is a single valued function of , which is given aswhere 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 , is given bywhere . For the deformed spot, the radius is given bywithwhere , is the equilibrium radius of the symmetric spot, and corresponds to the deviation. In the expansion of , corresponds to a -periodic deformation. According to the area conservation equation (9), . In addition, we exclude the terms of because the translational motion of the spot is incorporated in . Thus, we set . Using (24), the curvature and the normal component of the velocity are given up to the first order of the deviations asrespectively, where the overdot denotes the time derivative.
For a general function , the Fourier transformation and its reverse transformation in two dimensions are defined asrespectively, where for brevity. For an isolated domain that forms a single loop without crossing, the step function is transformed by the Fourier transformation asSubstituting (24) into (30), we obtain the expansion of up to the first order of aswhere , is the angle between and the -axis, and is a Bessel function of the first kind of order .
4.2. Equation of Motion of a Spot
We derive the equation of motion of a spot from (14). is given as functions of and by (17), in which we used the notation . When the velocity of the spot is small, the deformation of the circle is small. For this situation, we derive the deviations of and from their stationary values as a power series of . Using (18), we first expand and up to the fourth-order time derivatives asrespectively, wherewith
In the above expansions, we omitted the term proportional to in the expressions of and . In the derivation of and , the terms proportional to , , , , , and were omitted. In the derivation of and , all the terms of the derivatives of and 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 due to the orthogonality relation of trigonometric functions. We remark that we include up to the third-order of : and , 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 : and . 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 . The values of and at the interface are denoted as and , respectively. Using (34), is given by with the substitutions and , where . Similarly, is given using (35). When is small but finite, substituting and into with , (14) becomeswhere is used. We note that the Lagrange multiplier in (37) can be absorbed into the constant term in . In the limit , we can calculate and by using , which is given by (13) (the validation is given in Appendix B). Using these expressions, (37) gives the dependence of on and . The numerical results obtained by using (13) and (17) are shown in Figure 2(a). When , monotonically decreases with . However, for large , 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 . The dependence of on when is large is shown in Figure 2(b). From (9), it can be seen that, for large , approximately satisfies , and, therefore, is proportional to the square root of (see the dashed curve in Figure 2(b)).
(a)
(b)
When the spot moves with an infinitesimal velocity, and deviate from and , respectively. By putting small deviations as and , we iteratively derive and 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 . 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 aswhere, for a function , implies . Up to the first order of and , (38) results in . For higher order corrections, a supplementary relation between and is necessary. In order to relate with , we consider the profile of in the radially symmetric function given by (13), where is a function of () denoted by . For this function, . Then, and are expanded around and for small asrespectively. In (39) and (40), implies . We make an ansatz that and are not independent but have a linear relation such that with , where the prime corresponds to the derivative with respect to . This assumption enables us to calculate as a function of . Substituting into (38), becomes a quadratic equation of aswhereAs a result of numerical calculations, we remark that and for small , but for large . By solving (41) in terms of for small , the solution is expanded up to the second order of asHere, it should be noted that although there are two solutions for (41), we chose one solution such that when the second-order term is neglected.
In order to obtain up to the third order of , given by (17) is expanded up to the third order of and . Using , becomes a cubic equation of aswhere is given byIn order to obtain the third-order term of in , we add a correction ) to given by (43) and substitute it into (44). On solving this equation in terms of , up to the third order of is obtained asWe 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 in (46), and finally the power series expansion of in terms of is obtained aswhere () is defined byIn the absence of chemotaxis (), (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 . For the left hand side, we putwhere and are neglected because these terms are higher than or equal to . 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, . In addition, , and, therefore, can be neglected up to the third order of v. As the translational motion of the spot is incorporated in , we chose in the expansion of . This results in up to the first order of the deformation. The terms and , 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:whereHere, . The other terms on the right hand side of (51) are defined similarly. In the expansion (51), we neglected the terms and , because these terms yielded only the terms and , respectively; however, these terms are set to zero. disappeared in the integral over owing to the orthogonality relation of trigonometric functions [15]. Owing to the same reasons, the terms , , and 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 and .
After practical calculations of and , we finally obtain the equation of motion of the spot aswhere is defined as . Here, is given by with the substitutions and . is similarly defined with . The coefficients , , and on the left hand side of (52) are given byrespectively, where we defined and asrespectively, where , and are integers and 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 , the motionless spot is stable; however, the motionless spot is destabilized into a moving spot for . When is positive, the third-order nonlinear coupling term suppresses the divergence for large . 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 .
(a)
(b)
In the derivation of (47), we made an ansatz for the relation between and such that with . Equation (55) suggests that depends explicitly on through . 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 .
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 , the component () of (52) becomeswhere and correspond to the component of and , respectively. Using the tensors, the right hand side of (58) is expressed in a simple form. The second rank tensor and the third rank tensor represent the elliptical deformation and the head-tail asymmetric deformation (-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 -periodic deformation.
We first introduce a second-rank tensor (, ) using as follows:where is a positive constant, which represents the radial deviation of the spot from , and is the angle between the long axis of the ellipse and the -axis. The tensor elements (59) are represented by using a normal vector along the long axis asThe tensor 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 and using asrespectively, where is a positive constant and is the angle between one of the long axes of the deformed spot and the -axis. In order to relate and with a tensor, we introduce the third-rank tensor (, , and or 2). The tensor elements are represented by using vectors (, and 3) aswhere the normal vectors are defined byThe tensor 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 and as , , andThe spot with head-tail asymmetric deformation is represented as
The terms and in (58) are expressed using tensors. The detailed calculations are given in Appendix C. The final form is as follows:whereThus, (58) is written using in the formwhere . and yield the term. This is due to the periodicity of the function in the integral over ; in the expansion of and , only terms contribute to the nonzero integral, resulting in the 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 . The deformation and velocity coupling term 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 tensor . In this subsection, we derive the time evolution equations of the tensors up to . 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 , an all the terms in (73) are of the order . Here and ; therefore . In the later calculations, we can confirm that and . In the derivation of (73), we omitted the terms that include higher order derivatives of time. We estimate these terms, for example, , , , , , , and in the expansion of and . From these estimates, the omission of those terms is justified.
and are composed of and , respectively, so that and . The coupling terms such as () are much smaller. Therefore, we linearize (14) in terms of the deformation . Then, the time evolution equation of becomesWe consider terms up to the fourth order of , where and are given byrespectively. Here, and () are defined by (34) as and , respectively. is similarly defined by (35). We expand , , , and as , , , and , respectively.
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, terms in and are neglected because these terms cause small order terms; and appear in the time evolution equation of and , respectively. These terms are and so that they do not contribute to the equation of motion of tensors. In addition, terms in and are neglected because these terms result in the second-order time derivative of ; these terms are and 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 , 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 iswhere , , , , , and . All the parameters in these expressions are given in Appendix E. The second term in (77) originates from and , the third term from and , the fourth term from and , and the fifth term from and . Thus, in order to consider the terms up to , and are necessary in the expansion equations (75) and (76). From the first and second terms in (77), we note that .
Next, we consider the time evolution equation of . Following procedures similar to those for , the time evolution equation of up to iswhere and . and are defined asrespectively. In the above equation, all the terms are . The terms of disappeared because of the orthogonality relation of trigonometric functions. The next higher order term is , which is not included in expansions and in (75) and (76). The second term in (78) originates from and , and the third term from and . From the first and second terms in (78), we note that .
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) intorespectively. We name (81)–(83) as the full tensor model [37]. In the above full tensor model, we chose parameters such that and put and in (81). The other rescaled parameters in (82) and (83) were , , , , , , , and .
We summarize the coefficients in (81)–(83) as follows. is a damping coefficient in the time evolution equation of ; although the motionless spot is stable for , the motionless spot is destabilized for and there appears a moving spot. and are damping coefficients in the time evolution equations of and , respectively. In this study, we fix as a positive constant and consider in the range . and are coefficients of the quadratic and fourth order of , respectively, in the time evolution equation of . The term can be absorbed into with the replacement . From the results, it can be seen that the term is enhanced (suppressed) when the sign of is the same as (different from) that of . is the coefficient of the cubic term of in the time evolution equation of . These terms cause deformations with increasing velocity. The other parameters are coupling coefficients of deformation and velocity. and are the coupling coefficients of and , respectively, and these terms influence the time evolution of and , respectively. and are the coupling coefficients of and , respectively. When , is decoupled from and . can be absorbed into with the replacement . Then, the term is suppressed (enhanced) when is positive (negative), with .
The dependence of , , , , , , and on is shown in Figure 4. In each figure, these parameters (shown on the vertical axis) are scaled by . We see that both and are negative, and , even at . Although is a monotonically increasing function of , is a monotonically decreasing function of . and are positive but is much smaller than .
(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 is large and is small, the effect of on is small in (82). In this case, we put , and the , , and system described by (81)–(83) is reduced to a and 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 rewrite , , and by introducing the following variables:where we choose , , . From (81), we obtain the time evolution equations of and asrespectively. From (82), we obtain the time evolution equations of and asrespectively. 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 and fix the parameters as , , , and .
6.1. 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 . In the case of , the stationary solution of (85) and (86) becomes with . In order to avoid the divergence of , we assume that . For , we consider the dependence of the stationary solution of , , , and on the signs of and . In the stationary state of (86) with , there are two stationary solutions of depending on the sign of and : (I) and are positive with and (II) and are negative with , where is an integer. Using (87)–(90), the stationary solution of case (I) iswhere satisfies 0 and is a stationary solution of , which is determined later. On the other hand, for case (II), the stationary solution isFrom these results, it is seen that has the same expression irrespective of the signs of and . As is positive, satisfying is chosen asFor both cases ((I) and (II)), the stationary solution of is given byIn both cases ((I) and (II)), by substituting into (85) in the stationary state, we obtain two stationary solutions asIn order to satisfy the condition that in the limits , , and , only exists. In the physical images of the moving spot in cases (I) and (II) with , 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
When is large in (89), is rapidly damped to zero. In addition, when 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 in (87) and (88). Using (86) and (88), we can derive the time evolution equation of as
For the stationary states , , we discuss the stability condition of the straight motion. From (96), the stationary solution satisfies or , and their corresponding stability conditions are or , respectively. Since and are negative in our RD system, as shown in Figure 4(a), we consider the case of as follows. In this case, the stability condition is calculated using the given by (92). The bifurcation of straight-circular motion, when , was already studied in [16], and the critical value for the straight motion is expressed asThe phase diagram in the - plane is shown by the solid curve in Figure 5(a). When , the motionless spot is stable. When is in the range , the spot moves straight. However, when , the straight motion become unstable and changes into circular motion. On the other hand, when and , the stability condition for straight motion isThe 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 (the solid curve), it can be seen that when , the parameter region for the circular motion is larger. However, when , the parameter region for the circular motion is smaller. We can explain this result using (87); the term essentially reduces the parameter of the damping term such that the damping term becomes . When is positive (negative), the term 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 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.
(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 . In this parameter region, we examine the stationary circular motion of the spot with a constant angular frequency and velocity . For this, we set and substitute , , , and into (85)–(88), and obtain the relations among , , , and . The calculation is straightforward; the final expressions are shown.
We first consider the case when , where up to terms are considered, and obtain the relations , , , , and . When the spot rotates with an angular frequency , satisfies 0, and this condition leads to [16]. The dependence of the critical velocity for the stable angular frequency on () is shown in Figure 6(a). In this case, there is a lower limit of velocity for the circular motion but no upper limit of velocity for the same.
(a)
(b)
(c)
(d)
(e)
Next, we consider the case when and , where up to terms are considered. After similar calculations to the case of , we obtain the relations , , , and . In this case, is a quartic function of . For the condition , there are two cases depending on the sign of . When , exists for , where are the solutions of , which are given by with . The dependence of on is shown in Figure 6(b); there is only a lower limit of velocity . This property is the same with the case when . On the other hand, when , exists in a certain range of ; under the condition . In this case, there is a lower limit and an upper limit of velocity, and , respectively. The gap between and is defined by . Let us consider two limiting cases. In the limit , and . That is, the circular motion cannot exist. On the other hand, for small and , we expand in terms of and and obtain and . Thus, in the limits and , , , and . This suggests that and there is no upper limit of velocity, which are the same with the case when . From the numerical results shown in Figures 4(a), 4(b), and 4(c), in our RD system; therefore, it turns out that exists only in a finite range, . For larger , the gap is smaller, as shown in Figure 6(c). The dependence of on is shown in Figure 6(d), which suggests that monotonically decreases with increasing . Since increases monotonically as increases (Figure 4(c)), decreases as the chemotaxis increases. When the condition is satisfied, and , which yields the critical value asWhen , , and, therefore, 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 and increase monotonically with increasing , and, therefore, the system diverges in the simulation. For this reason, is shown in the range .
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 and , 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 and , 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 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 , decreases and increases as increases in the range . 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 and increases .
(a)
(b)
From the above analyses on (85)–(88) with , 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 terms are considered, there will be an upper limit for the velocity if up to 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 .
6.4. Simulation of the Full Tensor Model
In this subsection, we examine the effect of on and . For this, we consider the case of . 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) , the region of the circular motion is larger (smaller). The phase diagrams in the - plane and - 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 on the time evolution equation of is similar to that of .
(a)
(b)
(c)
(d)
The phase diagram of spot motion in the - 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 . 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 on the stationary velocity and the corresponding radius of circular motion are shown in Figure 10. On comparing these figures, it can be seen that although and depend on and , the influences of on and are much smaller than those of .
(a)
(b)
(c)
(d)
The dependence of the critical velocity and on is shown in Figure 11. There exists a critical velocity for as shown in Figure 11(a). For larger , the gap is larger. This suggests that enhances the circular motion, which is consistent with the results shown in Figure 8. For positive , becomes larger than that in the reduced tensor model (solid curve) as shown in Figure 11(b). On the other hand, for negative , there is no . In our system, 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.
(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 of up to (see (77)). There were three terms of : , , and . The term can be absorbed into the term by replacing with . Then, with increasing velocity, the term changes the value of only . In our theoretical analysis, considering that the influence of is small, we set and examined the effect of on the stationary solution in Sections 6.2 and 6.3. When up to 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 terms, even in the absence of chemotaxis (), . The term 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 , 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, , 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 . The spatial dependence of the inhibitor before the leading interface and that behind the rear interface are given by () and (), respectively, where , , and 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 for large . For the spot moving along a circle with a given radius , the upper limit of velocity will approximately satisfy the relation : . 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, , 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, 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 were considered, and we found that the terms and 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 , and the spot does not undergo circular motion for .
In practical experiments, the candidates of autonomous systems including chemotaxis are E. 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.
Appendix
A. Velocity of the Interface in One Dimension
Here, we derive the velocity of the interface in one dimension, given by (16), by following [43]. We denote as the value of at the interface at time . Applying the stretching transformations and to (5), we obtainwhere denotes the function evaluated at , and the position of the interface in the stretching coordinate is . In order to consider the traveling front solution of (A.1), we introduce the moving coordinate . Then, (A.1) becomeswhere is the velocity of the front and . The boundary conditions at areand the boundary conditions at are
The stationary solution of (A.2) iswhere is a coefficient and . Under the boundary conditions (A.4), we get and . Using these relations, we obtainIn the main text, in (A.6) is denoted as , which is explicitly expressed as (16).
B. Explict Expressions of and in the Limit
We show that in the limit , and calculated using (13) agree with the ones calculated using (34) and (35) with the substitutions and , where , respectively.
and obtained directly from (13) arerespectively.
We first derive from the definition (34). Substituting and into (34), we obtainWe apply the formulas of the integrals ([44])to (B.3); we can confirm that (B.3) agrees with (B.1).
Next, we derive from the definition (35). Substituting and into (35), we obtainApplying formula (B.4) to (B.6), we can confirm that (B.6) agrees with (B.2).
C. Derivation of (68)
Here, we give the detailed derivation of (68). We first derive . Up to the first order of the deviations, is decomposed aswhere and are given in (25) and (31), respectively.
Substituting into , the () component is obtained asSubstituting and into , is obtained asSubstituting and into , is obtained asUsing the formula for the Bessel function , becomes and are the same as those of the result in [16].
Next, we derive . Up to the first order of the deviations, is decomposed asSubstituting and into , the () component is obtained asSubstituting into , is obtained asSubstituting and into , is obtained asSubstituting and into , is obtained asIn the main text, is expressed asand is expressed aswhere we used as the formula for the Bessel function.
D. Expansions of and in (75) and (76)
The expansions of and in (75) and (76) arewhere and all the coefficients in (D.1) and (D.2) are given in Appendix E. In order to calculate the time evolution of , is multiplied to both sides of eq. (74), and integration is performed over in the range . Using the periodicity of trigonometric functions, we obtain the time evolution equation of . As the tensors and are linear combinations of , we can derive their time evolution equations.
E. Coefficients in (D.1) and (D.2)
We give the coefficients in (D.1) and (D.2) below. By using and , they are expressed as
Conflicts of Interest
The author declares that there are no conflicts of interest.