# Modeling the Interplay of Oscillatory Synchronization and Aggregation via
Cell-Cell Adhesion

Tilmann Glimm, Daniel Gruszka  Department of Mathematics, Western Washington
University, Bellingham WA 98225

###### Abstract.

We present a model of systems of cells with intracellular oscillators
(’clocks’). This is motivated by examples from developmental biology and from
the behavior of organisms on the threshold to multicellularity. Cells undergo
random motion and adhere to each other. The adhesion strength between
neighbors depends on their clock phases in addition to a constant baseline
strength. The oscillators are linked via Kuramoto-type local interactions. The
model is an advection-diffusion partial differential equation with nonlocal
advection terms. We demonstrate that synchronized states correspond to Dirac-
delta measure solutions of a weak version of the equation. To analyze the
complex interplay of aggregation and synchronization, we then perform a linear
stability analysis of the incoherent, spatially uniform state. This lets us
classify possibly emerging patterns depending on model parameters. Combining
these results with numerical simulations, we determine a range of possible
far-from equilibrium patterns when baseline adhesion strength is zero: There
is aggregation into separate synchronized clusters with or without global
synchrony; global synchronization without aggregation; or unexpectedly a
“phase wave” pattern characterized by spatial gradients of clock phases. A 2D
Lattice-Gas Cellular Automaton model confirms and illustrates these results.

##  1\. Introduction

###  1.1. Overview

The spontaneous synchronization of coupled oscillators is a common occurrence
in physics, biology and social sciences that has captured the interest of
mathemtical modelers for decades [55, 40, 6]. Well-studied, prominent
phenomena are the emergence of synchrony in neural networks, i.e,. the
synchronous firing of nerve cells [45, 24, 52]. and the synchronization of the
flashing of fireflies (f.ex. [18, 41]; see also the excellent article [46]).
Synchronization of intracellular molecular oscillatory processes is also
ubiquitous in developmental biology [30, 5, 54, 12].

The fundamental mathematical model of synchronization is the paradigmatic
Kuramoto model [33, 49, 48, 56] which has been studied extensively, in
particular on general complex networks [42]. In general, synchronization of
spatially distributed oscillators on networks is the focus of very active
research (see e.g. [6, 15] and the references therein). Typically these
networks, while often having complex topology, are static and do not encompass
spatiotemporal rearrangements of the nodes [15, 42].

This is a very appropriate assumption for many applications, e.g. neural
networks, but it does not apply to aggregation processes in which the
individual oscillators are moving in space, and crucially, there is an
interplay between clocks and motion. The study of the feedback between motion
and clock synchronization in systems of moving oscillators has only relatively
recently come into focus. Among those is the work of Tanaka and coworkers on
systems of chemotactic oscillators [50, 29]. The most influential work in this
area is O’Keefe et al.’s seminal investigation of swarming oscillators [37],
which has stimulated a number of further papers [43, 36] and has fascinating
applications to e.g. the behavior of fleets of robots [3]. O’Keeffe et al.
studied systems of Kuramoto-like oscillators (“swarmalators”) represented by
mass points which attract or repel each other.

One important application of systems of interacting moving oscillators is
aggregation of biological cells mediated by cell-cell adhesion. In such
situations the swarmalator model of O’Keefe et al. has some drawbacks, most
prominently that spatial attraction between two oscillators is independent of
the distance between them. Their model leads to highly symmetrical patterns of
aggregation (a single disk or ring). In contrast, adhesive interactions
between cells are local, which, combined with random motion, typically leads
to irregularly shaped aggregation clusters.

We therefore propose and investigate here a model where oscillators are
coupled to adhesive forces, using models of cell-cell adhesion. Before
describing our model in more detail, we expand on our motivation:
investigating the aggregation of biological cells with internal biological
clocks.

###  1.2. Aggregation and Molecular Clocks: Dictyostelia and Myxobacteria

Two model organisms have been studied extensively: the slime mold
Dictyostelia, in particular the best-studied species Dictyostelium discoideum
[7, 53, 34] and the slime bacteria Myxobacteria [44, 58, 39, 51]. While these
two lineages are phylogenetically very distant (Dictyostelia are eukaryotes
and Myxobacteria prokaryotes), they display convergent evolution as being on
the threshold of multicellularity: Both can display solitary cell behavior
(though only observed in a laboratory setting for Myxobacteria [1]), where
individual cells behave independently of each other. But in conditions of
nutrient scarcity, both transition into a developmental stage characterized by
cell aggregation which culminates in the formation of multicellular structures
called fruiting bodies.

The details of the mechanisms of aggregation are quite distinct for the two
lineages and research of the two lineages has been almost completely
independent. In Dictyostelia, aggregation has been tied to a complex interplay
of production of the chemotactic agent cAMP by cells and the movement of the
same cells up the cAMP gradient. This is the basis of the celebrated Keller-
Segel model [32, 27]. In Myxobacteria, aggregation has been explained by
coordinated cell motion mediated by alignment of the rod-shaped cells and
synchronization of the reversal rate of neighboring cells [57, 28]. In a
recent article, Del Angel et al. [1] pointed out that despite these
differences, there are some very profound similarities in the mechanisms of
aggregation and directed migration between the two organisms. They argue that
in both instances, it is a complex interplay between generic processes, i.e.
mechanical or physicochemical processes that operate similarly in living and
nonliving material and agent-like processes unique to living systems. The
latter are mediated through networks of interaction with neighboring cells
through mechanical and chemical signals. A paradigm for this is cell-cell
adhesion interacting with the spontaneous emergence of synchronicity in
intracellular oscillators. In Dictyostelia, the molecular basis is the
periodic production of cAMP. (It is not clear if these oscillations are an
intrinsic property of individual cells or if they exist only as a population-
level phenomenon [22].) In Myxobacteria, the periodic reversal of directions
is controlled by the Frz and MglAB intracellular oscillators [23].

###  1.3. Model and Results

In our bottom-up model of oscillators coupled to adhesive forces (’sticky
clocks’), we use a PDE model for a cell density $R(t,x,\theta)$ that depends
on time $t$, position $x$ and clock value $\theta$. To model cell-cell
adhesion, we adapt the model of Armstrong et al. [2]; see also [16, 9, 10, 20,
21], which was set up to model cell sorting by differential adhesion [47, 19]
in a PDE framework. In our approach, adhesion is modeled by an advection-
diffusion equation with nonlocal advection and Kuramoto-like interactions
between oscillators via advection terms in the clock variable $\theta$.

There is a baseline adhesion strength given by a parameter $J_{0}$. An
additional differential adhesion term depends on clock phases and is encoded
by $J$. The coupling strength of oscillators is denoted by $K$. For $K>0$,
neighbors want to synchronize, for $K<0$, they want to anti-synchronize.
Positive $J$ means that cells’ adhesion to neighbors is strongest if they are
in the same clock phase; for $J<0$, the adhesion is strongest for cells in
opposite phases. (This choice of notation is motivated by O’Keeffe et al’s
swarmalator model [37], where $J$ encodes attraction between particles.)

In this model, spatial aggregation and synchronization can co-occur and,
depending on the parameters, can help or hinder each other. We perform a
linear stability of the spatial homogeneous, incoherent (unsynchronized)
state. This leads to a classification of the possible types of emerging
patterns (“pattern onsets”). The pattern onsets display either globally
synchronized, locally synchronized or incoherent behavior with or without
spatial aggregation. We then compare the linear pattern onsets with far-from
equilibrium patterns from numerical solutions in the case of zero baseline
adhesion ($J_{0}=0$). We find that the pattern onsets obtained by linear
stability analysis generally capture the far-from equilibrium behavior:
Positive clock coupling $K$ leads to synchronization, where the value of the
spatial adhesion strength $J$ controls whether there is spatial aggregation
and if the synchronization is global or only local within aggregated clusters.
For negative $K$ and negative $J$, the incoherent state is stable and no
aggregation occurs. We found a surprising result when $K$ is negative but
small in absolute value and $J$ is sufficiently large. In this case, small
diffuse aggregates form. Oscillators are not synchronized, but instead ther
are spatial gradients of clock phases, so that there is a correlation between
clock phase and spatial position (a “phase wave” pattern). These results are
summarized in Figure 5 with examples in one spatial dimension in Figure 6, and
with a two-dimensional Lattice-Gas Cellular Automaton implementation in Figure
9.

The structure of the paper is as follows: In section 2, we present the model,
which consists of a diffusion-advection equations with nonlocal terms, and
discuss its features. We formulate a weak version of the PDE in section 3,
which enables us to mathematical rigorously identify synchronized states as
Dirac-delta functions. In section 4, we perform a linear stability analysis,
which leads to a classification of pattern onsets. Finally, in section 5, we
investigate the far-from-equilibrium behavior of soultions by combining the
results of the linear stability analysis with numerical solutions of the PDE
(in one spatial dimension) and an implementation via a Lattice-Gas Cellular
Automaton as described above.

##  2\. The Model

We consider the density $R(t,x,\theta)$ of moving oscillators. Here $t$ is
time, $x$ is spatial location and $\theta$ is the phase of the internal clock.
The number $N$ of ocillators located in a region $U$ of $n-$dimensional space
with phases between $\theta_{1}$ and $\theta_{2}$ at time $t$ is thus

$N=\int_{U}\int_{\theta_{1}}^{\theta_{2}}R(t,x,\theta)\,d\theta\,d^{n}x.$

This setup is applicable e.g. to biological cells that have an internal
molecular clock. The model is described by the following conservation law for
the density $R(t,x,\theta)$:

$\displaystyle\frac{\partial R}{\partial
t}=D_{x}\nabla^{2}R-\nabla\cdot(v_{1}(t,x,\theta;R)R)+D_{\theta}\frac{\partial^{2}R}{\partial\theta^{2}}-\,\frac{\partial}{\partial\theta}(v_{2}(t,x,\theta;R)R).$
(2.1)

The domain for (2.1) is

$(x,\theta)\in[0,L]^{2}\times[0,2\pi]$

and we assume periodic boundary conditions, i.e. the density $R(t,x,\theta)$
is required to be $2\pi$-periodic in $\theta$ and $L-$periodic in $x$.
Equivalently, we can think of $R(t,x,\theta)$ as function on the domain
$[0,T)\times\Omega\times\mathbb{T}$. Here $T>0$ and
$\Omega={\mathbb{R}}^{n}/\sim$ is the quotient of Cartesian space
$\mathbb{R}^{n}$ by the equivalence relation that identifies the sides of a
hypercube (square or cube in two or three dimensions) with each other, i.e.
$(x,y)\sim(x+L,y)\sim(x,y+L)$. Similarly,
$\mathbb{T}=\mathbb{R}/2\pi\mathbb{Z}$ is a circle, namely the quotient
$\mathbb{R}$ by the equivalence relation $\theta\sim\theta+2\pi$.

The term $v_{1}(t,x,\theta,R)$ models the effective velocity of cells due to
cell-cell adhesion. The term $v_{2}(t,x,\theta,R)$ is the speed at which the
clock phases advances.

To set up the functional forms of these functions, we assume that neighboring
cells influence each others’ cell clocks in a fashion similar to the Kuramoto
model. Specifically, we introduce a parameter $K$ below that describes and
quantifies the strength of this influence. For $K>0$, neighboring cells seek
to synchronize – simply put, if my neighbor’s clock is ahead of mine, my own
clock speeds up a little. If $K<0$, this relation is reversed, and neighboring
cells tend to anti-synchronize, i.e. if my neighbor’s clock is ahead of mine,
my own clock tends to slow down.

Cells also move in the space. They undergo diffusion, and adhere to each
other. The strength of the adhesion between neighbors depends on the absolute
difference $|\Delta\theta|$ of their clock phases. This is quantified by a
second parameter $J$. For $J>0$, the adhesion strength is at a maximum at
$|\Delta\theta|=0\,(\rm{mod}2\pi)$ – in other words “like attracts like” and
cells with the same phases adhere particularly strongly. In the case $J<0$,
the adhesion strength is at a minimum when $|\Delta\theta|=0$ and maximal when
$|\Delta\theta|=\pi$, a situation best described by “opposites attract.” (See
Table 1 for a summary of the meaning of the parameters $J$ and $K$. )

