#  Navigation of  _C. elegans_ in three-dimensional media: roll maneuvers and
planar turns

Alejandro Bilbao  Department of Mechanical Engineering, Texas Tech University,
Lubbock, Texas, U.S.A.  Amar K. Patel  Department of Mechanical Engineering,
Texas Tech University, Lubbock, Texas, U.S.A.  Mizanur Rahman  Department of
Chemical Engineering, Texas Tech University, Lubbock, Texas, U.S.A.  Siva A.
Vanapalli  Department of Chemical Engineering, Texas Tech University, Lubbock,
Texas, U.S.A.  Jerzy Blawzdziewicz  Departments of Mechanical Engineering and
Physics, Texas Tech University, Lubbock, Texas, U.S.A.

###### Abstract

Free-living nematode Caenorhabditis elegans is a powerful genetic model,
essential for investigations ranging from behavior to neuroscience to aging,
and locomotion is a key observable used in these studies. However, despite the
fact that in its natural environment C. elegans moves in three-dimensional
(3D) complex media (decomposing organic matter and water), quantitative
investigations of its locomotion have been limited to two-dimensional (2D)
motion. Based on our recent quantitative analysis of 2D turning maneuvers
[Phys. Fluids 25, 081902 (2013)] we follow with the first quantitative
description of how C. elegans moves in 3D environments. We show that by
superposing body torsion and 2D undulations, a burrowing or swimming nematode
can rotate the undulation plane. A combination of these roll maneuvers and 2D
turns associated with variation of undulation-wave parameters allows the
nematode to explore 3D space. We apply our model to analyze 3D chemotaxis of
nematodes burrowing in a gel and swimming in water; we conclude that the
nematode can achieve efficient chemotaxis in different environments without
adjusting its sensory-motor response to chemical signals. Implications of our
findings for understanding of 3D neuromuscular control of nematode body are
discussed.

Caenorhabditis elegans, locomotion, hydrodynamic bead models, chemotaxis,
neuromuscular control, 3D environment

##  1 Introduction

Investigations of locomotion of a free-living nematode C. elegans (one of the
most important model organisms) can reveal a wealth of information on
physiology, muscle biology, and neural control of movement and behavior.
Efficient locomotion depends on muscular strength Ghanbari  _et al._ (2012);
Johari  _et al._ (2013); Etheridge  _et al._ (2015), neuromuscular
coordination Wen  _et al._ (2012); Stirman  _et al._ (2011), and information
processing Gray  _et al._ (2005); Kocabas  _et al._ (2012). Locomotory
readouts are thus extensively used to describe impacts of genetic mutations
Sznitman  _et al._ (2010a); Brown  _et al._ (2013); Restif  _et al._ (2014)
and evaluate pharmacological effects Artal-Sanz  _et al._ (2011); Omura  _et
al._ (2012).

In its natural environment (soil rich in decaying organic matter), C. elegans
interacts with a complex three-dimensional (3D) surrounding medium that
includes rainwater and soft organic materials. Most studies of nematode
locomotion, however, have focused on two-dimensional (2D) motion Cohen and
Boyle (2010); Fang-Yen  _et al._ (2010); Sznitman  _et al._ (2010b);
Padmanabhan  _et al._ (2012); Stephens  _et al._ (2008), and little effort has
been made to characterize worm motility in 3D environments that would mimic C.
elegans natural habitat.

Only recently was 3D burrowing motion of C. elegans in soft substrates
quantitatively imaged Kwon  _et al._ (2013) and feasibility of using 3D motion
to assess genetic mutations Kwon  _et al._ (2015) and to identify muscular
dystrophy Beron  _et al._ (2015) experimentally demonstrated. Such 3D methods
of scoring C. elegans locomotory behavior are promising; however, unleashing
the full potential of such techniques requires complementary modeling efforts
to elucidate the biomechanics and neuromuscular control of 3D locomotion.

Very few quantitative analyses of 3D locomotion of C. elegans have been
reported Deng and Xu (2014); Delmotte  _et al._ (2015); no investigations
present a realistic mathematical description of nematode motion. Our study
aims to fill this knowledge gap by providing quantitative understanding of the
geometry and mechanics of 3D maneuvers of C. elegans burrowing in gel-like
materials and swimming in Newtonian fluids.

Our modeling approach generalizes the recently proposed piecewise harmonic
curvature (PHC) model of 2D nematode gait Padmanabhan  _et al._ (2012). The
PHC model has been shown to give an accurate parameterization of body postures
(defined by the configuration of the body centerline Stephens  _et al._
(2008); Sznitman  _et al._ (2010c); Husson  _et al._ ; Berman  _et al._
(2013)) and of entire trails of nematodes crawling on agar without sidewise
slip; it also yields a convenient quantification of time-lapse sequences of
body postures in swimming.

In the present paper we expand the PHC gait model by including curve torsion
in the description of the nematode postures. We show that combining the PHC
model for the line curvature with a piecewise constant torsion (PCT)
assumption, we can generate three-dimensional trajectories that mimic paths of
burrowing nematodes. Using, in addition, our hydrodynamic models of swimming
dynamics Bilbao  _et al._ (2013), we can describe 3D motion of C. elegans in
bulk fluids. In particular, our gait model reproduces the rapid undulation-
plane reorientation observed in our experiments.

The explicit PHC representation of the nematode gait has been shown to be a
useful tool for quantitative investigations of swimming efficiency and 2D
turning maneuvers Bilbao  _et al._ (2013). The combined PHC/PCT gait model
will help analyze diverse aspects of 3D undulatory locomotion of C. elegans in
complex media, including chemotaxis and mechanisms of neuro-muscular actuation
of body movements. To illustrate the former application we present an analysis
of the efficiency of 3D chemotaxis of burrowing and swimming C. elegans.

Our paper is organized as follows. In Sec. 2 we develop geometrical gait
models for C. elegans moving in 2D and 3D. We also assess the effectiveness of
reorientation maneuvers of animals crawling on a flat 2D substrate and
burrowing in soft 3D matter. In Sec. 3 we present our hydrodynamic analysis of
reorientation maneuvers of C. elegans swimming in a viscous fluid. The
geometrical and hydrodynamic locomotion descriptions developed in Secs. 2 and
3 are combined in Sec. 4 with simple models of chemosensation and motion
control to describe chemotaxis of burrowing and swimming C. elegans in 3D.
Implications of our findings for the emerging area of 3D undulatory locomotion
are discussed in Sec. 5.

Figure 1:  (Colour online) Family of curves reproduced from harmonic curvature
(3) for normalized amplitude $A/q$ as labeled. The green and blue arrows are
the unit tangent and normal vectors, respectively. The body of a nematode
moving along the curves is represented by the thick line segments of length
$L$. Two normalized wavevectors are shown (a) $qL=5.5$ (typical swimming gait)
and (b) $qL=9$ (typical crawling gait). Depending on the wavevector, the
normalized amplitude $A/q=1$ may correspond to (a) C-shaped (top) and (b)
W-shaped body postures; $A/q=1.8$ produces an $\Omega$-shaped body posture.
Figure 2:  A crawling nematode performing three-mode planar turn with the
default-mode amplitude $A_{1}/q=1$ and turning-mode amplitude $A_{2}/q$ as
labeled, for $qs^{c}_{1}=3\pi/2$ and $q\Delta s^{c}=0.6$. The turn angle
$\theta^{\textrm{turn}}$ and mode-switching points $s^{c}_{1}$ and $s^{c}_{2}$
are defined in the top inset. Bottom inset shows the dependence of the turn
angle $\theta^{\textrm{turn}}$ (in degrees) on the turning-mode amplitude
$A_{2}/q$. The black dots correspond to the amplitudes of the trajectories
shown.

##  2 Geometrical description of the nematode gait

Caenorhabditis elegans crawling on or burrowing through a soft gel-like
substrate (such as agar) often move with no transverse slip with respect to
the environment. On a moist porous surface, slip may be hindered by
interfacial forces associated with water meniscus surrounding the nematode
body; in a 3D gel-like medium, transverse displacements are prevented by a
nonzero yield stress.

In the absence of transverse slip, all segments of the nematode body follow a
single trail, which is completely determined by the geometry of the crawling
or burrowing gait. Thus, for a given sequence of nematode movements, the
trajectory does not depend on the substrate mechanics. Here we focus on this
purely geometrical no-slip case.

In Sec. 2.1 we summarize our previously reported PHC model Padmanabhan  _et
al._ (2012) for describing body postures and trajectories of C. elegans moving
in two dimensions. In Sec. 2.2 we expand our model by adding a torsion
component to the 2D gait representation to achieve a 3D generalization of the
previously established PHC approach.

Based on our earlier analyses of C. elegans swimming in two dimensions
Padmanabhan  _et al._ (2012); Bilbao  _et al._ (2013), we assume that PHC
model and its 3D generalization provide a realistic representation for the
swimming gait. In Sec. 3 the results derived in Sec. 2 are used to develop a
comprehensive analysis that combines the gait geometry and motion mechanics to
fully describe 3D maneuvers of swimming C. elegans.

###  2.1 Description in 2D

The body of C. elegans is highly elongated, and therefore its postures can be
faithfully represented by the body centerline Padmanabhan  _et al._ (2012);
Stephens  _et al._ (2008); Kim  _et al._ (2011). In our analysis, the shape of
the centerline is described by the worm curvature function
$\kappa_{w}(s^{\prime},t)$, where $t$ denotes time and $s^{\prime}$ is the
arclength coordinate along the body. In what follows $s^{\prime}=0$ is the
tail and $s^{\prime}=L$ is the head coordinate value (where $L$ is the
nematode length). The body configuration is fully determined by providing the
curvature function $\kappa_{w}$ and the position and orientation of the
nematode head (or any other reference point along the body).

A curvilinear trail followed by the body of a nematode crawling with no
transverse slip can be described using the trail curvature $\kappa(s)$, where
$s$ is the arclength coordinate along the trail. Since at any given time $t$
the configuration of the nematode body coincides with a trail segment (see
Fig. 1), the worm and trail curvature functions are related

$\kappa_{w}(s^{\prime},t)=\kappa(s^{\prime}+v_{s}t),$ (1)

where $v_{s}$ is the velocity along the trail. Body postures and the trail
followed by a nematode are constructed from the curvature (1) by solving the
Serret–Frenet equations Stoker (1969)

$\displaystyle\frac{d{\hat{\boldsymbol{t}}}}{ds}$
$\displaystyle=\phantom{-}\kappa{\hat{\boldsymbol{n}}},$ (2a)
$\displaystyle\frac{d{\hat{\boldsymbol{n}}}}{ds}$
$\displaystyle=-\kappa{\hat{\boldsymbol{t}}},$ (2b)

where ${\hat{\boldsymbol{t}}}$ and ${\hat{\boldsymbol{n}}}$ are mutually
perpendicular unit tangent and normal vectors.

Observations of C. elegans swimming in water indicate that the propagating-
curvature-wave relation (1) applies not only to no-slip motion [in which case
Eq. (1) results from geometrical constraints], but is also preserved with good
accuracy when significant transverse slip occurs Padmanabhan  _et al._ (2012).
This behavior is a signature of the proprioceptive coupling Wen  _et al._
(2012) that allows the nematode to maintain a well organized gait in
environments with diverse mechanical properties. Based on the similarity
between the curvature wave form for crawling and swimming nematodes we
formulate a unified gait description for different media.

####  2.1.1 PHC representation of nematode gait

Caenorhabditis elegans produces a variety of body postures, ranging from
shallow undulations to highly curved $\Omega$ shapes and loop shapes.
Recently, we have demonstrated Padmanabhan  _et al._ (2012) that all these
body shapes can be accurately represented using the harmonic curvature
function

$\kappa(s)={A}\cos\left(qs+\phi\right),$ (3)

where $A$ is the wave amplitude, $q=2\pi/\lambda$ is the wavevector (with
$\lambda$ denoting the wavelength), and $\phi$ is the phase angle.

Combining Eqs. (1)–(3) yields a family of curves representative of the most
common nematode postures. For moderate values of the dimensionless amplitude
$A/q\approx 1$, typical crawling postures (W-shaped) and swimming postures
(C-shaped) are obtained, as illustrated in Fig. 1. Higher amplitudes
correspond to $\Omega$-shapes ($A/q\approx 2$) and loop shapes ($A/q\approx
3$, described in Padmanabhan  _et al._ (2012)).

