{'fraction_non_alphanumeric': 0.04597354474557284, 'fraction_numerical': 0.01490800982377526, 'mean_word_length': 4.448356807511737, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Wavelength selection of rippling patterns in myxobacteria

Rippling patterns of myxobacteria appear in starving colonies before they aggregate to form fruiting bodies. These periodic traveling cell density waves arise from the coordination of individual cell reversals, resulting from an internal clock regulating them, and from contact signaling during bacterial collisions. Here we revisit a mathematical model of rippling in myxobacteria due to Igoshin et al.\ [Proc. Natl. Acad. Sci. USA {\bf 98}, 14913 (2001) and Phys. Rev. E {\bf 70}, 041911 (2004)]. Bacteria in this model are phase oscillators with an extra internal phase through which they are coupled to a mean-field of oppositely moving bacteria. Previously, patterns for this model were obtained only by numerical methods and it was not possible to find their wavenumber analytically. We derive an evolution equation for the reversal point density that selects the pattern wavenumber in the weak signaling limit, show the validity of the selection rule by solving numerically the model equations and describe other stable patterns in the strong signaling limit. The nonlocal mean-field coupling tends to decohere and confine patterns. Under appropriate circumstances, it can annihilate the patterns leaving a constant density state via a nonequilibrium phase transition reminiscent of destruction of synchronization in the Kuramoto model.


I. INTRODUCTION
Self-organized patterns are ubiquitous in bacterial colonies [1]. Cooperative behavior is often produced by motion following chemical (chemotaxis) or adhesive (haptotaxis) gradients and modeled using reaction-diffusion equations [2]. Instead of using diffusion of chemicals to communicate cells and produce patterns, direct cell contacts yield density waves in myxobacteria such as Myxoccocus xanthus [3][4][5][6][7]. Myxobacteria are rod-shaped gramnegative bacteria, components of soil, that move by gliding in the direction of their long axis on solid surfaces, either as individuals or in swarms [8]. This type of motility of M. xanthus cells is controlled by two different motors concentrated at the cell poles: the A-motility system (adventurous) specific for individual cell motion and the S-motility system (social), for group cell motion. The A-engine is a 'pusher' and works by secretion and hydration of slime, a polyelectrolyte gel, whereas the S-engine is a 'puller' and operates through extension and retraction of type IV pili [9].
When nutrients are abundant, myxobacteria aggregate into multicellular swarms and spread outwards from the border of the colony. However they respond to starvation conditions by reversing their directions, moving towards the colony center and forming a multicellular fruiting body with nonmotile spores. Spores can survive for long periods of time and, under adequate conditions, germinate giving rise to motile vegetative cells [10]. Fruiting body formation is a multi-step process of alignment, rippling, streaming, and aggregation [11]. Overall models of the myxobacteria life cycle depending on food availability can be found in [6,12,13]. Before aggregation begins, bacterial collective behavior and intercellular communication create fascinating collective concentric and spiral traveling waves called ripples. During the rippling phase, bacteria move in concert in such a way that colliding waves appear to pass through one another. The resulting periodic patterns consist of equally spaced ridges of high cell density that appear to advance through the bacterial population as rhythmically traveling waves moving in opposite directions. In these counter-propagating waves, individual cells are parallel to the direction of the ripples, move back and forth and exchange developmental signals (C-signals) only when they collide [14][15][16]. Most experimentally observed rippling patterns can be characterized as counter-propagating traveling waves [3,17]. Unlike waves generated by reaction-diffusion instabilities, which annihilate on collision, myxobacteria waves appear to pass through one another unaffected and are therefore analogous to solitons [4].
There are discrete [12,[18][19][20][21] and continuum [3,4,[22][23][24] theoretical models of cell behavior and cell-to-cell interaction based on C-signaling. Intercellular communication is by direct cell contact, without any diffusible morphogen signaling. Cells migrating in opposite directions can come into contact with each other (collide) and exchange C-signals. Cell state changes periodically so that isolated cells reverse periodically their motion. These changes can be modeled by the concept of an internal clock [3,4,7]. The internal clock or pacemaker regulates swarming that is driven by growth [25]. In some models rippling cannot occur unless there is a refractory period after cell reversal during which a cell does not respond to C-signals from other cells and does not reverse [3,4,18]. Rippling may occur even without a refractory phase in generalized discrete models if the interaction strength and the degree of cooperativity are large enough [23].
In this paper, we consider Igoshin et al's continuum one-dimensional (1D) model for rippling [3,4], solve it numerically and interpret the results by a study of the small nonlinearity limit. An extension of this model to 2D can be found in [5] and related agent-based models in [26]. Previous work has shown that numerical solutions of the continuum model exhibit rippling [4]. Surprisingly, the analysis of the small nonlinearity limit does not provide a selection rule for the observed wavelength of ripples [4]. In contrast with these results, we show that the small nonlinearity limit does provide a description of the rippling instability and it supplies the critical wavenumber of the periodic patterns in terms of the refractory period. In fact, careful consideration of flux continuity at the different stages of one cycle of the internal clock of myxobacteria shows that the limiting equation for the cell density is of Fokker-Planck type but it contains an additional source term. The extra source term is the cell density times a nonlocal growth rate that vanishes as the refractory period tends to zero. The source term is key new element in our analysis, as it produces dissipation even in the absence of diffusion. Thus the role of the source term is similar to that of collision terms in the Boltzmann transport equation. The balance between advection in space and time and dissipation due to the "collision term" selects the wavelength of the rippling patterns that issue forth from the uniform stationary state. These patterns are periodic in time and space.
Numerical solutions of the full model equations show that the pattern frequency decreases as the strength of the nonlinearity increases.
It is interesting to contrast the behavior of the present myxobacteria model with the well-known synchronization of globally coupled phase oscillators described by the Kuramoto model [27,28]. The Kuramoto oscillators move on a circle with their own random natural frequency and their mutual interaction through a mean field tends to synchronize them.
Partial or complete synchronization is achieved for sufficiently strong coupling through a nonequilibrium phase transition. Igoshin et al's model describes phase oscillators with an internal clock whose angular speed changes according to their interaction with oppositely moving oscillators. Since the Igoshin oscillators move with a constant positive or negative velocity, patterns arising from appropriate initial conditions persist in the absence of interaction. In this case, mean-field interaction among oscillators may produce loss of rhythmicity resulting in the destruction of patterns. This change appears as a nonequilibrium phase transition at critical values of parameters that has been explicitly shown in the weak coupling limit.
The rest of the paper is as follows. We describe Igoshin et al's continuum model [4] in Section II and write the corresponding equations in nondimensional form. The weak signaling limit is examined in Section III, where the limiting Fokker-Planck type equation with source term is derived. The linear stability of its uniform stationary solution is analyzed in Section IV. In the absence of diffusion and for disturbances with frequency one (same as that of the signaling solution without nonlinear terms), we find that rippling patterns appear for disturbances with wavenumber less than one, whatever the refractory period. Increasing the wavenumber k may result in the cancellation of patterns, as the uniform stationary solution becomes linearly stable. Since diffusion has a stabilizing role for the uniform solution (producing a negative term proportional to k 2 in the real part of eigenvalues), unavoidable numerical noise should tend to annihilate the periodic patterns and be more effective as

II. MODEL
When starving, myxobacteria aggregate and form a fruiting body and eventually a spore.
Before aggregation begins, there appear periodic patterns of equally spaced high density bands that move as traveling waves (ripples). The one-dimensional (1D) model for ripples of Igoshin et al [3,4] is based on the following experimental observations: -Internal biochemical clock : cells are aligned parallel and glide along their long axis in one direction during part of their internal period and glide in the opposite direction during the other part of their period.
-Contact signaling: a cell collides with an oppositely moving cell and interchanges a signal (a C-protein). As a result, both cells reverse their motion. The collision frequency depends on the local cell density.
-Refractory period : after one collision there is a refractory period during which the cell does not reverse its motion even if it collides again.
Let n(x, φ, t) be the number density of bacteria per unit length x and per unit internal phase φ. Cells with 0 < φ < π move to the right with velocity v and cells with −π < φ < 0 move to the left (velocity −v). Cell division and death are negligible during the rippling phase. The governing equation is as follows: Here D x ≥ 0 is a space diffusivity, the density n(x, φ, t) is a 2π-periodic function of φ and it satisfies periodic boundary conditions at x = ±L or it decays to zero if L = ∞.
To model the angular flux J(x, φ, t), we make precise the above experimental observations [4]. The cell velocity is given by (2) and reversals occur at φ = 0, ±π. The internal clock of a bacterium advances with constant angular velocity ω but when the cell collides with another one moving in the opposite direction, both cells exchange a signal. The collision frequency, and therefore the signaling intensity, is proportional to the local cell density. The cells may respond positively to this signal by accelerating their angular speed from ω to ω + ǫωΩ(N ± ) (see below) depending on their internal phase. The coefficient ǫ measures the relative change in angular velocity from that during the refractory period, ω, to the bacterium velocity during the signaling period, thereby characterizing signaling strength.
After each reversal (at φ = 0, ±π), the cell enters a refractory period during which does not respond to collision signaling and does not reverse. If α > 0 is the angular duration of the refractory period, the angular flux is J = ω n(x, φ, t) for 0 < φ < α and for −π < φ < −π+α where D φ is a positive number and with r > 0. We shall use r = 4 [4]. In Ref. 4, it is explained that both D x and D φ are small.
The total density of bacteria at point x and time t is N + (x, t) + N − (x, t) and the density of Here n(x, φ = 0+, t) and n(x, φ = −π+, t) are, respectively, the density of left-to-right and of right-to-left reversal points in spacetime. Thus their sum, as in (7), is the density of all reversal points in spacetime.
The total number of myxobacterium cells should be independent of time. This means In (1), n and ∂ φ n are continuous and 2π-periodic in φ, so that N is independent of time if Here we have used the boundary conditions at x = ±L (including the case L = ∞) and the and so on. For the angular flux (3), (9) becomes If there is no refractory period so that α = 0, 2π-periodicity of n and ∂ φ n yields instead of (10).  (1), (3) and (8) and dropping hats in the results, we obtain the following equations: where (the dimensionless parameters are assumed to be of order unity): For D x = D φ = 0, (12) resembles a hyperbolic equation (or a system of two hyperbolic equations for oppositely moving bacteria). However, strictly speaking this system is only hyperbolic in the spatial dimension. The phase fluxes (resulting from the nonlinear interaction between oppositely moving bacteria) are nonlocal in the clock angle: they are defined as integrals over the whole angular domain implying that the resulting 2D system (in x and φ) is integro-differential. The role of the nonlocal advection as generator of dissipation will be shown by the analysis of the weakly nonlinear limit and by solving (12) using a high order accurate weighted essentially non-oscillatory (WENO) numerical method [29].

III. WEAK SIGNALING LIMIT
For ǫ = 0, (12) becomes whose solution is [4] Here f (x, t) is an arbitrary function, 2π-periodic in its second argument. f (x, t) and f (x − π, t + π) represent the densities of left-to-right and of right-to-left reversals, respectively.
Then, according to (7), is the reversal point density (RPD) in the weakly nonlinear and weak diffusion limit as ǫ → 0+. This limit is also called the weak signaling limit [4]. In (16), n(x, t, φ) is continuous and 2π-periodic in φ. The , is a particular solution of (12).
Igoshin et al have derived a Fokker-Planck equation for f in the weak signaling limit as ǫ → 0 by using physical arguments and also by singular perturbation methods (see Appendix in [4]). Their derivation missed the collision-type source term that we find in this paper.
Finding this term requires delving more deeply in the perturbation method, therefore we describe this method from scratch. In the limit as ǫ → 0, we seek a solution so that n 1 is 2π-periodic in φ. Note that the function f in (18) has the form (16) with an additional dependence upon the slow time τ . The equation for n 1 is with Continuity of the flux (4), , across angle boundaries and (17) yield From (18) and (19), we obtain along the characteristics the following equation for n 1 : In (26) we have defined Ignoring the initial condition for n 1 , the solution of (26) along the characteristics (25) is where f = f (x 0 , t 0 , τ ), the constants of integration C j are independent of φ, and To determine the constants C j in (28), we impose the jump conditions (21)-(23). Using The integral on the right hand side of (32) can be rewritten as We have used that f (x, t, τ ) is 2π-periodic in t and the change of variable π − ψ → ψ to simplify the integral in (33). The first term in C α−π can be similarly simplified, thereby The condition (24) ensuring 2π-periodicity in φ provides the sought equation for f : acts as an effective collision term that produces dissipation. As we shall see in the next section, the source term provides a mechanism for wave number and speed selection of the ripples.
No such mechanism was found in Ref. 4.

IV. CONSTANT SOLUTION IN THE WEAK SIGNALING LIMIT AND ITS LIN-EAR STABILITY
The constant function f =N/(2π) is an exact solution of (35) that coincides with the following exact piecewise constant solution of the full model (12): when ǫ = 0. In (39), p and q given by (40) have been calculated from N ± = αp + (π − α)q = N /2 and from the condition that the flux be continuous. Note that substitution of f =N /(2π) in (28) yields Let us see what we find by a linear stability analysis of f =N /(2π) as a solution of (35).
Inserting f =N /(2π) + ν(x, t, τ ) with ν ≪ 1 in (35) and keeping only terms that are linear in ν, we find Assuming ν = e ikx+ilt+στ , we obtain the following eigenvalues: where For real k and l, the real and imaginary parts of (43) are respectively. Except for the last three terms in (45), this is the same as (57) in Ref. 4 (witĥ N replaced byN /2) provided l is an integer. We have σ = 0 for k = l = 0 indicating that we As ǫ → 0, f (x, t, τ ) approaches (18), which is a 2π-periodic function of t. The corre-sponding frequency is l = 1. Figure 1 shows the contour plot of the neutral stability curve Re σ(α, k, 1) = 0 in (45) for frequency l = 1 in the absence of diffusion, D x = D φ = 0.
The real part of the eigenvalue (45), Re σ(α, k, 1), is positive at the middle region enclosing k = 0 in Fig. 1 and its sign changes each time a line of the neutral stability curve is crossed.
Thus, for unit frequency and |k| ≤ 1, we expect to see patterns in For other values of the frequency, the neutral stability curve qualitatively changes. For instance, let us consider wave trains traveling to the right, so that k = −l, L = l, K = 0, as in [4]. (45) becomes In absence of diffusion, D x = D φ = 0, (47) shows that Re σ = 0 for l = l c (α) such that 2 sin 2 (αl) = sin 2 (πl). This l c (α) provides the wavelength of the observed patterns. The contour plot of Re σ = 0 in Fig. 2 shows that the function l c (α) is multivalued and that its leftmost branch (corresponding to the lowest values of l > 0) decreases from α = π/ √ 2 at l = 0 to 0 at l = 1. For this parameter range, the constant solution is linearly stable if l < l c (α) and unstable for l > l c (α) (but smaller than values at the next branch of l c (α)). For for 0 ≤ α ≤ π, −6 ≤ l ≤ 6, and D x = D φ = 0. Meaning of "s" and "u" as in Fig. 1.
sufficiently small l, α smaller (larger) than a critical value corresponds to stable (unstable) constant solution. As the exponent r increases and Ω ′ (N /2) correspondingly decreases, the contribution of the second term in (47) becomes smaller and can be obliterated by diffusion or numerical noise even if the said term is positive.

V. NUMERICAL SOLUTIONS
In this section, we present numerical solutions of the full kinetic model equations (1) - (6) and compare them with the linear stability results for the Fokker-Planck type equation (35) derived in the weak signaling limit. To visualize the results, we obtain the total density, , t), and the RPD, n RPD (x, t), from the numerical solution.
We construct the numerical solution using uniform grids for both the spatial variable x Here C is positive constant such that C < 1 to ensure stability of the scheme. We use C = 0.4 in all our computations.
Since the first order scheme introduces more dissipation than needed, we have formulated a high order accurate version in space and time based on the fifth-order weighted essentially nonoscillatory reconstruction procedure (WENO5) for the spatial variables and a thirdorder Runge-Kutta method to evolve in time. WENO5 procedure was designed to get fifth order accuracy in space using a nonlinear convex combination of essentially nonoscillatory parabolas, (see [29] for details). The third order Runge-Kutta method for the integration in time is the one proposed by Shu and Osher in [30] and allows the maximum time stepsize dictated by the first order upwind scheme.
If ǫ = 0 the hyperbolic system (1) Igoshin et al [3] use a heuristic argument to suggest that k c = 1 and estimate ǫ = 3, α = π/5 from their measurements. Different values of ǫ give idea of the solution either close or far from the weak signaling limit: ǫ = 0.1 corresponds to the weak signaling limit and larger values, such as ǫ = 3 [3], or ǫ = 10, go beyond the Fokker-Planck type description of Section IV. To check the stability results of that section, we first solve the nondimensional equation (12) with the following initial condition that is periodic in x and φ: n(x, φ, t = 0) = 1 + 0.1 sin(kx) sin(φ) on the interval [−π, π] × [−π, π]. For a first numerical simulation, we consider the case of ǫ = 0.1 and a fixed α = π/2. In Figure 3, we show the time-space surface plots of the total density with k = 1, k = 4, and k = 6 from left to right. As expected from the neutral stability curve of Fig. 1, the simulations of Fig. 3 show that periodic patterns with angular frequency l = 1 are found for k = 1, whereas the uniform stationary solution is linearly stable if k = 4. Fig. 3 The RPD exhibits the same behavior as the total density and its contour plot (not shown) is quite similar to it. In Figure 4, we depict the RPD at x = 0 for the same values of the parameters. For unit wavenumber, periodic patterns are stable. The slight decrease of the maxima of the RPD observed in Fig. 4(a) is due to unavoidable dissipation due to numerical errors. For k = 4, the uniform stationary solution is stable according to the neutral stability curve of Fig. 1. Fig. 4(b) shows a pronounced evolution towards a constant. However, as we are not far from the bifurcation point (critical wavenumbers bounding the stability region for the uniform stationary solution are 3 and 5), the evolution of the RPD towards a constant value is slow. For k = 6, neutral stability predicts periodic patterns. However at such large wavenumber, the dissipation due to numerical noise is so large that the pattern disappears and the uniform stationary solution becomes stable, as shown in Fig. 4(c).  If the two bumps are closer, then the standing waves generate local maxima as in Fig. 6(a) and (b), for ǫ = 0 and 1.2, respectively. As ǫ (nonlinearity) increases, the reversals occur in smaller regions about x = ±π/2, as shown in Figs. 6(c) and (d). Clearly, for bumps that are close enough, the waves they issue reinforce the density at the points where they cross.
For an initial condition n(x, φ, 0) = 1.5 + 0.1 sin[k(x − φ)], α = π/5, the standing wave patterns as in Figure 4 of [31] that are reminiscent of our patterns in Fig. 7 although the 2D agent-based models the authors of [31] introduce to explain the experiments are more sophisticated than the 1D model we study in this paper. with an internal clock in the limit as the number of bacteria N → ∞. The bacteria themselves satisfy the following equationṡ Here x j (t) and φ j (t) move on circles and can be considered to take values on the intervals [−L, L] and [−π, π], respectively, ν s is a scaling parameter, and the delta functions are regularized in an appropriate way. In the limit as we take away this regularization and N → ∞, the densities (52) approach their continuum limits (5). A typical bacterium x j (t) moves counterclockwise on a circle of length 2L if its internal phase φ j (t) ∈ [−π, π] is positive, and clockwise otherwise. Its internal phase can be accelerated due to a mean-field interaction with opposite moving bacteria that collide with it. Appropriate initial conditions produce periodic (rhythmic) spatiotemporal patterns in the absence of interaction (ǫ = 0).
Interaction tends to confine rhythmicity to parts of the circle (−L, L), as shown in The behavior of this model can be compared with the well-known Kuramoto model of globally coupled phase oscillators [27,28]. In the Kuramoto model, the phases of free oscillators increase following their natural frequencies that are random (unsynchronized or