Specifically, these assumptions are modeled via the following functional terms
for $v_{1}$ and $v_{2}$:

$\displaystyle v_{1}(t,x,\theta;R)$
$\displaystyle=\int_{0}^{2\pi}\int_{D_{\rho}(0)}(J_{0}+J\cos(\theta-{\tilde{\theta}}))\,\sigma(R(t,x+\tilde{x},{\tilde{\theta}})\,)\frac{\tilde{x}}{|\tilde{x}|}d^{n}\tilde{x}\,d{\tilde{\theta}}$
(2.2) $\displaystyle v_{2}(t,x,\theta;R)$
$\displaystyle=\omega+K\,\int_{0}^{2\pi}\int_{D_{\rho}(0)}\sin({\tilde{\theta}}-\theta)\,R(t,x+\tilde{x},{\tilde{\theta}})\,d^{n}\tilde{x}\,d{\tilde{\theta}}$
(2.3)

Here $D_{\rho}(0)$ is a disk of radius $\rho>0$ centered at $0$ – a template
for the interaction neighborhood: Cells within a spatial distance of $\rho$
influence each others’ phases and adhere to each other. For the function
$\sigma$ in the fluxes, one can typically choose a linear form or a saturating
logistic expression [2, 20]:

$\displaystyle\sigma(R)=R\quad\text{or}\quad\sigma(R)=\frac{R_{\max}}{R_{\max}-1}R\max\left(1-\frac{R}{R_{\max}},0\right)\text{
for some }R_{\max}>0.$ (2.4)

Here on the right, $R_{\max}$ is a maximum cell density. We will concentrate
on the linear form $\sigma(R)=R$ in the analysis, but will use the logistic
form for numerical simulations. The factor is chosen so that the two terms are
equal for $R=1$. (The simulations in section 5 use a constant value $R=1$ with
small random noise for the initial conditions.)

The form of the adhesion velocity $v_{1}$ follows Armstrong et al.’s model
([2]). The parameter $J_{0}$ is a “baseline” adhesion strength. For $J=0$, the
above equations are in fact identical to Armstrong et al’s model. Our addition
is that the strength of attraction is modulated by the clock phase difference,
described by the term $J\cos(\theta-{\tilde{\theta}})$, following a similar
approach as in [37]. The clock speed $v_{2}(t,x,\theta,R)$ comes from the
Kuramoto model [33].

Note that by considering the temporal shift $\theta\to\theta-\omega t$, we can
assume without loss of generality that $\omega=0$.

There are also diffusive terms in $x$ and $\theta$. The random motion of cells
in space leads to diffusion with spatial diffusion coefficient $D_{x}$. There
is also diffusion in $\theta$–space, describing a small random component in
the clock speed with diffusion coefficient $D_{\theta}$. We will put
particular emphasis on the case $D_{\theta}=0$, which is the most natural
assumption (no diffusion in clock-space). However, a small value of
$D_{\theta}$ is advantageous for numerical analysis because it prevents blowup
of solutions. In any case, $D_{\theta}$ should be “small”, in the sense that
the dimensionless quantity $D_{\theta}\rho^{2}/D_{x}$ satisfies

$\frac{D_{\theta}\rho^{2}}{D_{x}}<<1$

(See section 4.3.) The equation (2.1) is complemented by an initial condition
is

$R(t=0,x,\theta)=R_{0}(x,\theta)$

where $R_{0}(x,\theta)$ is a nonnegative function.

We note that while the model (2.1) is defined for any values $J_{0},J,K$, not
all of them are meaningful from a modeling point of view. The model hinges on
the assumption that patterns are sufficiently coarse grained that a
description using cell density is appropriate. Speficially, meaningful values
of $J_{0},J$ and $K$ are those that lead to patterns with a characteristic
length scale not smaller than the cell interaction radius $\rho$, see Appendix
B for further discussion.

Condition | Description | Shorthand
---|---|---
$J>0$ | max adhesion in-phase ($\Delta\theta=0$) | “like attracts like”
$J<0$ | max adhesion out-of-phase ($\Delta\theta=\pm\pi$) | “opposites attract”
$K>0$ | clock speeds up if behind neighbor’s | “neighbors seek to synchronize”
$K<0$ | clock slows down if behind neighbor’s | “neighbors seek to anti-synchronize”
Table 1. Explanation of the parameters $J$ and $K$. Here $\Delta\theta$ is the
difference of a cell’s clock phase and its neighbor’s phase.

##  3\. Nonnegativity, weak versions and steady state solutions

It is useful to formulate weak versions of the equation (2.1). This enables us
to state in a rigorous way that a certain family of Dirac-delta distributions
are steady states. First, let us emphasize that $R(t,x,\theta)$ is a density
function with respect to spatial position $x$ and clock phase $\theta$. Not
surprisingly, solutions remain nonnegative for nonnegative initial conditions
$R_{0}(x,\theta)$. For positive clock diffusion $D_{\theta}>0$, this follows
from a maximum principle for strictly hyperbolic PDEs ([17]), but for
completeness, we include a proof that covers the case $D_{\theta}\geq 0$:

######  Theorem 3.1.

(Nonnegativity of solutions) Let $R(t,x,\theta)\in
C^{2}([0,T)\times\Omega\times\mathbb{T})$ be a strong solution of (2.1), where
$T>0$ and $D_{\theta}\geq 0$. Suppose the initial condition
$R(t=0,x,\theta)=R_{0}(x,\theta)$ is a nonnegative continuous function. Then
$R(t,x,\theta)\geq 0$ for all $0\leq t<T$ and
$x\in\Omega,\theta\in\mathbb{T}$.

###### Proof.

Write $w_{i}(t,x,\theta)=v_{i}(t,x,\theta;R(t,x,\theta))$ for $i=1,2$ For
$0<T_{1}<T$, let

$\lambda>\max_{0\leq t\leq T_{1},x\in\Omega,t\in\mathbb{T}}(|\nabla\cdot
w_{1}|+|\nabla\cdot w_{2}|)$

Consider for $\varepsilon>0$ the function

$R_{\varepsilon}(t,x,\theta)=R(t,x,\theta)+\varepsilon\cdot e^{\lambda t}.$

Then

$\displaystyle\frac{\partial R_{\varepsilon}}{\partial t}$
$\displaystyle=D_{x}\nabla^{2}R_{\varepsilon}+D_{\theta}\frac{\partial^{2}R_{\varepsilon}}{\partial\theta^{2}}-\nabla\cdot(w_{1}\,R_{\varepsilon})-\frac{\partial}{\partial\theta}(w_{2}\,R_{\varepsilon})+\varepsilon
e^{\lambda t}(\nabla\cdot w_{1}+\nabla\cdot w_{2}+\lambda)$
$\displaystyle>D_{x}\nabla^{2}R_{\varepsilon}+D_{\theta}\frac{\partial^{2}R_{\varepsilon}}{\partial\theta^{2}}-\nabla\cdot(w_{1}\,R_{\varepsilon})-\frac{\partial}{\partial\theta}(w_{2}\,R_{\varepsilon})$

The function $R_{\varepsilon}(t,x,\theta)>0$ cannot take on nonpositive values
for $\varepsilon>0$. Indeed, suppose that $R_{\varepsilon}$ attains the value
$0$ for the first time at time $t=t^{*}$ at a point $x=x^{*}$ and clock value
$\theta=\theta^{*}$. Then $\frac{\partial R_{\varepsilon}}{\partial
t}(t^{*},x^{*},\theta^{*})\leq 0$. Also, we have that at those values $\nabla
R=0$ and $\frac{\partial R}{\partial\theta}=0$, and, because it is a minumumn,
, $\nabla^{2}R_{\varepsilon}\geq 0,$ and
$\frac{\partial^{2}R}{\partial\theta^{2}}\geq 0$. But this gives

$\frac{\partial R_{\varepsilon}}{\partial t}(t^{*},x^{*},\theta*)\geq 0,$

by the above estimate, a contradiction. The positivity of $R_{\varepsilon}$
for all $\varepsilon>0$ immediately implies that $R(t,x,\theta)\geq 0$ for all
$t,x,\theta$. ∎

In keeping with the standard approach ([38]), we can formulate a weak version
of the problem: Say a function $R\in
L^{\infty}([0,T)\times\Omega\times\mathbb{T})$ is a solution of (2.1) in the
distribution sense if for all test functions $\varphi\in
C^{\infty}_{c}([0,T)\times\Omega\times\mathbb{T})$, we have

$\displaystyle-\int_{0}^{T}\int_{\Omega}\int_{0}^{2\pi}R(t,x,\theta)$
$\displaystyle\left(\frac{\partial\phi}{\partial
t}+v_{1}(t,x,\theta;R)\cdot\nabla\varphi+D_{x}\nabla^{2}\phi+v_{2}(t,x,\theta;R)\frac{\partial\phi}{\partial\theta}+D_{\theta}\frac{\partial^{2}\phi}{\partial\theta^{2}}\right)\,d\theta\,d^{n}x\,dt$
$\displaystyle=\int_{\Omega}\int_{0}^{2\pi}\varphi(t=0,x,\theta)\,R_{0}(x,\theta)\,d\theta\,d^{n}x,$

where $v_{1}(t,x,\theta;R)$ and $v_{2}(t,x,\theta;R)$ are given by (2.2) and
(2.3), respectively.

For an even more general concept of a solution that permits for Dirac-delta-
like distributions in the clock variable $\theta$, we adapt the approach of
Hillen et al. [26]. Let ${\mathcal{B}}(\mathbb{T})$ be the space of regular
signed real-valued Borel measures on $\mathbb{T}=\mathbb{R}/2\pi\mathbb{Z}$.
This is a Banach space if the absolute variation is taken as a norm, see [26]
for details. Let ${\mathbb{X}}=L^{\infty}(\Omega,\mathcal{B}(\mathbb{T}))$. So
if $R\in\mathbb{X}$ and $x\in\Omega$, then $R(x)$ is a measure on
$\mathbb{T}$. Say $R\in\mathbb{X}$ is a weak steady state if the following
equality holds for all test functions $\varphi\in C^{\infty}(X,\mathbb{T})$:

$\displaystyle
0=\int_{\Omega}\int_{\mathbb{T}}\left(v_{1}(x,\theta,x;R)\cdot\nabla\varphi+v_{2}(x,\theta;R)\frac{\partial{\varphi}}{\partial\theta}+D_{x}\nabla^{2}\varphi+D_{\theta}\frac{\partial^{2}\varphi}{\partial\theta^{2}}\right)\,R(x,d\theta)\,d^{n}x$
(3.1)

Here $R(x,d\theta)$ denotes integration with respect to the measure $R(x)$ and

