An age-structured continuum model for myxobacteria

TL;DR

This paper introduces a novel age-structured macroscopic model for myxobacteria movement, derived from microscopic rules, highlighting the importance of a diffusion term for accurate wave behavior and analyzing the role of refractory periods.

Contribution

It presents a new continuum model based on microscopic interactions, incorporating age structure and emphasizing the need for a diffusion term for realistic wave dynamics.

Findings

The macroscopic model aligns well with biological observations after adding diffusion.

Refractory periods are not essential for wave formation but are crucial for synchronization.

Higher order approximations are necessary to include diffusion in the model.

Abstract

Myxobacteria are social bacteria, that can glide in 2D and form counter-propagating, interacting waves. Here we present a novel age-structured, continuous macroscopic model for the movement of myxobacteria. The derivation is based on microscopic interaction rules that can be formulated as a particle-based model and set within the SOH (Self-Organized Hydrodynamics) framework. The strength of this combined approach is that microscopic knowledge or data can be incorporated easily into the particle model, whilst the continuous model allows for easy numerical analysis of the different effects. However we found that the derived macroscopic model lacks a diffusion term in the density equations, which is necessary to control the number of waves, indicating that a higher order approximation during the derivation is crucial. Upon ad-hoc addition of the diffusion term, we found very good agreement between the age-structured model and the biology. In particular we analyzed the influence of a refractory (insensitivity) period following a reversal of movement. Our analysis reveals that the refractory period is not necessary for wave formation, but essential to wave synchronization, indicating separate molecular mechanisms.

Figures

Tables

Table 1. Table 1: Biological Parameters

Name	Meaning	Value	Comment
$v_0$	bacterial speed	$9\mu m/$min	[39]
${\lambda }_m$	spontaneous reversal rate	$0.07/$min	[39, 43, 46]
${\lambda }_M$	maximal reversal rate	$1.5/$min	[39, 43, 46]
$\overline{\rho }$	inflection density	$0.2/\mu m^2$	fitting parameter
$T$	refractory period	2 min	estimates in [26, 7]
$D$	angular diffusion constant	$0.1$/min	$\}$ leads to $d_1=0.99$
$\nu$	alignment frequency	$100/$min
$2m_0$	total average density	$0.25/\mu m^2$	[43]

Table 2. Table 2: Numerical Parameters

Name	Meaning	Value
Simulation parameters of IBM
$L_x,L_y$	width and length of simulation domain	$L_x=L_y=250\mu m$
$N$	total number of bacteria	$2m_0{L}_x{L}_y=15625$
$\Delta t$	time step	$0.01$ min
$R$	interaction radius	$5\mu m$
Simulation parameters of 1D 2-age model
$L$	simulation domain	$250\mu m$
$\Delta x$	spatial step	$0.313\mu m$
$\Delta t$	time step for transport operator	$\Delta x/\left(d_1{v}_0\right)0.95=0.033$ min
$\delta$	diffusion constant	$0.8\mu m^2/$ min

Equations

$$\frac{\mathrm{d}{X}_i}{\mathrm{d}t}=v_{0}v({\Theta }_i(t)),$$

$$\displaystyle\,\mathrm{d}\Theta_{i}=-\nu{\rm Sign}\big{(}\cos(\Theta_{i}-\bar{\Theta}_{i})\big{)}\sin(\Theta_{i}-\bar{\Theta}_{i})\,\mathrm{d}t-4D\cos(\Theta_{i}-\bar{\Theta}_{i})\sin(\Theta_{i}-\bar{\Theta}_{i})\,\mathrm{d}t$$

$$+\sqrt{2D{\cos }^2({\Theta }_i-{\bar{\Theta }}_i)}\mathrm{d}{B}_t^i.$$

$$\displaystyle\left(\begin{array}[]{c}\cos(2\bar{\Theta}_{i})\\ \sin(2\bar{\Theta}_{i})\end{array}\right)=\frac{J_{i}}{|J_{i}|}\qquad\text{ with }\quad J_{i}=\sum_{k:|X_{k}-X_{i}|\leq R}\left(\begin{array}[]{c}\cos(2\Theta_{k})\\ \sin(2\Theta_{k})\end{array}\right).$$

$$s_i(t)=\frac{t-{\tau }_i}{T},$$

