Author's Accepted Manuscript
Competitive aggregation phase wave signals
dynamics
using
Hidetsugu Sakaguchi, Satomi Maeyama
www.elsevier.com/locate/yjtbi
PII: DOI: Reference:
S0022-5193(14)00360-9 http://dx.doi.org/10.1016/j.jtbi.2014.06.017 YJTBI7784
To appear in:
Journal of Theoretical Biology
Received date: 10 March 2014 Revised date: 30 April 2014 Accepted date: 12 June 2014 Cite this article as: Hidetsugu Sakaguchi, Satomi Maeyama, Competitive aggregation dynamics using phase wave signals, Journal of Theoretical Biology, http://dx.doi.org/10.1016/j.jtbi.2014.06.017 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting galley proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Competitive aggregation dynamics using phase wave signals Hidetsugu Sakaguchi and Satomi Maeyama Department of Applied Science for Electronics and Materials, Interdisciplinary Graduate School of Engineering Sciences, Kyushu University, Kasuga, Fukuoka 816-8580, Japan
Abstract Coupled equations of the phase equation and the equation of cell concentration n are proposed for competitive aggregation dynamics of slime mold in two dimensions. Phase waves are used as tactic signals of aggregation in this model. Several aggregation clusters are formed initially, and target patterns appear around the localized aggregation clusters. Owing to the competition among target patterns, the number of the localized aggregation clusters decreases, and finally one dominant localized pattern survives. If the phase equation is replaced with the complex Ginzburg-Landau equation, several spiral patterns appear, and n is localized near the center of the spiral patterns. After the competition among spiral patterns, one dominant spiral survives. Key words: Mathematical model,Slime mold, Aggregation dynamics, Phase waves
1
Introduction
Nonlinear nonequilibrium systems can generate various interesting patterns [1]. Target patterns and spiral patterns are typical patterns found in excitable or oscillatory media such as the BZ reaction [2,3]. Similar type of spiral pattern is observed in the early stage of the aggregation of the slime mold Dictyostelium discoideum [4]. The amoebae cells of Dictyostelium discoideum secrete cAMP and the concentration of cAMP exhibits the limit-cycle oscillation. Martiel and Goldbeter proposed a mathematical model for the limit-cycle oscillation of cAMP [5]. Tyson et al. proposed a mathematical model for spiral waves Email address: [email protected] (Hidetsugu Sakaguchi and Satomi Maeyama).
Preprint submitted to Journal of Theoretical Biology
17 June 2014
based on the model of Martiel and Goldbeter [6]. Amoebae cells tend to aggregate to the center of the spiral pattern. The amoebae cells are considered to have chemotaxis to cAMP. Keller and Segel proposed a mathematical model of the aggregation dynamics [7]. The model equations are coupled equations of the concentration of amoebae cells and the concentration of cAMP. Similar type equations have been studied by several authors [8–10]. Recently, Gregor et al. studied collective behaviors of the social amoebae, and found a transition from an excitable state to an oscillatory one by the increase of cell concentration [11]. The phenomenon was interpreted as an example of dynamical quorum sensing [12]. The dynamical quorum sensing is a dynamical transition which occurs by sensing the cell concentration of the same species. The limitcycle oscillation of extra-cellular cAMP was observed in their experiment. The social amoebae begin to aggregate when the period of the oscillation becomes shorter. Several local aggregates are created initially. The competition occurs among the local aggregates, and the number of local aggregates decreases in time, and finally, only a dominant aggregate survives. The chemical signal of cAMP is considered to play an important role in the aggregation dynamics. Generally, the phase dynamics obtained by the phase reduction method is useful for the analysis of various nonlinear phenomena related to the limitcycle oscillation such as mutual synchronization and phase waves [2,13]. There are some reports that the phase information of the limit-cycle oscillation is used for active motions in slime mold Physarum [14]. In spiral waves of the Dictyostelium discodeum, the concentration of cAMP changes periodically in space. A simple chemotaxis toward the higher concentration of cAMP is not always efficient for the aggregation dynamics, because the concentration of cAMP changes between high and low values in the wave pattern such as the spiral wave, and the aggregations trapped in the high concentration regions might move away from the center of the spiral pattern with the cAMP waves. We proposed a one-dimensional model equation for the aggregation dynamics using the phase wave information [15]. The model equations are coupled equations of the cell concentration and the phase of cAMP oscillation in one dimension. This is a modified model of the Keller-Segel model. That is, the concentration of cAMP in the Keller-Segel model is replaced with the phase variable of cAMP oscillation. Van Oss et al. studied a model combining the Martiel-Goldbeter model of the cAMP oscillation with a certain aggregation dynamics and succeeded in reproducing the motion toward the center of a spiral pattern [9]. In their model, the concentration gradient of cAMP is used for the aggregation signal. However, it is furthermore assumed that the cells respond to the cAMP gradients in the forward of cAMP waves, but do not respond to the cAMP gradients in the back of cAMP wave. Roughly speaking, the phase gradient and the cAMP gradient have the same direction in the forward of cAMP wave, and the directions of the two gradients are reversed in the back of cAMP wave. The mechanism of the aggregation toward the wave 2
source in their model is similar to our model in the forward of cAMP wave. On the other hand, in the back of cAMP wave, the cells do not respond to the cAMP gradient in the model of Van Oss et al., but the cells respond to the inverse direction of the cAMP gradient in our model. In our model, the cells move toward the wave source both in the forward and the back of cAMP wave. In other words, the cells are assumed to respond to the propagation of cAMP wave in our model. To the best of our knowledge, there is no evidence of direct detection of phase gradient or wave propagation now. Therefore, our model is a hypothetical model now. In the phototaxis, some organisms move toward the source of light. The gradient of light intensity might be detected, but there is a possibility that the direction of light wave is detected by the position in the organism on which light ray strikes. We consider that the propagation direction of the cAMP wave might be detected by some similar mechanism. We hope that it can be experimentally checked whether cells can detect the propagation of cAMP wave or not by observing the difference of aggregation dynamics between in the cAMP wave and in the stationary cAMP gradient. In the one-dimensional model, we found a kind of dynamical quorum sensing, i.e., a dynamical transition from an excitable state to an oscillatory state when the cell concentration increases during the aggregation process and goes beyond a critical value. In the numerical simulation of the model equation, several aggregation clusters appear initially, but only one dominant cluster survives after the competition among aggregation clusters. Furthermore, we showed that the aggregation dynamics is facilitated in the oscillatory regime, because phase waves propagate far away without decay, and cells move with a velocity proportional to the local phase gradient toward the source of the phase waves. That is, the phase wave information is very efficient for the aggregation dynamics. In this paper, we extend the model equation to a two-dimensional system, because the one-dimensional model is simple but the aggregation process in the slime mold occurs in two dimensions. We find some singularity in the two-dimensional model, which was not observed in the one-dimensional model. Furthermore, we study coupled equations of the cell concentration and the complex Ginzburg-Landau equation to generate spiral patterns, because the phase equation cannot generate a phase singularity point which exists in the center of a spiral.
2
Competitive aggregation dynamics using target waves
The two-dimensional model equations are written as ∂φ αn = − b sin φ + ν∇2 φ + g(∇φ)2 , ∂t 1 + γn 3
(1)
∂n = D∇2 n − C ∂t
∂ ∂x
∂φ ∂ n + ∂x ∂y
∂φ n ∂y
.
(2)
where φ is the phase of the oscillation of the chemotactic factor such as cAMP, n is the concentration of cells, and α, γ, b, ν, g, D and C are parameters. The model in this paper is a qualitative model to understand the the aggregation dynamics. However, there is some relation with experimental results. The first term αn/(1 + γn) implies that natural frequency increases with n but saturates for sufficiently large n owing to the denominator. This type of behavior was reported in the experiment of Gregor et al. [11] That is, when the cell density is low, there is no pulsation of cAMP. The pulsation of cAMP concentration starts when cells begin to aggregate and the cell concentration beyond a critical value. At first it is sporadic like every 15 to 20 min. The period of cAMP oscillation decreases with time and finally takes its minimum value of 6 min. The highest frequency 1/6 per minute corresponds to the parameter α/γ in our model. The parameter α is the proportionality coefficient of the relation of the frequency and the cell concentration n when n is sufficiently small. The second term −b sin φ in Eq. (1) controls the transition from an excitable state to an oscillatory state. The third term in Eq. (1) expresses the phase diffusion and the fourth term is a nonlinear term which appears by the standard phase reduction method explained in [2], and the coefficients ν and g can be estimated by the phase reduction method using experimental data of cAMP oscillation in principle. The first term of Eq. (2) represents the diffusion of n and the second term denotes the chemotaxis using the phase wave information, that is, the cells are assumed to move toward a source of phase waves in proportion to the phase gradient −∂φ/∂x. The total number of cells N = n(x, y, t)dxdy is conserved in the time evolution of Eq. (2). In the numerical simulation, the system size is assumed to be L × L, and periodic boundary conditions are imposed. There is a uniform and stationary solution φ(x, y, t) = φ0 and n(x, y, t) = n0 = N/L2 , which satisfies αn0 = b sin φ0 , 1 + γn0
(3)
if αn0 /(1 + γn0 ) < b. However, the uniform solution becomes unstable for D < αCn0 /{b cos φ0 (1 + γn0 )2 }, which can be shown by the linear stability analysis around the uniform and stationary solution [15]. When αn0 /(1 + γn0 ) > b, there is no stationary solution, and phase oscillation occurs. Even in the case of the oscillatory state, a uniformly oscillating state is unstable, and the clustering of cells is expected to occur. Figure 1 shows snapshot patterns of φ and n at (a) t = 300, (b) 400, (c) 500 and (d) 1500 for α = 0.25, γ = 0.5, b = 0.04, ν = 0.02, g = 0.025, D = 4
Fig. 1. (color online) Snapshot patterns of φ and n in Eqs. (1) and (2) at (a) t = 300, (b) 400, (c) 500 and (d) 1500 for α = 0.25, γ = 0.5, b = 0.04, ν = 0.02, g = 0.025, D = 0.03, C = 0.05 and n0 = 0.2.
0.03, C = 0.05, and n0 = 0.2. The system size is L × L = 20 × 20. The initial condition is φ(x, y) = 0 and n(x, y) = 0.2 + rn(x, y), where rn(x, y) is a random number between -0.1 and 0.1. n(x, y) > 0.22 in blue regions, and cos φ > 0 in green regions. For these parameters, αn0 /(1 + γn0 ) > b, and phase oscillation occurs for the uniform state. The uniform state is unstable and the cells tend to aggregate and form localized clusters. Target waves are emitted from the localized clusters. The frequency of target waves is higher for a localized cluster with larger concentration n(x, y). When the frequency of a target pattern is faster than that of the neighboring target pattern, the boundary between the two target patterns moves toward the slower target pattern. As a result, the target pattern with faster frequency invades the neighboring slower target pattern. When the boundary approaches the core region of the slower target pattern, the cells localized near the core region of the slower target pattern begin to move toward the core region of the target pattern with the faster frequency owing to the tactic signals of phase waves, and are swallowed into the faster target pattern. This process is analogous to a struggle of territorial expansion. The competition between target patterns occurs sequentially, and the number of target patterns decreases. Finally only one circular localized cluster appears at t = 1500. Tongue-like patterns of n observed at t = 400 and 500 imply that the cells in invaded localized clusters are being swallowed by the dominant neighboring cluster at x ∼ 11, y ∼ 16. The survived target pattern accompanying a localized cluster of n in the center is an attractor of the nonlinear equations Eqs. (1) and (2). It can be interpreted as a localized dissipative structure which appears as a result of the instability of the uniform state. A localized target pattern was previously found in the electrohydrodynamic convection of liquid crystals [16], although the formation mechanism is different. A one-dimensional model for the localized target pattern was proposed by one of the present authors [17], however, we did not succeed in modeling in two dimensions. The existence of a localized dissipative structure in a two-dimensional model equation is not a trivial extension of that in a one-dimensional system. 5
4 220 (x,L/2)
n(x,L/2)
3 2
210
1 200
0 5
7
9
11 x
13
15
0
5
10 x
15
20
Fig. 2. (color online) (a) Snapshot pattern of φ and n in Eqs. (1) and (2) for α = 0.25, γ = 1, b = 0, ν = 0.02, g = 0.025, D = 0.03, C = 0.05 and n0 = 0.2. (b) n(x, y) at the cross section y = L/2. (c) φ(x, y) at the cross section y = L/2.
Hereafter, we consider only the case of b = 0. In this case, Eqs. (1) and (2) become a model of purely oscillatory systems. Figure 2(a) is a stationary pattern of φ and n for α = 0.25, γ = 1, b = 0, ν = 0.02, g = 0.025, D = 0.03, C = 0.05 2 2 and N = 8. The initial conditions are φ = 0 and n(x, y) = Ae−x −y with A = 8/π. Figures 2(b) and (c) show n(x, y) and φ(x, y) at the cross-section y = L/2 = 10. n(x, y) is localized around (x, y) = (L/2, L/2), the phase profile of φ(x, y) shows that phase waves are emitted from the localized center of n, and the wavenumber k of the phase waves is represented by the phase gradient −∂φ/∂x. In case of b = 0, Eq. (2) can be rescaled with respect to the parameter γ. That is, stationary solutions n(x, y) and φ(x, y) at parameter γ can be expressed as n(x, y) = (1/γ)n0 (x/γ 1/2 , y/γ 1/2) and φ(x, y) = φ0 (x/γ 1/2 , y/γ 1/2 ) using stationary solutions n0 and φ0 at γ = 1. Owing to this relation, the peak amplitude of n increases with 1/γ and the width of the localized solution decreases with γ 1/2 as γ → 0 for the same total number N of cells. We have checked this relation by direct numerical simulations. This relation implies that the localized solution diverges at γ = 0. On the other hand, in the one-dimensional model corresponding to Eqs. (1) and (2), we obtained an exact solution of the form n(x) = B/ cosh2 (kx) at γ = 0 and b = 0 in [17]. The divergence of the solution at γ = 0 is characteristic of the two-dimensional model. It might be related to the collapse phenomenon in the two-dimensional nonlinear Schr¨odinger equation [18]. The singularity in the Keller-Segel type equations has been an important topic in applied mathematics. Herroro and Vel´azquez showed the chemotaxis collapse in the Keller-Segel type model in two dimensions [19]. The origin of the singularity at γ = 0 is clear in our two-dimensional model and the result might be interesting as a problem of applied mathematics. 6
3
Competitive aggregation dynamics using spiral waves
The coupled equations (1) and (2) exhibit interesting pattern dynamics, however, spiral patterns cannot appear. It is because there is a phase singularity point in the center of a spiral, [2,3] and the phase singularity point cannot be expressed by the phase equation. We can use the complex Ginzburg-Landau equation as a simple model to incorporate the effect of the amplitude of oscillation [1,2,20,21]. The complex Ginzburg-Landau equation has a form ∂W = (1 + ω)W − (1 + ic2 )|W |2W + ν∇2 W, ∂t
(4)
where ω, c1 and c2 are parameters. By the phase description, W is assumed to be W ∼ eiφ(x,y,t) . The phase variable obeys the phase equation [2]: ∂φ = ω − c2 + ν∇2 φ + c2 ν(∇φ)2 . ∂t
(5)
If c2 is assumed to be g/ν, and ω is assumed to be ω = c2 + α/(1 + γn), Eq. (1) with b = 0 is obtained. The phase gradients ∂φ/∂x and ∂φ/∂y can be approximated respectively by (1/2)(W ∗∂W/∂x − W ∂W ∗ /∂x) and (1/2)(W ∗∂W/∂y − W ∂W ∗ /∂y) owing to the relation W ∼ eiφ . Using these relations, we propose coupled equations of the complex Ginzburg-Landau equation and the diffusion equation of n instead of Eqs. (1) and (2):
αn ∂W = 1+i W + (1 + ic2 )|W |2 W + (1 + ic2 )∇2 W, (6) ∂t 1 + γn ∂n ∂W ∗ ∂W ∗ C ∂ ∂ 2 ∗ ∂W ∗ ∂W = D∇ n − −W −W W n + W n . ∂t 2 ∂x ∂x ∂x ∂y ∂y ∂y (7) We have performed numerical simulation of Eqs. (6) and (7). Figure 3(a) shows a snapshot pattern of W and n at t = 1000 for ν = 0.02, α = 0.25, γ = 1, D = 2 2 0.03 and C = 0.05. The initial conditions are n(x, y) = Ae−x −y with A = 8/π and W = X + iY = 1 where X and Y are real and imaginary parts of W . n > 0.22 in blue regions, and X > 0 in green regions. Figure 3(a) is similar to Fig. 2(a), because the phase description of Eq. (6) reduces to Eq. (1). If the 2 2 initial conditions are set to be n(x, y) = Ae−x −y with A = 8/π and W = eiθ where θ denotes the polar angle tan−1 (y/x) around the center (L/2, L/2), a spiral pattern appears and the cells aggregate toward the spiral core, however, the whole synchronization cannot be attained owing to the too fast pace of spiral core at ν = 0.02, and a new spiral pair is generated near the boundaries. 7
Fig. 3. (color online) (a) Snapshot pattern of W and n for Eqs. (6) and (7) at t = 1000 for ν = 0.02, α = 0.25, γ = 1, D = 0.03 and C = 0.05. The initial 2 2 conditions are n(x, y) = Ae−x −y with A = 8/π and W = 1. (b) Snapshot pattern at t = 1000 for ν = 0.02, α = 0.25, γ = 1, D = 0.03 and C = 0.05. The initial 2 2 conditions are set to be n(x, y) = Ae−x −y with A = 8/π and W = eiθ . Snapshot pattern at t = 1000 for ν = 0.04, α = 0.25, γ = 0.5, D = 0.03 and C = 0.05. The 2 2 initial conditions are set to be n(x, y) = Ae−x −y with A = 8/π and W = eiθ .
Fig. 4. (color online) Snapshot patterns of W and n at (a) t = 400 and (b) t = 1200 for ν = 0.04, α = 0.25, γ = 1, D = 0.03 and C = 0.05. The initial condition of n(x, y) = 8/L2 and W = X + iY where X and Y are random numbers between -0.05 and 0.05.
Fig. 5. (color online) Three snapshot patterns of W and n at (a) t = 400 (b) t = 800 and (c) t = 1500 for ν = 0.03, α = 0.25, γ = 0.5, D = 0.03 and C = 0.05. The initial conditions are the same as the one shown in Fig. 4.
The state including four spirals is rather stable and is maintained even at t = 1000 as shown in Fig. 3(b). Figure 3(c) shows a snapshot pattern of W and n at t = 1000 for ν = 0.04, α = 0.25, γ = 1, D = 0.03 and C = 0.05. The initial 2 2 condition is n(x, y) = Ae−x −y with A = 8/π and W = eiθ , which is the same 8
Fig. 6. (color online) Two snapshot patterns of W and n at (a) t = 2000 (b) t = 10000 for ν = 0.02, α = 0.25, γ = 0.5, D = 0.03 and C = 0.05. The initial conditions are the same as the one shown in Fig. 4.
as the one shown in Fig. 3(b). At ν = 0.04, a regular spiral pattern appears. There are at least two phase singularity points under the periodic boundary conditions. However, the central vortex point becomes dominant and a spiral pattern appears around the central vortex, because the cell concentration n is localized only near the central vortex. At this parameter, the spiral state is stable and the synchronization is attained in the whole system. The center of spiral is a phase singularity point and the phase cannot be defined there. However, the phase can be defined except for the spiral center, and the phase increases toward the center of spiral. It is therefore natural that cells aggregate toward the spiral center. Next, we have performed numerical simulations from an initial condition of n(x, y) = 8/L2 and W = X +iY where X and Y are random numbers between -0.05 and 0.05. The concentration n(x, y) is uniform in this initial condition. Figure 4 shows two snapshot patterns of W and n at (a) t = 400 and (b) t = 1200 for ν = 0.04, α = 0.25, γ = 1, D = 0.03 and C = 0.05. Many spirals are generated owing to the random initial condition of W . Only one spiral survives at t = 1200. In the single complex Ginzburg-Landau equation Eq. (4), many spirals coexist stably for the random initial condition, because the interaction between spiral cores is very weak. In the coupled equations Eqs. (6) and (7), the competition occurs among spiral patterns owing to the frequency difference induced by the cell concentration n(x, y). As a result of the competition, the winner spiral pattern swallows the cells located in the core region of the loser spiral, and the number of spirals decreases. Finally, only one spiral survives. The competitive dynamics among spiral patterns is characteristic in our system Eqs. (6) and (7). We furthermore study a relation of the aggregation dynamics of cells and the complex dynamics of spiral patterns. Figure 5 shows three snapshot patterns of W and n at (a) t = 400 (b) t = 800 and (c) t = 1500 for ν = 0.03, α = 0.25, γ = 0.5, D = 0.03 and C = 0.05. The initial conditions are the same as the one shown in Fig. 4. Many spirals are generated as in the case of Fig. 4, and only one spiral survives at t = 800. However, the one-spiral state is not stationary because of the too fast pace of the spiral core. A vortex pair is created at a point apart from 9
the center of the spiral, and then one vortex collides with the original vortex existing in the center of the spiral pattern, and the vortex pair disappear near the center of the spiral. Finally, a stable target pattern is created at t = 1500. It is because the attractive force of target pattern for the cell aggregation is weaker than that of spiral pattern, and the peak concentration of n is smaller than that for the spiral pattern. Therefore, the frequency of the core region of target pattern is not so high, and the synchronization is attained in the whole system. As a result, the target pattern is stably maintained. However, at a smaller value of ν, spirals breaks up and spiral chaos appears. Figure 6 shows two snapshot patterns at (a) t = 2000 and (b) t = 10000 at ν = 0.02 for the same other parameters α = 0.25, γ = 0.5, D = 0.03 and C = 0.05. As a result of the spiral chaos, the clustering toward the only one cluster cannot occur. Several aggregation clusters survive around spiral cores and the positions of aggregation clusters change randomly owing to the spiral chaos.
4
Summary
We have proposed two-dimensional coupled equations of the phase equation and the equation for cell concentration n to understand the efficient aggregation dynamics of the slime mold qualitatively. Initially many aggregation clusters appears, but the cluster number decreases owing to the competitive dynamics among target patterns, and finally only one localized aggregation cluster survives. The long-range competitive interaction is due to the nonlinear phase dynamics. The survived target pattern accompanying the localized aggregation cluster at the center can be interpreted as a dissipative structure in a nonlinear-nonequilibrium system. In the two-dimensional model for b = 0, the peak amplitude of cell concentration n increases when the parameter γ decreases to 0. The peak amplitude diverges and a singularity appears at γ = 0 in the two-dimensional model. It can be interpreted as a kind of chemotaxis collapse, which is a phenomenon not observed in the one-dimensional model. We have further generalized the model equations to coupled equations of the complex Ginzburg-Landau equation and the equation of cell concentration. Only one spiral pattern survives as a result of the competitive aggregation dynamics, and the cell concentration is localized near the spiral core. We have found more complex dynamical behaviors when the desynchronization occurs around the spiral core. For a strongly unstable case, the clustering toward the only one cluster cannot be realized owing to the breakup of spirals. Then, several clusters coexist and the positions of localized clusters of n change in time randomly owing to spiral chaos. 10
References
[1] Cross, M., Greedside, H.,2009. Pattern Formation and Dynamics in Nonequilibrium Systems (Cambrdge University Press, Cambridge). [2] Kuramoto, Y.,1984. Chemical Oscillations, Waves and Turbulence (SpringerVerlag, New York). [3] Winfree, A. T.,1980. The Geometry of Biological Time (Springer-Verlag, New York). [4] Gross,J., Peacey, M., Trevan, D.,1976. Signal emission and signal propagation during early aggregation in Dictyostelium-Discoideum. J. Cell Sci. 22, 645-656. [5] Martiel,J., Goldbeter, A., 1987. A model based on receptor desensitization for cyclic-AMP signaling in Dictyostelium cells. Biophys. J. 52, 807-828. [6] Tyson, J., Alexander, K., Manoranjan, V., Murray, J., 1987. Spiral waves of cyclic-AMP in a model of slime-mold aggregation. Physica D 34, 193-207. [7] Keller, E., Segel, L., 1970. Initiation of slime mold aggregation viewed as an instability. J. Theor. Biol. 26, 399-415. [8] Nanjundiah, V., 1973. Chemotaxis, signal relaying and aggregation morphology. J. Theor. Biol. 42, 63-105. [9] Van Oss, C., Panfilov, Hogeweg, Siegert, A., Weijer, C., 1996. Spatial pattern formation during aggregation of the slime mold Dictyostelium discoideum. J. Theor. Biol. 181, 203-213. [10] Hillen, T., Painter, K., 2009. A user’s guide to PDE models for chemotaxis. J. Math. Biol. 58, 183-217. [11] Gregor, T., Fujimoto, K., Masaki, N., Sawai, S., 2010. The onset of collective behavior in social amoebae. Science 328, 1021-1025. [12] Garcia-Ojalvo,J., Elowitz, M.,Strogatz, S., 2004. Modeling a synthetic multicellur clock:repressillators coupled by quorum sensing. PNAS 101,1095510960. [13] Winfree, A., 1967. Biological rhythms and behavior of populations of coupled oscillators. J. Theor. Phys. 16, 15-42. [14] Matsumoto, K., Ueda, T., Kobatake. Y., 1988. Reversal of themotaxis with oscillatory stimulation in the plasmodium of Physarum polycephalum. J. Theor. Phys. 131, 175-182. [15] Sakaguchi, H., Maeyama, S., 2013. Dynamical quorum sensing and clustering dynamics in a population of spatially distributed active rotators. Phys. Rev. E 87, 024901.
11
[16] Nasuno, S., Sano, M., Sawada, Y., 1989. Phase wave propagation in the rectangular convective structure of nematic liquid-crystal. J. Phys. Soc. Jpn. 58, 1875-1878. [17] Sakaguchi, H., 1992. Localized oscillation in a cellular pattern. Prog. Theor. Phys. 87, 1049-1053. [18] Kivshar, Y., Agrawal, G., 2003. Optical Solitons (Academic Press, San Diego). [19] Herrero, M., Vel´ azquez, J., 1996. Singularity patterns in a chemotaxis model. Math. Ann. 306, 583-623. [20] Aranson, I., Kramer, L., 2002. The world of the complex Ginzburg-Landau equation. Rev. Mod. Phys. 74, 99-143. [21] Sakaguchi, H., 1993. Phase dynamics and localized solutions to the GinzburgLandau type amplitude equations. Prog Theor. Phys. 89, 1123-1146.
12
Competitive aggregation phase wave signals
dynamics
using
Hidetsugu Sakaguchi, Satomi Maeyama
www.elsevier.com/locate/yjtbi
PII: DOI: Reference:
S0022-5193(14)00360-9 http://dx.doi.org/10.1016/j.jtbi.2014.06.017 YJTBI7784
To appear in:
Journal of Theoretical Biology
Received date: 10 March 2014 Revised date: 30 April 2014 Accepted date: 12 June 2014 Cite this article as: Hidetsugu Sakaguchi, Satomi Maeyama, Competitive aggregation dynamics using phase wave signals, Journal of Theoretical Biology, http://dx.doi.org/10.1016/j.jtbi.2014.06.017 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting galley proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Competitive aggregation dynamics using phase wave signals Hidetsugu Sakaguchi and Satomi Maeyama Department of Applied Science for Electronics and Materials, Interdisciplinary Graduate School of Engineering Sciences, Kyushu University, Kasuga, Fukuoka 816-8580, Japan
Abstract Coupled equations of the phase equation and the equation of cell concentration n are proposed for competitive aggregation dynamics of slime mold in two dimensions. Phase waves are used as tactic signals of aggregation in this model. Several aggregation clusters are formed initially, and target patterns appear around the localized aggregation clusters. Owing to the competition among target patterns, the number of the localized aggregation clusters decreases, and finally one dominant localized pattern survives. If the phase equation is replaced with the complex Ginzburg-Landau equation, several spiral patterns appear, and n is localized near the center of the spiral patterns. After the competition among spiral patterns, one dominant spiral survives. Key words: Mathematical model,Slime mold, Aggregation dynamics, Phase waves
1
Introduction
Nonlinear nonequilibrium systems can generate various interesting patterns [1]. Target patterns and spiral patterns are typical patterns found in excitable or oscillatory media such as the BZ reaction [2,3]. Similar type of spiral pattern is observed in the early stage of the aggregation of the slime mold Dictyostelium discoideum [4]. The amoebae cells of Dictyostelium discoideum secrete cAMP and the concentration of cAMP exhibits the limit-cycle oscillation. Martiel and Goldbeter proposed a mathematical model for the limit-cycle oscillation of cAMP [5]. Tyson et al. proposed a mathematical model for spiral waves Email address: [email protected] (Hidetsugu Sakaguchi and Satomi Maeyama).
Preprint submitted to Journal of Theoretical Biology
17 June 2014
based on the model of Martiel and Goldbeter [6]. Amoebae cells tend to aggregate to the center of the spiral pattern. The amoebae cells are considered to have chemotaxis to cAMP. Keller and Segel proposed a mathematical model of the aggregation dynamics [7]. The model equations are coupled equations of the concentration of amoebae cells and the concentration of cAMP. Similar type equations have been studied by several authors [8–10]. Recently, Gregor et al. studied collective behaviors of the social amoebae, and found a transition from an excitable state to an oscillatory one by the increase of cell concentration [11]. The phenomenon was interpreted as an example of dynamical quorum sensing [12]. The dynamical quorum sensing is a dynamical transition which occurs by sensing the cell concentration of the same species. The limitcycle oscillation of extra-cellular cAMP was observed in their experiment. The social amoebae begin to aggregate when the period of the oscillation becomes shorter. Several local aggregates are created initially. The competition occurs among the local aggregates, and the number of local aggregates decreases in time, and finally, only a dominant aggregate survives. The chemical signal of cAMP is considered to play an important role in the aggregation dynamics. Generally, the phase dynamics obtained by the phase reduction method is useful for the analysis of various nonlinear phenomena related to the limitcycle oscillation such as mutual synchronization and phase waves [2,13]. There are some reports that the phase information of the limit-cycle oscillation is used for active motions in slime mold Physarum [14]. In spiral waves of the Dictyostelium discodeum, the concentration of cAMP changes periodically in space. A simple chemotaxis toward the higher concentration of cAMP is not always efficient for the aggregation dynamics, because the concentration of cAMP changes between high and low values in the wave pattern such as the spiral wave, and the aggregations trapped in the high concentration regions might move away from the center of the spiral pattern with the cAMP waves. We proposed a one-dimensional model equation for the aggregation dynamics using the phase wave information [15]. The model equations are coupled equations of the cell concentration and the phase of cAMP oscillation in one dimension. This is a modified model of the Keller-Segel model. That is, the concentration of cAMP in the Keller-Segel model is replaced with the phase variable of cAMP oscillation. Van Oss et al. studied a model combining the Martiel-Goldbeter model of the cAMP oscillation with a certain aggregation dynamics and succeeded in reproducing the motion toward the center of a spiral pattern [9]. In their model, the concentration gradient of cAMP is used for the aggregation signal. However, it is furthermore assumed that the cells respond to the cAMP gradients in the forward of cAMP waves, but do not respond to the cAMP gradients in the back of cAMP wave. Roughly speaking, the phase gradient and the cAMP gradient have the same direction in the forward of cAMP wave, and the directions of the two gradients are reversed in the back of cAMP wave. The mechanism of the aggregation toward the wave 2
source in their model is similar to our model in the forward of cAMP wave. On the other hand, in the back of cAMP wave, the cells do not respond to the cAMP gradient in the model of Van Oss et al., but the cells respond to the inverse direction of the cAMP gradient in our model. In our model, the cells move toward the wave source both in the forward and the back of cAMP wave. In other words, the cells are assumed to respond to the propagation of cAMP wave in our model. To the best of our knowledge, there is no evidence of direct detection of phase gradient or wave propagation now. Therefore, our model is a hypothetical model now. In the phototaxis, some organisms move toward the source of light. The gradient of light intensity might be detected, but there is a possibility that the direction of light wave is detected by the position in the organism on which light ray strikes. We consider that the propagation direction of the cAMP wave might be detected by some similar mechanism. We hope that it can be experimentally checked whether cells can detect the propagation of cAMP wave or not by observing the difference of aggregation dynamics between in the cAMP wave and in the stationary cAMP gradient. In the one-dimensional model, we found a kind of dynamical quorum sensing, i.e., a dynamical transition from an excitable state to an oscillatory state when the cell concentration increases during the aggregation process and goes beyond a critical value. In the numerical simulation of the model equation, several aggregation clusters appear initially, but only one dominant cluster survives after the competition among aggregation clusters. Furthermore, we showed that the aggregation dynamics is facilitated in the oscillatory regime, because phase waves propagate far away without decay, and cells move with a velocity proportional to the local phase gradient toward the source of the phase waves. That is, the phase wave information is very efficient for the aggregation dynamics. In this paper, we extend the model equation to a two-dimensional system, because the one-dimensional model is simple but the aggregation process in the slime mold occurs in two dimensions. We find some singularity in the two-dimensional model, which was not observed in the one-dimensional model. Furthermore, we study coupled equations of the cell concentration and the complex Ginzburg-Landau equation to generate spiral patterns, because the phase equation cannot generate a phase singularity point which exists in the center of a spiral.
2
Competitive aggregation dynamics using target waves
The two-dimensional model equations are written as ∂φ αn = − b sin φ + ν∇2 φ + g(∇φ)2 , ∂t 1 + γn 3
(1)
∂n = D∇2 n − C ∂t
∂ ∂x
∂φ ∂ n + ∂x ∂y
∂φ n ∂y
.
(2)
where φ is the phase of the oscillation of the chemotactic factor such as cAMP, n is the concentration of cells, and α, γ, b, ν, g, D and C are parameters. The model in this paper is a qualitative model to understand the the aggregation dynamics. However, there is some relation with experimental results. The first term αn/(1 + γn) implies that natural frequency increases with n but saturates for sufficiently large n owing to the denominator. This type of behavior was reported in the experiment of Gregor et al. [11] That is, when the cell density is low, there is no pulsation of cAMP. The pulsation of cAMP concentration starts when cells begin to aggregate and the cell concentration beyond a critical value. At first it is sporadic like every 15 to 20 min. The period of cAMP oscillation decreases with time and finally takes its minimum value of 6 min. The highest frequency 1/6 per minute corresponds to the parameter α/γ in our model. The parameter α is the proportionality coefficient of the relation of the frequency and the cell concentration n when n is sufficiently small. The second term −b sin φ in Eq. (1) controls the transition from an excitable state to an oscillatory state. The third term in Eq. (1) expresses the phase diffusion and the fourth term is a nonlinear term which appears by the standard phase reduction method explained in [2], and the coefficients ν and g can be estimated by the phase reduction method using experimental data of cAMP oscillation in principle. The first term of Eq. (2) represents the diffusion of n and the second term denotes the chemotaxis using the phase wave information, that is, the cells are assumed to move toward a source of phase waves in proportion to the phase gradient −∂φ/∂x. The total number of cells N = n(x, y, t)dxdy is conserved in the time evolution of Eq. (2). In the numerical simulation, the system size is assumed to be L × L, and periodic boundary conditions are imposed. There is a uniform and stationary solution φ(x, y, t) = φ0 and n(x, y, t) = n0 = N/L2 , which satisfies αn0 = b sin φ0 , 1 + γn0
(3)
if αn0 /(1 + γn0 ) < b. However, the uniform solution becomes unstable for D < αCn0 /{b cos φ0 (1 + γn0 )2 }, which can be shown by the linear stability analysis around the uniform and stationary solution [15]. When αn0 /(1 + γn0 ) > b, there is no stationary solution, and phase oscillation occurs. Even in the case of the oscillatory state, a uniformly oscillating state is unstable, and the clustering of cells is expected to occur. Figure 1 shows snapshot patterns of φ and n at (a) t = 300, (b) 400, (c) 500 and (d) 1500 for α = 0.25, γ = 0.5, b = 0.04, ν = 0.02, g = 0.025, D = 4
Fig. 1. (color online) Snapshot patterns of φ and n in Eqs. (1) and (2) at (a) t = 300, (b) 400, (c) 500 and (d) 1500 for α = 0.25, γ = 0.5, b = 0.04, ν = 0.02, g = 0.025, D = 0.03, C = 0.05 and n0 = 0.2.
0.03, C = 0.05, and n0 = 0.2. The system size is L × L = 20 × 20. The initial condition is φ(x, y) = 0 and n(x, y) = 0.2 + rn(x, y), where rn(x, y) is a random number between -0.1 and 0.1. n(x, y) > 0.22 in blue regions, and cos φ > 0 in green regions. For these parameters, αn0 /(1 + γn0 ) > b, and phase oscillation occurs for the uniform state. The uniform state is unstable and the cells tend to aggregate and form localized clusters. Target waves are emitted from the localized clusters. The frequency of target waves is higher for a localized cluster with larger concentration n(x, y). When the frequency of a target pattern is faster than that of the neighboring target pattern, the boundary between the two target patterns moves toward the slower target pattern. As a result, the target pattern with faster frequency invades the neighboring slower target pattern. When the boundary approaches the core region of the slower target pattern, the cells localized near the core region of the slower target pattern begin to move toward the core region of the target pattern with the faster frequency owing to the tactic signals of phase waves, and are swallowed into the faster target pattern. This process is analogous to a struggle of territorial expansion. The competition between target patterns occurs sequentially, and the number of target patterns decreases. Finally only one circular localized cluster appears at t = 1500. Tongue-like patterns of n observed at t = 400 and 500 imply that the cells in invaded localized clusters are being swallowed by the dominant neighboring cluster at x ∼ 11, y ∼ 16. The survived target pattern accompanying a localized cluster of n in the center is an attractor of the nonlinear equations Eqs. (1) and (2). It can be interpreted as a localized dissipative structure which appears as a result of the instability of the uniform state. A localized target pattern was previously found in the electrohydrodynamic convection of liquid crystals [16], although the formation mechanism is different. A one-dimensional model for the localized target pattern was proposed by one of the present authors [17], however, we did not succeed in modeling in two dimensions. The existence of a localized dissipative structure in a two-dimensional model equation is not a trivial extension of that in a one-dimensional system. 5
4 220 (x,L/2)
n(x,L/2)
3 2
210
1 200
0 5
7
9
11 x
13
15
0
5
10 x
15
20
Fig. 2. (color online) (a) Snapshot pattern of φ and n in Eqs. (1) and (2) for α = 0.25, γ = 1, b = 0, ν = 0.02, g = 0.025, D = 0.03, C = 0.05 and n0 = 0.2. (b) n(x, y) at the cross section y = L/2. (c) φ(x, y) at the cross section y = L/2.
Hereafter, we consider only the case of b = 0. In this case, Eqs. (1) and (2) become a model of purely oscillatory systems. Figure 2(a) is a stationary pattern of φ and n for α = 0.25, γ = 1, b = 0, ν = 0.02, g = 0.025, D = 0.03, C = 0.05 2 2 and N = 8. The initial conditions are φ = 0 and n(x, y) = Ae−x −y with A = 8/π. Figures 2(b) and (c) show n(x, y) and φ(x, y) at the cross-section y = L/2 = 10. n(x, y) is localized around (x, y) = (L/2, L/2), the phase profile of φ(x, y) shows that phase waves are emitted from the localized center of n, and the wavenumber k of the phase waves is represented by the phase gradient −∂φ/∂x. In case of b = 0, Eq. (2) can be rescaled with respect to the parameter γ. That is, stationary solutions n(x, y) and φ(x, y) at parameter γ can be expressed as n(x, y) = (1/γ)n0 (x/γ 1/2 , y/γ 1/2) and φ(x, y) = φ0 (x/γ 1/2 , y/γ 1/2 ) using stationary solutions n0 and φ0 at γ = 1. Owing to this relation, the peak amplitude of n increases with 1/γ and the width of the localized solution decreases with γ 1/2 as γ → 0 for the same total number N of cells. We have checked this relation by direct numerical simulations. This relation implies that the localized solution diverges at γ = 0. On the other hand, in the one-dimensional model corresponding to Eqs. (1) and (2), we obtained an exact solution of the form n(x) = B/ cosh2 (kx) at γ = 0 and b = 0 in [17]. The divergence of the solution at γ = 0 is characteristic of the two-dimensional model. It might be related to the collapse phenomenon in the two-dimensional nonlinear Schr¨odinger equation [18]. The singularity in the Keller-Segel type equations has been an important topic in applied mathematics. Herroro and Vel´azquez showed the chemotaxis collapse in the Keller-Segel type model in two dimensions [19]. The origin of the singularity at γ = 0 is clear in our two-dimensional model and the result might be interesting as a problem of applied mathematics. 6
3
Competitive aggregation dynamics using spiral waves
The coupled equations (1) and (2) exhibit interesting pattern dynamics, however, spiral patterns cannot appear. It is because there is a phase singularity point in the center of a spiral, [2,3] and the phase singularity point cannot be expressed by the phase equation. We can use the complex Ginzburg-Landau equation as a simple model to incorporate the effect of the amplitude of oscillation [1,2,20,21]. The complex Ginzburg-Landau equation has a form ∂W = (1 + ω)W − (1 + ic2 )|W |2W + ν∇2 W, ∂t
(4)
where ω, c1 and c2 are parameters. By the phase description, W is assumed to be W ∼ eiφ(x,y,t) . The phase variable obeys the phase equation [2]: ∂φ = ω − c2 + ν∇2 φ + c2 ν(∇φ)2 . ∂t
(5)
If c2 is assumed to be g/ν, and ω is assumed to be ω = c2 + α/(1 + γn), Eq. (1) with b = 0 is obtained. The phase gradients ∂φ/∂x and ∂φ/∂y can be approximated respectively by (1/2)(W ∗∂W/∂x − W ∂W ∗ /∂x) and (1/2)(W ∗∂W/∂y − W ∂W ∗ /∂y) owing to the relation W ∼ eiφ . Using these relations, we propose coupled equations of the complex Ginzburg-Landau equation and the diffusion equation of n instead of Eqs. (1) and (2):
αn ∂W = 1+i W + (1 + ic2 )|W |2 W + (1 + ic2 )∇2 W, (6) ∂t 1 + γn ∂n ∂W ∗ ∂W ∗ C ∂ ∂ 2 ∗ ∂W ∗ ∂W = D∇ n − −W −W W n + W n . ∂t 2 ∂x ∂x ∂x ∂y ∂y ∂y (7) We have performed numerical simulation of Eqs. (6) and (7). Figure 3(a) shows a snapshot pattern of W and n at t = 1000 for ν = 0.02, α = 0.25, γ = 1, D = 2 2 0.03 and C = 0.05. The initial conditions are n(x, y) = Ae−x −y with A = 8/π and W = X + iY = 1 where X and Y are real and imaginary parts of W . n > 0.22 in blue regions, and X > 0 in green regions. Figure 3(a) is similar to Fig. 2(a), because the phase description of Eq. (6) reduces to Eq. (1). If the 2 2 initial conditions are set to be n(x, y) = Ae−x −y with A = 8/π and W = eiθ where θ denotes the polar angle tan−1 (y/x) around the center (L/2, L/2), a spiral pattern appears and the cells aggregate toward the spiral core, however, the whole synchronization cannot be attained owing to the too fast pace of spiral core at ν = 0.02, and a new spiral pair is generated near the boundaries. 7
Fig. 3. (color online) (a) Snapshot pattern of W and n for Eqs. (6) and (7) at t = 1000 for ν = 0.02, α = 0.25, γ = 1, D = 0.03 and C = 0.05. The initial 2 2 conditions are n(x, y) = Ae−x −y with A = 8/π and W = 1. (b) Snapshot pattern at t = 1000 for ν = 0.02, α = 0.25, γ = 1, D = 0.03 and C = 0.05. The initial 2 2 conditions are set to be n(x, y) = Ae−x −y with A = 8/π and W = eiθ . Snapshot pattern at t = 1000 for ν = 0.04, α = 0.25, γ = 0.5, D = 0.03 and C = 0.05. The 2 2 initial conditions are set to be n(x, y) = Ae−x −y with A = 8/π and W = eiθ .
Fig. 4. (color online) Snapshot patterns of W and n at (a) t = 400 and (b) t = 1200 for ν = 0.04, α = 0.25, γ = 1, D = 0.03 and C = 0.05. The initial condition of n(x, y) = 8/L2 and W = X + iY where X and Y are random numbers between -0.05 and 0.05.
Fig. 5. (color online) Three snapshot patterns of W and n at (a) t = 400 (b) t = 800 and (c) t = 1500 for ν = 0.03, α = 0.25, γ = 0.5, D = 0.03 and C = 0.05. The initial conditions are the same as the one shown in Fig. 4.
The state including four spirals is rather stable and is maintained even at t = 1000 as shown in Fig. 3(b). Figure 3(c) shows a snapshot pattern of W and n at t = 1000 for ν = 0.04, α = 0.25, γ = 1, D = 0.03 and C = 0.05. The initial 2 2 condition is n(x, y) = Ae−x −y with A = 8/π and W = eiθ , which is the same 8
Fig. 6. (color online) Two snapshot patterns of W and n at (a) t = 2000 (b) t = 10000 for ν = 0.02, α = 0.25, γ = 0.5, D = 0.03 and C = 0.05. The initial conditions are the same as the one shown in Fig. 4.
as the one shown in Fig. 3(b). At ν = 0.04, a regular spiral pattern appears. There are at least two phase singularity points under the periodic boundary conditions. However, the central vortex point becomes dominant and a spiral pattern appears around the central vortex, because the cell concentration n is localized only near the central vortex. At this parameter, the spiral state is stable and the synchronization is attained in the whole system. The center of spiral is a phase singularity point and the phase cannot be defined there. However, the phase can be defined except for the spiral center, and the phase increases toward the center of spiral. It is therefore natural that cells aggregate toward the spiral center. Next, we have performed numerical simulations from an initial condition of n(x, y) = 8/L2 and W = X +iY where X and Y are random numbers between -0.05 and 0.05. The concentration n(x, y) is uniform in this initial condition. Figure 4 shows two snapshot patterns of W and n at (a) t = 400 and (b) t = 1200 for ν = 0.04, α = 0.25, γ = 1, D = 0.03 and C = 0.05. Many spirals are generated owing to the random initial condition of W . Only one spiral survives at t = 1200. In the single complex Ginzburg-Landau equation Eq. (4), many spirals coexist stably for the random initial condition, because the interaction between spiral cores is very weak. In the coupled equations Eqs. (6) and (7), the competition occurs among spiral patterns owing to the frequency difference induced by the cell concentration n(x, y). As a result of the competition, the winner spiral pattern swallows the cells located in the core region of the loser spiral, and the number of spirals decreases. Finally, only one spiral survives. The competitive dynamics among spiral patterns is characteristic in our system Eqs. (6) and (7). We furthermore study a relation of the aggregation dynamics of cells and the complex dynamics of spiral patterns. Figure 5 shows three snapshot patterns of W and n at (a) t = 400 (b) t = 800 and (c) t = 1500 for ν = 0.03, α = 0.25, γ = 0.5, D = 0.03 and C = 0.05. The initial conditions are the same as the one shown in Fig. 4. Many spirals are generated as in the case of Fig. 4, and only one spiral survives at t = 800. However, the one-spiral state is not stationary because of the too fast pace of the spiral core. A vortex pair is created at a point apart from 9
the center of the spiral, and then one vortex collides with the original vortex existing in the center of the spiral pattern, and the vortex pair disappear near the center of the spiral. Finally, a stable target pattern is created at t = 1500. It is because the attractive force of target pattern for the cell aggregation is weaker than that of spiral pattern, and the peak concentration of n is smaller than that for the spiral pattern. Therefore, the frequency of the core region of target pattern is not so high, and the synchronization is attained in the whole system. As a result, the target pattern is stably maintained. However, at a smaller value of ν, spirals breaks up and spiral chaos appears. Figure 6 shows two snapshot patterns at (a) t = 2000 and (b) t = 10000 at ν = 0.02 for the same other parameters α = 0.25, γ = 0.5, D = 0.03 and C = 0.05. As a result of the spiral chaos, the clustering toward the only one cluster cannot occur. Several aggregation clusters survive around spiral cores and the positions of aggregation clusters change randomly owing to the spiral chaos.
4
Summary
We have proposed two-dimensional coupled equations of the phase equation and the equation for cell concentration n to understand the efficient aggregation dynamics of the slime mold qualitatively. Initially many aggregation clusters appears, but the cluster number decreases owing to the competitive dynamics among target patterns, and finally only one localized aggregation cluster survives. The long-range competitive interaction is due to the nonlinear phase dynamics. The survived target pattern accompanying the localized aggregation cluster at the center can be interpreted as a dissipative structure in a nonlinear-nonequilibrium system. In the two-dimensional model for b = 0, the peak amplitude of cell concentration n increases when the parameter γ decreases to 0. The peak amplitude diverges and a singularity appears at γ = 0 in the two-dimensional model. It can be interpreted as a kind of chemotaxis collapse, which is a phenomenon not observed in the one-dimensional model. We have further generalized the model equations to coupled equations of the complex Ginzburg-Landau equation and the equation of cell concentration. Only one spiral pattern survives as a result of the competitive aggregation dynamics, and the cell concentration is localized near the spiral core. We have found more complex dynamical behaviors when the desynchronization occurs around the spiral core. For a strongly unstable case, the clustering toward the only one cluster cannot be realized owing to the breakup of spirals. Then, several clusters coexist and the positions of localized clusters of n change in time randomly owing to spiral chaos. 10
References
[1] Cross, M., Greedside, H.,2009. Pattern Formation and Dynamics in Nonequilibrium Systems (Cambrdge University Press, Cambridge). [2] Kuramoto, Y.,1984. Chemical Oscillations, Waves and Turbulence (SpringerVerlag, New York). [3] Winfree, A. T.,1980. The Geometry of Biological Time (Springer-Verlag, New York). [4] Gross,J., Peacey, M., Trevan, D.,1976. Signal emission and signal propagation during early aggregation in Dictyostelium-Discoideum. J. Cell Sci. 22, 645-656. [5] Martiel,J., Goldbeter, A., 1987. A model based on receptor desensitization for cyclic-AMP signaling in Dictyostelium cells. Biophys. J. 52, 807-828. [6] Tyson, J., Alexander, K., Manoranjan, V., Murray, J., 1987. Spiral waves of cyclic-AMP in a model of slime-mold aggregation. Physica D 34, 193-207. [7] Keller, E., Segel, L., 1970. Initiation of slime mold aggregation viewed as an instability. J. Theor. Biol. 26, 399-415. [8] Nanjundiah, V., 1973. Chemotaxis, signal relaying and aggregation morphology. J. Theor. Biol. 42, 63-105. [9] Van Oss, C., Panfilov, Hogeweg, Siegert, A., Weijer, C., 1996. Spatial pattern formation during aggregation of the slime mold Dictyostelium discoideum. J. Theor. Biol. 181, 203-213. [10] Hillen, T., Painter, K., 2009. A user’s guide to PDE models for chemotaxis. J. Math. Biol. 58, 183-217. [11] Gregor, T., Fujimoto, K., Masaki, N., Sawai, S., 2010. The onset of collective behavior in social amoebae. Science 328, 1021-1025. [12] Garcia-Ojalvo,J., Elowitz, M.,Strogatz, S., 2004. Modeling a synthetic multicellur clock:repressillators coupled by quorum sensing. PNAS 101,1095510960. [13] Winfree, A., 1967. Biological rhythms and behavior of populations of coupled oscillators. J. Theor. Phys. 16, 15-42. [14] Matsumoto, K., Ueda, T., Kobatake. Y., 1988. Reversal of themotaxis with oscillatory stimulation in the plasmodium of Physarum polycephalum. J. Theor. Phys. 131, 175-182. [15] Sakaguchi, H., Maeyama, S., 2013. Dynamical quorum sensing and clustering dynamics in a population of spatially distributed active rotators. Phys. Rev. E 87, 024901.
11
[16] Nasuno, S., Sano, M., Sawada, Y., 1989. Phase wave propagation in the rectangular convective structure of nematic liquid-crystal. J. Phys. Soc. Jpn. 58, 1875-1878. [17] Sakaguchi, H., 1992. Localized oscillation in a cellular pattern. Prog. Theor. Phys. 87, 1049-1053. [18] Kivshar, Y., Agrawal, G., 2003. Optical Solitons (Academic Press, San Diego). [19] Herrero, M., Vel´ azquez, J., 1996. Singularity patterns in a chemotaxis model. Math. Ann. 306, 583-623. [20] Aranson, I., Kramer, L., 2002. The world of the complex Ginzburg-Landau equation. Rev. Mod. Phys. 74, 99-143. [21] Sakaguchi, H., 1993. Phase dynamics and localized solutions to the GinzburgLandau type amplitude equations. Prog Theor. Phys. 89, 1123-1146.
12