$\displaystyle v_{1}(x,\theta;R)$
$\displaystyle=\int_{D_{\rho}(0)}\int_{\mathbb{T}}(J_{0}+J\,\cos(\theta-{\tilde{\theta}})\,R(x+\tilde{x},d{\tilde{\theta}})\,\frac{\tilde{x}}{|\tilde{x}|}d^{n}\tilde{x}$
(3.2) $\displaystyle v_{2}(x,\theta;R)$
$\displaystyle=K\,\int_{D_{\rho}(0)}\int_{\mathbb{T}}\sin({\tilde{\theta}}-\theta)\,R(x+\tilde{x},d{\tilde{\theta}})\,d^{n}\tilde{x}$
(3.3)

(These $v_{1},v_{2}$ are Borel-measurable functions on
$\Omega\times\mathbb{T}$.) With these definitions, we can now state that we
have two classes of steady states.

######  Theorem 3.2.

  1. (1)

For any constant $\bar{R}>0$, the constant function $R\in
C^{\infty}~{}(\Omega~{}\times~{}\mathbb{T})$

$R(x,\theta)=\bar{R}=const$

is a strong steady state solution. We refer to it as the incoherent state.

  2. (2)

Suppose $D_{\theta}=0$. Let $\theta_{0}\in[0,2\pi)$ be any clock value.
Suppose $R_{0}\in C^{2}(\Omega)$ is a strong solution of

$\displaystyle D_{x}\nabla^{2}R_{0}$
$\displaystyle-\nabla\cdot(w(x;R_{0})R_{0})=0$ $\displaystyle{\text{ with
}}w(x;R_{0})=\left(J_{0}+J\right)\int_{D_{\rho}(0)}R_{0}(x+\tilde{x})\frac{\tilde{x}}{|\tilde{x}|}d^{n}\tilde{x}.$

Then the function $R\in\mathbb{X}$

$R(x)=R_{0}(x)\,\delta(\theta-\theta_{0})$

is a weak steady state. We refer to it as a (perfectly) synchronized solution.

Here $\mu=\delta(\theta-\theta_{0})$ denotes the point mass (Dirac) measure on
$\mathbb{T}$ concentrated at $\theta=\theta_{0}$, i.e. for any Borel set
$B\subseteq\mathbb{T}$, we have $\mu(B)=0$ if $\theta_{0}\notin B$ and
$\mu(B)=1$ if $\theta_{0}\in B.$

###### Proof.

This is a straightforward computation. For (2), note that

$v_{1}(x,\theta;R)=\left(J_{0}+J\cos(\theta-\theta_{0})\right)\int_{D_{\rho}(0)}R_{0}(x+\tilde{x})\frac{\tilde{x}}{|\tilde{x}|}d^{n}\tilde{x}$

and

$v_{2}(x,\theta;R)=K\sin(\theta-\theta_{0})\,\int_{\Omega}R_{0}(x)\,d^{n}x.$

Using these expressions, it is easy to verify that (3.1) holds. ∎

##  4\. Analysis of Phases

To analyze the behavior of equation (2.1) in dependence of the parameters $K$,
$J$ and $J_{0}$, we perform a linear stability analysis of the incoherent
steady state $R(x,\theta)=\bar{R}=const$. The solutions of the linearized
system give indications about the possible types of far-from-equilibrium
patterns that can form. We use this to classify the emerging patterns
(‘pattern onsets’). We then compare these linearized patterns with numerical
solutions of the original equation in Section 5.

###  4.1. Linear Stability Analysis

For the linear stability analysis of the incoherent steady state
$R(x,\theta)=\bar{R}=const$, we seek solution of the form

$R(t,x,\theta)=\bar{R}+c(t)e^{i(q\theta+k\cdot x)},$

where $q$ is the wave number of the clock phase $\theta$ and $k$ is the wave
vector of the spatial coordinate $x$. The $2\pi-$periodicity of
$R(t,\theta,x)$ in $\theta$ dictates that $q\in\mathbb{Z}$ takes on integer
values. Solutions with temporally increasing coefficients $c(t)$ correspond to
emerging patterns in space and in clock synchronization. For $|q|>0$, clocks
synchronize; for nonzero $k$, spatial patterns form.

Utilizing this form in the model (2.1) and neglecting terms of order
$c(t)^{2}$ and higher yields first-order equations for $c(t)$. The following
Proposition summarizes the result of these computations in one and two spatial
dimensions in dependence of the clock wave number $q$; see the appendix A for
details on the derivation:

######  Proposition 4.1.

Let $d=1$ or $2$ denote the spatial dimension.

  1. (1)

For $q=0$, the linearization is

$\frac{dc}{dt}=\frac{D_{x}}{\rho^{2}}\sigma_{1}(|\rho k|)\,c(t)$

with $\sigma_{1}(|\rho k|)=-|\rho
k|^{2}+2\frac{J_{0}\bar{R}\rho^{d+1}}{D_{x}}\,F_{d}(|\rho k|)$; see below for
the definition of the function $F_{d}$.

  2. (2)

For $|q|\neq 0,1$, the linearization is

$\frac{dc}{dt}=(-|k|^{2}D_{x}-q^{2}D_{\theta})\,c(t)$

  3. (3)

For $|q|=1$, the linearization is

$\frac{dc}{dt}=\frac{D_{x}}{\rho^{2}}\sigma_{3}(|\rho k|)\,c(t).$

with $\sigma_{3}(\rho k)=-|\rho
k|^{2}-\frac{\rho^{2}D_{\theta}}{D_{x}}+\frac{J\bar{R}\rho^{d+1}}{D_{x}}F_{d}(|\rho
k|)+\frac{\bar{R}K\rho^{d+2}}{D_{x}}G_{d}(|\rho k|)$; see below for the
definitions of $F_{d}$ and $G_{d}$.

Here for each possible spatial dimension $d=1,2$, the pair $F_{d}(k^{*})$ and
$G_{d}(k^{*})$ are two explicit functions defined via certain integrals. In
one spatial dimension, they have the form

$F_{1}(k^{*})=2\pi(1-\cos(k^{*})),\quad
G_{1}(k^{*})=2\pi\frac{\sin(k^{*})}{k^{*}}.$

In two spatial dimensions, the functions $F_{2}$ and $G_{2}$ can be expressed
explicitly via Bessel and Struve functions; see the Appendix A for their
explicit functional forms.

The formula for $F_{1}$ is the same as in the clockless adhesion model of
Armstrong et al. [2] and was given there. Plots of the graphs of $F_{2}$ and
$G_{2}$ are given in Figure 1.

Note that the temporal eigenvalue in case 2 is always negative. Therefore, in
the linear stability analysis, only two cases of dependence on the clock phase
$\theta$ need to be considered: Either uniform dependence (wavenumber $q=0$ in
$\theta$-space) or waves with one mode ($q=\pm 1$).

For the other two cases in Proposition 4.1, there is a maximum temporal growth
rate which corresponds to a certain mode in $\theta-$ and $x-$space. For $q=0$
(case 1), we denote the maximum growth rate by $\lambda_{1}$; it is up to a
constant factor the maximum of the function $k^{*}\mapsto\sigma_{1}(k^{*})$.
For $q=\pm 1$ (case 3), it is $\lambda_{2}$, the maximum (up to a constant
factor) of $k^{*}\mapsto\sigma_{3}(k^{*})$:

$\displaystyle\lambda_{1}$ $\displaystyle=\max_{k^{*}\geq
0}\frac{D_{x}}{\rho^{2}}\sigma_{1}(k^{*})$ $\displaystyle\lambda_{2}$
$\displaystyle=\max_{k^{*}\geq 0}\frac{D_{x}}{\rho^{2}}\sigma_{3}(k^{*})$

In addition, we denote by $k_{2}$ the value where the maximum of $\sigma_{3}$
is attained111We assume the smallest possible value of $k_{2}$ if the maximum
is attained at multiple points.:

$k_{2}={\rm{argmax}}_{k^{*}\geq 0}\sigma_{3}(k^{*})$

Figure 1. Graphs of the functions $F_{2}(k^{*})$ (red, solid) and
$G_{2}(k^{*})$ (blue, dashed); see Proposition 4.1.

###  4.2. Classification of Pattern Onsets

By the previous analysis, linearized solutions $R(t,\theta,x)$ are sums of the
functions of the form $e^{\lambda t}e^{ikx}$ and $e^{\lambda t}e^{\pm
i\theta+ik\cdot x}$ with temporal growth rate $\lambda$ and spatial wavevector
$k$. This indicates the onset of certain types of possible patterns. Steady
state patterns can be categorized by their degree of synchronization (globally
synchronized, locally synchronized or incoherent, i.e. unsynchronized) and
whether they are spatially uniform, i.e. without spatial variations (patterns)
in the total spatial density
$R_{\text{tot}}(t,x)=\int_{0}^{2\pi}R(t,x,\theta)d\theta$. Specifically, we
classify the linearized solutions into six groups - uniform incoherent (UI),
uniform globally synchronized (UGS), unifom locally synchronized (ULS),
aggregated incoherent (AI), aggregated locally synchronized (ALS) and
aggregated globally synchronized (AGS). For this, consider linearized
solutions of the form

$\displaystyle R$ $\displaystyle(t,x,\theta)=\bar{R}+\Delta
R(t,x,\theta)\text{, where}$ $\displaystyle\Delta
R(t,x,\theta)=C_{1}e^{\lambda_{1}t}\exp(ik_{1}\cdot
x)+e^{\lambda_{2}t}(C_{2,1}\exp(i\theta+ik_{2}\cdot
x)+C_{2,2}\exp(-i\theta+ik_{2}\cdot x)).$

Here $\lambda_{1}$ and $\lambda_{2}$ are temporal growth rates of the
corresponding modes. We assume that these modes are the fastest growing, i.e.
that $\lambda_{1}$ is the maximum growth rate for the modes with $q=0$ (case 1
in Proposition 4.1) with corresponding spatial wavevector $k_{1}$, and
$\lambda_{2}$ is the maximum growth rate for $q=\pm 1$ (case 3 in Proposition
4.1) with corresponding $k_{2}$. If the growth rate $\lambda_{i}$ is positive,
the corresponding mode will grow in time; for negative $\lambda_{i}$, it will
decay. This gives a classification of pattern onsets: Each combination of
growth rates gives another category; for $\lambda_{2}>0$, we also need to
distinguish the cases $k_{2}=0$ [no spatial pattern, i.e. spatially uniform
case] and $k_{2}\neq 0$. Considering only terms with positive temporal growth
rate, we then obtain the classification summarized in Table 2 and Figure 2.
Here we also list for each case an example for the functional form of the
linearized pattern $\Delta R(x,\theta)$, mostly for the purpose of plotting an
example for each pattern onset (Figure 2).

$\lambda_{1}$ | $\lambda_{2}$ | $k_{2}$ | example $\Delta R(x,\theta)$ | Name
---|---|---|---|---
$\lambda_{1}<0$ | $\lambda_{2}<0$ |  | $0$ | uniform incoherent (UI)
$\lambda_{1}>0$ | $\lambda_{2}<0$ |  | $-\cos(k_{1}x)$ | aggregated incoherent (AI)
$\lambda_{1}<0$ | $\lambda_{2}>0$ | $k_{2}=0$ | $-\cos(\theta)$ | unif. globally synchronized (UGS)
$\lambda_{1}<0$ | $\lambda_{2}>0$ | $k_{2}\neq 0$ | $-\cos(\theta)\cos(k_{2}x)$ | unif. locally synchronized (ULS)
$\lambda_{1}>0$ | $\lambda_{2}>0$ | $k_{2}=0$ | $-\cos(k_{1}x)+\cos(\theta)$ | aggr. globally synchronized (AGS)
$\lambda_{1}>0$ | $\lambda_{2}>0$ | $k_{2}\neq 0$ | $-\cos(k_{1}x)-\cos(\theta)\cos(k_{2}x)$ | aggr. locally synchronized (ALS)
Table 2. Classification of linearized solutions (pattern onsets); see Figure 2

Figure 2. Graphical representations of the six pattern onsets listed in Table
2 in one spatial dimension. For each pattern onset, the top panel shows a heat
map of the linearized patterns $\Delta R(x,\theta)$ in $x-\theta-$ space, the
bottom panel shows the corresponding pattern in the total spatial density
$\Delta R_{\rm{tot}}(x)=\int_{0}^{2\pi}\Delta R(x,\theta)\,d\theta$. For
presentation purposes, we scaled densities to take values between $-1$ and
$+1$. Here the length of the spatial interval is $L=10$. For AI, AGS and ALS,
we chose $k_{1}=\frac{3\pi}{5}$. Additionally for ULS, $k_{2}=\frac{3\pi}{5}$
and for ALS, $k_{2}=\frac{\pi}{5}$.

###  4.3. Phase Diagram of Pattern Onsets

The linearization result in Proposition 4.1 along with the classification of
possible pattern onsets (Table 2 and Figure 2) makes it possible to determine
a phase diagram of possible behavior of the system (2.1) according to the
parameters $J_{0},J$ and $K$. Define first the nondimensional parameters

$J^{*}=\frac{\bar{R}\rho^{d+1}}{D_{x}}J,\quad
J_{0}^{*}=\frac{\bar{R}\rho^{d+1}}{D_{x}}J_{0},\quad
K^{*}=\frac{\bar{R}\rho^{d+2}}{D_{x}}K$

With this notation, the function $\sigma_{1}$ and $\sigma_{3}$ can be written
as

$\displaystyle\sigma_{1}(k^{*})$
$\displaystyle=-(k^{*})^{2}+2J_{0}^{*}\,F_{d}(k^{*})$
$\displaystyle\sigma_{3}(k^{*})$
$\displaystyle=-(k^{*})^{2}-D_{\theta}^{*}+J^{*}F_{d}(k^{*})+K^{*}G_{d}(k^{*})$

where $D_{\theta}^{*}=\rho^{2}D_{\theta}/D_{x}$. We restrict ourselves to
parameter ranges that are consistent with the assumption that a meaningful
cell density can be defined; see Appendix B.

The classification of patterns in Table 2 is given in terms of the growth
rates $\lambda_{1}$ and $\lambda_{2}$, which up to constants are the maxima of
the functions $\sigma_{1}$ and $\sigma_{3}$. Thus which onset pattern occurs
only depends on the dimensionless values $J_{0}^{*},J^{*},K^{*}$ and
$D_{\theta}^{*}$. Note that we assume that $D_{\theta}^{*}\ll 1$, and in fact
we are most interested in the case $D_{x}^{*}=0$, i.e. the absence of
diffusion in the clock phase $\theta$. We concentrate on this case
$D_{\theta}^{*}=0$ in this section.

Pattern onsets depend on the signs of $\lambda_{1}$ and $\lambda_{2}$ and on
whether $k_{2}$ is zero or nonzero (Table 2). The following two propositions
give us conditions on the parameters $J^{*},K^{*}$ and $J^{*}_{0}$. First, the
sign of $\lambda_{1}$ is determined by $J_{0}$:

######  Proposition 4.2.

In each dimension $d=1,2$ there is a value
$J^{d}_{\rm{crit}}=\frac{1}{F^{\prime\prime}_{d}(0)}$ such that

$\lambda_{1}<0\text{ for }J_{0}^{*}<J^{d}_{\rm{crit}}\quad\text{ and
}\lambda_{1}>0\text{ for }J^{*}_{0}>J^{d}_{\rm{crit}}.$

One computes

$\displaystyle J^{1}_{\rm{crit}}$ $\displaystyle=\frac{1}{2\pi}\text{ in 1D}$
$\displaystyle J^{2}_{\rm{crit}}$ $\displaystyle=\frac{3}{2\pi^{2}}\text{ in
2D}$

###### Proof.

The sign of $\lambda_{1}$ is determined by whether $\sigma_{1}(k^{*})$ has a
positive global maximum for $k^{*}\neq 0$ or not. By inspection, we have
$\sigma_{1}(0)=\sigma_{1}^{\prime}(0)=0$. It is thus the sign of
$\sigma^{\prime\prime}(0)=-2+2J_{0}^{*}F^{\prime\prime}_{d}(0)$ which
determines the sign of $\lambda_{1}$. The above condition on $J^{*}_{0}$
follows immediately. ∎

In terms of the model, Proposition 4.2 simply means that in every dimension,
there is a critical value for $J^{*}_{0}$ above which spatial aggregation
occurs spontaneously from incoherent, spatially homogeneous initial
conditions.

The sign of $\lambda_{2}$ and whether the condition $|k_{2}|=0$ is satisfied
is determined by $J^{*}$ and $K^{*}$. The $K^{*}J^{*}$-parameter plane is
divided into three regions:

######  Proposition 4.3.

In each dimension $d=1,2$ there is a functions $f_{d}(K^{*})$ such that the
following holds. (See Figure 3.)

For any point $(K^{*},J^{*})$ in admissible parameter space (see Appendix B),
we have:

  * •

If $K^{*}<0$ and $J^{*}<f_{d}(K^{*})$, then $\lambda_{2}<0$.

  * •

If $K^{*}>0$ and $J^{*}<f_{d}(K^{*})$, then $\lambda_{2}>0$ and $k_{2}=0$.

  * •

If $J^{*}>f_{d}(K^{*})$, then $\lambda_{2}<0$ and $k_{2}>0$.

The function $f_{d}(k^{*})$ is piecewise defined as follows:

For $K^{*}\geq 0$,

$f_{d}(K^{*})=\frac{2-K^{*}G_{d}^{\prime\prime}(0)}{F_{d}^{\prime\prime}(0)}=\begin{cases}\frac{1}{\pi}+\frac{1}{3}K^{*}&\text{
for dimension }d=1\\\ \frac{3}{\pi^{2}}+\frac{3}{8}K^{*}&\text{ for dimension
}d=2\end{cases}$

For $K^{*}<0$, the graph of $f_{d}(K^{*})$ can be parameterized via

$\left(\begin{array}[]{c }K^{*}\\\
f_{d}(K^{*})\end{array}\right)=\frac{1}{G_{d}(\kappa)F_{d}^{\prime}(\kappa)-F_{d}(\kappa)G^{\prime}_{d}(\kappa)}\left(\begin{array}[]{c
c}\kappa^{2}F^{\prime}_{d}(\kappa)-2\kappa F_{d}(\kappa)\\\
-\kappa^{2}G^{\prime}_{d}(\kappa)+2\kappa
G_{d}(\kappa)\end{array}\right)\quad(\kappa>0).$

###### Proof.

Recall the definition that $\lambda_{2}$ is the (up to a constant) the maximum
of $\sigma_{3}(k^{*})=-(k^{*})^{2}+J^{*}F_{d}(k^{*})+K^{*}G_{d}(k^{*})$ for
$k^{*}\geq 0$ and $k_{2}$ is the smallest value at which it is attained. Note
that $\sigma_{3}(0)=K^{*}G_{d}(0)$ with $G_{d}(0)>0$.

Consider first the case $K^{*}>0$. Then $\lambda_{2}>0$. Moreover, by
inspection of $F_{d}$ and $G_{d}$, for $J^{*}\leq 0$, the maximum of
$\sigma_{3}(k^{*})$ occurs at $k^{*}=0$, so $k_{2}=0$ then. (See also Figure
1.) For small positive $J^{*}$, the maximum will still occur at $k^{*}=0$. As
$J^{*}$ increases for fixed $K^{*}>0$, there is some critical value
$J^{*}=f_{d}(K^{*})$ above which the positive maximum occurs at a value
$k^{*}>0$. The condition for this critical value is that
$\sigma_{3}^{\prime\prime}(0)=0$. This means

$\sigma_{3}^{\prime\prime}(0)=-2+J^{*}F^{\prime\prime}_{d}(0)+K^{*}G^{\prime\prime}_{d}(k^{*})=0\quad\text{
at }J^{*}=f_{d}(K^{*})$

This gives the formula for $f_{d}(K^{*})$ for $K^{*}>0.$

Now consider the case $K^{*}<0$. Then $\sigma_{3}(0)<0$. If $J^{*}<0$, then
the maximum of $\sigma_{3}(k^{*})$ is negative by our assumptions on the
permissible range of $J^{*}$ and $K^{*}$ (see Appendix B), so $\lambda_{2}<0$
then. For positive $J^{*}$, as $J^{*}$ increases, there is a critical value
$J^{*}=f_{d}(K^{*})$ above which a positive maximum occurs for some value
$k^{*}=\kappa>0$. The condition for this is

$\left\\{\begin{array}[]{l}*\sigma_{3}(\kappa)=-\kappa^{2}+J^{*}F_{d}(\kappa)+K^{*}G_{d}(\kappa)=0,\\\
\sigma^{\prime}_{3}(\kappa)=-2\kappa+J^{*}F^{\prime}_{d}(\kappa)+K^{*}G^{\prime}_{d}(\kappa)=0\end{array}\right\\}\quad\text{
at }J^{*}=f_{d}(K^{*}).$

Solving this for $(K^{*},J^{*})$ gives a paramaterization (with parameter
$\kappa>0$) for the graph of $J^{*}=f_{d}(K^{*})$ for $K^{*}<0$. ∎

Propositions 4.2 and 4.3 in conjunction with Table 2 now allow for a complete
classification of the onset of patterns in dependence of the dimensionless
parameters $J^{*}_{0}$, $J^{*}$ and $K^{*}.$ Figure 4 summarizes this
graphically.

Figure 3. Graphical representation of the result in Proposition 4.3. Note that
the images are slightly different in one and two spatial dimensions, although
the overall appearance is very similar. Figure 4. Phase diagrams of onsets of
patterning in $J^{*}K^{*}-$space. We show the case of $d=1$ spatial dimension;
the case $d=2$ is very similar (Figure 3). See Table 2 for the classification
of the different pattern onsets. Results are based on Propositions 4.2 and
4.3. Left diagram: Case $J_{0}^{*}<J^{1}_{\rm{crit}}.$ In this case, there is
no spatial aggregation and pattern onsets are nonconstant in clock space only.
The three possibilities are uniform incoherent (UI), uniform locally
synchronized (ULS) and uniform globally synchronized (UGS). Right diagram:
Case $J_{0}^{*}>J^{1}_{\rm{crit}}.$ In this case, spatial aggregation occurs
simultaneously with patterning in clock space only. The three possibilities
are aggregated incoherent (AI), aggregated locally synchronized (ALS) and
aggregated globally synchronized (AGS). For the reader’s convenience, we
include the example densities patterns in $x\theta$-space from Figure 2.

##  5\. Far-From Equilibrium Behavior and Numerical Simulations

In the previous section, we used linear stability analysis to classify the
onset of patterns on the oscillator aggregation model from a spatially
homogeneous, incoherent state (summarized in Figure 4). This does not
necessarily give information about the far-from equilibrium behavior of the
system. There is one exception to this: the UI (uniform incoherence) case.
Then all modes are linearly stable. As a consequence, no patterns can form in
physical space or clock phase space. (For initial conditions that are
homogeneous in the spatial coordinate $x$, i.e.,
$R(t=0,x,\theta)=R_{0}(\theta)$, it is even possible to write out an explicit
solution using Fourier analsyis and rigorously show that for $K<0$, the
solution converges to the incoherent state, see Appendix C.)

In general, a natural expectation is that pattern onsets will generally lead
to steady state patterns with similar behavior. Specifically, if spatial
aggregation occurs in the pattern onset, i.e., if the total spatial density
$\bar{R}(t,x)=\int_{0}^{2\pi}R(t,x,\theta)d\theta$ has spatial peaks and
valleys, this will lead to far-from equilibrium patterns with spatial
aggregation.

Similarly, it is a natural expectation that pattern onsets of the “globally
synchronized” variety (UGS, AGS) actually give rise to globally synchronized
far-from-equilibrium patterns. In other words, for zero clock diffusivity
($D_{\theta}=0$), one expects convergence to a steady state pattern of the
form $R(x,\theta)=R_{0}(x)\delta(\theta-\theta_{0})$ for some fixed
$\theta_{0}$. This means that all oscillators are perfectly synchronized and
the total spatial density is given by $R_{0}(x)$; see also Theorem 2. Note
then that a linear stability can be applied to $R_{0}(x)$. A formal argument
gives the same result as Proposition 4.1 (1) with the only change being that
we have to replace $J_{0}$ by $J_{0}+J$. By Proposition 4.2, this means that
now we expect patterns in the total spatial density for
$J^{*}_{0}+J^{*}>J^{d}_{\rm{crit}}$.

These arguments yield that for the pattern onsets UI and AGS, one expects
similar behavior of the far-from equilibrium patterns of uniform
(unaggregated) incoherent patterns and globally synchronized patterns with
spatial aggregations, respectively. In the case of UGS (uniform globally
synchronized) pattern onsets in $d$ dimensions, the above arguments suggest
that for $J^{*}_{0}+J^{*}<J^{d}_{\rm{crit}}$, the far-from-equilibrium pattern
should be globally synchronized and also spatially uniform (i.e., without
spatial aggregations). But for $J^{*}_{0}+J^{*}>J^{d}_{\rm{crit}}$, one
expects that even the UGS pattern onset will give rise to aggregations once
cells are sufficiently synchronized. Similar arguments apply to the ULS
pattern onset.

In this section, we compare the pattern onsets from linear stability analysis
with the far-from equilibrium patterns obtained by numerical simulations. We
concentrate on the case $J_{0}=0$. This means that cells with a
$\frac{\pi}{2}$-phase shift neither repel nor attract each other; in
particular, cells in opposite phases repel each other and there is no spatial
aggregation in the incoherent phase. We investigate the behavior in dependence
of the parameter $J^{*}$ and $K^{*}$. We consider two computational models:
The first is a numerical solution of the system (2.1) in one spatial
dimensions. The second is an implementation of a Lattice-Gas Cellular
Automaton (LGCA) in two spatial dimensions representing a discretized version
of the PDE model (2.1).

###  5.1. Numerical Solution of the PDE 2.1 in one spatial dimension

We used the finite element solver FreeFEM++ [25] to solve the PDE (2.1) in one
spatial dimension for $J_{0}=0$. Figure 6 summarizes the results with steady
state patterns shown for different values of $J^{*}$ and $K^{*}$ obtained from
initial conditions with random small perturbations from the incoherent
spatially uniform state $R_{0}(x,\theta)=1$.

To quantify the degree of synchronization of a distribution $R(t,x,\theta)$,
we used the $r$-parameter first introduced by Kuramoto ($\sigma$ in Equation
(5.4.13) in [33]). Here this is given by

$r=r(t)=\left|\int_{\Omega}\int_{0}^{2\pi}e^{i\theta}R(t,x,\theta)d\theta
d^{n}x\middle/\int_{\Omega}\int_{0}^{2\pi}R(t,x,\theta)d\theta d^{n}x\right|$
(5.1)

It represents the modulus of the mean value of $e^{i\theta}$. Note that $0\leq
r\leq 1$. In the completely incoherent case $R(t,x,\theta)=R(t,x),$ we have
$r=0$. In contrast, perfect synchronization
($R(t,x,\theta)=R_{c}(t,x)\delta(\theta-\theta_{0})$ for some
$\theta_{0}\in[0,2\pi)$) corresponds to $r=1$.

The phase diagram in Figure 5 summarizes the results using a combination of
the linear stability analysis of pattern onsets (Figure 4, left panel) and
example numerical computations for certain pairs of parameters $(J^{*},K^{*})$
denoted by Roman numerals I–VII (Figure 6). The parameter space splits up into
five regions separated by the curves $K^{*}=0$,
$J^{*}=J^{1}_{\rm{crit}}=1/2\pi$ ($K^{*}>0$) and $J^{*}=f_{1}(K^{*})$ (see
Proposition 4.3). For the UI onset case, solutions also display spatial
uniformity and clock incoherence (pairs I and II). The region of the UGS onset
case displays two different possible behaviors: In the region
$J^{*}<J^{1}_{\rm{crit}}$, we obtain global synchronization without spatial
aggregation (IV, V). In the slice of this region where
$J^{*}>J^{1}_{\rm{crit}}$, there is some spatial aggregation with a small
amplitude in addition to global synchronization in the nonlinear system even
though the linear pattern onset is spatially uniform (VI). The reason is that
once synchronized, the adhesion strength between cells becomes sufficiently
strong to cause aggregation.

The ULS onset case is the most interesting. Our results indicate that for
$K^{*}>0$, the far-from-equilibrium behavior is as suggested by the linear
stability analysis, but with additional subsequent spatial aggregation: Local,
but not global synchronization occurs, followed by spatial aggregation into
segregated clusters. These clusters are synchronized within themselves, but
not among each other (VII). But for $K^{*}<0$, our results indicate a new kind
of phase that is substantially different from all the linear patterns (III).
This is the case of what in [37] is referred to as a “phase wave”: Here
spatial position is correlated with clock phase, i.e. the mean clock $\theta$
increases linearly with position $x$. We also observe some spatial aggregation
along with the establishment of the phase wave pattern.

Figure 5. Phase diagram of far-from-equilibrium patterning in one spatial
dimension for Equation (2.1) in $J^{*}K^{*}-$space (with $J^{*}=0$). This is
based on the onset diagram in Fig 3 (left panel). The horizontal line on the
$K^{*}>0$-half plane corresponds to the critical value $J^{d}_{\rm{crit}}$ in
$d=1$ dimension. For example numerical computation for the parameters
indicated by Roman numerals, see Figure 6. Specifically, $(J^{*},K^{*})$
values are: I $(-1,-1)$, II $(0.4,-1)$, III $(1,-1)$, IV $(-1,1)$, V
$(0.05,1)$, VI $(0.4,1)$, VII $(1,1)$.  Figure 6. Steady state distributions
from numerical solutions of (2.1) with the parameters indicated in Figure 5.
For each simulation, the graph of the Kuramato $r-$parameter as a function of
time is shown. Parameters were $J_{0}=0$, $D_{x}=1,D_{\theta}=0.05,\rho=1$ on
the spatial interval $[0,10]$ with periodic boundary conditions. A logistic
function with $R_{{\max}}=30$ was chosen in the spatial flux term, see (2.4).

###  5.2. A Lattice-Gas Cellular Automaton model in two spatial dimensions

The solution of the PDE (2.1) is numerically challenging in the case $K>0$,
i.e., when neighboring clock tend to synchronize. In the case of zero clock
diffusion ($D_{\theta}=0$), this leads to solutions with Dirac-delta steady
state solutions (see Proposition 2). But even for small nonzero diffusion
coefficient $D_{\theta}$, solutions display sharp gradients along lines of
constant $\theta$; see Figure 6.

To investigate spatial patterns in two dimensions and to confirm the one-
dimensional results, we set up a Lattice-Gas Cellular Automaton model based on
models of cell-cell adhesion and aggrgeation ([8] and [13] (chapter 7)); see
Figure 7. The model is using a 2-dimensional $50\times 50$ hexagonal lattice
with periodic boundary conditions. Each hexagonal node has six “channels.”
Each channel corresponds to a direction of motion. Each particle occupies one
channel in a hexagon. Additionally, each particle has a clock associated with
it. The “state” of a node is the collection of occupied channels in it with
corresponding clock values, i.e., collections $\\{s_{i}\\}_{i=1,\ldots,6}$,
where $s_{i}$ is the occupancy number (0 or 1) of the $i$th channel, and
corresponding clock values $\\{\theta_{i}\\}_{i=1,\ldots,6}$.

The key idea of the model is that spatiotemporal evolution is modeled via two
steps: a propagation step and an interaction step. (See Figure 7.) In the
propagation step, cells move from one node to a neighboring one in the
direction given by their channel. In the interaction step, the state of each
node is updated. For this, each possible rearrangement of the particles in a
node (state ${\mathbf{s}}=\\{s_{i}\\}_{i=1,\ldots,6}$ with clock values
${\bm{\theta}}=\\{\theta_{i}\\}_{i=1,\ldots,6}$) is assigned an energy
$H({\mathbf{s}},{\bm{\theta}})$. A new state $({\mathbf{s}},{\bm{\theta}})$ is
then chosen with probability

$P({\mathbf{s}},{\bm{\theta}})=\frac{1}{Z}\exp(H({\mathbf{s}},{\bm{\theta}})),$

where $Z$ is a constant chosen such that the sum of all probabilities is
unity. In our case, the energy $H({\mathbf{s}},{\bm{\theta}})$ is calculated
by considering how each particle interacts with the neighboring particles.
Specifically,

$H({\mathbf{s}},{\bm{\theta}})=\sum_{\begin{subarray}{c}\text{channel }\\\
i=1,\ldots,6\end{subarray}}\sum_{\text{neighbors}}s_{i}\,\vec{c}_{i}\cdot\vec{p}_{\text{neighbor}}\,(J_{0}+J\cos(\theta_{\text{neighbor}}-\theta_{i}))$
(5.2)

Here $\vec{c}_{i}$ is the vector of direction corresponding to channel $i$ and
$\vec{p}_{\text{neighbor}}$ is a vector in the direction of a (first order)
neighbor’s lattice site; see Figure 8. The clock values of the particle in
channel $i$ and the neighbor are $\theta_{i}$ and $\theta_{\text{neighbor}}$,
respectively. The function $H({\mathbf{s}},{\bm{\theta}})$ encodes a flux
along the direction of highest cell density, modulated by a term depending on
the cell clock. (The same approach – minus oscillations – is behind Armstrong
et al’s PDE model of cell-cell adhesion [2] which we adapted for our PDE model
(2.1); see also [9].)

Finally, the clocks are updated at every LGCA step based on a discrete
Kuramoto-type model. That is, the change in clock value of the particle at
channel $i$ is given by

$\Delta\theta_{i}=\omega\Delta
t+K\frac{1}{N_{i}}\sum_{j=1}^{N_{i}}\sin(\theta_{j}-\theta_{i})$

Here the sum is over the first order neighbors of the particle in channel $i$
and the particles in the other channels of the same hexagonal node (total
$N_{i}$ particles). As in the PDE model (2.1), we can assume without loss of
generality that $\omega=0$, by shifting the frame in clock phase by the term
$\omega t$, where $t$ is the number of time steps.

As in the PDE case (Equation (5.1)), one can define an $r-$parameter which
quantifies the degree of local synchronization. In the LGCA model, this is
given by

$r=\left|\frac{1}{N}\sum_{k=1}^{N}e^{i\theta_{k}}\right|,$

where $N$ is the total number of particles. Again $r=1$ means complete
synchronization and $r=0$ complete incoherence. We also considered
$r_{\rm{local}}$, a local variation of $r$, where the sum is taken at each
node over the first-order neighbors only and then averaged over all nodes.
This quantifies the degree of local synchronization.

Note that the model contains parameters $J_{0},J$ and $K$ whose meaning is
completely analogous to that of the PDE model (2.1), and equally the
interpretations in Table 1 apply here.

One thus expects similar behavior as the displayed by PDE version (2.1)
displayed in Figure 5. (Note that the phase diagram is for one spatial
dimension, but as shown in Figure 3, the 2D case is substantially the same, at
least its linearization.) To test this, we conducted numerical simulations of
this LGCA model with $J_{0}=0$ and various values of $J$ and $K$. As
illustrated in Figure 9, this two dimensional LGCA model shows the same phases
as the PDE model with the possible exception of the “aggregated global
synchronized” phase. While the results of Figure 5 do not apply quantitatively
to the LGCA model, their qualitative behavior is equivalent: For negative $K$
and $J$, the LGCA model also displays incoherent clock phases and lack of
aggregation. For positive $K$ and negative $J$, global synchronization occurs
without spatial aggregation. For positive $J$ and $K$, we observed local
aggregation into separate clusters, with synchronization within, but not among
clusters. (Figure 5 predicts a range of values of positive, small $J$ and
large $K$ where aggregation occurs along global synchronization. In the LGCA,
we only observed separation into clusters where most, but not all clusters
were synchronized222We note that for large $K>0$, our method for updating the
clock phases appears to undergo a period-doubling bifurcation so that
neighboring cells do not completely synchronize even in the absence of
movement. This a numerical artifact which goes against the spirit of the
Kuramoto model and we restricted our computations to values of $K>0$ where
local synchronization does occur. .) Finally, for negative values of $K$ close
to 0 and large values of $J$, we observed a phase apparently equivalent to the
“aggregated phase wave” in Figure 5. Here, some spatial aggregation occured,
but the individual clusters did not synchronize. They were not completely
incoherent either though, as indicated by its local $r-$value
$r_{\rm{local}}$, which was much larger than for the incoherent phase; see
Figure 9. This indicates a partial sorting by phases within clusters.

Figure 7. Schematics of the Lattice-Gas Cellular Automaton model. Figure 8.
Sketch for Equation (5.2) Figure 9. Typical results for $J_{0}=0$ and
different values of $J$ and $K$ in the LGCA model. Note that qualitatively,
these phases correspond to the ones obtained with the PDE model for Figure 5.
Different phases are indicated by colors; see the color wheel for reference.
The images show configurations obtained from random initial conditions after
2000 iterations on a 50$\times$50 hexagonal grid with periodic boundary
conditions. There were $N=6,000$ particles in each simulation, corresponding
to a confluency (percentage of occupied channels overall) of $40\%$.

##  6\. Discussion

We have investigated the interplay of oscillators (cells with intracellular
’clocks’) that undergo random motion and adhere to each other. The strength of
adhesion depends on the relative clock phases and neighbors can sync or anti-
sync. A linear stability analysis allows for classification of pattern onsets.
We further investigated far-from equilibrium patterns for the case that cells
with a $\frac{\pi}{2}$-phase shift neither repel nor attract each other
($J_{0}=0$). We identified five different types of patterns: a completely
synchronized state without spatial aggregation; aggregation into clusters that
are internally synchronized, but out of sync with each other; an incoherent,
unaggregated state; and a phase wave state with diffuse clusters and spatial
gradients of clock phases. A fifth state shows aggregation into diffuse,
globally synchronized clusters. This behavior is found only in a small region
of parameter space and was only approximately achieved in an LGCA
implementation of our model.

A number of open issues present themselves: First, we investigated the far-
from-equilibrium behavior only in the natural case of zero baseline adhesion
$(J_{0}=0)$. Our pattern onset classification extends to positive and negative
$J_{0}$ as well. In particular, the case $J_{0}>J^{d}_{\rm{crit}}$ (so that
cells aggregate even in the incoherent state) promises other types of behavior
and an in-depth investigation could be fruitful. We also did not investigate
the spatial dynamics of synchronization in our model in detail, i.e., how
exactly synchronization patterns arise in time and space, which again is an
interesting avenue for further work.

There are also interesting analytical questions: Numerical simulations
indicate that for postive $K$ and negative $J$, solution converge to a Dirac-
delta point measure in clock-space and uniform distribution in space when the
clock diffusivity is zero ($D_{\theta}=0$). Can this be proved rigorously?
Even more general, do solutions always converge to function of the form
$R(x)\delta(\theta-\theta_{0})$ for $K>0$? For these question, a further
examination of the properties of solutions of the measure-theoretic
formulation sketched in section 3 would be helpful (e.g. proving existence and
uniqueness).

As discussed in the introduction, one of our much broader future goal is to
use this model and others of its type to investigate generic mechanisms of
aggregation, synchronization and pattern formation in cell populations. This
pertains to organisms on the threshold to multicellularity like Myxobacteria
[28], and Dictyostelia [34], and also to developmental processes like pattern
formation mediated by oscillatory signaling in the Notch-Delta circuit [31]).
Our hope is to explore broad general principles in the spirit of the
conceptual framework of Dynamical Pattern Modules (DPMs) [35]. Such DPMs are
units of gene products associated with physical processes they mobilize,
forming a proposed “pattern language” for pattern formation inm developmental
biology.