$$\displaystyle\left\{\begin{array}[]{ccl}\frac{\,\mathrm{d}s_{i}}{\,\mathrm{d}t}&=\frac{1}{T},&\text{ if }t>\tau_{i}\text{ and }t\text{ is not a reversal time,}\\ s_{i}&=0,&\text{ if }t\text{ is a reversal time.}\end{array}\right.$$

$$\phi(s)=\left\{\begin{array}[]{ll}0,&\text{ if }0\leq s\leq 1\text{ (refractory period)},\\ 1,&\text{ if }s>1\text{ (C-factor sensitive period)}.\end{array}\right.$$

$$\rho^{i}_{\pm}=\frac{1}{|\mathcal{B}_{i}|}\text{Card}\{k\in\mathcal{N}\big{|}X_{k}\in\mathcal{B}_{i}\text{ and }\pm\cos(\Theta_{k}-\bar{\Theta}_{i})\geq 0\},$$

$$\Lambda (\rho ,s)=\lambda (\rho )\varphi (s).$$

$$\text{Pr}\{{\Theta }_i(t+\Delta t)=$$

$$\displaystyle\Lambda(\rho^{i}_{-{\rm Sign}\cos(\Theta_{i}-\bar{\Theta}_{i})},s_{i}(t))\Delta t\,\text{exp}\big{(}-\Lambda(\rho^{i}_{-{\rm Sign}\cos(\Theta_{i}-\bar{\Theta}_{i})},s_{i}(t))\Delta t\big{)}.$$

$${\partial }_{t}f+v_0{\nabla }_x\cdot (v(\theta )f)+\frac{1}{T}{\partial }_{s}f=Q_{\text{al}}^R(f)+Q_{\text{rev}}^R(f).$$

$$Q_{\text{al}}^R(f)={\partial }_{\theta }\left[\nu \mathrm{Sign}\left(\cos (\theta -{\bar{\theta }}_f^R)\right)\sin (\theta -{\bar{\theta }}_f^R)f+D{\cos }^2(\theta -{\bar{\theta }}_f^R){\partial }_{\theta }f\right],$$

$$v(2{\bar{\theta }}_f^R(x,t))=\frac{J_f^R(x,t)}{|J_f^R(x,t)|},J_f^R(x,t)={\int }_0^{\infty }{\int }_{|x-y|\le R}{\int }_{-\pi }^{\pi }v(2\theta )f(y,\theta ,s,t)\mathrm{d}y\mathrm{d}\theta \mathrm{d}s.$$

$${\sigma }_f^{R,\pm }(x,t)=\frac{1}{\pi R^2}{\int }_0^{\infty }{\int }_{\pm \cos (\theta -{\bar{\theta }}_f^R)>0}{\int }_{|x-y|\le R}f(y,\theta ,s,t)\mathrm{d}y\mathrm{d}\theta \mathrm{d}s.$$

$$Q_{\text{rev}}^R(f)=-\left[\Lambda ({\sigma }_f^{R,-},s){\chi }_{\{\cos (\theta -{\bar{\theta }}_f^R)>0\}}+\Lambda ({\sigma }_f^{R,+},s){\chi }_{\{\cos (\theta -{\bar{\theta }}_f^R)<0\}}\right]f(x,\theta ,s,t),$$

$$\displaystyle f(x,\theta,0,t)=T\!\!\int_{0}^{\infty}\!\!\Big{(}$$

$$\displaystyle+\Lambda(\sigma_{f}^{R,+},s^{\prime})\chi_{\{\cos(\theta-\bar{\theta}^{R}_{f})>0\}}\Big{)}f(x,\theta+\pi,s^{\prime},t)\,\mathrm{d}s^{\prime}.$$

$$\hat{f}(\hat{x},\theta ,s,\hat{t})=\frac{f(x,\theta ,s,t)}{(1/x_0^{\prime }{)}^2}\text{and}{\hat{\sigma }}_{\hat{f}}^{\hat{R},\pm }(\hat{x},\hat{t})=\frac{{\sigma }_f^{R,\pm }(x,t)}{(1/x_0^{\prime }{)}^2}.$$

$${\partial }_{\hat{t}}\hat{f}+{\nabla }_{\hat{x}}\cdot (v(\theta )\hat{f})$$

$$\frac{1}{\epsilon }{\partial }_{\theta }\left[\mathrm{Sign}\left(\cos (\theta -{\bar{\theta }}_{\hat{f}}^{\hat{R}})\right)\sin (\theta -{\bar{\theta }}_{\hat{f}}^{\hat{R}})\hat{f}+d{\cos }^2(\theta -{\bar{\theta }}_{\hat{f}}^{\hat{R}}){\partial }_{\theta }\hat{f}\right]$$

$$-\left[\Lambda ({\hat{\sigma }}_{\hat{f}}^{\hat{R},-},s){\chi }_{\{\cos (\theta -{\bar{\theta }}_{\hat{f}}^{\hat{R}})>0\}}+\Lambda ({\sigma }_{\hat{f}}^{\hat{R},+},s){\chi }_{\{\cos (\theta -{\bar{\theta }}_{\hat{f}}^{\hat{R}})<0\}}\right]\hat{f},$$

$$\displaystyle\hat{f}(\hat{x},\theta,0,\hat{t})=\hat{T}\!\!\int_{0}^{\infty}\!\!\Big{(}$$

$$\displaystyle+\Lambda(\hat{\sigma}_{\hat{f}}^{\hat{R},+},s^{\prime})\chi_{\{\cos(\theta-\bar{\theta}^{\hat{R}}_{\hat{f}})>0\}}\Big{)}\hat{f}(\hat{x},\theta+\pi,s^{\prime},\hat{t})\,\mathrm{d}s^{\prime}.$$

$$\displaystyle\varepsilon\big{(}\partial_{t}f^{\varepsilon}+\nabla_{x}\cdot(v(\theta)f^{\varepsilon})+\frac{1}{T}\partial_{s}f^{\varepsilon}\big{)}=Q_{\text{al}}(f^{\varepsilon})+\varepsilon Q_{\text{rev}}(f^{\varepsilon}),$$

$$Q_{\text{al}}(f)={\partial }_{\theta }\left[\mathrm{Sign}\left(\cos (\theta -{\bar{\theta }}_f)\right)\sin (\theta -{\bar{\theta }}_f)f+d{\cos }^2(\theta -{\bar{\theta }}_f){\partial }_{\theta }f\right],$$

$$Q_{\text{rev}}(f)=$$

$$f(x,\theta ,0,t)=$$

$$\displaystyle+\Lambda(\sigma_{f}^{+}(x),s^{\prime})\chi_{\{\cos(\theta-\bar{\theta}_{f})>0\}}\Big{)}f(x,\theta+\pi,s^{\prime},t)\,\mathrm{d}s^{\prime}.$$

$$v(2{\bar{\theta }}_f(x,t))=\frac{J_f(x,t)}{|J_f(x,t)|},\text{with }{J}_f(x,t)={\int }_0^{\infty }{\int }_{-\pi }^{\pi }v(2\theta )f(x,\theta ,s,t)\mathrm{d}\theta \mathrm{d}s,$$

$${\sigma }_f^{\pm }(x,t)={\int }_0^{\infty }{\int }_{\pm \cos (\theta -{\bar{\theta }}_f)>0}f(x,\theta ,s,t)\mathrm{d}\theta \mathrm{d}s.$$

$$\left\{{\bar{f}}_{{\rho }_{+}(s),{\rho }_{-}(s),\bar{\theta }}(\theta )|{\rho }_{\pm }:[0,\infty )\to [0,\infty ),\bar{\theta }\in [0,\pi )\right\},$$

$${\bar{f}}_{{\rho }_{+},{\rho }_{-},\bar{\theta }}(\theta )=\{\begin{array}{l}{\rho }_{+}{M}_{\bar{\theta }}(\theta )\text{for}\cos (\theta -\bar{\theta })>0\\ {\rho }_{-}{M}_{\bar{\theta }}(\theta )\text{for}\cos (\theta -\bar{\theta })<0.\end{array}$$

$$M_{\bar{\theta }}(\theta )=\frac{1}{Z_d}\exp \left(-\frac{1}{d|\cos (\theta -\bar{\theta })|}\right),\theta \in [-\pi ,\pi )$$

Full text

1. Department of Mathematics, Imperial College London

London, SW7 2AZ, United Kingdom

pdegond@imperial.ac.uk

2. Courant Institute of Mathematical Sciences,

New York University, 251 Mercer Street,

New York, NY 10012, USA

angelika.manhart@cims.nyu.edu

3. Yau Mathematical Sciences Center,

Jin Chun Yuan West Building, Tsinghua University,

100084 Beijing, China

huiyu@tsinghua.edu.cn

**Key words: ** Self-propelled particles; nematic alignment; hydrodynamic limit; generalized collision invariant; diffusion correction; myxobacteria; wave formation; refractory period.

**AMS Subject classification: ** 35L60, 35K55, 35Q70, 82C05, 82C22, 82C70, 92D50

1 Introduction

Myxobacteria are a fascinating example for how simple cell-cell interaction rules can lead to emergent, collective behavior. These single-celled organisms have the ability to move on two dimensional surfaces and form large colonies. When swarming, the colony exists as a rather uniform mono- or multi-layer of densely packed cells with single cells occasionally venturing away from the main swarm body. Myxobacterial swarms are predatory, searching and killing prey as a collective, which is one reason why these bacteria have often been called social bacteria [21]. Upon meeting prey, but also under starvation conditions the cells enter a ripple phase, during which periodic density waves are formed [42]. When two waves traveling in opposite directions collide, the waves appear to pass through each other unaffectedly. However, by tracking individual bacteria [39, 46] it was discovered that most cells in the wave crests in fact reverse their direction of movement, showing that the density waves are actually being reflected off each other. Myxobacteria reverse without turning, by internally exchanging the lagging and the leading pole. Isolated bacteria reverse spontaneously (on average every 10 min), however, their reversal rate increases as a response to higher densities of other bacteria around them. Although the precise function of rippling is not known, it often serves as a prelude and also overlaps with an aggregation phase: in this developmental stage bacteria aggregate into several growing mounds which eventually rise out of the plane and form large three dimensional fruiting bodies. Both waves and aggregates are macroscopic structures with typical length scales of 100 $\mu m$, whereas individual bacteria are only a few microns long. Biologically this makes myxobacteria an interesting and suitable research object for understanding the development of multicellular cooperation, the basis of all complex life forms. Finally, the myxobacteria’s unique metabolites have also rendered them an attractive source for potential new drugs [38].

The various social and cooperative behavior in observed myxobacterial colonies raises questions about the mechanisms of cell-cell communication. The most important mechanism responsible for inducing both ripple formation and aggregation has been found to be C-signaling: The C-factor is a 17-kD protein associated with the cell surface. It has been shown that direct cell-cell contact is necessary for C-signaling and that the exchange is facilitated via end-to-end contacts[29]. Isolated cells exposed to purified C-factor show an increase in reversal frequency [39], suggesting that cell-to-cell contacts increase the probability for a cell to reverse.

In this work we try to shed light on some of the questions associated with ripple formation: What primary effect of C-signaling causes a uniformly spread swarm to start forming ripples? Are density-dependent changes in reversal frequency enough to explain the formation of opposing, periodic wave trains? One idea brought forward in [26] and inspired by D. discoideum is that of an insensitivity or refractory period. Based on the observation that there seems to be a minimum time of around 40 sec [46] between two reversals of the same bacterium, it is suggested that bacteria become insensitive to C-signaling immediately after they have reversed. Using mathematical modeling we show that a refractory period is not necessary for the formation of traveling waves as such, but rather for controlling the width of the waves as well as their wavelength. By analyzing the composition of waves in terms of insensitive and sensitive cells, we discover a possible mechanism how periodic waves are created and maintained in myxobacterial colonies.

While mutation experiments provide valuable insight, computational models offer a powerful alternative to test and analyze different mechanistic biological models. Detailed measurements and statistics about single cell behavior [39, 46] as well as mutation experiments provide the quantitative data necessary to formulate, parametrize and validate mathematical models. In many cases when modeling biological or physical systems one of the first modeling decisions is whether to use an individual- or particle-based model (IBM), in which the individual agents (in our case bacteria) interact by simple rules or to use a continuum model, in which the evolution of macroscopic quantities such as densities or mean directions are described by differential equations [32]. Advantages of IBMs are that they generally allow for an easy incorporation of biological knowledge or hypotheses and can deal with noise in a straightforward way. However the analysis of the model is often limited to running a large number of simulations and little mechanistic insight is gained. For (macroscopic) differential equations on the other hand a large analytical toolbox ranging from asymptotic methods to linear stability analysis and bifurcation theory is available, which can produce precise results about the parameter dependence of solution behavior, etc. However for biological systems it is often difficult to derive continuous models, in many cases ad hoc models are used in which some desired system behavior is already built into the derivation, thereby limiting the explanatory potential of the model.

The derivation of macroscopic models from IBMs of collective dynamics has been the subject of an intense literature with a particular focus on applications to biology. This derivation proceeds through an intermediate modelling level called the kinetic or mean-field model [4, 5, 8]. The derivation from kinetic to macroscopic models of collective dynamics faces the problem of the lack of conservation relations (such as the lack of momentum conservation, see e.g. the review in [45]). A recent breakthrough is the so-called generalized collision invariant concept [18, 23] subsequently developed in a variety of biological contexts (tissue self-organization [13], flocking [14]) as well as physical (micromagnetism [15]) or social (economics [16]) contexts. Other models of collective dynamics can be found in [1, 2, 3, 9, 11, 12, 22, 31, 33, 37, 44]. In most works within the SOH framework, the derivation involves an expansion of the collision operator in terms of $\varepsilon$ – the ratio between the microscopic and the macroscopic scale – and a (formal) limit of $\varepsilon\rightarrow 0$. One effect of this approximation is the absence of any mass diffusion. For the myxobacteria model we found that this approach leads to a discrepancy between the particle model, in which the number of waves are controlled and the diffusion-free continuous model, in which they are not. Hence the $\mathcal{O}(\varepsilon)$ terms describing small, but finite effects seem to be crucial. In this work we use the IBM to estimate the correct size of the mass diffusion term, however future work will include a formal derivation using the Chapman-Enskog expansion (applied to an SOH model in [20]).

In this paper, to model the refractory period between two reversals of myxobacteria, we use the concept of a local time that is reset to zero after each reversal. This idea is borrowed from similar ideas used in neuron dynamics (see e.g. [36]) itself being a variant of the renewal equation used to model the cell cycle [25].

The rest of the manuscript is organized as follows: In Sec. 2 we present a hierarchy of models starting with an Individual Based Model (IBM, Sec. 2.1) and systematically deriving a PDE-based macroscopic model. The final, macroscopic 2-age model is presented in Sec. 2.3. In Sec. 3 we perform detailed simulations: first we validate the macroscopic model by comparison with the IBM (Sec. 3.1-3.2), then we further analyze the properties of the 2-age model, in particular with respect to experimental predictions and wave formation. Concluding remarks are found in Sec. 4.

2 Model Presentation

We present a model hierarchy at the individual and macroscopic levels. We start with an individual-based model (Sec. 2.1) describing the position and velocity of each bacterium as well as its internal biochemical age variable, which can be interpreted as memory of the bacterium. Here age refers to the time passed since the last reversal. From there we systematically derive an age-structured continuous model for the macroscopic quantities, density and nematic mean direction (Sec. 2.2) where the age variable is still continuous. As a last step we discretize the age variable and assume only two groups or ages: being in the refractory period (i.e. being insensitive to C-signaling) or being sensitive to C-signaling. This results in a macroscopic 2-age model (Sec. 2.3), that forms the basis of the subsequent analysis.

Model Assumptions. Fig. 1 shows the main model ingredients and assumptions:

M1: Bacteria move in 2D with constant speed, in the direction of their orientation. This orientation is subject to random noise (Fig. 1.B).

M2: Bacteria align nematically with other bacteria within their immediate vicinity (Fig. 1.A).

M3: Bacteria can reverse their orientation. Their reversal rate is a function of the local density of opposing bacteria (Fig. 1.C1/C3).

M4: After a reversal, bacteria go through a refratory period of fixed length, denoted by $T$, during which they cannot reverse (Fig. 1.D/C2).

M2 can be interpreted as the effect of physical interactions between hard rods and models size exclusion effects and steric interactions. M3 is a consequence of the contact-dependence of C-signaling, which we assume to act only over very short distances.

2.1 The Individual Based Model (IBM)

We describe the movement of $N$ individual myxobacteria. For $i\in\mathcal{N}:=\{1,..,N\}$, the $i$-th bacterium at time $t>0$ is characterized by its center of mass $X_{i}(t)\in\mathbb{R}^{2}$, its orientation angle $\Theta_{i}(t)\in[-\pi,\pi)$ (defined modulo $2\pi$) and an age variable $s_{i}(t)\geq 0$, related to the time since the bacterium’s last reversal.

Movement and Alignment.

Bacteria move with constant speed $v_{0}>0$ in the direction $v(\Theta_{i}(t)):=(\cos\Theta_{i}(t),\sin\Theta_{i}(t))^{T}$. To model nematic alignment, we follow [17] and assume that the following stochastic differential equations govern the evolution of $X_{i}(t)$ and $\Theta_{i}(t)$:

[TABLE]

Remark 2.1

Note that the term $4D\cos(\Theta_{i}-\bar{\Theta}_{i})\sin(\Theta_{i}-\bar{\Theta}_{i})\,\mathrm{d}t$ presents a subtle difference to the model presented in [17] and arises from a different interpretation of the stochastic differential equations (SDEs): Using the usual Ito convention, this term is necessary to be consistent with both the numerical implementation presented in App. A.1 and the Fokker-Planck equation derived in Sec. 2.2.1. In [17] the SDEs were interpreted in the Backward Ito sense (also called isothermal convention [30]), and hence formulated without this term.

The parameter $\nu>0$ measures the alignment frequency to the local mean direction $\bar{\Theta}_{i}$. Bacteria either align with it (if $\cos(\Theta_{i}-\bar{\Theta}_{i})>0$) or against it (if $\cos(\Theta_{i}-\bar{\Theta}_{i})<0$). ${\rm Sign}$ is defined by ${\rm Sign(a)}=\pm 1$ for $\pm a>0$. The nematic mean direction $\bar{\Theta}_{i}$ can be understood as an average of lines going through each bacterium. It is defined by

[TABLE]

$J_{i}$ represents the nematic mean current and the parameter $R>0$ specifies the interaction range of the alignment. Since it is an average of lines (as opposed to an average of angles) $\bar{\Theta}_{i}$ has to be understood modulo $\pi$ and we always choose $\bar{\Theta}_{i}\in[0,\pi)$. A more detailed discussion can be found in [17]. The angular noise is modeled by the stochastic process $\,\mathrm{d}B_{t}^{i}$, describing independent Browinian motion of intensity $D\cos^{2}(\Theta_{i}-\bar{\Theta}_{i})$, where $D>0$. The term $\cos^{2}(\Theta_{i}-\bar{\Theta}_{i})$ aids the separation into two groups of bacteria traveling in opposite direction and is described in detail in [17]. This concludes the modeling of assumptions M1 and M2.

Reversals and Insensitivity.

To model assumptions M3 and M4 we start by noting that bacteria can reverse their orientation, which changes $\Theta_{i}$ to $\Theta_{i}+\pi$. The reversal frequency depends on physical contact during which the signaling molecule C-factor is exchanged. We assume that immediately after a reversal, bacteria go through a refractory period of length $T$, during which they are insensitive to C-signaling. We therefore endow each bacterium with an age variable $s_{i}(t)>0$ which measures the time elapsed since its last reversal, normalized by the refractory period $T$:

[TABLE]

where $\tau_{i}$ records the time of the last reversal for the $i$-th particle and $s_{i}(t)$ is between [math] and $\tau^{\prime}_{i}-\tau_{i}$, where $\tau^{\prime}_{i}$ is the instant of the next reversal. Therefore, the dynamics of $s_{i}$ between two instants of reversal is given by

[TABLE]

Eqs. (2.2) and (2.5) are equivalent. While the former gives an explicit expression of the age variable $s_{i}$, the latter will be used in the formulation of the kinetic equation (2.8) below. We model C-signaling sensitivity by a step function $\phi(s)\in\{0,1\}$:

[TABLE]

To model M3, we assume that an individual bacterium’s reversal rate is a function of the local density of bacteria oriented opposite to it. To that end, we define the local densities $\rho_{\pm}^{i}$ in the $i$-th bacterium’s neighborhood $\mathcal{B}_{i}(X)=\{X:|X-X_{i}|\leq R\}$ as

[TABLE]

where ‘Card’ is the cardinal of a set and $|\mathcal{B}_{i}|$ the area of $\mathcal{B}_{i}$. Since both are related to physical contact, we choose the interaction radii of density sensing and alignment to be equal, however in general they could be different. The subscripts $\pm$ indicate whether the density refers to bacteria moving with ($+$) or against ($-$) the $i$-th particle. The reversal frequency as a function of density of opposing bacteria is denoted by $\lambda(\rho)$. Since higher concentration of opposing bacteria have been observed to cause more frequent reversals, we assume $\lambda(\rho)$ to be an increasing function of $\rho$. More discussion is provided in Sec. 2.4.

The total reversal function $\Lambda(\rho,s)$ takes into account both of the above factors, namely the density and the refractory period, and is defined by

[TABLE]

Finally, the probability that the $i$-th particle has reversed between $t$ and $t+\Delta t$ is modeled as a Poisson process:

[TABLE]

This completes the description of the Individual Based Model (IBM). Numerical results are presented in Sec. 3.1.

2.2 The Macroscopic Continuous-Age Model

The IBM presented in Sec. 2.1 consists of $4N$ coupled stochastic differential equations for typical bacterial colony sizes of $N\approx 10^{6}$. We therefore derive a macroscopic model that consists of only three partial differential equations. The derivation strategy uses the Self-Organized Hydrodynamics (SOH) framework, which allows a systematic derivation of hydrodynamic equations for particle systems that do not have enough conserved quantities, a common obstacle in biological systems. The derivation follows [17] and is described in more detail below in Sec. 2.2.1; however its structure can be summarized as follows. First a mean field model is derived, which leads to a Fokker-Planck equation for the 1-particle distribution function. A hydrodynamic scaling introduces a small parameter representing the difference in microscopic and macroscopic time and spatial scales. Taking this parameter to zero, one finds that the equilibrium distribution function is locally characterized by three quantities: The nematic mean direction $\bar{\theta}$, the density of particles aligned with it, $\rho_{+}$ and that of particles anti-aligned with it, $\rho_{-}$. The macroscopic model describes how these three quantities change in time and space and how the densities depend on the bacteria’s biochemical age.

2.2.1 Derivation of the Macroscopic Continuous-Age Model (2.19)

The particle model presented in Sec. 2.1 serves as the starting point for the derivation. We largely follow [17] and emphasize differences where appropriate.

Kinetic Equation.

Following the classical strategy for mean field models, presented e.g. in [18], we let the number of particles $N$ tend to infinity. Then the distribution function $f(x,\theta,s,t)$ satisfies the following Kolmogorov-Fokker-Planck type equation:

[TABLE]

We note in particular that the term $\frac{1}{T}\partial_{s}f$ comes from the age equation (2.5). We recall that $R$ denotes the interaction radius for both reversals and alignment. Here and in the following the superindex $R$ is used to emphasize the non-locality of the corresponding terms. The collision operator $Q^{R}_{\text{al}}$, caused by the alignment is almost identical to that derived in [17] and is given by

[TABLE]

where the nematic mean direction $\bar{\theta}^{R}_{f}$ is defined by

[TABLE]

Note that the definition of the mean nematic current $J_{f}^{R}$ requires integrating over all ages $s$, a small difference to the operator of the age-free model defined in [18], which reflects the fact that the biochemical age of a bacterium does not influence the alignment. The first term in $Q^{R}_{\text{al}}$ is a drift term in $\theta$ that moves the mass towards $\bar{\theta}^{R}_{f}$ and $\bar{\theta}^{R}_{f}-\pi$. The second term causes diffusion in $\theta$ with magnitude $D\cos^{2}(\theta-\bar{\theta}_{f}^{R})$.

The new operator $Q^{R}_{\text{rev}}(f)$, describes the reversals and is defined as follows. First we introduce two quantities $\sigma_{f}^{R,\pm}(x,t)$ representing the position-dependent, total densities of each group,

[TABLE]

Note that these quantities are densities, i.e. number of particles per unit area, hence the factor $\pi R^{2}$ at the denominator, representing the area of the set $\{y\,\,|\,\,|x-y|\leq R\}$. Then we have

[TABLE]

where $\chi_{\mathcal{S}}$ is the characteristic function on the set $\mathcal{S}$ and $\Lambda(\sigma,s)$ is defined in (2.7). Eq. (2.9) is nonzero only for $s>1$, since only then particles are sensitive to C-signaling and can reverse with a frequency that depends on the density of the opposing group. Note that (2.9) represents particles reversing away from their group. Indeed, bacteria having reversed only enter the balance at age $s=0$. Away from $s=0$ there is no influx of (turned) bacteria, which explains why in (2.9) the loss term appears alone. The contribution of bacteria having reversed is accounted for as a boundary conditions at age $s=0$ given by

[TABLE]

The integral reflects the fact that reversing particles of all ages that point in direction $\theta+\pi$, will add mass to the distribution function at angle $\theta$ and age $s=0$, since reversing resets the biochemical age to zero. An alternative formulation - which we discard here - would be to add a Dirac delta at $s=0$ multiplied by the right-hand side of (2.10) (divided by $T$) to the expression of (2.9). Both approaches are equivalent. In particular, the correctness of the mass balance can be checked by integration of (2.8) with respect to $s$ on $[0,\infty)$.

Scaling.

Analogous to [17], we perform the nondimensionalization and the hydrodynamic scaling in one step. On the microscopic scale the reference time and space units are given by $t_{0}=1/\nu$ and $x_{0}=v_{0}t_{0}$. The age variable $s$ remains unchanged since it is already dimensionless. The scaled diffusion constant is $d=D\,t_{0}$. On the macroscopic scale we use the coarse units $t_{0}^{\prime}=t_{0}/\varepsilon,x_{0}^{\prime}=x_{0}/\varepsilon$ where $\varepsilon>0$ is some small real number. Then the dimensionless macroscopic variables are $\hat{t}=\frac{t}{t_{0}^{\prime}},\hat{x}=\frac{x}{x_{0}^{\prime}}$. Further we set $\hat{R}=\frac{R}{x_{0}^{\prime}}$. The scaled distribution function $\hat{f}$ and densities $\hat{\sigma}_{f}^{\hat{R},\pm}$ are given by

[TABLE]

For the reversal term, we set $T=\hat{T}t_{0}^{\prime}$ and $\hat{\Lambda}(\hat{\sigma}_{\hat{f}}^{\hat{R},\pm},s)=\Lambda(\sigma^{R,\pm}_{f},s)t_{0}^{\prime}$. Note that $\bar{\theta}^{\hat{R}}_{\hat{f}}(\hat{x},\hat{t})=\bar{\theta}^{R}_{f}(x,t)$. The scaled equation reads

[TABLE]

where the second and third line represent the scaled alignment and reversal operators respectively. The scaled boundary condition at $s=0$ is given by

[TABLE]

At this point the definitions of the nematic mean direction and the densities $\hat{\sigma}_{\hat{f}}^{\hat{R},\pm}$ still involve space integrals, i.e. they are non-local. We assume purely local interactions for both alignment and reversals and therefore set $\hat{R}=\varepsilon r$ with $r=\mathcal{O}(1)$. Then Taylor expansion of $\bar{\theta}^{\varepsilon r}_{\hat{f}}$ and $\hat{\sigma}_{\hat{f}}^{\varepsilon r,\pm}$ around $\varepsilon=0$ shows that the functions can be approximated by the local-in-space functions $\bar{\theta}_{\hat{f}}$ and $\hat{\sigma}_{\hat{f}}^{\pm}$ respectively (see (2.13), (2.14)) with a remainder of order $\mathcal{O}(\varepsilon^{2})$. In the following we drop the hats for better readability and call the solution $f^{\varepsilon}$ to emphasize its dependence on $\varepsilon$.

Hyrodynamic Limit.

To derive the mean field equation, we need to find the solution $f^{\varepsilon}(x,\theta,s,t)$ as $\varepsilon\to 0$ in

[TABLE]

where

[TABLE]

and

[TABLE]

supplemented by the boundary condition at $s=0$

[TABLE]

The mean nematic direction is defined by

[TABLE]

and the local mass functions are

[TABLE]

Taking the hydrodynamic limit in the SOH framework [18] now involves two steps: (i) Characterizing the kernel of $Q_{\text{al}}(f)$ ($\theta$ dependence) and (ii) Using Generalized Collision Invariances (GCIs) to extract information about the $x$ and $t$ dependence from the transport and reversal terms.

Step (i): Since $Q_{\text{al}}(f)$ is similar to the collision operator analyzed in [17], but for the addition of the age variable, we simply state the result without proof:

Lemma 2.1

The kernel of $Q_{\text{al}}$ is given by

[TABLE]

where

[TABLE]

$M_{\bar{\theta}}(\theta)$* describes the Generalized von Mises (GVM) distribution defined by*

[TABLE]

The equilibria have the shape of two opposing maxima: One in direction $\bar{\theta}$ with mass $\rho_{+}(s)$ and one in direction $\bar{\theta}-\pi$ with mass $\rho_{-}(s)$ and are depicted in Fig. 2.

Step (ii): In a classical setting, one would at this point mutliply equation (2.11) with collision invariants (CIs) and integrate over all directions $\theta$ and all ages $s$. A CI is defined as a function $\Psi(\theta,s)$ on $[-\pi,\pi)\times[0,\infty)$, such that

[TABLE]

This would allow us to remove the term of order one from (2.11) and thereby yield three macroscopic equations that describe the constant $\bar{\theta}$ and the two functions $\rho_{\pm}(s)$ that characterize the equilibria defined in (2.15). However, the operator $Q_{\text{al}}$ has only one family of CIs: $\Psi(\theta,s)\equiv\psi_{0}(s)$, for an arbitrary function $\psi_{0}$ which represents mass conservation of particles of a given age $s$ (since alignment doesn’t modify the age). This necessitates the use of Generalized Collision Invariants (GCIs). Given an angle of lines $\bar{\theta}$ a GCI associated with $\bar{\theta}$ is a function $\Psi_{\bar{\theta}}(\theta,s)$ defined on $[-\pi,\pi)\times[0,\infty)$, such that

[TABLE]

where $Q_{\text{al},\bar{\theta}}$ is defined analogously to $Q_{\text{al}}$ in (2.12), but with $\bar{\theta}_{f}$ replaced by $\bar{\theta}$. A GCI is a more general concept than CIs. If enough families of GCIs can be found, Step (ii) of the derivation now requires integrating equation (2.11) against these GCIs associated to $\bar{\theta}_{f^{\varepsilon}}$. Similar to the classical CI-based approach, this removes the leading order singular term and allows us to derive the macroscopic equations. We omit the precise functional analytical setting as it can be easily extended from [17] and simply state the following result:

Lemma 2.2

Given an angle of lines $\bar{\theta}$, the space of GCIs of $Q_{\text{al}}$ associated with $\bar{\theta}$ is spanned by $\psi_{\bar{\theta}}^{\pm}(\theta,s):=\text{H}(\pm\cos{(\theta-\bar{\theta})})\psi_{0}(s)$ and $g_{\bar{\theta}}(\theta,s):=g(\theta-\bar{\theta})$ where H denotes the Heaviside function,

[TABLE]

which is extended to $[-\pi,\pi)$ by $g(-\theta)=-g(\theta)$ and $g(\pi-\theta)=-g(\theta)$, and $\psi_{0}(s)$ is an arbitrary functions of $s$.

Note that $\psi_{\bar{\theta}}^{+}(\theta,s)+\psi_{\bar{\theta}}^{-}(\theta,s)\equiv\psi_{0}(s)$ and we recover the CI associated with mass conservation. Now we proceed as explained, by integrating (2.11) against the three GCIs associated to $\bar{\theta}_{f^{\varepsilon}}$.

Proceeding similar to [17] finally yields

Proposition 2.3

Taking the (formal) limit $\varepsilon\rightarrow 0$ in (2.11) we obtain

[TABLE]

where $\bar{f}_{\rho_{+}(x,s,t),\rho_{-}(x,s,t),\bar{\theta}(x,t)}(\theta)$ is given by (2.15) and the macroscopic quantities $\rho_{\pm}(x,s,t)$ and $\bar{\theta}(x,t)$ are such that $\rho_{\pm}(x,s,t)\in[0,\infty)$ and $\bar{\theta}(x,t)$ is a real number defined modulo $\pi$, and fulfill

[TABLE]

supplemented by the boundary conditions

[TABLE]

where the coefficients $d_{1}$, $d_{2}$ and $\mu$ are given by

[TABLE]

$v(\bar{\theta})^{\perp}=(-\sin(\bar{\theta}),\cos(\bar{\theta}))^{T}$* and $\langle\phi\rangle_{M}$ represents the average with respect to $M(\theta)=M_{0}(\theta)$:*

[TABLE]

Finally we remove the non-dimensionalization and revert system (2.16)-(2.17) back to physical units yielding (2.19)-(2.20) (below). The next section summarizes the model.

2.2.2 The macroscopic continuous-age model.

We denote the spatial variable by $x\in\mathbb{R}^{2}$, the age variable by $s\in[0,\infty)$ and time by $t>0$. Then $\bar{\theta}(x,t)\in[0,\pi)$ describes the local nematic mean direction, which is independent of $s$. We recall that $v(\bar{\theta}):=\left(\cos(\bar{\theta}),\sin(\bar{\theta})\right)^{T}$ and define $v^{\perp}$ as its left-oriented orthogonal. We denote by $\rho_{+}(x,s,t)$ and $\rho_{-}(x,s,t)$ the local densities of bacteria of age $s$ that are transported in the directions $v(\bar{\theta})$ and $-v(\bar{\theta})$ respectively. $\bar{\theta}(x,t)$, $\rho_{\pm}(x,s,t)$ fulfill the following system of equations:

[TABLE]

supplemented by boundary conditions at the age $s=0$,

[TABLE]

where the coefficients $d_{1}$, $d_{2}$ and $\mu$ are given by (2.18), $\Lambda(\sigma,s)$ is defined in (2.7) and $\sigma_{\pm}(x,t)$ are the local masses of the two opposing groups:

[TABLE]

Fig. 3(a) illustrates the reversal- and age-related dynamics of (2.19)-(2.20): For $s\in[0,1]$, bacteria are in the refractory period and can not reverse, ensured by $\Lambda(\sigma_{\mp},s)=0$. For $s>1$, bacteria enter the C-factor sensitive period and the reversal frequency is governed by the reversal function $\Lambda(\sigma_{\mp},s)=\lambda(\sigma_{\mp})$ (using (2.6) and (2.7)). Assuming that there are in fact only two essential age states for the age variable $s$ allows to simplify (2.19) to remove the age $s$ as independent variable. This is done in the following section.

2.3 The 2-Age Model

To arrive at an easy-to-handle, yet powerful macroscopic model, we perform one last simplification step: we assume only two age group of bacteria, those who can, and those who cannot reverse. The main difference to the continuous-age model is that the aging itself is now described as a simple reaction term with rate $1/T$. Fig. 3(b) depicts the corresponding reaction schematic. The time evolution of the nematic mean direction remains unchanged. The derivation is given below in Sec. 2.3.1.

2.3.1 Derivation of the 2-age model

An age-discretized macro-model with $K+1$ age groups.

As a starting point we use the system of equations (2.19) together with the boundary conditions (2.20), i.e. the macroscopic myxobacteria model with a continuous age variable. To derive the corresponding discrete age system we discretize the age variable $s$ by $s_{k}=k\Delta s$ for $k=0,\ldots,K$, yielding $K+1$ age groups defined by

[TABLE]

To get equations for $\tilde{\rho}_{\pm}^{k}(x,t)$ we use a forward difference discretization of $\partial_{s}\rho_{\pm}$. Since the equation for $\bar{\theta}$ is independent of $s$, we only need to consider the two density equations. The only noteworthy point is that the largest age group i.e. $\tilde{\rho}^{K}_{\pm}$ can lose particles only by reversing, not by aging. The corresponding system for a discrete age system with $K+1$ age groups is

[TABLE]

(and analogously for $\tilde{\rho}_{-}^{k}$). The boundary conditions are included by using “virtual” age groups defined by discretizing the integrals in (2.20), yielding

[TABLE]

The total group densities $\tilde{\sigma}_{\pm}(x,t)$ are defined by

[TABLE]

The 2-age macro-model.

To obtain the 2-age model (3.25), we define

[TABLE]

in (2.3.1) and assume that $\phi(s_{0})=0$ and $\phi(s_{k})=1$ for $k\geq 1$. This yields a closed system for $(\rho^{0}_{\pm},\rho^{1}_{\pm})$. For better readability we drop the tilde on the expressions for the total densities $\sigma_{\pm}$.

[TABLE]

Setting $\Delta s=1$ we get the final system (2.22), below. In Sec. 3 we demonstrate that this 2-age model is a good approximation of the full age-dependent dynamics and sufficient to reproduce and explain almost all experimentally observed features of myxobacteria.

2.3.2 The 2-age model

The final 2-age model derived in Sec. 2.3.1 consists of five equations for the four densities $\rho_{\pm}^{0,1}(x,t)$ and the nematic mean direction $\bar{\theta}(x,t)$. Subindices $\pm$ denote densities of bacteria moving with and against the nematic mean direction respectively and superindices $0/1$ denote bacteria in a non-reversible and reversible state respectively.

[TABLE]

[TABLE]

where the local masses are given by $\sigma_{\pm}=\rho_{\pm}^{0}+\rho_{\pm}^{1}$. The constants $d_{1}$, $d_{2}$ and $\mu$ are defined analogously as for (2.19). Note that this model does not contain any spatial diffusion, which we will discuss further in Sec. 3.2.

2.4 The Reversal Frequency $\lambda(\rho)$

To complete the models presented in Secs. 2.1-2.3, the density dependence of the reversal frequency $\lambda(\rho)$ has to be specified based on the available information from experiments. Firstly, in most experiments the number of reversals increase with the density of opposing bacteria, i.e. $\lambda^{\prime}(\rho)\geq 0$. Next, in the absence of other bacteria isolated myxobacteria still reverse, i.e. $\lambda(0)=:\lambda_{m}>0$. Spontaneous reversal rates between $0.07$ and $0.09$ reversals per minute have been reported [39, 43, 46]. Thirdly, there seems to be an upper limit as to how frequent reversals can be, which confirms the biological intuition that the rearrangement of the internal movement machinery takes some time. In [46] a maximal rate of $1.5$ reversals per minute have been observed. Finally, in [39] experiments were performed in which the reversal rate of isolated bacteria was measured in response to externally adding the signaling molecule C-factor, which was thought to communicate the density information. At low concentration of C-factor the reversal rate remained the same, while with increasing concentration a growth in reversal rate was observed, which plateaued for very high concentrations of C-factor. Using these pieces of information we assume a sigmoid shape of $\lambda(\rho)$. As a convenient representation we use a $C^{1}$ and piecewise smooth function:

[TABLE]

Note that $\lambda(\rho)$ is parameterized by three quantities, the spontaneous reversal rate $\lambda_{m}$, the maximal reversal rate $\lambda_{M}$ and the inflection density $\bar{\rho}$, at which $\lambda(\rho)$ grows the fastest.

3 Numerical Analysis - Comparison to Experiments

Depending on the precise experimental set-up, various bacterial speeds $v_{0}$ were observed [43, 39], ranging from $2.7$ to $11\,\mu m$/min. We do not have reliable biological data on the refractory period $T$ and the inflection density $\bar{\rho}$. In [26] it was suggested that the refractory period must be less than 40 sec, while in [7] times around 5 min were suggested. We use $T=2$ min, but note that our analysis shows wave synchronization also for any value larger that $\approx 1$ min. For $\bar{\rho}$ we estimated that it should be of the order of half of the initial mean density to have an effect. A list of all parameters can be found in Tab. 1.

3.1 The particle model in 2D

As a first test of the model, we simulate the full 2D particle model described in Sec. 2.1. Details about the numerical method as well as simulation parameters can be found in App. A.1. Both the initial positions $X_{i}\in\mathbb{R}^{2}$ and the initial orientations $\Theta_{i}\in[-\pi,\pi)$ follow a uniform random distribution. The age variable was initialized with a uniform random distribution on $[0,1]$, i.e. all bacteria are assumed to have reversed before the start of the simulation, however also other choices will lead to the same behavior.

The particle model shows ripple formation.

Within about one hour an almost spatially constant nematic mean direction is established (in this simulation it is $\bar{\Theta}=120^{\circ}$) and all bacteria are either aligned or anti-aligned with it, with small deviations caused by the noise. Just like in experiments, macroscopic traveling density bands develop, in which bacteria travel in the same direction as the band itself and normal to its longitudinal axis. Fig. 4(a) depicts the simulation outcome at time $t=200$ min (see also SI video 1 for the whole time history of this simulation). To distinguish between bacteria aligned or anti-aligned with the (global) nematic mean direction, they are shown in red and blue respectively and in the following we will sometimes refer to them as right- and left-moving bacteria. The global ordering indicates that the nematic alignment quickly drives the system to a quasi 1D situation. To further analyze what happens along the nematic mean direction, we calculate the densities of the right- and left-moving bacteria within a thin strip in the simulation domain (Fig. 4(b) upper and middle rows). Additionally we examine the composition of each wave in terms of C-signaling sensitive and insensitive bacteria: the middle rows in Fig. 4(b) clearly show the density waves and indicates that the wave composition is different before (box A) and after (box B) a wave crest collision. This will be examined further below.

Individual bacteria reverse upon crest collisions.

From biological experiments it is known that when two waves meet, bacteria in the crests typically reverse their direction of movement. The macroscopic, experimental observations, i.e. that counter-propagating waves move with approximately the bacterial speed, persist over time and travel through each other seemingly unaffected, can be easily verified in Fig. 5, which shows a space-time plot of the total (1D) densities in the rectangular strip marked in Fig. 4(a) for 170 min $\leq t\leq 200$ min. To examine the behavior of individual bacteria in this macroscopic context, Fig. 5(a) also shows the space-time path (blue) of two individual bacteria (marked in Fig. 4(a)). Consistent with experiments, bacteria mostly reverse upon crest collision and rarely in between. This shows that the waves are mostly reflected off each other, confirming the accordion-like behavior known biologically.

3.2 Adding spatial diffusion in the 2-age model.

The above results show that the particle model provides a faithful approximation of the biological reality. Since particle models come with a high computational cost and are inherently hard to analyze analytically, we want to use the much simpler 2-age model to gain biological insight into wave formation. Motivated by the observation that alignment leads to a global nematic mean direction, we assume $\bar{\theta}\equiv 0$ in the 2-age model, which amounts to setting $v(\bar{\theta})=(1,0)^{T}$ and omitting the equation for $\bar{\theta}$. The remaining four equations constitute a non-linearly coupled system of transport-reaction equations with no spatial diffusion.

No wavelength control without spatial diffusion.

Simulating this system with initial conditions, randomly perturbed around the space-homogeneous steady state (3.26), we found that it produces counter-propagating density waves, just like its IBM counter-part. However, the wavelength can be arbitrarily small and is only limited by the diffusivity of the numerical scheme or the spatial resolution (see App. 11). The waves in the IBM have a controlled wavelength, indicating that spatial diffusion is inherent in the IBM, but is “lost” during the derivation of the SOH. Rigorous derivation of the spatial diffusion term will be a topic for future work and will involve performing a Chapman-Enskog expansion [10] as detailed e.g. in [20]. The fact that we need a higher order approximation, is most likely a consequence of the fact that for this system the microscopic scale (size of one bacterium, $\approx 5\,\mu m$) is small, but not negligibly small, compared to the size of the developing waves (order of 50-100$\mu m$). In this work we add a diffusion term of size $\delta$ to (2.22) and use the IBM to estimate this otherwise difficult to measure parameter by fitting the produced wavelengths. In all of the following simulations we use the following diffusion-corrected version of (2.22):

[TABLE]

The diffusion corrected 2-age model is a good approximation of the particle model.

We simulate (3.25) with randomly perturbed, constant initial conditions of equal mean density as for Fig. 4 (see caption of Fig. 4 for details). We used $\delta=0.8\,\mu m^{2}/$min, fitted to produce the correct wavelength. In the 2-age model, opposing traveling waves also emerge (discussed in more detail in Sec. 3.3). Fig. 4(b) compares the 1D densities calculated from the IBM to those of the 1D 2-age model: Crest and trough widths match closely. Further one can observe that the composition of the crests in terms of refractory and non-refractory bacteria match very well prior and after crest collisions (boxes A and B in the middle and lower rows). Finally also the macroscopic wave behavior is very similar to that observed in the space-time plots shown in Fig. 5. The remainder of the section is therefore devoted to analyzing the diffusion corrected 1D 2-age model (3.25).

3.3 The wildtype: emergence of waves

We want to understand wave formation, shape and behavior in more detail using the 1D 2-age model (3.25). All simulations are performed with periodic boundary conditions on an interval of length $L=250\,\mu m$ using the wildtype parameters listed in Table 1 (unless stated differently). The numerical method is discussed in App. A.1. Note that system (3.25) conserves the total mass and we define $2\,m_{0}$ as the average total density, which is constant in time:

[TABLE]

As initial conditions we determine spatially uniform steady-state solutions by finding values which make the reaction terms on the right-hand side of (3.25) zero, yielding

[TABLE]

and perturb them with a uniform random distribution. These steady-state solutions reflect the fractions of non-refractory and refractory cells in the absence of spatial patterning. Large reversal rates, i.e. large values of $\lambda(m_{0})$ will increase the fraction of refractory cells, because cells will spend less time on average in a non-refractory state. Small refractory periods $T$ on the other hand decrease them, as they will shortly become sensitive to C-signaling again.

In Fig. 6 a time series of one simulation is shown. After about one hour bands of oppositely moving ripples start to develop and are fully established after two hours, after which their general shape and speed do not change anymore. The ripple crests move with a speed close to the bacterial speed ($9\,\mu m$/min) to the left and right respectively. The density ratio between crests and troughs is about $10$, which corresponds to the experimental values found e.g. in [39]. Upon collisions of two such waves, the total bacteria densities (Fig. 6E-H) double, as described also in experiments. The shapes of the individual waves of left- and right-moving bacteria themselves seem to be almost unaffected by the collision, with only slight deformations. However, when inspecting the composition of the ripple crests in terms of refractory and non-refractory bacteria during a collision, one can observe two distinct phases: a collision phase and a reconstitution phase (see Fig. 7): while before the collision the fraction of refractory bacteria is low, it increases rapidly in the collision phase, indicating a high number of reversals taking place. In the reconstitution phase that follows, this fraction decreases again and resumes its original value. During this phase, the C-signal insensitive bacteria go through their refractory period and “age” back into C-signaling sensitive bacteria. To estimate the fraction of bacteria that reverse during a collision, we compare the total number of reversing bacteria in one wave in the course of one collision to the original number of bacteria present in that wave before the collision. Fig. 7 (left, solid line) shows that more than half of all bacteria reverse.

3.4 The influence of the refractory period $T$

A memory-free model still produces traveling waves.

A crucial part of the presented model is the introduction of the refractory period $T$. What is its influence on the bacterial behavior predicted by the 2-age model? In order to assess this, we investigate what changes in the absence of a refractory period, i.e. for $T\rightarrow 0$. The corresponding model can be interpreted as a memory-free model, as bacteria retain no information about their previous reversals (see explanation below). Mathematically this can be realized by taking the limit $T\rightarrow 0$, i.e. bacteria are susceptible to C-signaling immediately after they reverse. System (3.25) then reduces to $\rho_{\pm}^{0}\equiv 0$ and $\rho_{\pm}^{1}=\sigma_{\pm}$ fulfilling (compare with Fig. 3(c) for the reaction diagram)

[TABLE]

[TABLE]

Remark 3.1

The formally correct way to derive this system requires two steps: Firstly, system (3.25) needs to be non-dimensionalized by introducing a suitable reference timescale $\bar{t}$. A natural choice is $\bar{t}=\frac{L}{d_{1}v_{0}}\approx 27$ min, i.e. the time it would take a non-reversing bacteria to cross the domain. The above limit then amounts to saying that $\varepsilon=T/\bar{t}$ is small compared to this typical time scale. Secondly, one needs to Taylor expand the functions $\rho_{\pm}^{0}$ and $\rho_{\pm}^{1}$ in terms of $\varepsilon$. One then finds that to the first order $\rho_{\pm}^{0}$ are zero and $\rho_{\pm}^{1}$ fulfill system (3.27).

The diffusion free version of this system has already been described in [17], where a memory-free myxobacteria model without an internal age variable $s$ was derived. Note that without the reaction term (i.e. $\lambda\equiv 0$), the system describes the macroscopic limit of purely nematic interactions, a phenomena of great interest in physics and studied in various works [19, 35, 37]. Assuming a constant nematic direction $\bar{\theta}$, the system reduces to two coupled transport-reaction equations; equations of this type were examined in the context of pattern formation and aggregation in biological systems e.g. in [22, 31] (see also the discussion below).

We simulate the memory-free model (3.27) in one space dimension (assuming $\bar{\theta}\equiv 0$), with the same parameters as for the 2-age model (Tab. 1) and observe that the constant steady state again destabilizes under randomly perturbed initial conditions and traveling bumps occur. However, their widths vary greatly and do not seem to be controlled by the dynamics at all. In fact in [31] it was commented that without diffusion the system seems to converge to piecewise constant traveling waves, traveling precisely at speed $d_{1}v_{0}$, which fulfill (3.27) in a weak sense.

The refractory period causes wave synchronization.

The results from the memory-free model suggest that while the refractory period might not be necessary for the formation of opposing traveling waves, it seems to be responsible for controlling the width of the individual traveling crests and synchronizes the waves by controlling the wavelength. To examine this further we systematically vary $T$ in the 2-age model and examine the resulting effects. In Fig. 8 it can be observed that there are two parameter regimes. For very small refractory periods, the system behaves similarly to the memory-free, limiting system (3.27), i.e. while traveling waves occur, the crest width and wavelength varies greatly, indicating a lack of synchronization. After reaching a critical value of $T\approx 1$ min the waves become more synchronized and the crests share a common width. In this regime larger refractory periods lead to wider crests with larger wavelengths. How can we interpret these results?

The refractory period imposes a time-delay between successive collisions.

The reason for the lack of synchronization in the case of small $T$ lies in the fact that bacteria retain no memory of previous reversals, hence two waves cannot affect each other. To understand how the refractory period enhances wave synchronization, it is instructive to look again at Fig. 7: as described above, during the reconstitution phase after a collision, the number of non-refractory cells grows back to its equilibrium value. This growth depends on $T$ (the smaller $T$, the faster the growth). If now the crest meets another crest before sufficiently many bacteria have regained their sensitivity to the C-factor, the wave will not be fully reflected off the oncoming wave. Fig. 9 depicts the result of a simulation, in which one right-moving wave meets two left-moving waves traveling with a short wavelength for both the 2-age model (A-D) and the memory-free model (F-I). Fig. 9E and J depict how many bacteria reverse on average in each wave per minute (i.e. to be precise, if $x_{a}$ and $x_{b}$ mark the beginning and end positions of a wave, Fig. 9E/J depict $\frac{1}{x_{b}-x_{a}}\int_{x_{a}}^{x_{b}}\lambda(\sigma_{\mp})\rho_{\pm}^{1}\,\mathrm{d}x$ in both models). At the first collision the same (high) number of bacteria reverse in both waves, the waves are reflected off each other and the number of reversing bacteria is similar for both the 2-age model and the memory-free model. In the second collision for the 2-age model, however, the right-moving wave contains much fewer non-refractory cells (i.e. $\rho_{+}^{1}$ is lower) than the second left-moving wave, hence it is only partially reflected. The overall effect is a significant reduction of the second left-moving wave. For the memory-free case, the second collision resembles exactly the first collision, hence the second wave is unchanged.

3.5 Parameter-dependence of wave formation

Two necessary wave conditions.

In many experimental set-ups, one of the main output is whether or not myxobacteria colonies form waves. We therefore examine what parameter combinations lead to wave formation. The previous section suggests that the length of the refractory period $T$ has no influence on the appearance of waves. We therefore perform a rough parameter scan over the shape parameters of the reversal function $\lambda(\rho)$, i.e. the spontaneous reversal rate $\lambda_{m}$, the maximal reversal rate $\lambda_{M}$ and the inflection density $\bar{\rho}$. Rigorous mathematical stability analysis of system (3.25) will give more insight into the precise stability regions also for other shapes of $\lambda(\rho)$ and will be the subject of future work. However, several conclusions can already be drawn. We find two wave formation conditions:

Condition A: The maximal reversal rate needs to be large enough compared to the spontaneous reversal rate ($\lambda_{M}>3\lambda_{m}$).

Condition B: The inflection density and the average total density need to be of similar order ($\bar{\rho}\approx m_{0}$).

How can we links this to experiments?

Remark 3.2

A non-dimensionalization of the system shows that its behavior depends on four dimensionless quantities

[TABLE]

The first quantity compares the time it takes a just reversed bacterium to age back into a reversible state, to the time it takes a (non-reversing) bacterium to cross the domain. The next two quantities can be interpreted as the expected number of reversals during one domain crossing for reversal rates $\lambda_{m}$ and $\lambda_{M}$ in the absence of a refractory period. The last quantity compares the initial average density to the inflection density, i.e. the density of opposing bacteria, where the reversal function $\lambda$ is the most sensitive to changes in density. Examining $\lambda_{m}$, $\lambda_{M}$ and $\bar{\rho}$ therefore amounts to analyzing the last three quantities.

Mutation experiments: hypo- and hyper-reversing bacteria.

To assess how the reversal behavior impacts the ripple formation, Sager and Kaiser have used M. xanthus strains, that have an insertion mutation in the frzCD gene, which has been shown to impact the reversal probability [39]. Isolated individuals of the hypo-reverses change direction on average $\approx$10 times less frequent and hyper-reversers $\approx$5 times more often as compared to the wildtype. Both strains have lost the ability to form ripples and it is suggested that the mutation affects the spontaneous reversal rate. Wave condition (A) suggests that if the spontaneous reversal rate $\lambda_{m}$ is too large compared to the maximal reversal rate, no waves will form, in agreement with the hyper-reverser experiments. Our wave formation conditions do not explain the absence of waves in hypo-reversers, however measuring the reversal rate $\lambda(\rho)$ experimentally could shed light on this question.

Dilution experiments: changing the fraction of C-signal competent cells.

Another way to demonstrate the influence of C-signaling on the ripple formation was a second set of experiments performed by Kaiser and Sager [39]: two strains of myxobacteria were used, a C-signaling competent wildtype strain $csgA^{+}$ and a mutant strain $csgA^{-}$, which can not produce C-signal, but can respond to it. By themselves ensembles of $csgA^{-}$ cells do not have the ability to form ripples. In their experiments Kaiser and Sager changed the fraction of $csgA^{+}$ and $csgA^{-}$ cells and measured the changes in ripple wavelength, speed and width. To simulate this situation using our 2-age model, we note that, since $csgA^{-}$ cells react to C-signals in the same manner as the wildtype, the fraction of $csgA^{+}$ and $csgA^{-}$ will be constant everywhere. This was confirmed in [39]. Let the fraction of C-signaling competent $csgA^{+}$ cells be $q\in[0,1]$. At some density of opposing bacteria $\sigma$, therefore only $q\sigma$ cells contribute to the amount of signal produced, i.e. $\lambda(\sigma_{\pm})$ in system (3.25) is replaced by $\lambda(q\sigma_{\pm})$. Substituting this into (2.24), the expression for $\lambda$, we see that this in fact simply changes $\bar{\rho}$ to $\bar{\rho}/q$. In view of wave condition (B), this suggests that if $q$ is too small, this would inhibit wave formation. This is in agreement with the experiments described in [39], although it should be noted that for very small fractions of $csgA^{+}$ cells, for which waves were still present in the experiment, the corresponding large values of $\bar{\rho}$ do not produce waves anymore in our model. We further study how varying the inflection density $\bar{\rho}$ affects the wave shape (within the range that supports ripple formation). Fig. 10 depicts the impact of this variation on various wave characteristics. Wave speed is not affected and stays close to the individual bacterial speed (Fig. 10A). For regular waves, crest width increases, as does ripple wavelength (Fig. 10C/D), both of which are in agreement with the findings in [39]. Our model also predicts a non-monotonous dependence of the ratio of the number of cells in the crest to the number of cells in the trough as a function of $\bar{\rho}$ (Fig. 10B), which could easily be addressed experimentally.

4 Discussion

We presented a novel, continuous, age-structured macroscopic model of myxobacteria, derived systematically from an individual-based model. The main assumptions of the model are nematic alignment, a density-dependent reversal function and a refractory period of fixed length, which introduces a memory effect.

In excellent agreement with experimental data on myxobacteria, simulations of the full IBM show the development of periodic waves, traveling in opposing directions and being reflected upon collision. We performed an in-depth numerical investigation of the one-dimensional macro-model for the case of two sensitivity/age groups: refractory cells, incapable of reversing, and non-refractory cells, that can reverse and are sensitive to C-signaling. A main result of our analysis is that the refractory period is not responsible for wave formation, but for wave synchronization. This is because it introduces a memory effect, that controls wavelengths. The existence of a refractory period is known for example for D. discoideum [41]. The idea of a refractory period for myxobacteria has originally been brought forward in the model presented in [26], where it is suggested that myxobacteria also go through an insensitivity period following a reversal. Assuming that the length of the refractory period is unaffected by the density of opposing bacteria, we show that the memory effect introduced by a fixed refractory period is sufficient to explain ripple synchronization. We discovered two wave formation conditions that are consistent with experimental results: the maximal reversal rate needs to be large enough compared to the spontaneous reversal rate and the average density of the myxobacteria colony needs to be close to $2\bar{\rho}$ with $\bar{\rho}$ being the inflection density at which the reversal function reacts the most sensitively to density changes. This predicts that both very high and very low densities will inhibit wave formation in myxobacterial colonies.

A strength of the Self-Organized Hydrodynamics (SOH) models is that they are directly derived from the corresponding IBM by the method of Generalized Collision Invariants (GCI). This allows for easy and transparent inclusion of assumptions, such as the density-dependent reversal rate and the refractory period. The simulations presented in Sec. 3.2 show good agreement between the IBM and diffusion corrected 2-age model, however, to fully validate the correspondence between the macroscopic model and the IBM, a formal derivation of the mass diffusion term [10, 20] will be necessary. The lack of control over the number of waves in absence of diffusion shows that for phenomena where the ratio between microscopic and macroscopic scales is not very small, higher order terms can be crucial and hence the additional effort for deriving them is justified.

In this work we concentrated on wave formation, where both experimental results and particle simulations suggest that the dynamics can be studied in one space dimension. However, several other macroscopic patterns are known for myxobacteria, most notably the formation of large aggregates. The combination of simulation results of both the IBM and the two-dimensional 2-age model as well as analytical results on the latter might shed light on what parameters cause cells to switch between a uniform state to ripple formation and aggregation. Several aspects will need to be considered: could the reversal frequency depend on both the densities of the opposing group as well as on the aligned one? In [31] such cases were analyzed for a memory-free model and simulations also showed wave formation. Further the authors in [31] also noted that small changes in the reversal function $\lambda(\rho)$ can cause the system to switch from ripples to aggregation. It will be interesting to extend the results to the age-structured model. In [28] experiments were described in which C-signaling and thereby local densities affect the gliding speed of bacteria. Such density-dependent parameters can easily be incorporated into SOH models and significantly impact the dynamics [23].

A large number of IBMs exist for self-propelled particles such as myxobacteria [43] and detailed numerical and statistical analysis of their properties have significantly contributed to the understanding of emergent phenomena and phase transitions. While they allow for direct comparison between the experimental data about the behavior of the individuals, they are limited in terms of insights into macroscopic phenomena. For macro-systems such as the 2-age model on the other hand, stability or asymptotic analysis can be performed which can elucidate precise parameter dependencies and long-term behavior. The macro-system of [26] triggered a number of works examining wavelength determination [6, 27], demonstrating the potential insight gained through analytical methods. For the macro-system presented here, a rigorous analytical treatment will be the subject of future work. Several works deal with the linkage between IBMs and meso- or macro-models, e.g. both in [7] and [24] where continuous models are derived from particle models. The methods presented there as well as the GCI method offer a powerful option of combining the strengths of both particle and continuous based methods.

Our work together with the results in [7, 26, 43] strongly suggests the existence of a refractory period. Several aspects of our findings can be addressed experimentally: our main results is that wave formation and wave synchronization are independent phenomena, which would suggest separate molecular mechanisms. We therefore predict that mutants that form non-synchronized traveling waves of various width have a density-dependent reversal frequency, but no refractory period. As to the length of the refractory period, in [26] experimental data argued for refractory periods of under 1 min. Our model argues that $T\approx 1$ min presents the lower limit for the formation of synchronized waves. Larger refractory periods also lead to periodic waves, but it takes the system much longer to produce steadily moving wave. Hence the value $T\approx 1$ min is the fastest way to make synchronized waves, presenting a possible, evolutionary answer as to why this particular refractory period length evolved.

Appendix A Appendix

A.1 Numerical Methods

Particle Model.

The IBM of Sec. 2.1 was simulated in a square domain with periodic boundary conditions, using the method described in [34].

2-Age Model.

We simulated the 2-age model (3.25) in one space dimension with periodic boundary conditions and $\bar{\theta}\equiv 0$. The transport and reaction terms were implemented using operator splitting with explicit upwind or downwind (for the $+$ and $-$ family respectively) finite differences for the transport term, an implicit finite difference scheme for the diffusion term and an explicit treatment of the reaction term. We also tested implementing the reaction term with an explicit Runge-Kutta (4,5) formula using the ode45 solver of Matlab and a Lax-Friedrichs Scheme for the transport term, both with no significant gains in accuracy. Tab. 2 shows the parameters used for both simulation set-ups.

A.2 Controlling the number of waves.

We want to asses whether (2.22) encodes information about the number of waves naturally created. We use as initial conditions small perturbations from (3.26), however using completely random data can be problematic due to its non-smootheness and because the result could be influenced by the spatial discretization. To avoid these problems, but mimic randomness, we use sums of Fourier-modes as perturbations:

[TABLE]

where $a_{k},b_{k}$ are chosen randomly in $[-1,1]$. Finally we normalize the resulting perturbations so that the maximal amplitude doesn’t exceed $\varepsilon=0.05$. For the time discretization we used $\Delta t=\Delta x/(d_{1}v_{0})\times$CFL, where CFL is the Courant-Friedrichs-Lewy number. The up-/down winding scheme for the transport terms implies that we need CFL$\leq 1$ for stability. Without the reaction terms, the transport equations could be solved numerically exactly, if CFL=1. Then the numerical diffusion introduced by the discretization of the transport term increases for smaller CFL.

Arbitrary small wavelengths in the absence of diffusion.

To examine the effects of numerical diffusion and spatial resolution on the wavelength, we vary 1.) the maximal number of Fourier modes $K$ used in the initial conditions, 2.) the spatial discretization $\Delta x$ and 3.) the CFL condition. For each choice of $K$ we use the same initial conditions, i.e. the same random numbers $a_{k},b_{k}$. Each simulation is run until time = $300$ min and the number of waves is determined as the average number of peaks of $\sigma_{+}$ and $\sigma_{-}$. Fig. 11A shows the outcome. For large $K$, for each fixed CFL number, a finer grid leads to more waves. On the other hand for a fixed grid size, the closer CFL is to $1$, the more waves we observe. This indicates that the upper bound of the number of waves is caused by the numerical diffusion and the spatial resolution, rather than by the equations themselves. However, the minimal number of waves seems to be encoded in the equations, as it is largely independent of the spatial discretization and CFL number. The issue of wavelength selection for a similar system was also discussed in [40].