According to the PHC model Padmanabhan  _et al._ (2012), trails of crawling
nematodes can be accurately represented by the harmonic function (3) with
abruptly changing wave parameters $A$, $q$, and $\phi$ between harmonic-
curvature intervals. Such changes generate either small fluctuations in the
direction of motion or significant body reorientation, depending on the
magnitude of the parameter variation and the mode switching point.

####  2.1.2 Planar turns

To investigate the turning ability of C. elegans, we consider three-mode PHC
maneuvers. More complex trajectories can be obtained by combining three-mode
turns. For simplicity, we focus on maneuvers in which only the amplitude of
the harmonic curvature varies,

$\kappa(s)=\left\\{\begin{matrix}{A_{1}}\cos\left(qs\right),\quad
s<s^{c}_{1},\phantom{--.}\\\ {A_{2}}\cos\left(qs\right),\quad
s^{c}_{1}<s<s^{c}_{2},\\\ {A_{1}}\cos\left(qs\right),\quad
s^{c}_{2}<s,\phantom{--.}\\\ \end{matrix}\right.$ (4)

where $A_{1}$ is the base crawling amplitude (with $A_{1}/q\approx 1$),
$A_{2}>A_{1}$ is the turning amplitude, and $s^{c}_{1}$ and $s^{c}_{2}$ are
the mode-switching points. Planar turns can also be obtained by tuning the
wavevector and phase angle; an analysis of such turns is presented in the
Electronic Supplementary Material (ESM).

Figure 2 shows an example of a crawling nematode executing a three-mode planar
turn for fixed switching points $s^{c}_{1}$ and $s^{c}_{2}$ and several values
of the normalized turning mode amplitudes $A_{2}/q$. The dependence of the
turn angle $\theta^{\textrm{turn}}$ (defined in the top inset of Fig. 2) on
the turn amplitude $A_{2}$ is depicted in the bottom inset of Fig. 2.

Figure 3:  (Colour online). Curves defined by the harmonic curvature
(normalized amplitude $A/q=1$) and constant torsion (normalized magnitude
$\tau_{0}L$ as labeled); worm contours (thick line segments) correspond to
$qL=9$. The green, blue and red arrows are the unit tangent, normal, and
binormal vectors, respectively.  Figure 4: Roll rotation rate per undulation
period $\omega^{\textrm{roll}}$ (in degrees) vs normalized body torsion
$\tau_{0}/q$ for one-mode torsional roll with harmonic-curvature amplitude
$A/q$ as labeled. Inset shows the dependence of the line slope
$\omega^{\textrm{roll}}q/\tau_{0}$ on the amplitude $A/q$.  Figure 5: Nematode
reorientation in 3D. The nematode (a) undulates in the initial plane of
motion; (b) executes a torsional roll maneuver; and (c) performs planar turn
in the new plane of motion.

###  2.2 Description in 3D

Caenorhabditis elegans, which crawls and burrows in rotting vegetation and
swims in rainwater, performs 2D and 3D maneuvers Kwon  _et al._ (2013);
Voorhies  _et al._ (2005); Frezal and Felix (2015) to find optimal
environmental conditions. Existing literature focuses nearly exclusively on 2D
motion, and 3D models of nematode gait and locomotion mechanics are
unavailable. In this section we expand our PHC model to incorporate 3D
nematode postures and trails into the gait representation.

The proposed description is motivated by recent observations of 3D body
postures in crawling Deng and Xu (2014), burrowing Kwon  _et al._ (2013), and
swimming. A crawling C. elegans, for example, can produce torsional
deformation of its body, resulting in detachment of some body sections from
the underlying surface, as shown in the images of crawling worms in Ref. Deng
and Xu (2014). Swimming nematodes undergo a similar out-of-plane body
deformations, accompanied by rapid rotation of the undulation plane (see Sec.
3). Three-dimensional body actuation seems also to occur for nictating
nematodes (see supplemental movies in Lee  _et al._ (2012)).

To incorporate out-of-plane body postures into the PHC model, we use the
complete 3D set of Serret–Frenet equations,

$\displaystyle\frac{d{\hat{\boldsymbol{t}}}}{ds}$
$\displaystyle=\phantom{-}\kappa{\hat{\boldsymbol{n}}},$ (5a)
$\displaystyle\frac{d{\hat{\boldsymbol{n}}}}{ds}$
$\displaystyle=-\kappa{\hat{\boldsymbol{t}}}+\tau{\hat{\boldsymbol{b}}},$ (5b)
$\displaystyle\frac{d{\hat{\boldsymbol{b}}}}{ds}$
$\displaystyle=-\tau{\hat{\boldsymbol{n}}},$ (5c)

which includes both the curvature and torsion components. Here $\tau=\tau(s)$
is the torsion of the curve defining the 3D nematode trail,
${\hat{\boldsymbol{t}}}$ is the tangent unit vector, ${\hat{\boldsymbol{n}}}$
is the normal unit vector in the local undulation plane, and
${\hat{\boldsymbol{b}}}={\hat{\boldsymbol{t}}}\times{\hat{\boldsymbol{n}}}$ is
the binormal unit vector (normal to the local undulation plane). The curve
torsion $\tau$, which corresponds to the rotation of the Serret–Frenet basis
$({\hat{\boldsymbol{t}}},{\hat{\boldsymbol{n}}},{\hat{\boldsymbol{b}}})$ about
the tangent direction ${\hat{\boldsymbol{t}}}$, results in out-of-plane
deformations of the curve.

####  2.2.1 Constant-torsion assumption

Quantitative experimental information regarding out-of-plane body postures of
C. elegans is very limited. However, a qualitative examination of nematode
images suggests that key features of 3D nematode gait can be reproduced by
combining the harmonic curvature model (3) with the constant-torsion
assumption,

$\tau(s)=\tau_{0}.$ (6)

The solution of the Serret–Frenet equations (5) with the curvature and torsion
given by equations (3) and (6) is depicted in Fig. 3. The results indicate
that the general shape of the body wave is preserved in the presence of
nonzero torsion, but the undulation plane rotates with the angular velocity
proportional to the torsion amplitude $\tau_{0}$. This behavior is in a
qualitative agreement with the crawling nematode image shown in Ref. Deng and
Xu (2014). The gait described by Eqs. (3) and (6) will be called a single-mode
torsional roll.

Figure 4 shows the torsional-roll rotation rate per undulation period,

$\omega^{\textrm{roll}}=\cos^{-1}\left[{\hat{\boldsymbol{b}}}(0)\boldsymbol{\cdot}{\hat{\boldsymbol{b}}}(2\pi/q)\right].$
(7)

The rotation rate (7) is determined from the angle between the binormal
vectors at $s=0$ and $s=2\pi/q$; it depends linearly on the dimensionless
torsion $\tau_{0}/q$. According to the results in Fig. 4, the slope of the
lines $\omega^{\textrm{roll}}$ vs $\tau_{0}/q$ decreases with the increasing
dimensionless curvature-wave amplitude $A/q$. The decrease stems from the fact
that only the projection of the torsion vector
$\boldsymbol{\tau}=\tau{\hat{\boldsymbol{t}}}$ onto the average direction of
motion contributes to the roll.

####  2.2.2 Three-mode roll maneuvers

In the case of variable torsion, we assume that the torsion wave
$\tau_{w}(s^{\prime},t)$ propagates backward along the nematode body,

$\tau_{w}(s^{\prime},t)=\tau(s^{\prime}+v_{s}t),$ (8)

similar to the curvature wave (1). For no-slip motion (e.g., burrowing in a
gel), relation (8) is required by the constraints that all segments of the
nematode body follow a single trail. By analogy with our 2D results
Padmanabhan  _et al._ (2012), we assume that the form (8) also applies to
swimming.

According to our qualitative observations of nematode motion in water, C.
elegans survey their environment producing roll reorientation maneuvers at
irregular intervals. To replicate such motions in burrowing and swimming, we
introduce a three-mode roll, where the nematode performs one-mode harmonic
curvature undulations (1) superposed with piecewise constant torsion

$\tau(s)=\left\\{\begin{matrix}{0},\quad s<s^{\tau}_{1},\phantom{--.}\\\
{\tau_{0}},\quad s^{\tau}_{1}<s<s^{\tau}_{2},\\\ {0},\quad
s^{\tau}_{2}<s,\phantom{--.}\\\ \end{matrix}\right.$ (9)

(where $s^{\tau}_{1}$ and $s^{\tau}_{2}$ are the mode switching points, and
$\Delta s^{\tau}=s^{\tau}_{2}-s^{\tau}_{1}$ is the roll-mode length). By
combining torsional roll maneuvers with three-mode planar turns, a nematode is
capable of exploring 3D space (see Fig. 5 and Video S1 in ESM).

Figure 6:  (Colour online) A sequence of snapshots of a swimming C-shaped
nematode ($qL=5.5$) performing a three-mode planar turn with the default-mode
amplitude $A_{1}/q=1$ and turning mode amplitude as labeled, for
$qs^{c}_{1}=3\pi/2$ and $q\Delta s^{c}=0.6$. The gray image represents the
starting position for all trajectories. The inset shows the dependence of the
turn angle $\theta^{\textrm{turn}}$ (in degrees) on the turning-mode amplitude
$A_{2}/q$. The dots correspond to the snapshots shown.  Figure 7:  Planar turn
angle $\theta^{\textrm{turn}}$ (in degrees) vs normalized turning-mode length
$\Delta s^{c}/L$ for (a) C-shaped ($qL=5.5$) and (b) W-shaped ($qL=9$)
nematodes (see insets) and turning-mode amplitude $A_{2}/q$ as labeled. Turns
initiated at $qs^{c}_{1}=0$ (solid lines) and $\frac{3}{2}\pi$ (dashed). The
normalized amplitude of the default forward mode is $A_{1}/q=1$.

##  3 Nematode maneuverability in fluids

The PHC and PHC/PCT gait representations developed in Sec. 2 will now be
combined with a hydrodynamic bead-based model to describe body reorienting
maneuvers performed by swimming C. elegans. To account for hydrodynamic
forces, the nematode body is approximated by an active chain of
hydrodynamically coupled spheres arranged along the body centerline described
by Eqs. (4) and (9) (see simulation images shown in Fig. 6).

The hydrodynamic forces and torques acting on the chain are calculated in the
Stokes-flow limit using the generalized Rotne–Prager approximation Wajnryb
_et al._ (2013). This approximation is much faster (although less accurate)
than the Hydromultipole algorithm Cichocki  _et al._ (1994) used in our
previous study Bilbao  _et al._ (2013).

Our methods for calculating the hydrodynamic friction and mobility tensors are
described in Ref. Bilbao  _et al._ (2013) and in ESM. We note that in an
independent study Berman et al. Berman  _et al._ (2013) also used a
hydrodynamic bead-based model to investigate undulatory swimming of a
nematode.

###  3.1 Swimming in 2D

To assess 2D reorientation ability of swimming C. elegans, we consider
trajectories of a nematode performing three-mode PHC turns, and evaluate the
angle $\theta^{\textrm{turn}}$ between the initial and final swimming
directions. Similar to our analysis of no-slip turns (Sec. 2.1.2), the
nematode initially swims using the default forward-locomotion mode with the
amplitude $A/q=1$, switches to a higher-amplitude turning mode, and returns to
the default mode. Due to hydrodynamic effects, the turn angle in swimming
depends not only on the curvature-wave and torsion-wave parameters, but also
on the nematode length $L$.

A comparison of the numerical simulations depicted in Fig. 6 with the
corresponding results for a nematode crawling without sidewise slip (Fig. 2)
indicates that a swimming nematode produces shallower turns than its crawling
counterpart. Moreover, due to the nonlinear dependence of the hydrodynamic
interactions on the body posture, in swimming the relation between the turn
angle $\theta^{\textrm{turn}}$ and the amplitude of the turning mode $A_{2}/q$
is nonlinear (in contrast to the linear dependence for no-slip motion).

Figure 7 provides further details regarding the dependence of the turn angle
$\theta^{\textrm{turn}}$ on the turning-mode parameters. The results, for
different values of the turning mode amplitude, are shown vs turning-mode
length $\Delta s^{c}/L$ for two wavevector values: $qL=5.5$ and $qL=9$ (i.e.,
for C-shaped and W-shaped gaits). In the scaling used in Fig. 7, i.e., with
the mode length $\Delta s^{c}$ normalized by the worm length $L$, the results
differ significantly between the two values of $qL$. We note, however, that
the dependence of the turn angle on the gait parameter $qL$ is much weaker
when the $s$ coordinate is normalized by the wavevector $q$ Bilbao  _et al._
(2013), which indicates that the nematode length has only a moderate effect on
the turn dynamics.

According to the results shown in Figs. 6 and 7, a swimming nematode can
control the turn angle by varying the turning-mode amplitude $A$ as well as
the turn initiation point $s^{c}_{1}$ and turning-mode length $\Delta s^{c}$.
Similar to the crawling case, a swimming nematode can fully explore 2D space
by performing a sequence of elementary turns. For a qualitative comparison of
a 2D multi-mode turn of a swimming wild-type nematode with a corresponding
numerical simulation using the PHC model see video S2 in ESM.

Figure 8:  Time-lapse images of a swimming C. elegans switching the plane of
undulations.  Figure 9:  (a) Non-planar body posture of a swimming C. elegans
and (b) qualitatively matching calculated body shape with constant torsion
$\tau_{0}L=0.9$ and harmonic curvature of the normalized amplitude $A/q=1$ and
wavevector $qL=5.5$.  Figure 10: Stereoscopic views of a C-shaped nematode
($qL=5.5$, $A/q=1$) performing a three-mode swimming roll maneuver (projection
planes as labeled). At times $t_{1}$ and $t_{3}$ the nematode performs planar
undulations; at $t_{2}$ it performs a roll with normalized torsion
$\tau_{0}L=1.5$ and length $q\Delta s^{\tau}=3.24$. The roll results in a
change of the undulation plane by more than $90^{\circ}$.  Figure 11:  Roll
rotation rate per undulation period $\omega^{\textrm{roll}}$ (in degrees) vs
normalized body torsion $\tau_{0}L$, for normalized curvature-wave amplitude
$A/q=0.8,1,1.2,1.4,1.6,1.8$ (from above) for (a) C-shaped (q$L=5.5$) and (b)
W-shaped (q$L=9$) swimming nematodes.  Figure 12: Roll reorientation angle
$\theta^{\textrm{roll}}$ vs normalized torsional mode length $q\Delta
s^{\tau}$. Panels (a) and (c) show results for C-shaped gait ($qL=5.5$), and
panels (b) and (d) show results for W-shaped ($qL=9$) gait. Normalized
magnitude of body torsion $\tau_{0}L$ as labeled. The initial mode-switching
point occurs when the curvature at the head position is at a maximum [panels
(a) and (b)] and at a minimum [panels (c) and (d)].  Figure 13:  Roll rotation
rate per undulation period $\omega^{\textrm{roll}}$ (in degrees) vs normalized
magnitude of body torsion $\tau_{0}/q$ for nematode performing one-mode rolls
for swimming (solid lines) and burrowing (dashed lines). Results correspond to
$A/q=1$ for C-shaped and W-shaped locomotive gaits, as labeled.  Figure 14:
(a) Front and (b) side views of a trajectory reproduced from sinusoidal
curvature with $A/q=1$ and constant torsion $\tau/q=0.15$. The first period of
undulation is shown in gray and the second in black. The symbol $\odot$
denotes the direction of propagation of the curve (out of the page) and the
green arrows represent the unit tangent vector ${\hat{\boldsymbol{t}}}$ (in
the direction of nematode motion). In the presence of a fluid, the nematode
moving along the curve experiences hydrodynamic resistance force (represented
by orange arrows) oriented in the same direction as the geometrical rotation
of the undulation plane.

###  3.2 Swimming in 3D

Observations of C. elegans swimming in bulk fluids show that the animal is
able to actively explore 3D environments. A swimming nematode can perform a
series of planar turns and is also capable of rapidly changing the plane of
undulations (see Fig. 8 and Video S3 in ESM).

Video images of C. elegans undergoing undulation-plane reorientation indicate
that during such maneuvers the nematode body postures are distorted out of the
undulation plane. We find that these nonplanar postures resemble the
torsional-roll contours defined by curvature and torsion relations (3) and (6)
(see the contours depicted in Fig. 3).

Figure 9 shows a direct comparison between a video image of a swimming
nematode and a calculated torsional-roll posture with a curvature wavevector
$qL=5.5$ (C-shaped gait) and constant torsion $\tau_{0}L=0.9$. We note that
the out-of-plane deformation of the nematode body is quite small for this
torsion value; our hydrodynamic analysis shows, however, that even a moderate
torsional deformation may have a significant dynamic effect, producing rolls
that are significantly stronger than the corresponding rolls in no-slip
burrowing.

To elucidate the effect of the torsional roll in swimming, we present results
of our hydrodynamic simulations for single-mode and three-mode roll maneuvers,
described by Eq. (6) and (9), respectively. The overall undulation-plane
rotation rate for a swimming nematode includes geometrical and hydrodynamic
contributions. The geometrical component results from the rotation of the
Serret–Frenet frame
$({\hat{\boldsymbol{t}}},{\hat{\boldsymbol{n}}},{\hat{\boldsymbol{b}}})$
according to equations (5); the hydrodynamic component is associated with the
rigid-body rotation of the whole nematode body, as required by the zero net
hydrodynamic force and torque conditions.

Figure 10 shows stereoscopic snapshots of the nematode configuration at the
beginning, during, and at the end of a simulated three-mode roll. Another
example of a three-mode roll maneuver is presented in Video S4 in ESM. In both
cases the nematode undergoes a roll rotation of approximately $90^{\circ}$.
Simulated nematode motion in the Video S4 closely resembles the behavior of a
wild-type nematode seen in Video S3.

Quantitative results regarding the effectiveness of roll reorientation
maneuvers in fluids are presented in Figs. 11–13. The roll rotation rate per
undulation period (7) during a single-mode roll is plotted in Fig. 11, and the
total roll reorientation angle

$\theta^{\textrm{roll}}=\cos^{-1}\left[{\hat{\boldsymbol{b}}}(s^{\tau}_{1})\boldsymbol{\cdot}{\hat{\boldsymbol{b}}}(s^{\tau}_{2})\right]$
(10)

for a nematode performing a three-mode roll is plotted in Fig. 12.

The roll rotation rate $\omega^{\textrm{roll}}$ depicted in Fig. 11 is shown
vs the dimensionless torsion amplitude $\tau_{0}L$, for several values of the
normalized curvature amplitude $A/q$, for C-shaped and W-shaped gaits. Note
that in swimming the roll rotation rate depends on the nematode length whereas
in no-slip burrowing $\omega^{\textrm{roll}}$ depends only on the torsion and
curvature-wave parameters. Hence a different normalization is used in Figs. 11
(swimming) and 4 (burrowing).

In swimming, similar to the results for the no-slip roll (Fig. 4), the roll
rotation rate is an increasing function of the normalized curvature-wave
amplitude $A/q$ and body torsion $\tau_{0}L$, but the dependence on the
torsion is nonliner due to hydrodynamic effects. The results for C-shaped and
W-shaped gaits have a similar magnitude for given values of $A/q$ and
$\tau_{0}L$.

The total reorientation angle $\theta^{\textrm{roll}}$ during a three-mode
roll maneuver is plotted in Fig. 12 vs the dimensionless torsional mode length
$q\Delta s^{\tau}$, normalized by the wavevector of the curvature undulations.
The results are similar for the C-shaped and W-shaped gaits and are only
moderately affected by the initial mode switching point. The magnitude of the
reorientation angle $\theta^{\textrm{roll}}$ is consistent with the single-
mode roll rotation rate presented in Fig. 11.

Our experimental observations of C. elegans swimming maneuvers indicate that
the nematode typically needs one to two periods of motion to roll by
$90^{\circ}$, and rolls by $90^{\circ}$ in less than one period are not
unusual. According to the results shown in Figs. 11 and 12, a $90^{\circ}$
roll corresponds to the body torsion of $\tau_{0}L\approx 1$ for the roll mode
length ranging from half to one undulation period. Shorter modes require
larger torsion (see, e.g. Video S4 in ESM).

####  3.2.1 Hydrodynamic enhancement in swimming rolls

To compare the swimming roll rotation rate $\omega^{\textrm{roll}}$ shown in
Fig. 11 with the corresponding purely geometrical no-slip result (Fig. 4), the
swimming and no-slip simulation data are replotted in Fig. 13 using the same
scaling of the torsion $\tau_{0}/q$ for both data sets. The results indicate
that hydrodynamic forces in swimming can enhance the rotation rate by a factor
ranging from approximately two (C-shaped gait) to even three (W-shaped gait).
While this hydrodynamic enhancement may seem surprising, its origin can be
explained by a careful examination of hydrodynamic forces acting on a nematode
performing a torsional roll in the absence of sidewise slip.

The projection of a torsional no-slip trail onto the plane normal to the
overall direction of motion is depicted in Fig. 14. The graph shows that the
plane of undulations gradually rotates anticlockwise, whereas the nematode
moves clockwise around the axis of the trajectory. In fluids, this clockwise
motion produces the oppositely oriented hydrodynamic resistance, which causes
hydrodynamically induced nematode rotation in the anticlockwise direction.
Since the geometrical and hydrodynamic contributions to
$\omega^{\textrm{roll}}$ have the same sign, they add up, resulting in the
enhanced reorientation.

We believe that swimming C. elegans use the above hydrodynamic enhancement
mechanism in the observed rapid roll maneuvers (see Video S3 in ESM). Our
numerical simulations presented in Video S4 in ESM reproduce this behavior.

Figure 15:  (Colour online) (a) Top view of the simulated chemotaxis assay in 3D. Nematodes are initialized with a random orientation at the green symbol $\times$ located midway between the position of a Gaussian chemoattractant concentration peak (right) and a spherical control region represented by the left red circle. The right blue circle represents a spherical test region. (b) Two views of a simulated path followed by a chemotaxing C. elegans in a 3D fluid.  Figure 16:  Chemotaxis index CI vs time for burrowing and swimming nematodes in 3D. The sensing parameters for the burrowing and swimming cases are the same, but the undulation frequency of the swimmer is four times larger than the undulation frequency of the burrower.  Parameter | Name | Value
---|---|---
$\alpha_{s}^{\prime}$ | Chemoattractant sensitivity parameter | $5\,\textrm{s}/\textrm{mm}$
$r_{s}^{0}$ | Turn initiation rate | $10\textrm{ turns}/\textrm{min}$
$r_{s}^{H}$ | Upper turn-rate saturation limit | $15\textrm{ turns}/\textrm{min}$
$r_{s}^{L}$ | Lower turn-rate saturation limit | $5\textrm{ turns}/\textrm{min}$
$r_{\tau}^{0}$ | Roll initiation rate | $37\textrm{ rolls}/\textrm{min}$
$\epsilon$ | Memory time constant | $1.6\,\textrm{s}^{-1}$
$A_{1}/q$ | Forward mode normalized amplitude | 1
$f_{b}$ | Burrowing gait frequency | $0.44\,\textrm{s}^{-1}$
$f_{s}$ | Swimming gait frequency | $1.75\,\textrm{s}^{-1}$
$q_{b}L$ | Burrowing gait normalized wavevector | $9$
$q_{s}L$ | Swimming gait normalized wavevector | $5.5$
$L$ | Worm length | $1\,\textrm{mm}$
Table 1:  System parameters for 3D simulations of chemotaxis of burrowing and
swimming C. elegans. The chemoattractant sensitivity parameter
$\alpha_{s}^{\prime}=\alpha_{s}/\nabla C|_{r=\sigma}$ is normalized by the
maximum chemical concentration gradient located at $r=\sigma$. The memory time
constant $\epsilon$ corresponds to an approximately $5\,\textrm{s}$ half-decay
response time of sensory neurons to the step increase of chemoattractant
concentration Suzuki  _et al._ (2008).

##  4 Application to Chemotaxis

According to the locomotion model developed in Secs. 2 and 3, C. elegans
navigates 3D environments by combining forward motion, planar turns, and
torsional rolls (see Videos S1 and S5 in ESM for simulations of burrowing and
swimming behavior). Here we show that by modulating the rate of the above
elementary maneuvers in response to a chemoattractant concentration variation,
C. elegans can efficiently chemotax in 3D.

It has been demonstrated that in 2D environments C. elegans employs two
chemotaxis mechanisms: klinotaxis (gradual turn) where the nematode actively
reorients itself towards the concentration gradient direction Iino and Yoshida
(2009), and klinokinesis (biased random walk), where it adjusts time intervals
of motion in favorable and unfavorable directions to achieve a net
displacement towards the higher concentration Pierce-Shimomura  _et al._
(1999).

While we expect both strategies to be effective in 3D, in the present paper we
focus on the biased random walk (BRW) mechanism. A more comprehensive study of
nematode chemotaxis in complex environment will be described elsewhere.

###  4.1 Chemical sensing model

A moving C. elegans assesses the spatial gradients of the chemoattractant
concentration by monitoring its time variation $C(t)$ at the current position
of chemical receptors (located in the nematode head segment) Lockery (2011).
To model this process we use an approach that has been successfully applied to
describe bacterial chemotaxis Celani and Vergassola (2010).

Accordingly, we assume that the rate of turns is modulated by the chemical
signal

$Q(t)=\int_{0}^{t}C(t-t^{\prime})M(t^{\prime})dt^{\prime},$ (11)

where $M(t)$ is the memory function that describes the response of sensory
neurons to chemoattractant. Assuming that the nematode responds to changes of
the concentration of chemoattractant but not to the absolute concentration
value Lockery (2011), we impose the ideal adaptation condition Celani and
Vergassola (2010)

$\int_{0}^{\infty}M(t)dt=0.$ (12)

We also assume the normalization

$\int_{0}^{\infty}tM(t)dt=1.$ (13)

In our numerical calculation we use the memory function defined by the
expression

$M(t)=\epsilon^{2}\mathrm{e}^{\epsilon t}\left[(\epsilon
t)-\frac{1}{2}(\epsilon t)^{2}\right],$ (14)

where $\epsilon^{-1}$ is the decay time Celani and Vergassola (2010). The
memory function (14) is positive for $0<t<\epsilon^{-1}$ and negative for
$t>\epsilon^{-1}$; the positive and negative regions have equal integrals. The
convolution integral (11) filters out short-time fluctuations of the
chemoattractant concentration and provides information whether the
concentration is increasing ($Q>0$) or decreasing ($Q<0$) on the timescale
$\epsilon^{-1}$.

###  4.2 Biased random walk

A nematode performing BRW generates a trajectory drift towards a peak of
chemoattractant concentration by controlling the frequency of random turns in
response to sensory signals. Namely, to achieve longer intervals of motion
towards the peak, the turn rate is reduced (increased) when the animal moves
towards (away from) the peak. In our chemotaxis model, the chemical signal $Q$
controls the frequency of elementary three-mode turns.

The nematode burrows or swims using the default harmonic-curvature locomotion
mode with $A/q\approx 1$, and switches randomly to a turning mode at a turn
initiation rate $r_{s}(t)$, which depends on the chemotaxis signal (11). The
nematode also performs roll maneuvers initiated at a rate $r_{\tau}(t)$.

To mimic the experimentally observed sigmoidal character of the turn response
to the chemoattractant concentration variation Pierce-Shimomura  _et al._
(1999), we assume that the turn initiation rate $r_{s}$ has a linear-response
region and a saturation cutoff,

$r_{s}(t)=\left\\{\begin{matrix}r_{s}^{L};\hskip 11.38109ptf(Q)\leq
r_{s}^{L}\\\ f(Q);\hskip 11.38109ptr_{s}^{L}<f(Q)<r_{s}^{H},\\\
r_{s}^{H};\hskip 11.38109ptr_{s}^{H}\leq f(Q)\\\ \end{matrix}\right.$ (15)

where

$f(Q)=r_{s}^{0}\left[1-\alpha_{s}Q(t)\right].$ (16)

Here $r_{s}^{0}$ is the turn rate in the absence of chemoattractant,
$\alpha_{s}$ is a chemoattractant sensitivity parameter, and $r_{s}^{H}$ and
$r_{s}^{L}$ are the saturation limits. Since there are no experimental data
regarding the effect of chemoattractant on roll maneuvers, we assume that the
roll initiation frequency is constant and independent of chemical sensing,

$r_{\tau}(t)=r_{\tau}^{0}.$ (17)

In Sec. 4.4 we show that assumptions (15) and (17) are sufficient to produce
efficient 3D chemotaxis.

###  4.3 The geometry of simulated chemotaxis assay

In our 3D numerical simulations of chemotaxis of burrowing and swimming
nematodes [see Fig. 15(a)], the simulated worms are placed at a designated
starting point at the edge of a Gaussian chemoattractant concentration
distribution. The chemotaxis efficiency is determined by comparing the number
of nematodes that reach a spherical test region at the chemoattractant
concentration peak with the number of those that reach a symmetrically placed
control region in the low-concentration domain. Our 3D arrangement mimics a
typical experimental setup for investigating 2D chemotaxis on agar surface
Bargmann and Horvitz (1991); Nishino  _et al._ (2015).

The chemotaxis efficiency is quantified using the chemotaxis index

$\textit{CI}(t)=\frac{n_{t}(t)-n_{c}(t)}{n_{t}(t)+n_{c}(t)},$ (18)

where $n_{t}(t)$ and $n_{c}(t)$ denote the number of nematodes that reached
the test and control regions, respectively, during the simulation interval
$[0,t]$. The simulation of a nematode trajectory is terminated when the
nematode enters either the test or control region (in 2D laboratory
experiments these regions usually contain anesthetic to immobilize the
nematodes Nishino  _et al._ (2015)). A sample trajectory of a swimming C.
elegans reaching the test region is shown in Fig. 15(b).

In our simulations we consider a system with the following parameters. The
spatial concentration of chemoattractant is given by a 3D Gaussian
distribution

$C\sim\exp\left(-r^{2}/2\sigma^{2}\right)$ (19)

with the peak position at a distance $W=36L$ from the starting point and the
width $\sigma=0.4W$. The diameter of the test and control regions is $D=10L$.
For the worm length $L=1\,\textrm{mm}$, the starting point is thus at a
distance $3.6\,\textrm{cm}$ from the chemoattractant concentration maximum,
and the test and control regions have a $1\,\textrm{cm}$ diameter (similar to
the geometry of 2D experiments in a Petri dish).

###  4.4 Simulation results

We have performed simulations for two cases: (a) burrowing W-shaped nematodes
(no transverse slip) and (b) swimming C-shaped nematodes. The nematode gait
parameters and parameters of the chemotaxis model are listed in Table 1 (see
Sec. 6.2 for more details). The gait parameter values are based on
observations of crawling and swimming C. elegans. We use the same
chemosensation and turn-rate parameters in the burrowing and swimming
simulations to isolate the effect of the locomotion mechanics on the
chemotaxis efficiency.

Chemotaxis index (18), shown in Fig. 16, is significantly greater than zero,
indicating a clear trajectory bias towards the chemoattractant concentration
peak. At long times the value of CI for both burrowing and swimming is
approximately the same. However, the burrowing nematodes are able to reach the
chemoattractant peak more quickly, due to the higher burrowing velocity (twice
the swimming velocity for the gait parameters considered).

The chemotaxis index for burrowing and swimming nematodes is similar due to
two compensating effects: (a) for swimming nematodes the chemical signal (11)
is weaker because, as a result of their lower velocity, swimming worms detect
a smaller time variation of chemoattractant concentration; while (b) the
number of turns per unit trajectory length is larger for swimming nematodes
because of their lower velocity. More turns implies a stronger bias towards
chemoattractant, which compensates for a weaker chemical signal.

We expect that similar compensating factors occur in other environments as
well. Thus, the nematode may achieve efficient chemotaxis without adjusting
turn-rate control parameters in response to a variation in mechanical
properties of the medium in which it moves; a more complex control system that
relies on mechanical feedback is unnecessary.

The values of the 2D and 3D chemotaxis index do not differ either, but in 3D
many more nematodes escape without reaching either the test or control region;
this effect should be taken into account in design of a testing platform for
future 3D chemotaxis experiments.

##  5 Summary and discussion

To provide mathematical tools for analysis of 3D nematode maneuvers, we have
developed a gait model that combines piecewise-harmonic-curvature (PHC) and
piecewise-constant-torsion (PCT) body-shape assumptions. We have demonstrated
that nematode reorientation in an arbitrary direction can be achieved by a mix
of 2D turns (driven by changes of curvature-wave parameters) and 3D roll
maneuvers (associated with nonzero body torsion).

We have applied our model to investigate 3D motion of nematodes burrowing in
gels and swimming in Newtonian fluids. We examined the effectiveness of the
turns and rolls; our modeling of torsional roll maneuvers reveals a strong
enhancement of the roll-rotation rate due to hydrodynamic forces. This
hydrodynamic effect may explain the observed rapid reorientation of the
undulation plane of swimming C. elegans (see video S3 in ESM).

We have used our gait representation to study nematode chemotaxis in gels and
fluids. We have shown that the nematode does not need to adjust its sensory-
motor apparatus to effectively respond to chemical stimuli in environments
with different mechanical properties. Chemotaxis of C. elegans in 3D
environments was experimentally demonstrated Pierce-Shimomura  _et al._
(2008); Beron  _et al._ (2015), but, to our knowledge, the formulation
presented here is the first quantitative analysis that takes into account the
gait geometry and swimming mechanics.

By supplying an analytical description of elementary body movements, our
PHC/PCT gait model is applicable to analyze nematode behavior in 3D complex
media, including inhomogeneous materials and non-Newtonian fluids. We expect
that our results may also shed light on other subtle aspects of neuromuscular
control of 3D locomotion of C. elegans.

Body movements of C. elegans are actuated by coordinated action of four
quadrants of body muscles, which are arranged into double rows that span the
entire body length Altun and Hall (2009). While this muscle anatomy is
compatible with 3D body actuation, it is less clear how the neural system can
activate differential muscle contractions needed to generate 3D body
deformations out of the dorsoventral undulation plane.

The head and neck body muscles are innervated by the nerve ring motor neurons,
which synapse onto cells in individual quadrants (and thus are capable to
actuate each quadrant individually). In contrast, the remaining body muscles
are innervated only by the ventral-cord motor neurons, which synapse onto
muscles of either two ventral or two dorsal quadrants, consistent with 2D
actuation White  _et al._ (1986); Altun and Hall (2009); Riddle  _et al._
(1997). Based on this topology of neural connections, it is generally assumed
that only the head and neck segments can perform 3D motion.

However, the stereoscopic experiments by Kwon et al. Kwon  _et al._ (2013)
clearly show that burrowing nematodes can adopt a variety of 3D postures, in
which the entire body undergoes out-of-plane deformations. In the swimming
case the evidence is less conclusive because stereoscopic images are not
available. Nonetheless, two-dimensional images, such as the one shown in Fig.
9, agree qualitatively with our 3D body posture model.

Our hydrodynamic calculations of the roll dynamics further support a
conclusion that the entire nematode body participates in generating the roll
motion. In particular, we find that torsional wave propagating throughout the
entire body length can produce a 90∘ roll in less than one undulation period.
In contrast, a stroke in which only the head and neck segments are displaced
out of the undulation plane produces weaker rolls than the ones observed in
our experiments (see the results of our additional simulations in Fig. S2 in
ESM).

We hypothesize that the generation of a nonzero torsion in the posterior part
of the nematode body involves a 3D proprioceptive feedback, which enables
propagation of the torsional wave to the body region innervated only by the
ventral-cord motor neurons. The existence of 2D proprioceptive feedback in C.
elegans and its important role in generation of a planar undulation wave have
been demonstrated Wen  _et al._ (2012). We thus conjecture that C. elegans may
employ 3D proprioception to generate well controlled out-of-plane body
movements. (An alternative hypothesis would be that the worm may control its
3D body movements using direct communication between muscle cells via gap
junctions.)

Current evidence is not sufficient to rule out a passive, purely mechanical
origin of propagation of 3D deformations along a portion of the nematode body
outside the head and neck regions. In particular, the body shapes observed
during burrowing Kwon  _et al._ (2013) can, in principle, be generated by 3D
actuation of the head and neck muscles, with the rest of the body passively
sliding along the tunnel carved by the nematode head. Such motion, however,
would be inefficient.

According to Wen et al. Wen  _et al._ (2012), “[t]he cellular economy of the
C. elegans wiring diagram implies that individual neurons may have high levels
of complexity” (for example, B-type motor neurons are able to transduce
proprioception). We conclude that it is thus conceivable that functional
complexity on a single-neuron level allows symmetrically connected neurons to
produce differential muscle excitation.

Nature created a four-quadrant muscle arrangement that is anatomically capable
of generating 3D body deformations. Therefore, our question is: why should
there not be a compatible neuromuscular-control system that would allow the
nematode to take advantage of this (existing) muscle structure? Further
analyses of 3D locomotion of C. elegans will help determine if such a control
system exists.

##  6 Methods

###  6.1 Imaging

Wild type C. elegans was used to record the 3D swimming in our experiments.
Worms were cultured in NGM plate (Nematode Growth Medium) seeded with bacteria
E. coli OP50 at 20∘C. A liquid pool was created for swimming assay by adding
5–7 mL of M9 buffer into a 60 mm Petri dish. The pool height was approximately
2 mm. The surface of the Petri dish was treated with 5 % Pluoronic F127 for 5
min before adding M9 buffer. Individual young adults were manually transferred
into a liquid pool of M9 buffer using a worm pick. Swimming episodes were
recorded at 5 frames per second using an SVSI camera and Zeiss stemi 2000-C
stereo microscope imaging system at 1.6X magnification with a pixel resolution
of 0.21 pixels/$\mu$m. All imaging was carried out in a food-free environment
at 20∘C.

###  6.2 Chemotaxis simulation details

In our simulations, the planar turns and rolls are initiated randomly at the
rates $r_{s}$ and $r_{\tau}$, provided that a reorientation maneuver of the
same kind is not already taking place. The turns and rolls are performed
independently.

The default gait parameters (used between reorientation maneuvers) have fixed
values listed in Table 1. The parameters of the turn and roll modes are chosen
randomly from Gaussian distributions. The average parameter values (indicated
by overbar) and standard deviations $\sigma$ are: normalized turning-mode
amplitude $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{}_{2}/q=1.67$ and $\sigma_{A_{2}/q}=0.28$; turning-mode length
$\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{}^{c}/L=0.5$ and $\sigma_{\Delta s^{c}/L}=0.167$; normalized torsion
$\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{}_{0}L=2$ and $\sigma_{\tau_{0}L}=1$; and torsional-mode length
$\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{}^{\tau}/L=0.5$ and $\sigma_{\Delta s^{\tau}/L}=0.1$. The parameter
distributions are the same for burrowing and swimming.

The turning rates listed in Table 1 correspond to the directional persistence
length of burrowing trajectories $l_{p}\approx 12L$ (which is of the same
order as the persistence length observed in crawling). We define the
persistence length $l_{p}$ as half-decay length of the directional correlation
function. The rate of roll maneuvers was chosen such that on average the worm
spends half the time in a torsionless state.

##### Authors’ contributions.

AB wrote the manuscript, developed locomotion and hydrodynamic codes, and
performed locomotion and chemotaxis simulations. AP provided chemotaxis
formulations and codes and performed some chemotaxis simulations. MR provided
experimental images of the worm. SAV provided experimental images and helped
develop conceptual framework and draft the manuscript. JB conceived the study
and wrote the manuscript.

##### Funding.

AB would like to acknowledge financial support form NSF Grant CBET 1059745 and
TTU Doctoral Dissertation Completion Fellowship. JB was partially supported
from NSF Grant CBET 1603627\. SAV and MR acknowledge partial support from
Grants NIH 1R21AG050503-01 and NASA NNX15AL16G.

##### Acknowledgments.

We would like to thank our visiting undergraduate researcher, Jose Montoya,
for his help in developing chemotaxis codes and acknowledge Dr. Szewczyk and
Dr. Driscoll for useful discussions.

##### Competing interests.

We have no competing interests.

## References

  * Ghanbari  _et al._ (2012) A. Ghanbari, V. Nock, S. Johari, R. Blaikie, X. Chen, and W. Wang, “A micropillar-based on-chip system for continuous force measurement of C. elegans,” J. Micromech. Microeng. 22, 095009 (2012).
  * Johari  _et al._ (2013) S. Johari, V. Nock, M. M. Alkaisi, and W. Wang, “On-chip analysis of C. elegans muscular forces and locomotion patterns in microstructured environments,” Lab. Chip. 13, 1699–1707 (2013).
  * Etheridge  _et al._ (2015) T. Etheridge, M. Rahman, C. J. Gaffney, D. Shaw, F. Shephard, J. Magudia, D. E. Solomon, T. Milne, J. Blawzdziewicz, D. Constantin-Teodosiu, P. L. Greenhaff, S. A. Vanapalli, and N. J. Szewczyk, “The integrin-adhesome is required to maintain muscle structure, mitochondrial ATP production, and movement forces in Caenorhabditis elegans,” FASEB J 29, 1235–1246 (2015).
  * Wen  _et al._ (2012) Q. Wen, M. D. Po, E. Hulme, S. Chen, X. Liu, S. W. Kwok, M. Gershow, A. M. Leifer, V. Butler, C. Fang-Yen, T. Kawano, W. R. Schafer, G. Whitesides, M. Wyart, D. B. Chklovskii, M. Zhen, and A. D. T. Samuel, “Proprioceptive Coupling within Motor Neurons Drives C. elegans Forward Locomotion,” Neuron 76, 750–761 (2012).
  * Stirman  _et al._ (2011) J. N. Stirman, M. M. Crane, S. J. Husson, S. Wabnig, C. Schultheis, A. Gottschalk, and H. Lu, “Real-time multimodal optical control of neurons and muscles in freely behaving Caenorhabditis elegans,” Nat. Methods. 8, 153–159 (2011).
  * Gray  _et al._ (2005) J. M. Gray, J. J. Hill, and C. I. Bargmann, “A circuit for navigation in Caenorhabditis elegans,” Proc. Natl. Acad. Sci. U.S.A. 102, 3184–3191 (2005).
  * Kocabas  _et al._ (2012) A. Kocabas, C.-H. Shen, Z. V. Guo, and S. Ramanathan, “Controlling interneuron activity in Caenorhabditis elegans to evoke chemotactic behaviour,” Nature 490, 273–277 (2012).
  * Sznitman  _et al._ (2010a) J. Sznitman, P. K. Purohit, P. Krajacic, T. Lamitina, and P. E. Arratia, “Material properties of Caenorhabditis elegans swimming at low Reynolds number,” Biophys. J. 98, 617–626 (2010a).
  * Brown  _et al._ (2013) A. E. X. Brown, E. I. Yemini, L. J. Grundy, T. Jucikas, and William R. Schafer, “A dictionary of behavioral motifs reveals clusters of genes affecting Caenorhabditis elegans locomotion,” Proc. Natl. Acad. Sci. U.S.A. 110, 791–796 (2013).
  * Restif  _et al._ (2014) C. Restif, C. Ibanez-Ventoso, M. M. Vora, S. Guo, D. Metaxas, and M. Driscoll, “CeleST: Computer Vision Software for Quantitative Analysis of C. elegans Swim Behavior Reveals Novel Features of Locomotion,” PLoS Comput. Biol. 10, 1–12 (2014).
  * Artal-Sanz  _et al._ (2011) M. Artal-Sanz, L. de Jong, and N. Tavernarakis, “Caenorhabditis elegans: A versatile platform for drug discovery,” J. Biotechnol. 8, 1405–1418 (2011).
  * Omura  _et al._ (2012) D. T. Omura, D. A. Clark, A. D. T. Samuel, and H. R. Horvitz, “Dopamine signaling is essential for precise rates of locomotion by C. elegans,” PLoS ONE 7, e38649 (2012).
  * Cohen and Boyle (2010) N. Cohen and J. H. Boyle, “Swimming at low Reynolds number: a beginners guide to undulatory locomotion,” Contemp. Phys. 51, 103–123 (2010).
  * Fang-Yen  _et al._ (2010) C. Fang-Yen, M. Wyart, J. Xie, R. Kawai, T. Kodger, S. Chen, Q. Wen, and A. D. T. Samuel, “Biomechanical analysis of gait adaptation in the nematode Caenorhabditis elegans,” Proc. Natl. Acad. Sci. U.S.A. 107, 20323–8 (2010).
  * Sznitman  _et al._ (2010b) J. Sznitman, X. Shen, P. K. Purohit, and P. E. Arratia, “The effects of fluid viscosity on the kinematics and material properties of C. elegans swimming at low Reynolds number,” Exp. Mech. 50, 1303–1311 (2010b).
  * Padmanabhan  _et al._ (2012) V. Padmanabhan, Z. S. Khan, D. E. Solomon, A. Armstrong, K. P. Rumbaugh, S. A. Vanapalli, and J. Blawzdziewicz, “Locomotion of C. elegans: a piecewise-harmonic curvature representation of nematode behavior,” PLoS ONE 7, e40121 (2012).
  * Stephens  _et al._ (2008) G. J. Stephens, B. Johnson-Kerner, W. Bialek, and W. S. Ryu, “Dimensionality and dynamics in the Behavior of C. elegans,” PLoS Comput. Biol. 4, e1000028 (2008).
  * Kwon  _et al._ (2013) N. Kwon, J. Pyo, S. Lee, and J. Je, “3-D Worm Tracker for Freely Moving C. elegans,” PLoS ONE 8, e57484 (2013).
  * Kwon  _et al._ (2015) N. Kwon, A. B. Hwang, Y.-J. You, S.-J. V. Lee, and J. H. Je, “Dissection of C. elegans behavioral genetics in 3-D environments,” Sci. Rep. 5, 09564 (2015).
  * Beron  _et al._ (2015) C. Beron, A. G. Vidal-Gadea, J. Cohn, A. Parikh, G. Hwang, and J. T. Pierce-Shimomura, “The burrowing behavior of the nematode Caenorhabditis elegans: a new assay for the study of neuromuscular disorders,” G2B 14, 357–368 (2015).
  * Deng and Xu (2014) X. Deng and J. Xu, “A 3d undulatory locomotion model inspired by C. elegans through DNN approach,” Neurocomputing 131, 248–264 (2014).
  * Delmotte  _et al._ (2015) B. Delmotte, E. Climent, and F. Plouraboue, “A general formulation of bead models applied to flexible fibers and active filaments at low Reynolds number,” J. Comput. Phys. 286, 14–37 (2015).
  * Sznitman  _et al._ (2010c) J. Sznitman, X. Shen, R. Sznitman, and P. E. Arratia, “Propulsive force measurements and flow behavior of undulatory swimmers at low Reynolds number,” Phys. Fluids 22, 12901 (2010c).
  * (24) S. J. Husson, W. S. Costa, C. Schmitt, and A. Gottschalk, “Keeping track of worm trackers (September 10, 2012),”  _The C. elegans Research Community_ , WormBook, doi/10.1895/wormbook.1.156.1, http://www.wormbook.org.
  * Berman  _et al._ (2013) R. S. Berman, O. Kenneth, J. Sznitman, and A. M. Leshansky, “Undulatory locomotion of finite filaments: lessons from C. elegans,” New J. Phys. 15, 075022 (2013).
  * Bilbao  _et al._ (2013) A. Bilbao, E. Wajnryb, S. A. Vanapalli, and J. Blawzdziewicz, “Nematode Locomotion in Unconfined and Confined Fluids,” Phys. Fluids. 25, 081902 (2013).
  * Kim  _et al._ (2011) D. Kim, S. Park, L. Mahadevan, and J. H. Shin, “The shallow turn of a worm,” J. Exp. Biol. 214, 1554–1559 (2011).
  * Stoker (1969) J. Stoker,  _Differential Geometry_ (John Wiley and Sons, New York, 1969).
  * Voorhies  _et al._ (2005) W Voorhies, J. Fuchs, and S. Thomas, “The longevity of Caenorhabditis elegans in soil,” Biol. Lett. , 247–249 (2005).
  * Frezal and Felix (2015) L. Frezal and M. Felix, “The natural history of model organisms: C. elegans outside the petri dish,” eLife 4, e05849 (2015).
  * Lee  _et al._ (2012) H. Lee, M.-K. Choi, D. Lee, H.-S. Kim, H. Hwang, H. Kim, S. Park, Y.-K. Paik, and J. Lee, “Nictation, a dispersal behavior of the nematode Caenorhabditis elegans, is regulated by IL2 neurons,” Nature Neurosci 15, 107–112 (2012).
  * Wajnryb  _et al._ (2013) E. Wajnryb, K. A. Mizerski, P. J. Zuk, and P. Szymczak, “Generalization of the Rotne-Prager-Yamakawa mobility and shear disturbance tensors,” J. Fluid Mech. 731, R3 (2013).
  * Cichocki  _et al._ (1994) B. Cichocki, B. U. Felderhof, K. Hinsen, E. Wajnryb, and J. Bławzdziewicz, “Friction and mobility of many spheres in Stokes flow,” J. Chem. Phys. 100, 3780–3790 (1994).
  * Suzuki  _et al._ (2008) H. Suzuki, T. R. Thiele, S. Faumont, M. Ezcurra, S. R. Lockery, and W. R. Schafer, “Functional asymmetry in Caenorhabditis elegans taste neurons and its computational role in chemotaxis,” Nature 454, 114–118 (2008).
  * Iino and Yoshida (2009) Y. Iino and K. Yoshida, “Parallel Use of Two Behavioral Mechanisms for Chemotaxis in Caenorhabditis elegans,” J. Neurosci. 29, 5370–5380 (2009).
  * Pierce-Shimomura  _et al._ (1999) J. T. Pierce-Shimomura, T. M. Morse, and S. R. Lockery, “The fundamental role of pirouettes in Caenorhabditis elegans chemotaxis,” J. Neurosci. 19, 9557–9569 (1999).
  * Lockery (2011) S. R. Lockery, “The computational worm: spatial orientation and its neuronal basis in C. elegans,” Curr Opin Neurobiol 21, 782–790 (2011).
  * Celani and Vergassola (2010) A. Celani and M. Vergassola, “Bacterial strategies for chemotaxis response,” Proc. Natl. Acad. Sci. U.S.A. 107, 1391–1396 (2010).
  * Bargmann and Horvitz (1991) C. I. Bargmann and H. R. Horvitz, “Chemosensory neurons with overlapping functions direct chemotaxis to multiple chemicals in c. elegans,” Neuron. 7, 729–742 (1991).
  * Nishino  _et al._ (2015) A. Nishino, F. Sato, K. Ito, and T. Matsuura, “Extension of the established period of diacetyl adaptation by oxygen intermediates in the nematode Caenorhabditis elegans,” Comp. Biochem. Physiol. A: Physiol. 184, 156–162 (2015).
  * Pierce-Shimomura  _et al._ (2008) J. T. Pierce-Shimomura, B. L. Chen, J. J. Mun, R. Ho, R. Sarkis, and S. L. McIntire, “Genetic analysis of crawling and swimming locomotory patterns in C. elegans,” Proc. Natl. Acad. Sci. U.S.A. 105, 20982–20987 (2008).
  * Altun and Hall (2009) Z. F. Altun and D. H. Hall, “Muscle system, somatic muscle.” In WormAtlas (2009).
  * White  _et al._ (1986) J. G. White, E. Southgate, J. N. Thomson, and S. Brenner, “The structure of the nervous system of the nematode Caenorhabditis elegans,” Philos. Trans. R. Soc. London, Ser. B 314, 1–340 (1986).
  * Riddle  _et al._ (1997) D. L. Riddle, T. Blumenthal, B. J. Meyer, and J. R. Priess,  _C. elegans II. 2nd edition._ (Spring Harbor, 1997) section I, Introduction: the Neural Circuit For Locomotion.


# Vertex Reconstruction with MaskFormers

Samuel Van Stroud1, Nikita Pond1, Max Hart1, Jackson Barr1,2,
Sébastien Rettie1,3, Gabriel Facini1, and Tim Scanlon1
1University College London, 2DESY, 3CERN

###### Abstract

In high-energy particle collisions, secondary decays can be reconstructed as
displaced vertices using the measured trajectories of charged particles. Such
vertices are useful in identifying and studying jets originating from $b$\- or
$c$-hadrons, which is a key component of the physics programs of modern
collider experiments. While machine learning has become mainstream in particle
physics, most applications are on an per-object basis, for example the
prediction of class labels or the regression of object properties. However,
vertex reconstruction is a many-to-many problem, in which a set of input
tracks must be grouped into a second variable length set of vertices. In this
work, we propose a fully learned approach to reconstruct secondary vertices
inside jets based on recent advancements in object detection from computer
vision. We demonstrate and discuss the advantages of this approach, in
particular its ability to estimate the properties of any number of vertices,
and conclude that the same methodology could be applicable to other
reconstruction tasks in particle physics.

_Keywords_ Particle Physics $\cdot$ Reconstruction $\cdot$ Machine Learning
$\cdot$ Object Detection $\cdot$ Vertexing $\cdot$ Flavour Tagging

##  1 Introduction

Facilities like the Large Hadron Collider (LHC) at CERN [1] have deepened our
understanding of fundamental physics through high-energy proton–proton ($pp$)
collisions. These collisions, observed in detectors like ATLAS [2] and CMS
[3], often result in the creation of particles which decay into conical sprays
of particles, known as jets. A key process in categorising collisions is jet
tagging, especially flavour tagging, which is the identification of jets
originating from heavy flavour hadrons ($b$\- and $c$-hadrons) [4, 5].
Identification mostly relies on charged particle trajectories, or tracks,
reconstructed using dedicated algorithms. Decays from $b\text{-hadrons}$
typically produce four to five charged particles which must be distinguished
from the roughly $10$–$100$ total particles within jets at the energies
typically probed at the LHC.

A vertex is a point where particles interact or decay. Traditionally, vertex
reconstruction happens in two stages. Firstly, vertex finding groups tracks
which have a compatible spatial origin. Secondly, vertex fitting deduces
properties of the vertex from the set of compatible tracks. Weakly decaying
heavy flavour hadrons have non-negligible lifetimes [6], resulting in
secondary decays which may be displaced from the primary $pp$ collision. These
decay points can be reconstructed as secondary vertices. Heavy flavour jets
typically contain one to three such vertices, while light-jets typically
contain displaced vertices at a lower rate and with a lower mass. As a
consequence, the reconstruction of vertices within jets is a key component of
flavour tagging.

Machine learning (ML) has been used in particle physics for several decades
[7]. In recent years, the application of deep learning has become widespread
and has led to dramatic performance improvements. Existing state-of-the-art
flavour tagging algorithms now fully embrace these methods and offer
exceptional jet tagging performance. However, they either ignore the task of
vertex reconstruction altogether [8, 9], or perform a limited form of vertex
finding [10, 11], and no vertex fitting.

The main contribution of this article is to enable the full reconstruction
(both finding and fitting) of multiple vertices inside a jet using a learned
approach. In this article, we focus on reconstructing secondary vertices
produced by heavy flavour hadron decays, but the approach could easily be
extended to additionally reconstruct other vertices inside jets. To achieve
this, we propose a novel approach based on MaskFormer [12], which was
originally designed for image segmentation and is based on the scalable
transformer architecture [13]. Our implementation of this model processes
unordered sets of input tracks and outputs a variable number of fully
characterised vertices. Each output vertex is characterised by a predicted
class label, a series of fitted properties, and an assignment of tracks, the
latter of which is provided in the form of a series of binary masks. In
addition, the model is fused with the existing state-of-the-art ML-based
flavour tagging algorithms similar to those developed by ATLAS [10]. By doing
this, we show that the capacity to fully reconstruct heavy flavour vertices
using MaskFormer improves the model’s ability to identify jets originating
from the decay of heavy hadrons.

The article is structured as follows. Section 2 contains an overview of
related work. The model architecture is described in Section 3, while Section
4 provides details about the simulated samples used and training specifics.
Results are shown in Section 5 and conclusions are given in Section 6.

##  2 Related Work

Previous work in both high-energy physics and computer vision lays a robust
foundation for our contributions. We draw from the latest advancements in
object detection and segmentation in computer vision to address the prevailing
challenges in vertex reconstruction in high-energy physics, paving the way for
more powerful physics analyses.

#### High Energy Physics

High-energy physics experiments, such as those conducted at ATLAS and CMS,
heavily rely on vertex reconstruction algorithms for analysing collision
events [14]. Traditional secondary vertex reconstruction algorithms [15, 16,
17] focus on iteratively clustering tracks around potential common spatial
points of origin, and then fitting the vertex to extract its properties. While
effective, these approaches demand extensive manual optimisation for new
environments and may not perform optimally across wide regions of phase space,
for example across the wide range of jet transverse momenta
(${{{{p_{\mathrm{T}}}}}}$) produced in modern collider experiments. Moreover,
their need for manual fine-tuning presents operational challenges.

The integration of ML techniques mark a significant shift in this domain.
Existing ML-based vertexing approaches use edge classification on input graphs
to predict whether each pair of tracks in the jet belong to a common vertex
[18, 10]. Compatible track pairs are used to form vertices by extracting the
connected components of the resulting graph in a post-processing step. As the
formation of multi-track vertices is not directly learned, this approach may
be suboptimal. Additionally, these models cannot be easily used to estimate
the properties (e.g. position and momentum) of an unknown number of vertices.
Other recent work [19] focuses on performing a differentiable vertex fit on a
single inclusive secondary vertex, but the procedure is currently unable to
reconstruct multiple vertices. We overcome these limitations by directly
predicting the presence and properties of an arbitrary number of vertices in a
single step using the MaskFormer architecture.

#### Computer Vision

In computer vision, object detection is the task of identifying objects within
an image and determining their properties. Object detection combines semantic
segmentation, the assignment of a semantic class label to each pixel, and
instance segmentation, which differentiates between individual instances
within the same semantic class. Techniques in this field have evolved to
predict object classes alongside either bounding boxes [20, 21, 22] and/or
pixel masks [23, 12]. The use of bounding boxes relies on the dense grid
structure of images and the implicit ordering of their pixels. However, per-
pixel classifiers can readily be extended for sparse, unordered collections of
inputs.

The Segment Anything model [24] is a state-of-the-art object detector based on
a MaskFormer architecture [12, 25]. This architecture takes in an image and
outputs a series of detected objects. Each output object consists of a class
label and a binary pixel mask. The object classes are semantic labels for
objects in images (e.g. “dog” or “cat”), while the pixel masks are used to
associate a set of pixels in the image to each labelled object. Multiple
instances of the same class are resolved as different masks with common class
labels. This concept of resolving multiple object instances of the same class
is particularly relevant to our work.

Our research draws inspiration from these advancements in computer vision. We
analogise a grid of image pixels to an unordered set of tracks within a jet,
and image-depicted objects to heavy flavour hadron decay vertices. This
correspondence allows us to adapt computer vision techniques for particle
physics, enabling us to predict both the presence and properties of multiple
vertices in a novel and effective manner.

##  3 Model Architecture

Our approach is based on the MaskFormer architecture as implemented in Refs.
[12, 25, 24]. We include modifications which tailor the model to handle sparse
inputs and additional regression tasks. In this approach, input tracks are
first processed by a transformer encoder backbone. In the object decoder a set
of object queries, each representing a possible output vertex, exchange
information with the tracks and each other. Finally, several task heads are
used to produce the model outputs. For each object query, these are: class
predictions query, vertex properties (via regressions) and masks, the latter
of which provide the assignment of tracks to vertices. Additionally, the model
is configured to output a prediction of the jet flavour. The architecture is
depicted in Fig. 1 and described in more detail below.

Figure 1:  Overview of the MaskFormer network architecture. $M$ input elements
are fed into an initial transformer encoder. The object decoder then takes a
set of $N$ object queries and iteratively updates them with information from
the input elements and other object queries. Finally, three task heads are
used to: categorise each object as being from one of $C+1$ object classes
(including a null class), $R$ regression targets, and $N\times M$ binary masks
which provide the assignment of input tracks to output vertices. The model can
also be configured to predict a global class label, which is not shown here
for simplicity.

###  3.1 Input Encoder

The first stage of the model is an initial dense embedding layer which
produces a $d=256$ dimension vector representation for each of the $M$ input
tracks. This is followed by a sinusoidal positional encoding [13] of track
impact parameters and angular displacements from the jet axis. The positional
encoding allows the self-attention mechanism to readily identify nearby tracks
without the computational overhead of adding explicit edge features as
suggested in Ref. [9], and was found to improve performance. The input
embeddings are then fed into a 4-layer self-attention transformer encoder
network [13], which performs feature extraction and outputs a second embedding
for each track. This transformer encoder has a model dimension $d$ and a feed-
forward update width of $2d$. ReLU activations are used throughout the model.

###  3.2 Object Decoder

The object decoder computes an explicit representation of a set of object
queries, which represent possible output vertices. At the start of the
decoding stage, object queries are initialised as a set of $N=5$ learned
vectors of dimension $d$. The number of object queries is chosen as the
maximum number of target vertices per jet in the training sample, but could
easily be adapted for other applications. In a given decoder layer, each
object query is updated via information exchange with the updated track
representations from the encoder (via cross-attention) and other object
queries (via self-attention), enabling the model to efficiently aggregate
information about the entire jet. For the cross-attention operations, we use
the MaskAttention operator Ref. [25] to sparsify the attention matrix based on
intermediate mask predictions. This operator enables object queries to attend
only to relevant tracks and improves performance. During the decoding stage,
the track representations are also updated by attending to the object queries.
We use a three layer decoder. The outputs of the final layer are the set of
updated track representations and the set of updated object queries, which are
used to predict the presence and properties of any vertices, along with their
track masks, as described in the following.

###  3.3 Task Heads

A series of task-specific dense networks are used to produce the various model
outputs. Each network has three hidden layers composed of 256, 256, and 128
nodes respectively. More details are provided about each type of output below.

#### Vertex Class

Vertex class labels are either $b\text{-hadron}$, $c\text{-hadron}$, or null,
with the latter indicating that no target vertex is present. If the model
predicts an object query as null, other outputs for that object are ignored.
The other object queries are then interpreted as output vertices, enabling the
model to reconstruct a variable number of vertices. Meanwhile, the separation
of valid (i.e. not null) vertices by flavour allows a full understanding of
the decay process inside the jet. For example, when considering jets produced
from a single hadron, $c\text{-hadrons}$ inside $b$-jets can be assumed to
originate from $b\text{-hadrons}$.

#### Masks

The assignment of input tracks to output vertices is given by a series of $N$
1D masks, each of length $M$. In order to obtain the masks, a dense network is
first used to predict a per-vertex mask token which has the same dimension,
$d$ as the track representation. The masks for each vertex are then computed
as binarised dot products between the mask token and the updated track
embeddings produced by the decoder. A value of 1 in the $i$th position of mask
$j$ indicates that the $i$th input track is assigned to vertex $j$, while 0
indicates no assignment.

#### Vertex Properties

Vertex fitting is implemented as a dense, multi-target regression network with
three output nodes; one for each regression target. We use simulation-level
properties of the heavy flavour hadrons to define these targets. While there
is flexibility to add predictions for a wide range of properties, three have
been chosen for this study: the hadron’s transverse momentum
(${{{{p_{\mathrm{T}}}}}}$), the radial flight length of the hadron with
respect to the beamline (${{{{L_{xy}}}}}$), and the angular displacement
between the jet axis and the hadron (${{{{dR}}}}$). Since each regression
target is always positive, a ReLU activation is applied to the output layer.
The same network is used to predict properties for each object query.

#### Jet Flavour Prediction

The model outputs a probability associated with each possible jet class:
$b$-jet, $c$-jet, or light-jet. In order to perform jet flavour
classification, the track and vertex representations are pooled independently
using attention pooling. The resulting pooled representations are concatenated
and fed into a fourth and final dense network task which outputs the jet
flavour probabilities.

###  3.4 Losses

The loss function is made up of an aggregate loss term for the vertexing tasks
and a categorical cross entropy (CE) loss for the jet flavour prediction.

$L=L_{\text{Vertexing}}+L_{\text{Jet Flavour}}$ (1)

The vertexing loss $L_{\text{Vertexing}}$ is comprised of several terms. The
first of these is a categorical CE loss for the vertex class labels. For the
masks, a linear combination of a binary CE loss and a dice loss [26] are
applied. Mean absolute error losses are applied for the regression targets.
The vertexing loss function is therefore as follows:

$L_{\text{Vertexing}}=0.5L_{\text{Class}}+\underbrace{12L_{\text{BCE}}+8L_{\text{Dice}}}_{L_{\text{Mask}}}+2L_{\text{Regression}},$
(2)

where weights have been applied to each term to produce an acceptable trade-
off between the performance of the different tasks. The optimal bipartite
matching between predictions and targets is used to define the vertexing loss
following the procedure described in Ref. [22].

##  4 Dataset & Experimental Setup

###  4.1 Dataset

Samples of hadronically decaying $t\overline{t}$ events from simulated $pp$
collisions are generated using Pythia 8.307 [27] at a centre-of-mass energy of
$\sqrt{s}=13$ TeV with an average of $50$ interactions per bunch crossing. The
detector response is simulated with Delphes 3.5.0 [28] using a modified
version of the “ATLAS Pileup” configuration. The modifications are intended to
make the simulation more realistic, and are listed below. Firstly, the track
reconstruction efficiency parameterisation was modified to more closely match
current ATLAS tracking software performance as measured in Ref. [29].
Secondly, parameterised impact parameter smearing is applied based on impact
parameter resolutions obtained by ATLAS after the addition of the IBL [30].
Finally, smearing is applied to track $\eta$ and $\phi$, based on the ATLAS
Inner Detector Technical Design Report [31].

Jets are required to satisfy $20<{{{{p_{\mathrm{T}}}}}}<250$ GeV and
$|\eta|<2.5$, and are assigned a flavour based on the presence of heavy
flavour hadrons with ${{{{p_{\mathrm{T}}}}}}>5$ GeV found within $dR<0.4$ of
the jet. A jet is labelled as a $b$-jet if one or more $b\text{-hadrons}$ are
present, a $c\text{-jet}$ if one or more $c\text{-hadrons}$ and no
$b\text{-hadrons}$ are present, and a light-jet otherwise. Resampling is
applied to match the jet ${{{{p_{\mathrm{T}}}}}}$ and $\eta$ distributions of
each jet flavour. Associated heavy flavour hadrons ($b$\- and $c$-hadrons)
with ${{{{p_{\mathrm{T}}}}}}>5$ are used as targets for the vertex
reconstruction tasks.

The full dataset used for training and evaluation has been published at [32].

###  4.2 Training Setup

#### Inputs

The model takes as inputs the 50 highest ${{{{p_{\mathrm{T}}}}}}$ tracks for
each jet. For each track, the following variables are used as inputs: the
signed impact parameters ($d_{0}$ and $z_{0}$) where the sign is positive for
tracks located in front of the primary vertex and negative for those behind,
the fraction of the jet ${{{{p_{\mathrm{T}}}}}}$ which is carried by the
track, and track $\eta$ and $\phi$ relative to the jet axis. The jet
${{{{p_{\mathrm{T}}}}}}$ and $\eta$ are also concatenated to each track. A
total of $13.5$ million jets are used for training, equally split between the
three jet flavours. Two further sets of $1.35$ million jets are used for each
of validation and testing.

#### Training

Models are trained on 3 NVIDIA A100 GPUs, with a combined batch size of 6000
jets. The AdamW [33] optimiser is used with a weight decay of $1\times
10^{-5}$ and a OneCycle learning rate schedule [34] with an initial warm-up
learning rate of $1\times 10^{-7}$, peaking at $5\times 10^{-4}$ once $5\%$ of
training steps have passed, before decreasing to a final value of $1\times
10^{-5}$. Models are trained for 50 epochs and evaluated at the epoch with the
lowest validation loss.

##  5 Results

We evaluate the performance of MaskFormer in three areas: vertex finding,
vertex fitting, and flavour-tagging. For a vertex finding and flavour-tagging
benchmark, we trained an EdgeClassifier model described in Section 2. The
benchmark EdgeClassifier model is trained on the same input jets and has a
similar architecture to the MaskFormer model (without the object decoder). The
width of the transformer encoder for the EdgeClassifier is doubled to
$d^{\prime}=512$ compared to the MaskFormer such that both models have 6.3
million trainable parameters.

###  5.1 Vertex Finding

Vertex finding is the task of correctly assigning tracks to vertices. In order
to measure successful vertex finding, we define two matching criteria; perfect
and loose, which are based on the recall and purity of tracks assigned to a
given predicted vertex. The perfect matching requires predicted vertices to
have a recall and purity of 100%, while the loose criteria requires a minimum
recall of 50% and minimum purity of 50%, where recall and purity are defined
below.

$\text{Recall}=\frac{N_{\text{Correctly Associated Tracks}}}{N_{\text{Tracks
in Target Vertex}}}~{},~{}~{}~{}~{}\text{Purity}=\frac{N_{\text{Correctly
Associated Tracks}}}{N_{\text{All Associated Tracks}}}$ (3)

As the EdgeClassifier reconstructs vertices from compatible pairs of tracks,
only vertices with 2 or more tracks associated at truth level can be
identified. While the MaskFormer has no such intrinsic limitation, we consider
only vertices with at least two tracks associated at truth level to facilitate
a comparison of the two approaches. In addition, hadrons with
${{{{L_{xy}}}}}<1$ mm at simulation level are excluded from the evaluation, as
these are not sufficiently displaced from the primary $pp$ interaction point
to be resolved.

The heavy-flavour vertex reconstruction efficiency of the two models as a
function of hadron ${{{{L_{xy}}}}}$ are shown in Fig. 2 for both the perfect
and loose matching criteria. The results are shown for all vertices within
heavy flavour jets: these include $b\text{-hadron}$ and $b\rightarrow$
$c\text{-hadron}$ vertices (inside $b$-jets) and $c\text{-hadron}$ vertices
(inside $c\text{-jets}$). For the perfect matching criteria, MaskFormer shows
improvements of approximately $10\%$ or more with respect to the
EdgeClassifier over the full hadron ${{{{L_{xy}}}}}$ range. In addition,
MaskFormer is able to find more than $75\%$ of vertices with
${{{{L_{xy}}}}}>5$ mm with the loose matching criteria. MaskFormer also
displays a reduced dependence on ${{{{L_{xy}}}}}$.

Figure 2:  Heavy-flavour vertex reconstruction efficiency as a function of
hadron ${{{{L_{xy}}}}}$ for vertices in $b$\- and $c$-jets (HF-jets) for the
perfect match (left) and loose match (right). Only vertices with at least two
tracks at truth level are considered. The shaded region represents binomial
error bands.

To study the fake rate (i.e. the rate at which the models predict a vertex
when none exists), we examine the overall heavy flavour vertex reconstruction
rate in light-jets. As light-jets by definition contain zero heavy flavour
hadrons, any reconstructed vertices can be considered to be false positives.
The fake rate for each model is compatible within the statistical uncertainty
and broadly constant as a function of jet ${{{{p_{\mathrm{T}}}}}}$, and was
found to be approximately $5\%$ percent within the events studied. The vertex
finding efficiency and fake rate can be traded-off by tuning the thresholds
placed on the model outputs.

MaskFormer has good vertex finding performance but a reduction is observed as
the multiplicity of heavy flavour hadrons in the jet increases above two
vertices inside a single jet. Vertices within jets with containing single
hadron (typically $c\text{-jets}$ or $b$-jets without a b→cb→c decay) are
reconstructed with a loose match efficiency of $79\%$ and $77\%$ in jets with
exactly 2 hadrons (typically $b$-jets with one $c\text{-hadron}$). Meanwhile
three or more hadrons (typically $b$-jets with two $c\text{-hadrons}$) are
reconstructed with $60\%$ efficiency.

MaskFormer can also assign the correct track to a hadronic vertex with only
one single reconstructed decay product with an efficiency of 48%. Meanwhile
the rate to reconstruct single track vertices in light-jets is 0.055. This can
be compared to an overall loose match efficiency (rate in light-jets) of 74%
(0.048) for vertices with two or more reconstructed tracks at truth level.

###  5.2 Vertex Fitting

Vertex fitting is the estimation of the properties of a vertex from its
constituent tracks. MaskFormer fits the properties of all vertices in the jet
simultaneously via a series of regressions to simulation-level quantities of
heavy flavour hadrons, as described in Section 3.3. As discussed, this is a
key improvement over the EdgeClassifier, which is unable to fit an unknown
number of vertices in-model.

The normalised residuals of the MaskFormer vertex fit are shown in Fig. 3 for
vertex ${{{{p_{\mathrm{T}}}}}}$, ${{{{L_{xy}}}}}$, and ${{{{dR}}}}$, as a
function of the respective target. The results are shown for all reconstructed
vertices without applying quality requirements on the assigned tracks.
Overall, MaskFormer displays a promising ability to recover the properties of
hadronic vertices inside jets. Normalised residuals are largest at low values
of the target distribution, but are otherwise generally consistent in the bulk
of the distributions. For each target, and especially for
${{{{p_{\mathrm{T}}}}}}$, the model tends to overestimate below the mean and
vice versa. Further work is required to reduce this bias towards predicting
the central values of each distribution.

Figure 3:  Normalised fit residuals
$(x_{\text{Pred}}-x_{\text{True}})/x_{\text{True}}$ for
${{{{p_{\mathrm{T}}}}}}$ (left), ${{{{L_{xy}}}}}$ (middle), and ${{{{dR}}}}$
(right) as a function of the respective target properties for all
reconstructed hadrons. The mean-normalised residual is shown by the red line,
while the box and whiskers demarcate the region containing $68\%$ and $95\%$
of the data respectively.

Normalised residuals are also shown for ${{{{p_{\mathrm{T}}}}}}$,
${{{{L_{xy}}}}}$, and ${{{{dR}}}}$ as a function of the number of correctly
assigned tracks in the predicted vertex in Fig. 4. For all variables, the
accuracy of the prediction has a significant positive correlation with the
number of correctly assigned tracks, and the residuals are largest when no
tracks are correctly assigned. The number of correctly assigned tracks
required for $95\%$ of vertices to have predictions within $50\%$ of their
target value is 5 for the prediction of ${{{{p_{\mathrm{T}}}}}}$, but only 2
for the prediction of ${{{{L_{xy}}}}}$, reflecting the fact that it is easier
to recover the hadron ${{{{L_{xy}}}}}$ with fewer tracks.

Figure 4:  Normalised fit residuals
$(x_{\text{Pred}}-x_{\text{True}})/x_{\text{True}}$ for
${{{{p_{\mathrm{T}}}}}}$ (left), ${{{{L_{xy}}}}}$ (middle), and $\Delta
R(\text{Jet, Hadron})$ as a function of the number of correctly associated
tracks to the hadronic vertex for all reconstructed hadrons.

###  5.3 Flavour-tagging

Previous approaches in flavour tagging have shown that the inclusion of
information about secondary vertices improves flavour identification
performance [35, 17]. As discussed in Section 1, a more recent shift towards
directly learning the presence of vertices has been taking place [10]. In this
work, we complete the transition with the full in-model reconstruction of an
arbitrary number of secondary vertices combined with a prediction of the jet
flavour.

In Fig. 5, the flavour-tagging performance of the EdgeClassifier baseline and
MaskFormer models is shown. Both models classify jets using a discriminant
following the procedure used in Ref. [10]. A significant improvement in
background rejection is observed for the MaskFormer, with $50\%$ ($15\%$)
improvement in light-jet rejection ($c$-jet rejection) observed for a $b$-jet
selection efficiency of $70\%$. This result suggests that the MaskFormer
architecture is able to more efficiently extract relevant information from the
training data. Future work will include examining this comparison for more
realistic simulations and increased training statistics.

Figure 5:  Background rejection as a function of $b$-jet selection efficiency
for light-jets (solid) and $c\text{-jets}$ (dashed). The bottom panels show
the ratio of rejection between the MaskFormer and EdgeClassifier. Binomial
errors are shown as shaded bands.

##  6 Conclusion

In this work, we train a neural network with a MaskFormer-based architecture
to fully reconstruct secondary vertices from heavy flavour hadron decays. This
architecture is able to find and fit an arbitrary and unknown number of
secondary vertices within jets, making it an example of full many-to-many
reconstruction. When comparing the vertex finding performance to an existing
state-of-the-art approach, MaskFormer is found to significantly improve vertex
finding efficiency at a comparable fake rate. Crucially, MaskFormer also
enables the prediction of vertex properties via a class label prediction and
set of regression tasks. This allows for a unified approach to vertex
reconstruction, with vertex finding and fitting optimised via a single learned
model, allowing it to be versatile and easily re-optimised to new
environments. The MaskFormer model also demonstrates improved flavour tagging
performance when compared to state-of-the-art edge classification approaches
which only perform vertex finding.

Further work will involve deploying this vertex reconstruction tool in a more
complex environment, such as the intricate topology of boosted large-radius
jets. The model could also be adapted to reconstruct primary $pp$ interaction
vertices and other secondary vertices not from heavy flavour hadron decay. We
hypothesise that the MaskFormer architecture will prove to be extremely
flexible and could be repurposed for a wide range of many-to-many
reconstruction tasks in particle physics, for example track reconstruction or
the analysis of full events.

## References

  * [1] Lyndon Evans and Philip Bryant  “LHC Machine”  In _JINST_ 3, 2008, pp. S08001  DOI: 10.1088/1748-0221/3/08/S08001
  * [2] ATLAS Collaboration  “The ATLAS Experiment at the CERN Large Hadron Collider”  In _JINST_ 3, 2008, pp. S08003  DOI: 10.1088/1748-0221/3/08/S08003
  * [3] CMS Collaboration  “The CMS Experiment at the CERN LHC”  In _JINST_ 3, 2008, pp. S08004  DOI: 10.1088/1748-0221/3/08/S08004
  * [4] ATLAS Collaboration  “ATLAS flavour-tagging algorithms for the LHC Run 2 pp collision dataset”  In _Eur. Phys. J. C_ 83.7, 2023, pp. 681  DOI: 10.1140/epjc/s10052-023-11699-1
  * [5] CMS Collaboration  “Identification of heavy-flavour jets with the CMS detector in pp collisions at 13 TeV”  In _JINST_ 13.05, 2018, pp. P05011  DOI: 10.1088/1748-0221/13/05/P05011
  * [6] P.A. Zyla  “Review of Particle Physics”  In _Prog. Theor. Exp. Phys._ 2020, 2020, pp. 083C01  DOI: 10.1093/ptep/ptaa104
  * [7] B. Brandl et al.  “Multivariate analysis methods to tag b-quark events at LEP/SLC”  In _Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment_ 324.1, 1993, pp. 307–316  DOI: https://doi.org/10.1016/0168-9002(93)90993-R
  * [8] Huilin Qu and Loukas Gouskos  “ParticleNet: Jet Tagging via Particle Clouds”  In _Phys. Rev. D_ 101.5, 2020, pp. 056019  DOI: 10.1103/PhysRevD.101.056019
  * [9] Huilin Qu, Congqiao Li and Sitian Qian  “Particle Transformer for Jet Tagging”, 2022  arXiv:2202.03772 [hep-ph]
  * [10] ATLAS Collaboration  “Graph Neural Network Jet Flavour Tagging with the ATLAS Detector”, ATL-PHYS-PUB-2022-027, 2022  URL: https://cds.cern.ch/record/2811135
  * [11] ATLAS Collaboration  “Transformer Neural Networks for Identifying Boosted Higgs Bosons decaying into $b\bar{b}$ and $c\bar{c}$ in ATLAS”, ATL-PHYS-PUB-2023-021, 2023  URL: https://cds.cern.ch/record/2866601
  * [12] Bowen Cheng, Alexander G. Schwing and Alexander Kirillov  “Per-Pixel Classification is Not All You Need for Semantic Segmentation”, 2021  arXiv:2107.06278 [cs.CV]
  * [13] Ashish Vaswani et al.  “Attention Is All You Need”, 2023  arXiv:1706.03762 [cs.CL]
  * [14] Rudolf Frühwirth and Are Strandlie  “Pattern Recognition, Tracking and Vertex Reconstruction in Particle Physics Experiments”  Springer, 2021  DOI: 10.1007/978-3-030-65771-0
  * [15] ATLAS Collaboration  “Expected performance of the ATLAS $b$-tagging algorithms in Run-2”, ATL-PHYS-PUB-2015-022, 2015  URL: https://cds.cern.ch/record/2037697
  * [16] ATLAS Collaboration  “Topological $b$-hadron decay reconstruction and identification of $b$-jets with the JetFitter package in the ATLAS experiment at the LHC”, ATL-PHYS-PUB-2018-025, 2018  URL: https://cds.cern.ch/record/2645405
  * [17] CMS Collaboration  “Identification of heavy-flavour jets with the CMS detector in $pp$ collisions at $13\,\text{TeV}$”  In _JINST_ 13, 2018, pp. P05011  DOI: 10.1088/1748-0221/13/05/P05011
  * [18] Jonathan Shlomi et al.  “Secondary vertex finding in jets with neural networks”  In _Eur. Phys. J. C_ 81.6, 2021, pp. 540  DOI: 10.1140/epjc/s10052-021-09342-y
  * [19] Rachel E.. Smith et al.  “Differentiable Vertex Fitting for Jet Flavour Tagging”, 2023  arXiv:2310.12804 [hep-ex]
  * [20] Joseph Redmon, Santosh Divvala, Ross Girshick and Ali Farhadi  “You Only Look Once: Unified, Real-Time Object Detection”, 2016  arXiv:1506.02640 [cs.CV]
  * [21] Shaoqing Ren, Kaiming He, Ross Girshick and Jian Sun  “Faster R-CNN: Towards Real-Time Object Detection with Region Proposal Networks”, 2016  arXiv:1506.01497 [cs.CV]
  * [22] Nicolas Carion et al.  “End-to-End Object Detection with Transformers”, 2020  arXiv:2005.12872 [cs.CV]
  * [23] Kaiming He, Georgia Gkioxari, Piotr Dollár and Ross Girshick  “Mask R-CNN”, 2018  arXiv:1703.06870 [cs.CV]
  * [24] Alexander Kirillov et al.  “Segment Anything”, 2023  arXiv:2304.02643 [cs.CV]
  * [25] Bowen Cheng et al.  “Masked-attention Mask Transformer for Universal Image Segmentation”, 2022  arXiv:2112.01527 [cs.CV]
  * [26] Carole H. Sudre et al.  “Generalised Dice Overlap as a Deep Learning Loss Function for Highly Unbalanced Segmentations”  In _Lecture Notes in Computer Science_ Springer International Publishing, 2017, pp. 240–248  DOI: 10.1007/978-3-319-67558-9˙28
  * [27] Christian Bierlich  “A comprehensive guide to the physics and usage of PYTHIA 8.3”  In _SciPost Phys. Codeb._ 2022, 2022, pp. 8  DOI: 10.21468/SciPostPhysCodeb.8
  * [28] J. Favereau et al.  “DELPHES 3: a modular framework for fast simulation of a generic collider experiment”  In _Journal of High Energy Physics_ 2014.2 Springer ScienceBusiness Media LLC, 2014  DOI: 10.1007/jhep02(2014)057
  * [29] ATLAS Collaboration  “Software Performance of the ATLAS Track Reconstruction for LHC Run 3”, 2023  arXiv:2308.09471 [hep-ex]
  * [30] ATLAS Collaboration  “ATLAS Insertable B-Layer: Technical Design Report”, 2010  URL: https://cds.cern.ch/record/1291633
  * [31] ATLAS Collaboration  “ATLAS Inner Detector: Technical Design Report, Volume 1”, 1997  URL: https://cds.cern.ch/record/331063
  * [32] N. Pond  “Top quark pair events for heavy flavour tagging and vertexing at the LHC”, 2023  DOI: 10.5281/zenodo.10371998
  * [33] Ilya Loshchilov and Frank Hutter  “Decoupled Weight Decay Regularization”, 2019  arXiv:1711.05101 [cs.LG]
  * [34] Leslie N. Smith and Nicholay Topin  “Super-Convergence: Very Fast Training of Neural Networks Using Large Learning Rates”, 2018  arXiv:1708.07120 [cs.LG]
  * [35] ATLAS Collaboration  “ATLAS b-jet identification performance and efficiency measurement with $t{\bar{t}}$ events in pp collisions at $\sqrt{s}=13$ TeV”  In _Eur. Phys. J. C_ 79.11, 2019, pp. 970  DOI: 10.1140/epjc/s10052-019-7450-8