Finally, there are some interesting sociological interpretations of our model
in the context of polarization of political opinions and the formation of echo
chambers aided by social media [11]. (Mathematical models of these phenomena
typically use network models, e.g. [4].) Our model’s behavior in the
aggregated, locally synchronized state (clusters that are synchronized within
themselves, but not between each other) shows clear formal similarities to
such echo chambers. Our results indicate that this is a natural outcome if
neighbors tend to synchronize and “like attracts like” ($J>0$). One way to
achieve global synchronization without formation of such echo chambers is then
to consider negative $J$; the case of “opposites attract.” This suggest that a
way out of the societal polarization is to get people with different views in
contact – a truism maybe, but this certainly suggests that the model has the
potential for less trivial sociological investigations.

### Acknowledgements

The authors acknowledge the financial support of the John Templeton Foundation
(#62220). The opinions expressed in this paper are those of the authors and
not those of the John Templeton Foundation.

## References

  * [1] J. A. Arias Del Angel, V. Nanjundiah, M. Benítez, and S. A. Newman, Interplay of mesoscale physics and agent-like behaviors in the parallel evolution of aggregative multicellularity, EvoDevo, 11 (2020), pp. 1–18.
  * [2] N. J. Armstrong, K. J. Painter, and J. A. Sherratt, A continuum approach to modelling cell-cell adhesion, Journal of Theoretical Biology, 243 (2006), pp. 98–113.
  * [3] A. Barciś and C. Bettstetter, Sandsbots: Robots that sync and swarm, IEEE Access, 8 (2020), pp. 218752–218764.
  * [4] F. Baumann, P. Lorenz-Spreen, I. M. Sokolov, and M. Starnini, Modeling echo chambers and polarization dynamics in social networks, Physical Review Letters, 124 (2020), p. 048301.
  * [5] R. Bhat, T. Glimm, M. Linde-Medina, C. Cui, and S. A. Newman, Synchronization of hes1 oscillations coordinates and refines condensation formation and patterning of the avian limb skeleton, Mechanisms of Development, 156 (2019), pp. 41–54.
  * [6] S. Boccaletti, A. N. Pisarchik, C. I. d. Genio, and A. Amann, Synchronization: From Coupled Systems to Complex Networks, Cambridge University Press, Mar. 2018.  Google-Books-ID: y8tMDwAAQBAJ.
  * [7] J. T. Bonner, The social amoebae, in The Social Amoebae, Princeton University Press, 2008.
  * [8] H. J. Bussemaker, Analysis of a pattern-forming lattice-gas automaton: mean-field theory and beyond, Physical Review E, 53 (1996), p. 1644.
  * [9] A. Buttenschön, T. Hillen, A. Gerisch, and K. J. Painter, A space-jump derivation for non-local models of cell–cell adhesion and non-local chemotaxis, Journal of Mathematical Biology, 76 (2018), pp. 429–456.
  * [10] L. Chen, K. Painter, C. Surulescu, and A. Zhigun, Mathematical models for cell migration: a non-local perspective, Philosophical Transactions of the Royal Society B, 375 (2020), p. 20190379.
  * [11] M. Cinelli, G. De Francisci Morales, A. Galeazzi, W. Quattrociocchi, and M. Starnini, The echo chamber effect on social media, Proceedings of the National Academy of Sciences, 118 (2021), p. e2023301118.
  * [12] V. E. Deneke, A. Melbinger, M. Vergassola, and S. Di Talia, Waves of cdk1 activity in s phase synchronize the cell cycle in drosophila embryos, Developmental cell, 38 (2016), pp. 399–412.
  * [13] A. Deutsch, S. Dormann, A. Deutsch, and S. Dormann, Cellular automata, Cellular Automaton Modeling of Biological Pattern Formation: Characterization, Examples, and Analysis, (2017), pp. 65–111.
  * [14] NIST Digital Library of Mathematical Functions.  https://dlmf.nist.gov/, Release 1.1.9 of 2023-03-15.  F. W. J. Olver, A. B. Olde Daalhuis, D. W. Lozier, B. I. Schneider, R. F. Boisvert, C. W. Clark, B. R. Miller, B. V. Saunders, H. S. Cohl, and M. A. McClain, eds.
  * [15] F. Dörfler and F. Bullo, Synchronization in complex networks of phase oscillators: A survey, Automatica, 50 (2014), pp. 1539–1564.
  * [16] J. Dyson, S. A. Gourley, and G. F. Webb, A non-local evolution equation model of cell–cell adhesion in higher dimensional space, Journal of biological dynamics, 7 (2013), pp. 68–87.
  * [17] L. C. Evans, Partial differential equations, American Mathematical Society, Providence, R.I., 2010.
  * [18] L. F. Faust, Natural History and Flash Repertoire of the Synchronous Firefly Photinus carolinus (Coleoptera: Lampyridae) in the Great Smoky Mountains National Park, Florida Entomologist, 93 (2010), pp. 208–217.  Publisher: Florida Entomological Society.
  * [19] J. A. Glazier and F. Graner, Simulation of the differential adhesion driven rearrangement of biological cells, Physical Review E, 47 (1993), p. 2128.
  * [20] T. Glimm, R. Bhat, and S. A. Newman, Modeling the morphodynamic galectin patterning network of the developing avian limb skeleton, Journal of Theoretical Biology, 346 (2014), pp. 86–108.
  * [21] T. Glimm and J. Zhang, Numerical approach to a nonlocal advection-reaction-diffusion model of cartilage pattern formation, Mathematical and Computational Applications, 25 (2020), p. 36.
  * [22] T. Gregor, K. Fujimoto, N. Masaki, and S. Sawai, The onset of collective behavior in social amoebae, Science, 328 (2010), pp. 1021–1025.
  * [23] M. Guzzo, S. M. Murray, E. Martineau, S. Lhospice, G. Baronian, L. My, Y. Zhang, L. Espinosa, R. Vincentelli, B. P. Bratton, et al., A gated relaxation oscillator mediated by frzx controls morphogenetic movements in myxococcus xanthus, Nature microbiology, 3 (2018), pp. 948–959.
  * [24] D. Hansel, G. Mato, and C. Meunier, Synchrony in excitatory neural networks, Neural Computation, 7 (1995), pp. 307–337.
  * [25] F. Hecht, New development in freefem++, J. Numer. Math., 20 (2012), pp. 251–265.
  * [26] T. Hillen, P. Hinow, and Z.-A. Wang, Mathematical analysis of a kinetic model for cell movement in network tissues, Discrete Contin. Dyn. Syst. Ser. B, 14 (2010), pp. 1055–1080.
  * [27] D. HORSTMANN, From 1070 until present: the Keller-Segel model in chemotaxis and its consequences, Jahresber. Dtsch. Math.-Ver., 105 (2003), pp. 103–165.
  * [28] O. A. Igoshin, A. Goldbeter, D. Kaiser, and G. Oster, A biochemical oscillator explains several aspects of myxococcus xanthus behavior during development, Proceedings of the National Academy of Sciences, 101 (2004), pp. 15760–15765.
  * [29] M. Iwasa, K. Iida, and D. Tanaka, Hierarchical cluster structures in a one-dimensional swarm oscillator model, Physical Review E, 81 (2010), p. 046220.
  * [30] Y.-J. Jiang, B. L. Aerne, L. Smithers, C. Haddon, D. Ish-Horowicz, and J. Lewis, Notch signalling and the synchronization of the somite segmentation clock, Nature, 408 (2000), pp. 475–479.
  * [31] R. Kageyama, H. Shimojo, and A. Isomura, Oscillatory control of notch signaling in development, Molecular Mechanisms of Notch Signaling, (2018), pp. 265–277.
  * [32] E. F. Keller and L. A. Segel, Initiation of slime mold aggregation viewed as an instability, Journal of theoretical biology, 26 (1970), pp. 399–415.
  * [33] Y. KURAMOTO, Chemical oscillations, waves and turbulence, Synergetics, 19 (1984).
  * [34] A. F. Marée and P. Hogeweg, How amoeboids self-organize into a fruiting body: multicellular coordination in dictyostelium discoideum, Proceedings of the National Academy of Sciences, 98 (2001), pp. 3879–3883.
  * [35] S. A. Newman and R. Bhat, Dynamical patterning modules: a” pattern language” for development and evolution of multicellular form, International Journal of Developmental Biology, 53 (2009), pp. 693–705.
  * [36] K. O’Keeffe, S. Ceron, and K. Petersen, Collective behavior of swarmalators on a ring, Physical Review E, 105 (2022), p. 014211.
  * [37] K. P. O’Keeffe, H. Hong, and S. H. Strogatz, Oscillators that sync and swarm, Nature Communications, 8 (2017), p. 1504.  Number: 1 Publisher: Nature Publishing Group.
  * [38] B. Perthame, Transport Equations in Biology, Frontiers in Mathematics, Birkhäuser Basel, 2006.
  * [39] F. Peruani, J. Starruß, V. Jakovljevic, L. Søgaard-Andersen, A. Deutsch, and M. Bär, Collective motion and nonequilibrium cluster formation in colonies of gliding bacteria, Physical review letters, 108 (2012), p. 098102.
  * [40] A. Pikovsky, M. Rosenblum, and J. Kurths, Synchronization: A Universal Concept in Nonlinear Sciences, Cambridge University Press, Cambridge, 1st edition ed., May 2003.
  * [41] G. M. Ramírez-Ávila, J. Kurths, S. Depickère, and J.-L. Deneubourg, Modeling Fireflies Synchronization, in A Mathematical Modeling Approach from Nonlinear Dynamics to Complex Systems, E. E. N. Macau, ed., Nonlinear Systems and Complexity, Springer International Publishing, Cham, 2019, pp. 131–156.
  * [42] F. A. Rodrigues, T. K. D. Peron, P. Ji, and J. Kurths, The Kuramoto model in complex networks, Physics Reports, 610 (2016), pp. 1–98.
  * [43] G. K. Sar, S. N. Chowdhury, M. Perc, and D. Ghosh, Swarmalators under competitive time-varying phase interactions, New Journal of Physics, 24 (2022), p. 043004.
  * [44] L. J. Shimkets and D. Kaiser, Induction of coordinated movement of myxococcus xanthus cells, Journal of Bacteriology, 152 (1982), pp. 451–461.
  * [45] W. Singer, Neuronal synchrony: a versatile code for the definition of relations?, Neuron, 24 (1999), pp. 49–65, 111–125.
  * [46] J. Sokol, How Do Fireflies Flash in Sync? Studies Suggest a New Answer., Quanta Magazine, (2022).
  * [47] M. S. Steinberg, Does differential adhesion govern self-assembly processes in histogenesis? equilibrium configurations and the emergence of a hierarchy among populations of embryonic cells, Journal of Experimental Zoology, 173 (1970), pp. 395–433.
  * [48] S. H. Strogatz, From Kuramoto to Crawford: exploring the onset of synchronization in populations of coupled oscillators, Physica D: Nonlinear Phenomena, 143 (2000), pp. 1–20.
  * [49] S. H. Strogatz and R. E. Mirollo, Stability of incoherence in a population of coupled oscillators, Journal of Statistical Physics, 63 (1991), pp. 613–635.
  * [50] D. Tanaka, General chemotactic model of oscillators, Physical review letters, 99 (2007), p. 134103.
  * [51] S. Thutupalli, M. Sun, F. Bunyak, K. Palaniappan, and J. W. Shaevitz, Directional reversals enable myxococcus xanthus cells to produce collective one-dimensional streams during fruiting-body formation, Journal of The Royal Society Interface, 12 (2015), p. 20150049.
  * [52] P. Uhlhaas, G. Pipa, B. Lima, L. Melloni, S. Neuenschwander, D. Nikolić, and W. Singer, Neural synchrony in cortical networks: history, concept and current status, Frontiers in Integrative Neuroscience, 3 (2009).
  * [53] C. Van Oss, A. V. Panfilov, P. Hogeweg, F. Siegert, and C. J. Weijer, Spatial pattern formation during aggregation of the slime mould dictyostelium discoideum, Journal ofTheoretical bBology, 181 (1996), pp. 203–213.
  * [54] O. F. Venzin and A. C. Oates, What are you synching about? emerging complexity of notch signaling in the segmentation clock, Developmental biology, 460 (2020), pp. 40–54.
  * [55] A. T. Winfree, The timing of biological clocks, vol. 19, Macmillan, 1987\.
  * [56] D. Witthaut and M. Timme, Kuramoto dynamics in hamiltonian systems, Physical Review E, 90 (2014), p. 032917.
  * [57] Y. Wu, A. D. Kaiser, Y. Jiang, and M. S. Alber, Periodic reversal of direction allows myxobacteria to swarm, Proceedings of the National Academy of Sciences, 106 (2009), pp. 1222–1227.
  * [58] D. R. Zusman, A. E. Scott, Z. Yang, and J. R. Kirby, Chemosensory pathways, motility and development in myxococcus xanthus, Nature Reviews Microbiology, 5 (2007), pp. 862–872.

##  Appendix A Linear Stability Analysis

We summarize the computation for the linear stability analysis here whose
results are stated in Proposition 4.1. We only cover the case of $d=2$ spatial
dimensions. The simpler case of one spatial dimension ($d=1$) is a
straightforward modification and we leave out the details that yield the
corresponding functional forms of $F_{1}$ and $G_{1}$.

Consider solutions of the form

$R(t,x,\theta)=\bar{R}+c(t)e^{i(q\theta+k\cdot x)},$ (A.1)

to (2.1). Here $\bar{R}$ is a constant solution, i.e. an incoherent solution
where all clock phases have equal spatially constant density. The vector
$k\in{\mathbb{R}}^{2}$ is the spatial wavevector. Periodic boundary conditions
in the clock phase $\theta$ force that $q\in\mathbb{Z}$. Below, we compute the
terms $\nabla\cdot(v_{1}(t,\theta,x;R)R)$ and
$\frac{\partial}{\partial\theta}(v_{2}(t,\theta,x;R)R)$ explicitly:

Term $\mathbf{\nabla\cdot(v_{1}(t,\theta,x;R)R)}$

One computes that

$\displaystyle v_{1}(t,x,\theta;R)$
$\displaystyle=\int_{0}^{2\pi}\int_{D_{\rho}(0)}(J_{0}+J\cos(\theta-{\tilde{\theta}}))\,\left(\bar{R}+c(t)e^{i(q{\tilde{\theta}}+k\cdot(x+\tilde{x}))}\right)\,\frac{\tilde{x}}{|\tilde{x}|}d^{2}\tilde{x}\,d{\tilde{\theta}}$
(A.2)
$\displaystyle=e^{ikx}c(t)\int_{D_{\rho}(0)}\int_{0}^{2\pi}(J_{0}+J\cos(\theta-{\tilde{\theta}}))\,e^{iq{\tilde{\theta}}}\,d{\tilde{\theta}}\,e^{ik\cdot\tilde{x}}\,\frac{\tilde{x}}{|\tilde{x}|}d^{2}\tilde{x}$
(A.3) $\displaystyle=\begin{cases}\pi c(t)Je^{i(q\theta+k\cdot
x)}\int_{D_{\rho}(0)}e^{ik\cdot\tilde{x}}\frac{\tilde{x}}{|\tilde{x}|}d^{2}\tilde{x}&\text{
for }q=\pm 1\\\ 2\pi c(t)J_{0}e^{ik\cdot
x}\int_{D_{\rho}(0)}e^{ik\cdot\tilde{x}}\frac{\tilde{x}}{|\tilde{x}|}d^{2}\tilde{x}&\text{
for }q=0\\\ 0&\text{ otherwise}\end{cases}$ (A.4)

Here the identity
$\int_{0}^{2\pi}e^{iq{\tilde{\theta}}}\cos(\theta-{\tilde{\theta}})d{\tilde{\theta}}=\begin{cases}\pi
e^{iq\theta}&q=\pm 1\\\ 0&\text{otherwise}\end{cases}$ was used. To evaluate
the integral in the above formula for $v_{1}(t,\theta,x;R)$ in the case $k\neq
0$, rewrite the integral in polar coordinates
$\tilde{x}=r\cos\phi\,\frac{k}{|k|}+r\sin\phi\,k^{\perp}$, where $k^{\perp}$
is a unit vector perpendicular to $k$. We then obtain

$\displaystyle\int_{D_{\rho}(0)}e^{ik\cdot\tilde{x}}\frac{\tilde{x}}{|\tilde{x}|}d^{2}\tilde{x}$
$\displaystyle=\int_{0}^{2\pi}\int_{0}^{\rho}e^{i|k|r\cos\phi}r(\cos\phi\,\frac{k}{|k|}+\sin\phi\,k^{\perp})dr\,d\phi$
$\displaystyle=i\,\rho^{2}f(|k|,\rho)\frac{k}{|k|},$

where
$f(|k|,\rho)=\frac{1}{\rho^{2}}\int_{0}^{\rho}\int_{0}^{2\pi}\sin(|k|r\cos\phi)r\cos\phi\,d\phi\,dr=\int_{0}^{1}\int_{0}^{2\pi}\sin(\rho|k|r\cos\phi)r\cos\phi\,d\phi\,dr$
is a real-valued function. (The real part of the integral evaluate to zero, as
can be seen by symmetry considerations and periodicity of sine and cosine.)
Note two things: First, $f(|k|=0,\rho)=0$, so the functional form of $f$ also
applies to the case $k=0$ in (A.4). Second, $f(|k|,\rho)$ actually only
depends on the (dimensionless) product $k^{*}=\rho|k|$. This gives

$\displaystyle\nabla\cdot(v_{1}(t,\theta,x;R)R)=\begin{cases}-\pi
J\bar{R}|k|\rho^{2}f(|k|,\rho)e^{i(q\theta+k\cdot
x)}c(t)+\mathcal{O}(c(t)^{2})&\text{ for }q=\pm 1\\\ -2\pi
J_{0}\bar{R}|k|\rho^{2}f(|k|,\rho)e^{ik\cdot
x}c(t)+\mathcal{O}(c(t)^{2})&\text{ for }q=0\\\ 0&\text{
otherwise}\end{cases}$

Term $\frac{\partial}{\partial\theta}(v_{2}(t,\theta,x;R)R)$

One computes that

$\displaystyle v_{2}(t,x,\theta;R)$
$\displaystyle=K\int_{0}^{2\pi}\int_{D_{\rho}(0)}\sin({\tilde{\theta}}-\theta))\,\left(\bar{R}+c(t)e^{i(q{\tilde{\theta}}+k\cdot(x+\tilde{x}))}\right)\,d^{2}\tilde{x}\,d{\tilde{\theta}}$
$\displaystyle=K\,c(t)e^{ik\cdot
x}\int_{D_{\rho}(0)}\int_{0}^{2\pi}\sin({\tilde{\theta}}-\theta))e^{iq{\tilde{\theta}}}d{\tilde{\theta}}\,e^{ik\cdot\tilde{x}}\,d^{2}\tilde{x}$
$\displaystyle=\begin{cases}i\pi\,q\,Kc(t)e^{i(q\theta+k\cdot
x)}\int_{D_{\rho}(0)}e^{ik\cdot\tilde{x}}d^{2}\tilde{x}&\text{ for }q=\pm 1\\\
0&\text{ otherwise}\end{cases}$