To test the influence of diffusion wavelengths, we repeat the same simulations, using the diffusion-corrected model (3.25). Fig. 11B shows that now the behavior is consistent for different $\Delta x$ and the maximal number of waves is controlled. Since in the IBM we observe regular waves of relatively fixed wavelengths, we conclude that a small amount of density diffusion is inherent in the IBM and leads to a control in wavelength. This diffusion is lost during the derivation of the SOH model, but is necessary to avoid very small wavelengths (see discussion in Sec. 3.2).

Acknowledgements

PD acknowledges support by the Engineering and Physical Sciences Research Council (EPSRC) under grant no. EP/M006883/1, by the Royal Society and the Wolfson Foundation through a Royal Society Wolfson Research Merit Award no. WM130048 and by the National Science Foundation (NSF) under grant no. RNMS11-07444 (KI-Net). PD is on leave from CNRS, Institut de Mathématiques de Toulouse, France. AM was supported by the Austrian Science Fund (FWF) through the doctoral school Dissipation and Dispersion in Nonlinear PDEs (project W1245) as well as the Vienna Science and Technology Fund (WWTF) (project LS13/029). HY acknowledges the support by Division of Mathematical Sciences [KI-Net NSF RNMS grant number 1107444]; DFG Cluster of Excellence Production technologies for high-wage countries [grant number DFG STE2063/1-1], [grant number HE5386/13,14,15-1]. HY and AM gratefully acknowledge the hospitality of the Department of Mathematics, Imperial College London, where part of this research was conducted.

Data availability

No new data were collected in the course of this research.

Supplmentary information

Link to video on figshare