Table of Contents Author Guidelines Submit a Manuscript
Advances in Mathematical Physics
Volume 2018, Article ID 6152961, 24 pages
https://doi.org/10.1155/2018/6152961
Research Article

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

Department of Complex and Intelligent Systems, School of Systems Information Science, Future University Hakodate, Hakodate 041-8655, Japan

Correspondence should be addressed to Satoshi Kawaguchi; pj.ca.nuf@ihsotas

Received 7 December 2017; Accepted 29 March 2018; Published 21 May 2018

Academic Editor: Christos Volos

Copyright © 2018 Satoshi Kawaguchi. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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 [46]. 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 [912]. 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 [2427], 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.

Figure 1: Boundary layer and interface. The solid and dotted curves correspond to and along a radial direction, respectively, that is, and . In the finite ,   has a boundary layer around with the thickness of , where is determined by the condition .   does not have a boundary layer but has a sharp interface at .

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)).

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

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 .

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

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