Here the identity
$\int_{0}^{2\pi}e^{iq{\tilde{\theta}}}\sin({\tilde{\theta}}-\theta)d{\tilde{\theta}}=\begin{cases}i\,q\pi
e^{iq\theta}&q=\pm 1\\\ 0&\text{otherwise}\end{cases}$ was used. To evaluate
the integral in the formula for $v_{2}(t,x,\theta;R)$, switch again to polar
coordinates $x=r\cos\phi\,\frac{k}{|k|},y+r\sin\phi\,k^{\perp}$. This gives

$\displaystyle\int_{D_{\rho}(0)}e^{ik\cdot\tilde{x}}d^{2}\tilde{x}$
$\displaystyle=\int_{0}^{2\pi}\int_{0}^{\rho}e^{i|k|r\cos\phi}rdr\,d\phi=\rho^{2}\,g(|k|,\rho).$

Here
$g(|k|,\rho)=\frac{1}{\rho^{2}}\int_{0}^{2\pi}\int_{0}^{\rho}\cos({|k|r\cos\phi})rdr\,d\phi=\int_{0}^{2\pi}\int_{0}^{1}\cos({\rho|k|r\cos\phi})rdr\,d\phi$
is a real-valued function. Note that $g(|k|,\rho)$ actually only depends on
the (dimensionless) product $k^{*}=\rho|k|$. This gives

