of 46
1
Wildebeest Herds on Rolling Hills:
Flocking on Arbitrary Curved Surfaces
Christina L. Hueschen
1
,
, Alexander R. Dunn
1
, Rob Phillips
2
,
3
,
1 Department of Chemical Engineering, Stanford University, Palo Alto, California, 94305
2 Department of Physics, California Institute of Technology, Pasadena, California, 91125
3 Division of Biology and Biological Engineering, California Institute of Technology,
Pasadena, California, 91125
Correspondence: chueschen@gmail.com, phillips@pboc.caltech.edu
Contents
Abstract
2
1 Introduction
2
2 Flocking Theory for Arbitrary Curved Surfaces
3
2.1 A Minimal Toner-Tu Theory in the Plane . . . . . . . . . . . . . . . . . . . . . . . . . . .
3
2.2 Formulating the Theory for Arbitrary Curved Surfaces . . . . . . . . . . . . . . . . . . . .
4
3 Formulating the Surface Toner-Tu Equations for Numerical Implementation
9
4 Parameters and Dimensionless Ratios for the Theory
12
4.1 Parameter Choices for Wildebeest Herds . . . . . . . . . . . . . . . . . . . . . . . . . . . .
12
4.2 Dimensionless Representation of the Theory . . . . . . . . . . . . . . . . . . . . . . . . . .
14
5 Solving Toner-Tu on Highly Symmetric Geometries
15
5.1 Solving the Toner-Tu Equations in the Plane . . . . . . . . . . . . . . . . . . . . . . . . .
15
5.2 Formulating the Toner-Tu Equations for Parameterized Surfaces . . . . . . . . . . . . . .
22
5.3 Solving the Toner-Tu Equations on a Cylindrical Surface . . . . . . . . . . . . . . . . . . .
23
5.4 Solving the Toner-Tu Equations on a Spherical Surface . . . . . . . . . . . . . . . . . . . .
25
6 Wildebeest on Spheres and Hills: Complex Flocking Dynamics
26
6.1 Toner-Tu Dynamics on the Sphere . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
27
6.2 Herding over a Gaussian Hill . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
29
6.3 Herding on an Undulating Island Landscape . . . . . . . . . . . . . . . . . . . . . . . . . .
29
7 Conclusion
30
8 Appendix
32
8.1 Calculus on Surfaces Via Projection Operators . . . . . . . . . . . . . . . . . . . . . . . .
32
8.2 Implementation in COMSOL Multiphysics . . . . . . . . . . . . . . . . . . . . . . . . . . .
35
8.3 Analytical Solutions for Toner-Tu on the Cylinder . . . . . . . . . . . . . . . . . . . . . .
37
8.4 Analytical Solution for Toner-Tu on the Sphere . . . . . . . . . . . . . . . . . . . . . . . .
39
Video Legends
43
Acknowledgements
43
References
44
.
CC-BY-NC-ND 4.0 International license
available under a
(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprint
this version posted June 25, 2022.
;
https://doi.org/10.1101/2022.06.22.497052
doi:
bioRxiv preprint
2
Abstract
The collective behavior of active agents, whether herds of wildebeest or microscopic actin filaments
propelled by molecular motors, is an exciting frontier in biological and soft matter physics. Almost three
decades ago, Toner and Tu developed a hydrodynamic theory of the collective action of flocks, or herds,
that helped launch the modern field of active matter. One challenge faced when applying continuum active
matter theories to living phenomena is the complex geometric structure of biological environments. Both
macroscopic and microscopic herds move on asymmetric curved surfaces, like undulating grass plains
or the surface layers of cells or embryos, which can render problems analytically intractable. In this
work, we present a formulation of the Toner-Tu flocking theory that uses the finite element method to
solve the governing equations on arbitrary curved surfaces. First, we test the developed formalism and
its numerical implementation in channel flow with scattering obstacles and on cylindrical and spherical
surfaces, comparing our results to analytical solutions. We then progress to surfaces with arbitrary
curvature, moving beyond previously accessible problems to explore herding behavior on a variety of
landscapes. Our approach allows the investigation of transients and dynamic solutions not revealed by
analytic methods. It also enables versatile incorporation of new geometries and boundary conditions and
efficient sweeps of parameter space. Looking forward, the work presented here lays the groundwork for a
dialogue between Toner-Tu theory and data on collective motion in biologically-relevant geometries, from
drone footage of migrating animal herds to movies of microscopic cytoskeletal flows within cells.
1. Introduction
The beautiful collective motions of flocking or herding animals have mesmerized the human mind
for millennia [1] and inspired the modern field of study of active matter [2, 3, 4]. In recent decades,
active matter theories have been used to describe collective motion in living systems across nearly a
billion-fold difference in scales [5], from starling flocks wheeling over Rome at 10 m/s [6], to 50
μ
m/s
flows of cytoplasm in cm-size internodal cells of the algae
Chara
[7], to 0.1
μ
m/s flows of actin filaments
and myosin molecules in 50
μ
m developing
C. elegans
worm embryos [8]. In these macroscopic and
microscopic contexts, beautiful self-organization phenomena often take place on surfaces. The collective
motion of flocking sheep or kilometer-wide migrating wildebeest herds are influenced by the hills, valleys,
or canyons of the landscape on which they travel. Similarly, the flows of actin and myosin molecules
mentioned above occur in thin layers at the surface of the embryo or cell, and are thus constrained by
its surface topology and shape. While modern techniques such as dynamical drone imaging [9] and high-
resolution fluorescence microscopy [10, 11] enable experimental measurement of these animal, cell, or
molecule kinematics, a dialogue between measurement and theory requires predictions that incorporate
the complex shapes of the real-world surfaces on which they move.
In this work, we present a general curved-surface formulation and numerical implementation of a
minimalistic continuum theory of flocking or herding active matter, with the hope that it will prove
useful to others interested in exploring continuum theory predictions on complex geometries. For active
matter systems such as bacterial swarms [12], active colloidal fluids [13], self-propelled rods [14], or purified
cytoskeletal networks [15], studying the contribution of engineered confinement geometries to emergent
patterns has proven a fruitful path [16, 17, 18, 19, 20]. Understanding the contribution of
biological
geometries to pattern formation is an exciting direction of growth [21, 22, 23, 24, 25, 26]. Indeed, our
own interest in solving active matter theories on arbitrary curved surfaces was initially inspired by a need
to predict emergent actin polarity patterns at the surface of single-celled
Toxoplasma gondii
parasites,
whose gliding motility is driven and directed by this surface actin layer [27]. These cells have a beautiful
but complex shape which lacks the symmetries that license traditional analytic approaches.
We focus on a classic continuum active matter model originally developed by John Toner and Yuhai
Tu [3, 4, 28], inspired by the work of T ́amas Vicsek
et al.
[2], to describe the collective behavior of flocking
or herding animals. The Toner-Tu theory helped launch the modern field of active matter [29] and can
.
CC-BY-NC-ND 4.0 International license
available under a
(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprint
this version posted June 25, 2022.
;
https://doi.org/10.1101/2022.06.22.497052
doi:
bioRxiv preprint
3
describe collections of dry, polar, self-propelled agents at any length scale. We develop a general curved
surface framework for the Toner-Tu theory and implement it in the finite element setting, enabling its
convenient use on complex surfaces. While flocking theory - and in particular the herding of wildebeest
- serves as our example in this study, the general curved surface formulation and finite element method
(FEM) approach presented here may prove useful for any continuum active matter theory. There is
a rich and beautiful literature associated with formulating hydrodynamic equations on arbitrary curved
surfaces [30, 31, 32, 33, 34, 35, 36]. It has been said of these formulations, “The complexity of the equations
may explain why they are so often written but never solved for arbitrary surfaces” [34]. Because of our
interest in collective motions on biological surfaces, we could not afford to avoid solving these equations
on such surfaces. We use the finite element method, which can be flexibly adjusted to an arbitrary choice
of geometry and permits an exploration of active, self-organized solutions predicted by the Toner-Tu
theory. With our general surface formulation and finite element implementation in hand, we test our
framework and approach on cylindrical and spherical geometries for which analytic solutions exist as well.
Satisfied by the agreement between our numerical and analytic results, we explore flocking phenomena
on a broad collection of surfaces.
The remainder of the paper is organized as follows. In Section 2, we show how the flocking theory of
Toner and Tu can be recast in a general surface form to describe flocking motions on arbitrary curved
surfaces. In Section 3, we construct a finite element surface formulation that respects the low symmetry
of realistic curved geometries while permitting numerical analysis of the dynamics. In Section 4, we
turn to the great wildebeest herds for inspiration and perform simple parameter estimates for Toner-Tu
wildebeest herding. We then explore the dimensionless ratios that appear when recasting the Toner-
Tu equations in dimensionless form. In Section 5, with the full surface formulation and its numerical
implementation in hand, we explore scaling relationships and changes of dynamical state that arise for
herds in channels with scattering obstacles, and we use cylindrical and spherical surfaces to compare
our finite element method results to corresponding analytic solutions. Our approach also allows the
exploration of transients and dynamic solutions not revealed by analytic methods. Finally, in Section 6,
we playfully use the curved-space formalism and its finite element implementation for case studies of
wildebeest herding on landscapes of rolling hills, and we note that this work serves as a theoretical
foundation to solve flocking theories on real-world curved surfaces of biological interest.
2. Flocking Theory for Arbitrary Curved Surfaces
2.1. A Minimal Toner-Tu Theory in the Plane
The first step in the development of field theories of active matter involves selecting the relevant
field variables. Following the classic work of Toner and Tu [3, 4] and inspired by observations of herds
of land animals like wildebeest, we consider a two-dimensional density field
ρ
(
r
,t
) and a corresponding
two-dimensional velocity field
v
(
r
,t
), which captures both the orientation and speed of a polar agent.
Throughout this study, wildebeest herds will serve as inspiration and example, although in a sense the
term “wildebeest” is our shorthand for a “self-propelled polar agent” at any length scale. We also note
that the Toner-Tu model used here describes
dry
systems, in which momentum is not conserved and
hydrodynamic coupling between the agents is negligible relative to frictional drag. With variables
ρ
and
v
in hand, our next step is to write the partial differential equations that describe the spatiotemporal
evolution of those field variables. The first governing equation is the continuity equation given by
∂ρ
∂t
+
(
ρv
i
)
∂x
i
= 0
,
(1)
which allows modulations in the density field but enforces mass conservation, forbidding wildebeest birth
or death in the midst of herding phenomena. Note that we are using the Einstein summation convention
that tells us to sum over all repeated indices; for example,
a
·
b
=
3
i
=1
a
i
b
i
=
a
i
b
i
. The second governing
.
CC-BY-NC-ND 4.0 International license
available under a
(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprint
this version posted June 25, 2022.
;
https://doi.org/10.1101/2022.06.22.497052
doi:
bioRxiv preprint
4
equation is a minimal representation of the dynamics of the
v
field offered by Toner and Tu [3] and is
given by
∂v
i
∂t
= [
α
(
ρ
ρ
c
)
βv
j
v
j
]
v
i
σ
∂ρ
∂x
i
+
D
2
v
i
λv
j
∂v
i
∂x
j
,
(2)
where
ρ
c
is the critical density above which the herd moves coherently, the ratio of
α
(
ρ
ρ
c
) and
β
sets the wildebeest mean-field speed
v
pref
=
α
(
ρ
ρ
c
)
, the term
σ ∂ρ/∂x
i
provides an effective
pressure,
D
tunes wildebeest alignment and speed matching with neighbors, and
λ
tunes velocity self-
advection. The intuition for each of these terms is sketched for a one-dimensional wildebeest herd in
Figure 1. We can conceptualize each term as an “update rule” that computes incremental changes in the
velocity vector at each point in space and at every step in time. The first term on the right hand side,
the preferred speed term, pushes the velocity magnitude toward the characteristic speed of a wildebeest,
v
pref
. The second is based on an equation of state that relates density and pressure. This pressure term
punishes gradients in density, adjusting velocity to flatten the density field. Given mass conservation in
a finite system, the pressure term alone would push the density field to a uniform steady-state
ρ
=
ρ
0
.
The third, the neighbor coupling term, provides a smoothing or diffusion of velocity (orientation and
speed) that reflects coordination between nearby wildebeest. The final term has analogy to the gradient
component of the material time derivative in the Navier-Stokes equations. In essence, the velocity field
advects itself; wildebeest move along in a direction dictated by their orientation and velocity, and they
bring that orientation and velocity with them. In the case of pure velocity self-advection,
λ
= 1, but this
is not necessarily the case for active flocks [37, 38]. In the case of wildebeest herds, we can conceptualize
this effect by considering that
λ
= 1 +
ξ
, where
ξ
reflects a behavioral response to gradients in velocity;
for example, wildebeest may resist running quickly into a steep gradient of decreasing velocity, and may
slow down.
2.2. Formulating the Theory for Arbitrary Curved Surfaces
Inspired by a desire to make contact with phenomena of the natural world, like wildebeest navigating
an undulating landscape, sheep flocks crossing hilly pastures, and the surfaces flows of flocking actin we
study in reference [27], we sought to solve the minimal Toner-Tu theory presented above on complex and
asymmetric surface geometries. We consider here flocking on non-deformable surfaces, but we refer the
interested reader to earlier work on surface hydrodynamics in the completely general case in which the
surface itself can evolve over time [32, 33, 34, 36]. Predicting flocking or herding behavior on arbitrary
surfaces requires us to reformulate the theory in a more general way that accounts for curvature. Instead
of basing our formulation on a parameterized surface and the intrinsic differential geometry tools that
this approach licenses, as done in beautiful earlier work [39], from the outset we have in mind arbitrary
surfaces that can be represented using finite element meshes and described by a field of local normal
vectors,
n
.
Our choice and use of this extrinsic differential geometry and finite element method approach was
aided by the work of many, including refs. [33, 34, 35, 40, 41, 42, 43, 44, 45] and chapter 3 of the
supplemental material for [46]. The finite element setting permits us the convenience of carrying out the
mathematics in the full three-dimensional setting of
R
3
, while using our knowledge of the normal vectors
everywhere on the surface of interest to project our governing equations onto the surface. While velocity
is described by the 3-dimensional vector
v
, both
v
and the scalar
ρ
are defined only on the surface. At
every point on the surface, derivatives evaluated in the usual
R
3
way are projected onto the tangent plane
using the local normal. For an insightful description of this extrinsic differential geometry approach to
handling curved surfaces and its mathematical equivalence to the intrinsic differential geometry strategy,
we recommend chapters 22 and 23 of Needham [47]. Central to the extrinsic geometry approach is the
projection operator, defined as
projection onto tangent plane =
P
=
I
n
n
,
(3)
.
CC-BY-NC-ND 4.0 International license
available under a
(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprint
this version posted June 25, 2022.
;
https://doi.org/10.1101/2022.06.22.497052
doi:
bioRxiv preprint
5
-0.8
-0.4
0
0.4
0.8
-0.2
-0.1
0
0.1
0.2
x (m)
0
2
4
6
8
10
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0
2
4
6
8
–0.6
–0.4
–0.2
0
0
0.2
0.4
0.6
0.8
1
1.2
1.4
v
pref
v
m
s
(
(
v
m
s
(
(
v
m
s
(
(
m
s
2
(
(
(
α (
ρ – ρ
c
) –
β
|
v
|
2
)
v
(
α (
ρ – ρ
c
) –
β
|
v
|
2
)
v
m
s
2
(
(
σ
ρ
σ
ρ
m
s
2
(
(
D∇
2
v
2
v
m
s
2
(
(
λ
v
v
preferred speed
preferred speed
pressure
neighbor
coupling
pressure
02468
10
x (m)
02468
10
x (m)
02468
10
neighbor
coupling
advection
–∇
v
v
v
∂t
=
+
D∇
2
v
λ
v
v
advection
(A)
(B)
(C)
(D)
Figure 1.
Physics of the velocity field in the Toner-Tu theory. The different terms in the Toner-Tu analysis are
represented graphically for a one-dimensional herd. We can think of each term as providing an update to the current
velocity. In each case, an example velocity profile (blue) or density profile (inset in B) is shown alongside the corresponding
velocity update. (A) The preferred speed term increases the velocity of wildebeest that are moving too slow and decreases
the velocity of wildebeest that are moving too fast. (B) The pressure term punishes gradients in density, adjusting
velocity in order to flatten these gradients. Given mass conservation in a finite system, the pressure term alone would lead
to a steady state of uniform density
ρ
=
ρ
0
. (C) The neighbor coupling term captures the velocity adjustment made by a
wildebeest to better match its neighbors, smoothing out differences in velocity. A given wildebeest (middle black dot)
adjusts its velocity by an amount represented by the green arrow: an average of the difference between its velocity and its
two neighbors’ velocities. This averaging of two differences is mathematically analogous to taking a local second derivative
or local Laplacian,
2
v
. (D) Finally, the advection term ensures that the filament velocity field is swept along according
to its own velocities, just as fluid velocity is self-advected in the Navier-Stokes equations. Mathematically, the velocity
adjustment needed to ensure velocity advection is a function of the spatial gradient of velocity (
−∇
v
, how mismatched in
velocity nearby wildebeest are) and the velocity itself (
v
, how fast that mismatch is carried along by the herd).
.
CC-BY-NC-ND 4.0 International license
available under a
(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprint
this version posted June 25, 2022.
;
https://doi.org/10.1101/2022.06.22.497052
doi:
bioRxiv preprint
6
n
v
v
v
v
= (
I
n
×
n
)
v
=
v
– (
v
n
)
n
=
v
v
projection
operator,
P
Figure 2.
Illustration of the projection operator essential to the general surface formulation and finite element treatment
of the Toner-Tu equations on curved surfaces. The surface of interest is represented by a collection of nodes (vertices of
triangles) and a field of surface normal vectors,
n
(black arrows). An arbitrary vector
v
is projected onto the tangent
plane (shown in gold), permitting a decomposition of the form
v
=
v
+
v
.
where
I
is the identity matrix and
n
n
is the outer product of the surface normal vector as shown in
Figure 2. To make sense of this expression mathematically, we recall that the outer product
n
n
is
defined through its action on a vector
v
as
(
n
n
)
v
=
n
(
n
·
v
)
.
(4)
We can write
v
=
v
+
v
, where
v
is the component of
v
in the tangent plane of the surface and
v
is normal to the surface. The action of
I
n
n
on a vector
v
is given by
(
I
n
n
)
v
=
v
n
(
n
·
v
) =
v
v
=
v
,
(5)
as illustrated in Figure 2. We can write the projection operator in component form as
P
ij
=
δ
ij
n
i
n
j
,
(6)
recalling that the normal vector to the surface is given by
n
= (
n
1
,n
2
,n
3
), or in full component form as
P
=
1
n
2
1
n
1
n
2
n
1
n
3
n
1
n
2
1
n
2
2
n
2
n
3
n
1
n
3
n
2
n
3
1
n
2
3
.
(7)
Using the projection operator allows us to perform calculations in the ordinary three-dimensional space
within which the surface of interest is embedded, but then to pick off only the pieces of the resulting
vectors that live within the surface.
In the remainder of this section, we examine each term in the minimal Toner-Tu equations and
translate it into its projected form. In Appendix Section 8.1, we define the projected surface form of
calculus operators as an additional reference for the reader. First, to modify the preferred speed term to
its curved-surface implementation, we note that it is now
v
, the in-plane velocity, that has a privileged
magnitude. This magnitude is imposed through the condition
preferred speed term = [
α
(
ρ
ρ
c
)
βv
j
v
j
]
v
i
.
(8)
If
|
v
|
is either larger or smaller than the privileged value
v
pref
=
α
(
ρ
ρ
c
)
, this term will adjust
the velocity towards that steady-state magnitude.
.
CC-BY-NC-ND 4.0 International license
available under a
(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprint
this version posted June 25, 2022.
;
https://doi.org/10.1101/2022.06.22.497052
doi:
bioRxiv preprint
7
Next, we consider the pressure term, whose physical origin comes from a model of pressure in powers
of density of the form introduced by Toner and Tu [3],
P
(
ρ
) =
n
=1
σ
n
(
ρ
ρ
0
)
n
(9)
where
ρ
0
is the mean density. In the minimal Toner-Tu theory we adopt here, only the first order term
in that expansion is kept, with the notational simplification that
σ
1
=
σ
, resulting in
pressure term =
∂P
∂x
i
=
σ
∂ρ
∂x
i
.
(10)
We note that while
ρ
0
is not explicitly present in the gradient of the first order pressure term, for the
case considered here of a finite surface on which total density is conserved, this pressure term effectively
maintains a density range centered around the mean density
ρ
0
established by our choice of initial
condition. The curved-space version of the pressure term requires projecting the full 3D gradient of the
density onto the tangent plane, using the projection operator defined in eqn. 3. We follow Jankuhn
et
al.
[34] in introducing the notation
Γ
for the projected surface gradient operator, where the subscript
Γ indicates that the gradient is evaluated on the surface of interest. Mathematically, this amounts to
computing
Γ
ρ
= (
I
n
n
)
ρ,
(11)
where
is the ordinary, three-dimensional gradient in Cartesian coordinates. This can be rewritten in
component form as
(
Γ
ρ
)
i
= [(
I
n
n
)
ρ
]
i
= (
δ
ij
n
i
n
j
)
∂ρ
∂x
j
=
∂ρ
∂x
i
n
i
n
j
∂ρ
∂x
j
.
(12)
The next term in the Toner-Tu equations that we consider in its curved-space format is the advection
term in eqn. 2, namely,
advection term =
λv
j
∂v
i
∂x
j
.
(13)
In direct notation, the curved-space version of this term has the form
advection term =
λ
(
v
·∇
Γ
)
v
,
(14)
which involves the curved-space version of the velocity gradient tensor. In Appendix Section 8.1, we
describe how to compute this tensor, which stated simply is
Γ
v
=
P
(
v
)
P
,
(15)
where
P
is the projection operator defined in eqn. 3. In indicial notation, this leads to the result
(
Γ
v
)
ij
= [
P
(
v
)
P
]
ij
=
P
ik
∂v
k
∂x
l
P
lj
.
(16)
Invoking the definition of the projection operator from eqn. 6, this expression simplifies to the result
(
Γ
v
)
ij
= [
P
(
v
)
P
]
ij
=
(
∂v
i
∂x
j
n
l
n
j
∂v
i
∂x
l
)
n
i
n
k
(
∂v
k
∂x
j
n
l
n
j
∂v
k
∂x
l
)
.
(17)
Thus, the curved-space advection term can be written as
λ
[
P
(
v
)
P
]
ij
v
j
=
λ
[(
∂v
i
∂x
j
n
l
n
j
∂v
i
∂x
l
)
n
i
n
k
(
∂v
k
∂x
j
n
l
n
j
∂v
k
∂x
l
)]
v
j
.
(18)
.
CC-BY-NC-ND 4.0 International license
available under a
(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprint
this version posted June 25, 2022.
;
https://doi.org/10.1101/2022.06.22.497052
doi:
bioRxiv preprint
8
With these examples of tangent-plane calculus established, we now turn to the most tricky of the terms
in the Toner-Tu equations, the surface-projected version of the neighbor coupling term,
D
2
v
. First, we
recall that the Laplacian of a vector field in normal 3D Cartesian space is defined as
2
v
=
∇·
(
v
).
Further, we note that
v
is itself a tensor. As shown by Jankuhn
et al.
[34], the surface-projected version
of the Laplacian term (see their eqn. 3.16 for the surface Navier-Stokes equations for the tangential velocity
on a stationary surface), is therefore given by
2
v
︷︷
flat Toner-Tu
=
div
(
v
)
︷︷
in the plane
−→
P
div
Γ
(
G
(
v
))
︷︷
curved surface
.
(19)
To begin to unpack this expression, we note that
P
is the projection operator defined in eqn. 3 and that
the tensor
G
is the surface velocity gradient,
G
(
v
) =
Γ
v
,
(20)
already presented in eqn. 17 and repeated here in component form as
G
ij
=
(
∂v
i
∂x
j
n
l
n
j
∂v
i
∂x
l
)
n
i
n
k
(
∂v
k
∂x
j
n
l
n
j
∂v
k
∂x
l
)
.
(21)
We note that Jankuhn
et al.
[34] use a symmetrized version of the velocity gradient since they are deriving
the surface versions of the Navier-Stokes equations and are thus taking the divergence of a stress. For the
Toner-Tu case of interest here, the “psychological viscosity” that comes from neighboring wildebeest in
the herd comparing their velocities is equivalent to the divergence of the velocity gradient itself. We then
invoke the definition of the surface-projected divergence of a tensor
G
presented in Appendix Section 8.1,
giving rise to a vector of the form
div
Γ
(
G
)
l
=
∂G
lj
∂x
j
n
j
n
k
∂G
lj
∂x
k
.
(22)
We also introduce the shorthand notation for this divergence as
(
∂G
lj
∂x
j
)
Γ
=
∂G
lj
∂x
j
n
j
n
k
∂G
lj
∂x
k
(23)
to simplify some of the complex expressions to follow. Finally, we can write the
i
th
component of the
surface version of the Toner-Tu term
D
2
v
in indicial notation as
D
[
P
div
Γ
(
G
(
v
))]
i
=
DP
il
(
∂G
lj
∂x
j
n
j
n
k
∂G
lj
∂x
k
)
.
(24)
We now assemble the results of eqns. 8, 12, 18 and 24 to construct a complete curved-surface
formulation of the minimal Toner-Tu equations. Specifically, we have
∂v
i
∂t
= [
α
(
ρ
ρ
c
)
βv
j
v
j
]
v
i
σ
(
∂ρ
∂x
i
n
i
n
j
∂ρ
∂x
j
)
+
DP
il
(
∂G
lj
∂x
j
n
j
n
k
∂G
lj
∂x
k
)
λ
[
P
(
v
)
P
]
ij
v
j
,
(25)
which can be streamlined to the alternative form
∂v
i
∂t
= [
α
(
ρ
ρ
c
)
βv
j
v
j
]
v
i
σ
(
∂ρ
∂x
i
)
Γ
+
DP
il
(
∂G
lj
∂x
j
)
Γ
λv
j
[(
∂v
i
∂x
j
)
Γ
n
i
n
k
(
∂v
k
∂x
j
)
Γ
]
.
(26)
.
CC-BY-NC-ND 4.0 International license
available under a
(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprint
this version posted June 25, 2022.
;
https://doi.org/10.1101/2022.06.22.497052
doi:
bioRxiv preprint
9
Similarly, the curved-space formulation of the governing equation for density can be written as
∂ρ
∂t
=
(
ρv
i
)
∂x
i
+
n
i
n
j
(
ρv
i
)
∂x
j
=
(
∂ρv
i
∂x
i
)
Γ
.
(27)
We now have a complete formulation of our minimal Toner-Tu equations for the general surface context,
requiring only a description of the surface in the language of normal vectors. We next turn to the
implementation of this general surface formulation in a fashion consistent with finite element treatments
on arbitrary surfaces.
3. Formulating the Surface Toner-Tu Equations for Numerical Implementation
Numerically solving active matter equations on complex surfaces relevant to the living world presents
a practical challenge. It precludes the differential geometric formalism used to describe parameterized
surfaces, which would lead to equations featuring covariant derivatives (e.g., eqn. 81 in Section 5.2). In
the finite element setting, the surface of interest is represented by a collection of nodes and corresponding
surface normals, as shown in Figure 2. Using these surface normals, we perform surface projections of
the full 3-space derivatives following Jankuhn
et al.
[34]. In the previous section, we presented a formal
statement in eqn. 25 for handling the minimal Toner-Tu model on an arbitrary surface, using these surface
normals. We now recast those curved-surface equations once more, in a fashion consonant with a finite
element method solver. To formulate the equations in an expression convenient for the finite element
method, we aim to rewrite the Toner-Tu equations in the form
v
∂t
+
∇·
J
(
v
)
=
f
(
v
)
.
(28)
Here, the whole formulation comes down to the definitions of
J
(
v
)
and
f
(
v
)
. The flux-like quantity
J
(
v
)
is
a 3
×
3 matrix defined such that
∇·
J
(
v
)
= (
∂x
1
,
∂x
2
,
∂x
3
)
J
11
J
12
J
13
J
21
J
22
J
23
J
31
J
32
J
33
,
(29)
recalling that the divergence of a 2nd rank tensor is a vector. In indicial notation, this can be written as
(
∇·
J
(
v
)
)
i
=
∂J
(
v
)
ji
∂x
j
.
(30)
We need to define the components of the tensor
J
(
v
)
such that they yield the correct Toner-Tu terms.
As we will see below,
J
(
v
)
is not symmetric despite its superficial resemblance to a stress tensor. To set
notation and to make sure that the strategy is clear, we begin by demonstrating how to implement the
flat-space version of the Toner-Tu equations in a form consistent with eqn. 2. If we define the force term
f
(
v
)
as
f
(
v
)
i
= [
α
(
ρ
ρ
c
)
βv
j
v
j
]
v
i
σ
∂ρ
∂x
1
λv
j
∂v
i
∂x
j
,
(31)
then the remaining terms are captured if we define the tensor
J
(
v
)
as
J
(
v
)
=
D
∂v
1
∂x
1
D
∂v
2
∂x
1
D
∂v
3
∂x
1
D
∂v
1
∂x
2
D
∂v
2
∂x
2
D
∂v
3
∂x
2
D
∂v
1
∂x
3
D
∂v
2
∂x
3
D
∂v
3
∂x
3
.
(32)
.
CC-BY-NC-ND 4.0 International license
available under a
(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprint
this version posted June 25, 2022.
;
https://doi.org/10.1101/2022.06.22.497052
doi:
bioRxiv preprint
10
Considering the 1-component of velocity, we see that
(
∇·
J
(
v
)
)
1
=
∂J
j
1
∂x
j
=
∂J
11
∂x
1
+
∂J
21
∂x
2
+
∂J
31
∂x
3
.
(33)
Plugging in the components of
J
(
v
)
gives the result
(
∇·
J
(
v
)
)
1
=
∂x
1
(
D
∂v
1
∂x
1
)
+
∂x
2
(
D
∂v
1
∂x
2
)
+
∂x
3
(
D
∂v
1
∂x
3
)
.
(34)
Evaluating these derivatives leads to the result
(
∇·
J
(
v
)
)
1
=
D
2
v
1
∂x
2
1
D
2
v
1
∂x
2
2
D
2
v
1
∂x
2
3
,
(35)
as expected from the original Toner-Tu equations. Combining the contributions from eqns. 31 and 35
in the form
v
/∂t
+
∇·
J
(
v
)
=
f
(
v
)
, we recover eqn. 2 precisely as we set out to do. By defining the
quantities
J
(
v
)
and
f
(
v
)
, we have successfully reframed the flat-space Toner-Tu equations in a format
that will be conveniently implemented in the finite element setting. We now need to tackle the more
demanding formulation for an arbitrary curved surface.
Abstractly, our finite element version of the curved-space Toner-Tu equations is written as
v
∂t
=
div
Γ
J
(
v
)
+
f
(
v
)
,
(36)
where div
Γ
is the surface projected version of the divergence introduced in eqn. 22. We now introduce
the definition
J
(
v
)
ji
=
DP
il
G
lj
,
(37)
where we recall
P
from eqn. 3 and
G
from eqn. 20. Given this definition, we can attempt to write our
Toner-Tu equations as
∂v
i
∂t
=
D
(
(
P
il
G
lj
)
∂x
j
)
Γ
+
f
(
v
)
i
,
(38)
where we use the shorthand notation introduced in eqn. 23. This can be rewritten as
∂v
i
∂t
=
DP
il
(
∂G
lj
∂x
j
)
Γ
︷︷
P
div
Γ
(
Γ
v
)
+
D
(
∂P
il
∂x
j
)
Γ
G
lj
︷︷
unwanted term
+
f
(
v
)
i
.
(39)
Unfortunately, as written, these equations contain an extra term that is not present in the Toner-Tu
formulation. To remove this unwanted extra term, we must introduce a new term in the force
f
(
v
)
that
subtracts off the unwanted term,
f
fict
i
=
D
(
∂P
il
∂x
j
)
Γ
G
lj
.
(40)
We view this as a mathematical trick that allows us to use a formalism convenient for finite element
method analysis and have not sought to find a “physical interpretation” of this force in the way that
ideas such as the Coriolis force arises in mechanics. In this case, our
real
Toner-Tu equations can be
written as
∂v
i
∂t
=
DP
il
(
∂G
lj
∂x
j
)
Γ
︷︷
P
div
Γ
(
Γ
v
)
+
D
(
∂P
il
∂x
j
)
Γ
G
lj
︷︷
unwanted term
+
f
(
v
)
i
+
f
fict
i
.
(41)
.
CC-BY-NC-ND 4.0 International license
available under a
(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprint
this version posted June 25, 2022.
;
https://doi.org/10.1101/2022.06.22.497052
doi:
bioRxiv preprint
11
We now have precisely the equations we want, cast in the form we will use in the finite element setting.
As we saw above, the force term
f
(
v
)
needs to account for those terms that are not present in
∇·
J
(
v
)
and to subtract off terms that are present in
∇·
J
(
v
)
but unwanted. To that end, the 1-component of
f
(
v
)
takes the form
f
(
v
)
1
=
(
α
(
ρ
ρ
c
)
β
((
v
1
)
2
+ (
v
2
)
2
+ (
v
3
)
2
)
)
v
1
σ
(
∂ρ
∂x
1
)
Γ
λ
(
v
1
∂v
1
∂x
1
+
v
2
∂v
1
∂x
2
+
v
3
∂v
1
∂x
3
v
1
n
1
n
1
∂v
1
∂x
1
v
1
n
1
n
2
∂v
2
∂x
1
v
1
n
1
n
3
∂v
3
∂x
1
v
2
n
1
n
1
∂v
1
∂x
2
v
2
n
1
n
2
∂v
2
∂x
2
−−
v
2
n
1
n
3
∂v
3
∂x
2
v
3
n
1
n
1
∂v
1
∂x
3
v
3
n
1
n
2
∂v
2
∂x
3
−−
v
3
n
1
n
3
∂v
3
∂x
3
)
D
((
∂P
11
∂x
1
)
Γ
G
11
+
(
∂P
11
∂x
2
)
Γ
G
12
+
(
∂P
11
∂x
3
)
Γ
G
13
)
D
((
∂P
12
∂x
1
)
Γ
G
21
+
(
∂P
12
∂x
2
)
Γ
G
22
+
(
∂P
12
∂x
3
)
Γ
G
23
)
D
((
∂P
13
∂x
1
)
Γ
G
31
+
(
∂P
13
∂x
2
)
Γ
G
32
+
(
∂P
13
∂x
3
)
Γ
G
33
)
.
(42)
Note that the
α
and
β
terms capture the preferred speed contribution, the
σ
term captures the pressure,
the terms multiplied by
λ
capture the advection contribution to the Toner-Tu equations, and the final
terms involving
P
and
G
subtract off the fictitious force. We can repeat a similar analysis for the 2- and
3- components of the force as
f
(
v
)
2
=
(
α
(
ρ
ρ
c
)
β
((
v
1
)
2
+ (
v
2
)
2
+ (
v
3
)
2
)
)
v
2
σ
(
∂ρ
∂x
2
)
Γ
λ
(
v
1
∂v
2
∂x
1
+
v
2
∂v
2
∂x
2
+
v
3
∂v
2
∂x
3
v
1
n
2
n
1
∂v
1
∂x
1
v
1
n
2
n
2
∂v
2
∂x
1
v
1
n
2
n
3
∂v
3
∂x
1
v
2
n
2
n
1
∂v
1
∂x
2
v
2
n
2
n
2
∂v
2
∂x
2
−−
v
2
n
2
n
3
∂v
3
∂x
2
v
3
n
2
n
1
∂v
1
∂x
3
v
3
n
2
n
2
∂v
2
∂x
3
−−
v
3
n
2
n
3
∂v
3
∂x
3
)
D
((
∂P
21
∂x
1
)
Γ
G
11
+
(
∂P
21
∂x
2
)
Γ
G
12
+
(
∂P
21
∂x
3
)
Γ
G
13
)
D
((
∂P
22
∂x
1
)
Γ
G
21
+
(
∂P
22
∂x
2
)
Γ
G
22
+
(
∂P
22
∂x
3
)
Γ
G
23
)
D
((
∂P
23
∂x
1
)
Γ
G
31
+
(
∂P
23
∂x
2
)
Γ
G
32
+
(
∂P
23
∂x
3
)
Γ
G
33
)
,
(43)
.
CC-BY-NC-ND 4.0 International license
available under a
(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprint
this version posted June 25, 2022.
;
https://doi.org/10.1101/2022.06.22.497052
doi:
bioRxiv preprint
12
and
f
(
v
)
3
=
(
α
(
ρ
ρ
c
)
β
((
v
1
)
2
+ (
v
2
)
2
+ (
v
3
)
2
))
v
3
σ
(
∂ρ
∂x
3
)
Γ
λ
(
v
1
∂v
3
∂x
1
+
v
2
∂v
3
∂x
2
+
v
3
∂v
3
∂x
3
v
1
n
3
n
1
∂v
1
∂x
1
v
1
n
3
n
2
∂v
2
∂x
1
v
1
n
3
n
3
∂v
3
∂x
1
v
2
n
3
n
1
∂v
1
∂x
2
v
2
n
3
n
2
∂v
2
∂x
2
−−
v
2
n
3
n
3
∂v
3
∂x
2
v
3
n
3
n
1
∂v
1
∂x
3
v
3
n
3
n
2
∂v
2
∂x
3
−−
v
3
n
3
n
3
∂v
3
∂x
3
)
D
((
∂P
31
∂x
1
)
Γ
G
11
+
(
∂P
31
∂x
2
)
Γ
G
12
+
(
∂P
31
∂x
3
)
Γ
G
13
)
D
((
∂P
32
∂x
1
)
Γ
G
21
+
(
∂P
32
∂x
2
)
Γ
G
22
+
(
∂P
32
∂x
3
)
Γ
G
23
)
D
((
∂P
33
∂x
1
)
Γ
G
31
+
(
∂P
33
∂x
2
)
Γ
G
32
+
(
∂P
33
∂x
3
)
Γ
G
33
)
.
(44)
By setting
v
/∂t
equal to the sum of
div
Γ
J
(
v
)
and the
f
(
v
)
terms, we have fully reproduced the
curved-space version of the Toner-Tu equations.
We express the continuity equation in similar form,
∂ρ
∂t
=
div
Γ
J
(
ρ
)
+
f
(
ρ
)
,
(45)
where
J
(
ρ
)
is now a vector and the scalar
f
(
ρ
)
can be thought of as a source term. We set div
Γ
J
(
ρ
)
equal
to zero and thus define
f
(
ρ
)
=
div
Γ
(
ρ
v
) =
(
(
ρv
l
)
∂x
l
n
l
n
k
(
ρv
l
)
∂x
k
)
=
(
∂ρv
1
∂x
1
)
Γ
(
∂ρv
2
∂x
2
)
Γ
(
∂ρv
3
∂x
3
)
Γ
.
(46)
Altogether, eqns. 32, 42, 43, 44, and 46 comprise a complete curved-surface FEM implementation
of the minimal Toner-Tu equations and make possible the numerical results presented in the remainder
of this work. We also refer the interested reader to Section 8.2 in the Appendix for details on our
implementation in the specific commercial finite element package COMSOL Multiphysics
R
©
.
4. Parameters and Dimensionless Ratios for the Theory
4.1. Parameter Choices for Wildebeest Herds
To put our finite element formulation into numerical action, we must of course adopt specific values
for the parameters
ρ
0
,
ρ
c
,
α
,
β
,
σ
,
D
, and
λ
that appear in our governing partial differential equations. In
this study, we focus on the macroscopic length scale of animal herds, although we consider the microscopic
activity of cytoskeletal proteins elsewhere [27]. While our goal here is to explore the phenomenology of
Toner-Tu flocks on curved surfaces in a general way, not to claim an understanding of specific animal
behavior, our parameter choices are loosely inspired by migrating wildebeest herds. By inspecting aerial
photographs of wildebeest herds, we estimated an average wildebeest density of
ρ
0
= 0
.
25 1
/
m
2
.
(47)
.
CC-BY-NC-ND 4.0 International license
available under a
(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprint
this version posted June 25, 2022.
;
https://doi.org/10.1101/2022.06.22.497052
doi:
bioRxiv preprint
13
For the critical density above which coordinated herding behavior occurs, we make the estimate
ρ
c
= 0
.
05 1
/
m
2
,
(48)
based roughly on the observation that for densities much higher than this, the photographed wildebeest
have an organized herd structure. We note that it would be very interesting to carefully measure these
parameters in the context of herds of wildebeest with a Toner-Tu framework in mind.
The remainder of the parameters are determined in the spirit of exploring how the different terms
compete to alter the density and velocity fields. Thus, parameters are chosen in order to make all the terms
comparable in magnitude. We provide here the actual values used in our finite element calculations, fully
cognizant that our parameter choices are at best approximate. We begin by estimating the magnitude of
∂v/∂t
. We picture wildebeests circling a small hill in the landscape, moving with a characteristic speed
of 1 m/s and taking 20 s to change direction completely. This scenario implies
∂v
∂t
1 m/s
(
1 m/s)
20 s
0
.
1 m/s
2
.
(49)
Using the strategy of balancing the magnitudes of the different terms, we estimate the coefficient for the
preferred speed terms by considering that
0
.
1 m/s
2
=
α
(
ρ
0
ρ
c
)
v
α
×
(0
.
25 1
/
m
2
0
.
05 1
/
m
2
)
×
1 m/s
,
(50)
which leads us to adopt
α
= 0
.
5 m
2
/s
.
(51)
For the parameter
β
, we choose a value that sets the correct mean field speed,
v
pref
, for wildebeests.
Thus, we adopt
β
=
α
(
ρ
ρ
c
)
v
2
pref
(
0
.
5 m
/
s
2
)(
0
.
2 1
/
m
2
)
(1 m
/
s)
2
= 0
.
1 s
/
m
2
.
(52)
We can also check independently that the magnitude of the
β
term is of order 0.1
m/s
2
. Indeed,
βv
3
0
.
1
s
m
2
×
(1 m/s)
3
= 0
.
1 m/s
2
.
(53)
To find the parameter
σ
, the coefficient of the pressure term, we consider the gradients in density seen at
the edge of a wildebeest herd and estimate that
ρ
drops from
ρ
0
= 25 1
/m
2
to 0 1
/m
2
over 25 m. Using
these numbers implies
0
.
1 m/s
2
=
σ
∂ρ
∂x
σ
0
.
25 1
/
m
2
25 m
(54)
which leads us to adopt
σ
= 10 m
4
/
s
2
.
(55)
We can now apply this thinking to make an estimate of the neighbor coupling coefficient
D
by using the
equality
0
.
1 m/s
2
=
D
2
v
∂x
2
D
1 m/s
(10 m)
2
(56)
which leads to the coefficient choice
D
= 10 m
2
/
s
.
(57)
To estimate the magnitude of
λ
, we imagine that the wildebeest will come to a full stop from a speed of
1 m/s over a distance of roughly 10 m, permitting us to make the correspondence
0
.
1 m/s
2
=
λv
∂v
∂x
λ
×
1 m/s
×
1 m/s
10 m
(58)
.
CC-BY-NC-ND 4.0 International license
available under a
(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprint
this version posted June 25, 2022.
;
https://doi.org/10.1101/2022.06.22.497052
doi:
bioRxiv preprint
14
which implies that
λ
= 1
.
(59)
This estimated set of parameters provides us with a complete description of the Toner-Tu herd in the
minimal model we seek to explore. Throughout the remainder of the paper, we consistently use these
parameter values. In later sections, in some cases we will go beyond this idealized parameter set to
perform parameter sweeps that permit different regimes of behavior to emerge. Unless noted otherwise,
time-dependent finite element method simulations were initialized with a uniform density field
ρ
=
ρ
0
and a disordered velocity field, with velocity orientations drawn randomly from a uniform distribution
of angles between 0 and 2
π
and with velocity magnitude
v
(0) = 1 m/s, our estimated characteristic
wildebeest speed. We note again that these parameters were chosen to highlight competition between
the different terms in the dynamics; they are far from the final word for describing real wildebeest herds.
4.2. Dimensionless Representation of the Theory
We next recast the Toner-Tu equations in dimensionless form, rendering the meaning and magnitude of
the terms more transparent. For the positional coordinate, we take
x
=
x/L
, where
L
is a characteristic
length scale in the problem. In this case, we conceptualize
L
as the length scale of the wildebeest herd.
We note that in the problems considered here, the scale of the herd is the same as the scale of the surface
landscape on which the herd moves. Similarly, we use a characteristic velocity scale
U
to define the
dimensionless velocity as
v
=
v/U
. The density can be rescaled by using the critical density
ρ
c
as our
scaling variable, resulting in
ρ
=
ρ/ρ
c
. Lastly, in light of the definitions above, we can determine a time
scale
L/U
, the time it takes for a given wildebeest to cross the entire surface landscape. Thus, we define
t
=
t/
(
L/U
).
Using the various definitions given above, we can now rewrite the Toner-Tu equations using the
dimensionless versions of
t
,
x
,
ρ
and
v
as
U
2
L
∂v
i
∂t
=
αρ
c
U
(
ρ
1)
v
i
βU
3
|
v
|
2
v
i
σρ
c
L
∂ρ
∂x
i
+
DU
L
2
2
v
i
λU
2
L
v
j
∂v
i
∂x
j
.
(60)
Dividing everything by
U
2
/L
results in five dimensionless parameters whose magnitudes provide a sense
of the relative contributions of the different terms, within the full dynamical equations of the form
∂v
i
∂t
=
αρ
c
L
U
(
ρ
1)
v
i
βLU
|
v
|
2
v
i
σρ
c
U
2
∂ρ
∂x
i
+
D
UL
2
v
i
λv
j
∂v
i
∂x
j
.
(61)
Using definitions of dimensionless variables given above, we find that the generalized continuity equation
takes the form
∂ρ
∂t
=
(
ρ
v
i
)
∂x
i
.
(62)
Next, we briefly turn to interpreting the dimensionless ratios that appear when casting the minimal
Toner-Tu theory in dimensionless form as exhibited in eqn. 61. The two components of the preferred
speed term have dimensionless parameters given by
preferred speed
α
term =
αρ
c
L
U
=
L/U
1
/
(
αρ
c
)
=
time for wildebeest to cross herd
time for speed to increase to
v
pref
(63)
and
preferred speed
β
term =
βLU
=
L/U
1
/
(
βU
2
)
=
time for wildebeest to cross herd
time for speed to decrease to
v
pref
.
(64)
Each of these terms provides intuition about how quickly the wildebeest will return to their steady-state
speed given some perturbation that disturbs them from that value. Recall that
L/U
(the “time for
.
CC-BY-NC-ND 4.0 International license
available under a
(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprint
this version posted June 25, 2022.
;
https://doi.org/10.1101/2022.06.22.497052
doi:
bioRxiv preprint