$\displaystyle\frac{\partial}{\partial\theta}(v_{2}(t,\theta,x;R)R)=\begin{cases}-\pi
K\bar{R}\rho^{2}g(|k|,\rho)e^{i(q\theta+k\cdot
x)}c(t)+\mathcal{O}(c(t)^{2})&\text{ for }q=\pm 1\\\ 0&\text{
otherwise}\end{cases}$

Complete linearized equation By inserting (A.1) into (2.1), using the above
computations and neglecting terms of second order and higher in $c(t)$, we
obtain the result of Proposition 4.1 in $d=2$ spatial dimensions:

######  Proposition A.1.

For $q=0$, the linearization is

$\frac{dc}{dt}=(-|k|^{2}D_{x}+2J_{0}\bar{R}\rho\,F_{2}(|\rho k|))\,c$

For $|q|\neq 0,1$, the linearization is

$\frac{dc}{dt}=(-|k|^{2}D_{x}-q^{2}D_{\theta})\,c$

For $|q|=1$, the linearization is

$\frac{dc}{dt}=(-|k|^{2}D_{x}-D_{\theta}+J_{0}\bar{R}\rho\,F_{2}(|\rho
k|)+\bar{R}K\rho^{2}\,G_{2}(|\rho k|))\,c.$

Explicitly, the two functions $F_{2}(k^{*})$ and $G_{2}(k^{*})$ for real
$k^{*}\geq 0$ are given by the following formulas, where we used Mathematica
to get explicit expressions for the integrals:

$\displaystyle F_{2}(k^{*})$ $\displaystyle=\pi k^{*}\,f(k^{*},1)=\pi
k^{*}\int_{0}^{1}\int_{0}^{2\pi}\sin(k^{*}r\cos\phi)r\cos\phi\,d\phi\,dr$
$\displaystyle=\pi^{3}({\mathcal{J}}_{1}(k^{*})\,H_{0}(k^{*})-\mathcal{J}_{0}(k^{*})\,H_{1}(k^{*}))$

and

$\displaystyle G_{2}(k^{*})$ $\displaystyle=\pi
g(k^{*},1)=\pi\int_{0}^{1}\int_{0}^{2\pi}\cos(k^{*}r\cos\phi)rd\phi dr$
$\displaystyle=2\pi^{2}\frac{1}{k^{*}}{\mathcal{J}}_{1}(k^{*}).$

Here $\mathcal{J}_{0},\mathcal{J}_{1}$ denotes Bessel functions of th efirst
and second kind, respectively, and $H_{0},H_{1}$ denote Struve functions [14].
Plots of $F_{2}$ and $G_{2}$ are shown in Fig 1.

##  Appendix B Constraints on $J_{0}^{*},J^{*},K^{*}$

As pointed out in the main tetxt, the PDE model (2.1) hinges on the assumption
that patterns have sufficiently large characteristic scales that a description
using cell density is appropriate. More specifically, we are only interested
in values of $J_{0},J$ and $K$ which lead to patterns with a characteristic
length scale comparable to the cell interaction radius $\rho$ or larger. In
this Appendix, we identify a region in parameter space which satisfies this
constraints.

For this, consider the linearization discussed in section 4.1. For a
linearization pattern (mode) with spatial wave vector $k$, a characteristic
length scale is the distance between peaks and valleys $L=\pi/|k|$. The
condition that $L$ does not exceed the interaction radius $\rho$ is thus
equivalent to $k^{*}=\rho k<\pi$. Proposition 4.1 shows that the growth rates
of different modes are governed by the functions $\sigma_{1}(k^{*})$ and
$\sigma_{3}(k^{*})$ given in $d$ dimensions by

$\displaystyle\sigma_{1}(k^{*})$
$\displaystyle=-(k^{*})^{2}+2J_{0}^{*}\,F_{d}(k^{*})$
$\displaystyle\sigma_{3}(k^{*})$
$\displaystyle=-(k^{*})^{2}-D_{\theta}^{*}+J^{*}F_{d}(k^{*})+K^{*}G_{d}(k^{*})$

where

$J^{*}=\frac{\bar{R}\rho^{d+1}}{D_{x}}J,\quad
J_{0}^{*}=\frac{\bar{R}\rho^{d+1}}{D_{x}}J_{0},\quad
K^{*}=\frac{\bar{R}\rho^{d+2}}{D_{x}}K$

The dominant modes are those that correspond to global maxima of those
functions.

We now use the following condition to constrain the parameters $J_{0},J$ and
$K$: We require

$\max_{k^{*}\geq\pi}(\sigma_{i}(k^{*}))<\max(0,\max_{k^{*}\geq
0}\sigma_{i}(k^{*}))\quad\text{for }i=1,3$ (B.1)

This somewhat cumbersome formula simply means that we require that the
dominant modes $k^{*}_{\rm{dom}}$for both $\sigma_{1}$ and $\sigma_{3}$ to
either be less than $\pi$ or have negative temporal growth rates. In other
words, a pattern with length scale less than the interaction length $\rho$
cannot be the fastest growing pattern.

Using the explicit forms of $F_{d}(k^{*})$ and $G_{d}(k^{*})$ given in section
4.1, one can derive with a combination of simple calculus arguments and
numerical investigation the following results. (In dimension 2, these are not
sharp.)

######  Proposition B.1.

  1. (1)

In $d=1$ spatial dimensions, the condition (B.1) is satified in the following
parameter range:

$-\pi<K^{*}<\infty,\quad-\infty<J_{0}^{*},J^{*}<\infty.$

  2. (2)

In $d=2$ spatial dimensions, the condition (B.1) is satifies in the following
parameter range:

$-1.4<K^{*}<\infty,\quad-7.3<J^{*}_{0}<3.1,\quad-3.65<J^{*}<1.9$

(For (1), note that $K^{*}>-\pi$ guarantees $\sigma_{3}(\pi)<0,$ which implies
that the condition (B.1) is satisfied. The result in (2) was obtained
numerically using Mathematica.)

##  Appendix C Solutions of the form $R(t,x,\theta)=R(t,\theta)$

We consider the special case that the cell density $R(t,x,\theta)$ in (2.1) is
constant with respect to the spatial position $x$:

$R(t,x,\theta)=R(t,\theta).$

Strictly speaking, this is of course only applicable in the case that the
initial conditions are perfectly spatially homogeneous. But the analysis is
more generally relevant in the case that all spatially dependent modes have
negative growth rates, see Propositions 4.2 and 4.3. The differential equation
for $R(t,\theta)$ takes the form

$\frac{\partial R}{\partial
t}=D_{\theta}\frac{\partial^{2}R}{\partial\theta^{2}}-\,\frac{\partial}{\partial\theta}(v_{2}(t,\theta;R)R).$
(C.1)

with

$v_{2}(t,\theta;R)=K\,|D_{\rho}|\int_{0}^{2\pi}\sin({\tilde{\theta}}-\theta)\,R(t,{\tilde{\theta}})\,d{\tilde{\theta}},$

where $|D_{\rho}|$ is the $n-$dimensional area of the ball of radius $\rho$
(i.e. the length of the interaction neighborhood in one spatial dimension or
its area in two dimensions.)

Consider the Fourier expansion of $R(t,\theta)$:

$R(t,\theta)=\sum_{j=-\infty}^{\infty}A_{j}(t)e^{ij\theta}.$ (C.2)

We consider real valued $R(t,\theta)$, so $A_{-j}(t)=\bar{A}_{j}(t)$. Note
that $A_{j}(t=0)$ are the Fourier coefficients of the initial condition
$R(t=0,\theta)=R_{0}(\theta)$.

We have the following proposition:

######  Proposition C.1.

The Fourier coefficients $A_{j}(t)$ satisfy the infinite system of
differential equations

$\frac{dA_{j}}{dt}=j|D_{\rho}|\pi
K(A_{1}(t)A_{j-1}(t)-\bar{A}_{1}(t)A_{j+1}(t))-j^{2}D_{\theta}A_{j}(t)\quad\quad{\text{
for }}j\in\mathbb{Z}$ (C.3)

In particular

$\frac{d}{dt}A_{0}(t)=0,$

so that
$A_{0}(t)=const=\frac{1}{2\pi}\int_{0}^{2\pi}R_{0}(\theta)d\theta=\bar{R}$

###### Proof.

Compute first $v_{2}(t,\theta;R)$ with the Fourier expansion C.2:

$\displaystyle v_{2}(t,\theta;R)$
$\displaystyle=\frac{K}{2i}\,|D_{\rho}|\int_{0}^{2\pi}(e^{i({\tilde{\theta}}-\theta)})-e^{-i({\tilde{\theta}}-\theta))})\,R(t,{\tilde{\theta}})\,d{\tilde{\theta}}$
$\displaystyle=\frac{K}{2i}\,|D_{\rho}|\sum_{j=-\infty}^{\infty}A_{j}(t)\left(e^{-i\theta}\int_{0}^{2\pi}e^{i(j+1){\tilde{\theta}}}\,d{\tilde{\theta}}-e^{i\theta}\int_{0}^{2\pi}e^{i(j-1){\tilde{\theta}}}d{\tilde{\theta}}\right)$
$\displaystyle=i\pi
K|D_{\rho}|\left(-A_{-1}(t)e^{-i\theta}+A_{1}(t)e^{i\theta}\right).$

Using this expression in the differential equation (C.1) yields

$\displaystyle\frac{\partial R}{\partial t}$
$\displaystyle=\sum_{j=-\infty}^{\infty}\left(K|D_{\rho}|\pi\left(-(j-1)A_{-1}(t)A_{j}(t)e^{i(j-1)\theta}+(j+1)A_{1}(t)A_{j}(t)e^{i(j+1)\theta}\right)-j^{2}D_{\theta}A_{j}(t)e^{ij\theta}\right)$

Comparing Fourier coefficients of the left and right hand sides gives the
system of equations for $A_{j}(t)$ in (C.3) ∎

######  Theorem C.2.

For $K<0$ and $D_{\theta}>0$, the incoherent state of the space-independent
problem (C.1) is globally stable, i.e. $R(t,\theta)$ converges to the constant
$\bar{R}=\frac{1}{2\pi}\int_{0}^{2\pi}R(t=0,\theta)d\theta$ as $t\to\infty$.

###### Proof.

Consider the function

$u(t)=\sum_{k=1}^{\infty}\frac{1}{k}|A_{k}|^{2}$

which converges absolutely by Parseval’s Theorem. One computes that

$\frac{1}{k}\frac{d}{dt}|A_{k}|^{2}=K|D_{\rho}|\pi\left(A_{1}A_{k-1}\bar{A}_{k}-\bar{A}_{1}A_{k+1}\bar{A}_{k}+\bar{A}_{1}\bar{A}_{k-1}A_{k}-A_{1}\bar{A}_{k+1}A_{k}\right)-2D_{\theta}k|A_{k}|^{2}.$

Using this, one finds that

$\displaystyle\frac{du}{dt}$
$\displaystyle=\sum_{k=1}^{N}\frac{1}{k}\frac{d}{dt}|A_{k}|^{2}$
$\displaystyle=K|D_{\rho}|\pi\left(2A_{0}|A_{1}|^{2}-\bar{A}_{1}A_{N+1}\bar{A}_{N}-A_{1}\bar{A}_{N+1}A_{N}\right)-2D_{\theta}k|A_{k}|^{2}$

due to a telescoping sum. Note that the remainder goes to zero due to the fact
that $\sum_{k=1}^{\infty}|A_{k}|^{2}$ converges. Therefore

$\frac{du}{dt}=2K|D_{\rho}|\pi
A_{0}|A_{1}|^{2}-2D_{\theta}\sum_{k=1}^{\infty}k|A_{k}|^{2}$

If $K<0$, we can estimate

$\frac{du}{dt}\leq-2D_{\theta}u(t).$

Provided $D_{\theta}>0$, Gronwall’s lemma implies $u(t)\to 0$ as $t\to\infty$.
This means that $|A_{k}(t)|^{2}\to 0$ as $t\to\infty$ for all $k=1,2,\ldots$.
Thus

$R(t,\theta)\to A_{0}=\frac{1}{2\pi}\int_{0}^{2\pi}R(t=0,\theta)d\theta.$

∎

We strongly expect the result to still hold for the case of zero clock
diffusion ($D_{\theta}=0$). In this case

$\frac{du}{dt}=2K|D_{\rho}|\pi A_{0}|A_{1}|^{2}$

This means that $u(t)=\sum_{k=1}^{\infty}\frac{1}{k}|A_{k}|^{2}$ is strictly
decreasing, as $A_{0}>0$ since we assume nonnegative initial conditions. We
conjecture that in fact $u(t)\to 0$ as $t\to\infty$. This would then imply the
statement of Theorem C.2.