\section{Introduction}

Investigations of locomotion of a free-living nematode {\it 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
\cite{%
  Ghanbari-Nock-Johari-Blaikie-Chen-Wang:2012,%
  Johari-Nock-Alkaisi-Wang:2013,%
  Etheridge-Rahman-Gaffney-Shaw-Shephard-Magudia-Solomon-Milne-Blawzdziewicz-Constantin-Teodosiu-Greenhaff-Vanapalli-Szewczyk:2015},
neuromuscular coordination
\cite{%
  Wen-Po-Hulme-Chen-Liu-Kwok-Gershow-Leifer-Butler-Fang-Yen-Kawano-Schafer-Whitesides-Wyart-Chklovskii-Zhen-Samuel:2012,%
  Stirman-Crane-Husson-Wabnig-Schultheis-Gottschalk-Lu:2011},
and information processing
\cite{%
  Gray-Hill-Bargmann:2005,%
  Kocabas-Shen-Guo-Ramanathan:2012%
}.
Locomotory readouts are thus extensively used to describe impacts of
genetic mutations
\cite
   Sznitman-Purohit-Krajacic-Lamitina-Arratia:2010,%
   Brown-Yemini-Grundy-Jucikas-Schafer:2013,%
   Restif-Ibanez_Ventoso-Vora-Guo-Metaxas-Driscoll:2014%
}
and evaluate pharmacological effects
\cite{Artal_Sanz-deJong-Tavernarakis:2006,%
  Omura-Clark-Samuel-Horvitz:2012}.

In its natural environment (soil rich in decaying organic matter),
{\it 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
\cite{Cohen-Boyle:2010,%
Fang_Yen-Wyart-Xie-Kawai-Kodger-Chen-Wen-Samuel:2010,%
Sznitman-Shen-Purohit-Arratia:2010,%
Padmanabhan-Khan-Solomon-Armstrong-Rumbaugh-Vanapalli-Blawzdziewicz:2012,%
Stephens-Johnson-Kerner-Bialek-Ryu:2008},
and little effort has been made to characterize worm motility in
3D\ environments that would mimic {\it C. elegans}\  natural habitat.

Only recently was 3D\ burrowing motion of {\it C. elegans}\ in soft
substrates quantitatively imaged
\cite{Kwon-Pyo-Lee-Je:2013}
and feasibility  of using 3D\  motion to assess genetic mutations
\cite{Kwon-Hwang-You-Lee-Je:2015}
and to identify muscular dystrophy
\cite{Beron-Vidal-Gadea-Cohn-Parikh-Hwang-Pierce-Shimomura:2015}
experimentally demonstrated.  Such 3D\ methods of scoring
{\it 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 {\it C. elegans}\ have been reported
\cite{Deng-Xu:2014,%
  Delmotte-Climent-Plouraboue: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 {\it 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
\cite{Padmanabhan-Khan-Solomon-Armstrong-Rumbaugh-Vanapalli-Blawzdziewicz:2012}.
The PHC model has been shown to give an accurate parameterization
of body postures (defined by the configuration of
the body centerline
\cite{%
  Stephens-Johnson-Kerner-Bialek-Ryu:2008,%
  Sznitman-Shen-Sznitman-Arratia:2010,%
  Husson-Costa-Schmitt-Gottschalk:2012,%
  Berman-Kenneth-Sznitman-Leshansky: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
\cite{Bilbao-Wajnryb-Vanapalli-Blawzdziewicz:2013}, we can describe 3D
motion of {\it 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
\cite{Bilbao-Wajnryb-Vanapalli-Blawzdziewicz:2013}.  The combined
PHC/PCT gait model will help analyze diverse aspects of
3D\ undulatory locomotion of {\it 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 {\it C. elegans}.

Our paper is organized as follows.  In Sec.\ \ref{crawling section} we
develop geometrical gait models for {\it 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.\ \ref{maneuverability in fluids} we present our
hydrodynamic analysis of reorientation maneuvers of
{\it C. elegans}\ swimming in a viscous fluid.  The geometrical and
hydrodynamic locomotion descriptions developed in Secs.\ \ref{crawling
  section} and \ref{maneuverability in fluids} are combined in
Sec.\ \ref{chemotaxis in 3d} with simple models of chemosensation and
motion control to describe chemotaxis of burrowing and swimming
{\it C. elegans}\ in 3D.  Implications of our findings for the emerging
area of 3D\ undulatory locomotion are discussed in
Sec.\ \ref{Summary and discussion}.

\figsgone{
\ \clearpage

\begin{figure}
\includegraphics[width=0.34\textwidth]{Fig1.eps}
\caption
{ ({\it Colour online}) Family of curves reproduced from harmonic
curvature \eqref{single mode curvature} 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) $\wavevectorL=5.5$ (typical
swimming gait) and (b) $\wavevectorL=9$ (typical
crawling gait).  Depending on the wavevector, the
normalized amplitude $A/q=1$ may correspond to
(a) \textit{C}-shaped\ (top) and (b) \textit{W}-shaped\ body postures;
$A/q=1.8$ produces an $\Omega$-shaped body
posture.
\label{curves with no torsion}
}
\end{figure}

\begin{figure}
\includegraphics[width=0.45\textwidth]{Fig2.eps}
\caption{ 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 $q
s^\curvatural_1=3 \pi /2$ and $q \Delta
s^\curvatural=0.6$.  The turn angle $\theta^\turn$ and
mode-switching points $s^\curvatural_1$ and $s^\curvatural_2$
are defined in the top inset.  Bottom inset shows the dependence
of the turn angle $\theta^\turn$ (in degrees) on the turning-mode
amplitude $A_2 / q$. The black dots correspond
to the amplitudes of the trajectories shown.
\label{slipless turn trajectory}
}
\end{figure}
}
\section{Geometrical description of the nematode gait}
\label{crawling section}

\sloppypar {\it 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.\ \ref{Shape and
  motion in two-dimensional space} we summarize our previously
reported PHC model
\cite{Padmanabhan-Khan-Solomon-Armstrong-Rumbaugh-Vanapalli-Blawzdziewicz:2012}
for describing body postures and trajectories of {\it C. elegans}\ moving in
two dimensions.  In Sec.\ \ref{crawl 3d description} 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 {\it C. elegans}\ swimming in two dimensions
\cite{Padmanabhan-Khan-Solomon-Armstrong-Rumbaugh-Vanapalli-Blawzdziewicz:2012,%
Bilbao-Wajnryb-Vanapalli-Blawzdziewicz:2013},
we assume that PHC model and its 3D\ generalization provide a
realistic representation for the swimming gait.  In
Sec.\ \ref{maneuverability in fluids} the results derived in
Sec.\ \ref{crawling section} are used to develop a comprehensive
analysis that combines the gait geometry and motion mechanics to fully
describe 3D\ maneuvers of swimming {\it C. elegans}.

\subsection{Description in 2D}
\label{Shape and motion in two-dimensional space}

The body of {\it C. elegans}\ is highly elongated, and therefore its postures
can be faithfully represented by the body centerline
\cite{Padmanabhan-Khan-Solomon-Armstrong-Rumbaugh-Vanapalli-Blawzdziewicz:2012,
  Stephens-Johnson-Kerner-Bialek-Ryu:2008,
  Kim-Park-Mahadevan-Shin:2011}.
In our analysis, the shape of the centerline is described by the worm
curvature function $\curvature_w(\wormcoord^\prime,t)$, where
$t$ denotes time and $\wormcoord^\prime$ is the arclength coordinate
along the body.  In what follows $\wormcoord^\prime=0$ is the tail and
$\wormcoord^\prime=L$ is the head coordinate value (where
$L$ is the nematode length).  The body configuration is
fully determined by providing the curvature function $\curvature_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.\ \ref{curves with no torsion}), the worm and trail curvature
functions are related
\eqgone{
\begin{equation}
  \curvature_w(\wormcoord^\prime,t)=
  \kappa(\wormcoord^\prime+v_s
  t),
  \label{curvature equation 1}
\end{equation}
}
where $v_s$ is the velocity along the trail.  Body postures and
the trail followed by a nematode are constructed from the curvature
\eqref{curvature equation 1} by solving the Serret--Frenet equations
\cite{Stoker:1969}
\eqgone{
   \begin{subequations}
      \label{2D serret-frenet equations}
      \begin{alignat}{2}
         \frac{d {\hat{\boldsymbol{t}}}}{d s}&=\phantom{-}
            \kappa {\hat{\boldsymbol{n}}}, \\
         \frac{d {\hat{\boldsymbol{n}}}}{d s}&=-\kappa {\hat{\boldsymbol{t}}},
      \end{alignat}
   \end{subequations}
}
where ${\hat{\boldsymbol{t}}}$ and ${\hat{\boldsymbol{n}}}$ are mutually perpendicular unit tangent
and normal vectors.

Observations of {\it C. elegans}\ swimming in water indicate that the
propagating-curvature-wave relation \eqref{curvature equation 1}
applies not only to no-slip motion [in which case Eq.\ \eqref{curvature
    equation 1} results from geometrical constraints], but  is also
  preserved with good accuracy when significant transverse slip occurs
\cite{Padmanabhan-Khan-Solomon-Armstrong-Rumbaugh-Vanapalli-Blawzdziewicz:2012}.
This behavior is a signature of the proprioceptive coupling
\cite{Wen-Po-Hulme-Chen-Liu-Kwok-Gershow-Leifer-Butler-Fang-Yen-Kawano-Schafer-Whitesides-Wyart-Chklovskii-Zhen-Samuel: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.

\gone{
In Sec.\ \ref{Harmonic-curvature representation
  of nematode gait} we describe a simple analytical model for the
nematode curvature wave \eqref{curvature equation 1}.

\subsubsection{PHC representation of nematode gait}
\label{Harmonic-curvature representation of nematode gait}

{\it Caenorhabditis elegans}\ produces a variety of body postures, ranging from shallow
undulations to highly curved $\Omega$ shapes and loop
shapes. Recently, we have demonstrated
\cite{Padmanabhan-Khan-Solomon-Armstrong-Rumbaugh-Vanapalli-Blawzdziewicz:2012}
that all these body shapes can be accurately represented using the
harmonic curvature function
\eqgone{
\begin{equation}
  \kappa(s)=
            {A} \cos \left(
            q s + \phi \right),
            \label{single mode curvature}
\end{equation}
}
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.\ \eqref{curvature equation 1}--\eqref{single mode
  curvature} yields a family of curves representative of the most
common nematode postures.  For moderate values of the dimensionless
amplitude $A/q\approx1$, typical crawling postures (\textit{W}-shaped)
and swimming postures (\textit{C}-shaped) are obtained, as illustrated in
Fig.\ \ref{curves with no torsion}.
Higher amplitudes
correspond to $\Omega$-shapes ($A/q\approx2$) and loop shapes
($A/q\approx3$,  described in
\cite{Padmanabhan-Khan-Solomon-Armstrong-Rumbaugh-Vanapalli-Blawzdziewicz:2012}).

According to the PHC model
\cite{Padmanabhan-Khan-Solomon-Armstrong-Rumbaugh-Vanapalli-Blawzdziewicz:2012},
trails of crawling nematodes can be accurately represented by the
harmonic function \eqref{single mode curvature} 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.

\subsubsection{Planar turns}
\label{Planar turns}

To investigate the turning ability of {\it 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,
\eqgone{
\begin{equation}
     \kappa(s)= \left\{
     \begin{matrix}
        {A_1} \cos \left(
           q s \right), \quad
           s < s^\curvatural_1, \phantom{--.} \\
        {A_2} \cos \left(
           q s \right), \quad
             s^\curvatural_1 < s < s^\curvatural_2, \\
          {A_1} \cos \left(
             q s \right), \quad
             s^\curvatural_2 < s, \phantom{--.} \\
       \end{matrix}
       \right.
\label{three mode curvature}
\end{equation}
}
where $A_1$ is the base crawling amplitude (with
$A_1/q\approx1$), $A_2 > A_1$ is
the turning amplitude, and $s^\curvatural_1$ and $s^\curvatural_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 \ref{slipless turn trajectory} shows an example of a crawling
nematode executing a three-mode planar turn for fixed switching points
$s^\curvatural_1$ and $s^\curvatural_2$ and several values of the
normalized turning mode amplitudes $A_2/q$.  The
dependence of the turn angle $\theta^\turn$ (defined in the top
inset of Fig.  \ref{slipless turn trajectory}) on the turn amplitude
$A_2$ is depicted in the bottom inset of Fig.
\ref{slipless turn trajectory}.

\figsgone{
\ \clearpage

\begin{figure}
\includegraphics[width=0.34\textwidth]{Fig3.eps}
\caption
{ ({\it 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
$q L=9$.  The green, blue and red arrows are
the unit tangent, normal, and binormal vectors, respectively.
\label{curves with and no torsion}
}
\end{figure}

\begin{figure}
\includegraphics[width=0.46\textwidthF]{Fig4.eps}
\caption
{Roll rotation rate per undulation period $\omega^\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 $\turnvelBellyq/\tau_0$ on the
amplitude $A/q$.
\label{turn rate off plane dry}
}
\end{figure}

\begin{figure}
\includegraphics[width=0.45\textwidth]{Fig5.eps}
\caption
{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.
\label{planes of motion exp}
   }
\end{figure}

\ \clearpage
}

\subsection{Description in 3D}
   \label{crawl 3d description}

{\it Caenorhabditis elegans}, which crawls and burrows in rotting
vegetation and swims in rainwater, performs 2D\ and
3D\ maneuvers
\cite{Kwon-Pyo-Lee-Je:2013,%
  Voorhies-Fuchs-Thomas:2005,%
  Frezal-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
\cite{Deng-Xu:2014},
burrowing
\cite{Kwon-Pyo-Lee-Je:2013},
and swimming.  A crawling {\it 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.\ \cite{Deng-Xu:2014}.  Swimming nematodes
undergo a similar out-of-plane body deformations, accompanied by
rapid rotation of the undulation plane (see Sec.\
\ref{maneuverability in fluids}).  Three-dimensional body actuation
seems also to occur for nictating nematodes (see supplemental movies
in
\cite{Lee-Choi-Lee-Kim-Hwang-Kim-Park-Paik-Lee:2012}).

To incorporate out-of-plane body postures into the PHC model, we use the
complete 3D\ set of Serret--Frenet equations,
\eqgone{
\begin{subequations}
      \label{full serret-frenet equations}
      \begin{alignat}{3}
         \frac{d {\hat{\boldsymbol{t}}}}{d s}&=\phantom{-}
            \kappa {\hat{\boldsymbol{n}}}, \\
         \frac{d {\hat{\boldsymbol{n}}}}{d s}&=-\kappa {\hat{\boldsymbol{t}}}+
            \tau {\hat{\boldsymbol{b}}}, \\
         \frac{d {\hat{\boldsymbol{b}}}}{d s}&=-\tau {\hat{\boldsymbol{n}}},
      \end{alignat}
\end{subequations}
}
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.

\subsubsection{Constant-torsion assumption}

Quantitative experimental information regarding out-of-plane body
postures of {\it 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 \eqref{single mode curvature} with the
constant-torsion assumption,
\eqgone{ \begin{equation} \tau(s)=\tau_0.  \label{constant
torsion} \end{equation} }

The solution of the Serret--Frenet equations \eqref{full
serret-frenet equations} with the curvature and torsion given by
equations \eqref{single mode curvature} and \eqref{constant torsion}
is depicted in Fig.\ \ref{curves with and no torsion}.  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.\ \cite{Deng-Xu:2014}.  The gait
described by Eqs.\ \eqref{single mode curvature} and \eqref{constant
torsion} will be called a single-mode torsional roll.

Figure \ref{turn rate off plane dry} shows the torsional-roll rotation
rate per undulation  period,
\eqgone{
\begin{equation}
  \omega^\roll= \cos^{-1}
  \left[
    {\hat{\boldsymbol{b}}}(0)\boldsymbol{\cdot}{\hat{\boldsymbol{b}}}(2 \pi / q)
    \right].
  \label{definition of turn velocity}
\end{equation}
}
The rotation rate \eqref{definition of turn velocity} 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.\ \ref{turn
  rate off plane dry}, the slope of the lines $\omega^\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.

\subsubsection{Three-mode roll maneuvers}

In the case of variable torsion, we assume that the
torsion wave $\tau_w(\wormcoord^\prime,t)$ propagates backward
along the nematode body,
\eqgone{
\begin{equation}
      \tau_w (\wormcoord^\prime,t)=
         \tau(\wormcoord^\prime+v_s t),
      \label{torsion equation 1}
\end{equation}
}
similar to the curvature wave \eqref{curvature equation 1}.  For
no-slip motion (e.g., burrowing in a gel), relation \eqref{torsion
  equation 1} is required by the constraints that all segments
of the nematode body follow a single trail.  By analogy with our
2D\ results
\cite{Padmanabhan-Khan-Solomon-Armstrong-Rumbaugh-Vanapalli-Blawzdziewicz:2012},
we assume that the form \eqref{torsion equation 1} also applies to
swimming.

According to our qualitative observations of nematode motion in
water, {\it 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
\eqref{curvature equation 1} superposed with piecewise constant
torsion
\eqgone{
\begin{equation}
  \tau(s)=
  \left\{
  \begin{matrix}
    {0} , \quad
    s < s^\torsional_1, \phantom{--.} \\
               {\tau_0} , \quad
               s^\torsional_1 < s < s^\torsional_2, \\
                              {0} , \quad
                              s^\torsional_2 < s, \phantom{--.} \\
  \end{matrix}
  \right.
  \label{three mode torsional turn}
\end{equation}
}
(where $s^\torsional_1$ and $s^\torsional_2$ are the mode
switching points, and
$\Deltas^\torsional=s^\torsional_2-s^\torsional_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.\ \ref{planes of motion exp} and Video
S1 in ESM).

\begin{figure}
\includegraphics[width=0.45\textwidth]{Fig6.eps}
\caption
{ ({\it Colour online}) A sequence of snapshots of a swimming
\textit{C}-shaped\ nematode ($q L=5.5$) performing a
three-mode planar turn with the default-mode amplitude
$A_1/q=1$ and turning mode amplitude as labeled,
for $q s^\curvatural_1=3 \pi /2$ and $q
\Delta s^\curvatural=0.6$.  The gray image represents the
starting position for all trajectories.  The inset shows the
dependence of the turn angle $\theta^\turn$ (in degrees) on the
turning-mode amplitude $A_2 / q$. The dots
correspond to the snapshots shown.
\label{on plane turn swimming}
}
\end{figure}

\begin{figure}
\includegraphics[width=0.46\textwidthF]{Fig7.eps}
\caption
{ Planar turn angle $\theta^\turn$ (in degrees) vs normalized
turning-mode length $\Delta s^\curvatural / L$ for (a)
\textit{C}-shaped\ ($q L=5.5$) and (b)
\textit{W}-shaped\ ($q L=9$) nematodes (see insets)
and turning-mode amplitude $A_2 / q$ as
labeled. Turns initiated at $q s^\curvatural_1=0$
(solid lines) and $\frac{3}{2} \pi$ (dashed). The normalized
amplitude of the default forward mode is $A_1 /
q=1$.
\label{planar turning performance of swimming worm rpy}
}
\end{figure}

\section{Nematode maneuverability in fluids}
   \label{maneuverability in fluids}

The PHC and PHC/PCT gait representations developed in
Sec.\ \ref{crawling section} will now be combined with a hydrodynamic
bead-based model to describe body reorienting maneuvers performed by
swimming {\it 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.\ \eqref{three mode curvature} and
\eqref{three mode torsional turn} (see simulation images shown in
Fig.\ \ref{on plane turn swimming}).

The hydrodynamic forces and
torques acting on the chain are calculated in the Stokes-flow limit
using the generalized Rotne--Prager approximation
\cite{Wajnryb-Mizerski-Zuk-Szymczak:2013}.
This approximation is much faster (although less accurate) than the
\textsc{Hydromultipole} algorithm
\cite{Cichocki-Felderhof-Hinsen-Wajnryb-Blawzdziewicz:1994}
used in our previous study
\cite{Bilbao-Wajnryb-Vanapalli-Blawzdziewicz:2013}.

Our methods for calculating the hydrodynamic friction and mobility
tensors are described in
Ref.\ \cite{Bilbao-Wajnryb-Vanapalli-Blawzdziewicz:2013} and in ESM.
We note that in an independent study Berman \textit{et al.}
\cite{Berman-Kenneth-Sznitman-Leshansky:2013}
also used a hydrodynamic bead-based model to investigate undulatory
swimming of a nematode.

\subsection{Swimming in 2D}
\label{Swimming in 2D}

To assess 2D\ reorientation ability of swimming {\it C. elegans}, we
consider trajectories of a nematode performing three-mode PHC turns,
and evaluate the angle $\theta^\turn$ between the initial and
final swimming directions.  Similar to our analysis of no-slip turns
(Sec.\ \ref{Planar turns}), 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.\ \ref{on
  plane turn swimming} with the corresponding results for a nematode
crawling without sidewise slip (Fig.\ \ref{slipless turn trajectory})
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^\turn$ and the amplitude
of the turning mode $A_2/q$ is nonlinear (in
contrast to the linear dependence for no-slip motion).

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

According to the results shown in Figs.\ \ref{on plane turn swimming}
and \ref{planar turning performance of swimming worm rpy}, a swimming
nematode can control the turn angle by varying the turning-mode
amplitude $A$ as well as the turn initiation point $s^\curvatural_1$ and turning-mode length $\Deltas^\curvatural$.
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.

\figsgone{
\ \clearpage

\begin{figure}
\includegraphics[width=0.45\textwidth]{Fig8.eps}
\caption
{ Time-lapse images of a swimming {\it C. elegans}\ switching the plane of
undulations.
\label{real celegans turning 3d}
}
\end{figure}

\begin{figure}
\includegraphics[width=0.33\textwidth]{Fig9.eps}
\caption
{ (a) Non-planar body posture of a swimming {\it C. elegans}\ and (b)
qualitatively matching calculated body shape with constant
torsion $\tau_0L=0.9$ and harmonic curvature
of the normalized amplitude $A/q=1$ and
wavevector $\wavevectorL=5.5$.
\label{real versus simulated worms 3d turn comparison}
}
\end{figure}

\begin{figure}
\includegraphics[width=0.45\textwidth]{Fig10.eps}
\caption
{Stereoscopic views of a \textit{C}-shaped\ nematode
($\wavevectorL=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_0L=1.5$ and length
$q\Deltas^\torsional=3.24$.  The roll results in
a change of the undulation plane by more than $90^\circ$.
\label{off plane turn stereoscopic}
}
\end{figure}

\begin{figure}
\includegraphics[width=0.46\textwidthF]{Fig11.eps}
\caption
{ Roll rotation rate per undulation period $\omega^\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) \textit{C}-shaped\ (q
$L=5.5$) and (b) \textit{W}-shaped\ (q $L=9$)
swimming nematodes.  \label{turning performance of swimming worm rpy}
}
\end{figure}

\begin{figure*}
\includegraphics[width=\textwidth]{Fig12.eps}
\caption
{Roll reorientation angle $\theta^\roll$ vs normalized
torsional mode length $q\Delta s^\torsional$.
Panels (a) and (c) show results for \textit{C}-shaped\ gait
($\wavevectorL=5.5$), and panels (b) and (d) show
results for \textit{W}-shaped\ ($q L=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)].
\label{off plane turns of swimming worm rpy}
}
\end{figure*}

\begin{figure}
\includegraphics[width=0.46\textwidthF]{Fig13.eps}
\caption { Roll rotation rate per undulation period $\omega^\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 \textit{C}-shaped\ and
\textit{W}-shaped\ locomotive gaits, as labeled.
\label{comparison crawling swimming turn} }
\end{figure}

\begin{figure}[t]
\includegraphics[width=0.45\textwidth]{Fig14.eps}
\caption
{ (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.
\label{explanation for drywet comparison}
}
\end{figure}

\ \clearpage
}

\subsection{Swimming in 3D}
\label{offplane turns swimming}

Observations of {\it 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.\ \ref{real
  celegans turning 3d} and Video S3 in ESM).

Video images of {\it 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 \eqref{single mode curvature} and \eqref{constant
  torsion} (see the contours depicted in Fig.\ \ref{curves with and no
  torsion}).

Figure \ref{real versus simulated worms 3d turn comparison} shows a
direct comparison between a video image of a swimming nematode and a
calculated torsional-roll posture with a curvature wavevector
$\wavevectorL=5.5$ (\textit{C}-shaped\ gait) and constant torsion
$\tau_0L=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.\ \eqref{constant torsion} and
\eqref{three mode torsional turn}, 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 \eqref{full
  serret-frenet equations}; 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 \ref{off plane turn stereoscopic} 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.\ \ref{turning performance
  of swimming worm rpy}--\ref{comparison crawling swimming turn}.  The
roll rotation rate per undulation period \eqref{definition of turn
  velocity} during a single-mode roll is plotted in Fig.\ \ref{turning
  performance of swimming worm rpy}, and the total roll reorientation
angle
\eqgone{
\begin{equation}
  \theta^\roll=\cos^{-1}
  \left[
    {\hat{\boldsymbol{b}}}(s^\torsional_1) \boldsymbol{\cdot}
    {\hat{\boldsymbol{b}}}(s^\torsional_2)
    \right]
  \label{definition of offplane turn angle}
\end{equation}
}
for a nematode performing a three-mode roll is plotted in
Fig.\ \ref{off plane turns of swimming worm rpy}.

The roll rotation rate $\omega^\roll$ depicted in Fig.\ \ref{turning
  performance of swimming worm rpy} is shown vs the dimensionless
torsion amplitude $\tau_0L$, for several values of the
normalized curvature amplitude $A / q$, for
\textit{C}-shaped\ and \textit{W}-shaped\ gaits.  Note that in swimming the roll
rotation rate depends on the nematode length whereas in no-slip
burrowing $\omega^\roll$ depends only on the torsion and curvature-wave
parameters.  Hence a different normalization is used in
Figs.\ \ref{turning performance of swimming worm rpy} (swimming) and
\ref{turn rate off plane dry} (burrowing).

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

The total reorientation angle $\theta^\roll$ during a three-mode
roll maneuver is plotted in Fig.\ \ref{off plane turns of swimming
  worm rpy} vs the dimensionless torsional mode length
$q\Deltas^\torsional$, normalized by the wavevector of
the curvature undulations.  The results are similar for the
\textit{C}-shaped\ and \textit{W}-shaped\ gaits and are only moderately affected by the
initial mode switching point.  The magnitude of the reorientation
angle $\theta^\roll$ is consistent with the single-mode roll
rotation rate presented in Fig.\ \ref{turning performance of swimming
  worm rpy}.

Our experimental observations of {\it 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.\ \ref{turning performance of swimming worm rpy} and \ref{off
  plane turns of swimming worm rpy}, a $90^\circ$ roll corresponds to
the body torsion of $\tau_0L\approx1$ for the roll mode
length ranging from half to one undulation period.  Shorter modes
require larger torsion (see, e.g.\ Video S4 in ESM).

\subsubsection{Hydrodynamic enhancement in swimming rolls}

To compare the swimming roll rotation rate $\omega^\roll$ shown in
Fig.\ \ref{turning performance of swimming worm rpy} with the
corresponding purely geometrical no-slip result
(Fig.\ \ref{turn rate off plane dry}), the swimming and no-slip
simulation data are replotted in Fig.\ \ref{comparison crawling
  swimming turn} 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 (\textit{C}-shaped\ gait) to even
three (\textit{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.\ \ref{explanation
  for drywet comparison}.  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^\roll$ have the same sign, they add up,
resulting in the enhanced reorientation.

We believe that swimming {\it 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.

\figsgone{
  \ \clearpage

\begin{figure}
\includegraphics[width=0.45\textwidth]{Fig15.eps}
\caption
{ (\textit{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 {\it C. elegans}\ in a 3D\ fluid.
\label{chemotaxis apparatus}
}
\end{figure}

\begin{figure}
\includegraphics[width=0.46\textwidthF]{Fig16.eps}
\caption
{ Chemotaxis index $\textit{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.
\label{chemotaxis saturation}
}
\end{figure}

\begingroup
  \begin{table*}
  \begin{ruledtabular}
   \begin{tabular}{ l l l }
   Parameter & Name & Value \\
   \hline
   $\alpha_s^\prime$ & Chemoattractant sensitivity parameter& $5\,\textrm{s}/\textrm{mm}$ \\
   $r_s^0$ & Turn initiation rate&$10\textrm{ turns}/\textrm{min}$ \\
   $r_s^{\high}$ &Upper turn-rate saturation limit& $15\textrm{ turns}/\textrm{min}$ \\
   $r_s^{\low}$ &Lower turn-rate saturation limit& $5\textrm{ turns}/\textrm{min}$ \\
   $r_\torsional^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_bL$ & Burrowing gait normalized wavevector & $9$ \\
   $q_sL$ & Swimming gait normalized wavevector & $5.5$ \\
   $L$ & Worm length & $1\,\textrm{mm}$\\
   \end{tabular}
   \end{ruledtabular}
   \caption
   {
      \label{chemosensing parameters}
      System parameters for 3D\ simulations of chemotaxis of
      burrowing and swimming {\it 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
      \cite{Suzuki-Thiele-Faumont-Ezcurra-Lockery-Schafer:2008}.}
  \end{table*}
  \endgroup

     \ \clearpage
}

\section{Application to Chemotaxis}
   \label{chemotaxis in 3d}

According to the locomotion model developed in Secs.\ \ref{crawling
  section} and \ref{maneuverability in fluids}, {\it 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,
{\it C. elegans}\ can efficiently chemotax in 3D.

It has been demonstrated that in 2D\ environments {\it C. elegans}\ employs
two chemotaxis mechanisms: klinotaxis (gradual turn) where the
nematode actively reorients itself towards the concentration gradient
direction \cite{Iino-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 \cite{Pierce-Shimomura-Morse-Lockery: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.

\subsection{Chemical sensing model}

\label{chemical sensing}

A moving {\it 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) \cite{Lockery:2011}.  To
model this process we use an approach that has been successfully
applied to describe bacterial chemotaxis
\cite{Celani-Vergassola:2010}.

Accordingly, we assume that the rate of turns is modulated by the
chemical signal
\eqgone{
\begin{equation}
  Q(t)=\int_0^t
  C(t-\chronos')
  M(\chronos') d \chronos',
  \label{chemo signal}
\end{equation}
}
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 \cite{Lockery:2011}, we
impose the ideal adaptation condition \cite{Celani-Vergassola:2010}
\eqgone{
\begin{equation}
  \int_0^\infty M(t) d t=0
  \label{chemo norm 1}.
\end{equation}
}
We also assume the normalization
\eqgone{
\begin{equation}
  \label{normalization of M}
\int_0^\infty tM(t) d t=1.
\end{equation}
}

In our numerical calculation we use the memory function defined by the
expression
\eqgone{
\begin{equation}
  M(t)=\epsilon^2 \mathrm{e}^{\epsilon t}
  \left[
    (\epsilon t)-\frac{1}{2}(\epsilon t)^2
    \right],
  \label{response function}
\end{equation}
}
where $\epsilon^{-1}$ is the decay time \cite{Celani-Vergassola:2010}.
The memory function \eqref{response function} is positive for
$0<t<\epsilon^{-1}$ and negative for $t>\epsilon^{-1}$;
the positive and negative regions have equal integrals.
The convolution integral \eqref{chemo signal} 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}$.

\subsection{Biased random walk}

\label{biased random motion}

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\approx1$, and switches randomly to a turning
mode at a turn initiation rate $r_s(t)$, which depends on
the chemotaxis signal \eqref{chemo signal}.  The nematode also
performs roll maneuvers initiated at a rate $r_\torsional(t)$.

To mimic the experimentally observed sigmoidal character of the turn
response to the chemoattractant concentration variation
\cite{Pierce-Shimomura-Morse-Lockery:1999}, we assume that the turn
initiation rate $r_s$ has a linear-response region and a saturation
cutoff,
\eqgone{
\begin{equation}
  r_s(t)=
  \left\{
  \begin{matrix}
    \hspace{-12mm} r_s^{\low}; \hspace{4mm} f(Q) \le
    r_s^{\low} \\
    f(Q); \hspace{4mm} r_s^{\low} <
    f(Q) < r_s^{\high}, \\
    \hspace{-11mm} r_s^{\high};
    \hspace{4mm} r_s^{\high} \le f(Q) \\
  \end{matrix}
  \right.
  \label{rate of turn}
\end{equation}
}
where
\eqgone{
\begin{equation}
  f(Q)=
  r_s^0 \left[
    1 - \alpha_s Q(t)
    \right].
  \label{rate of turn mid}
\end{equation}
}
Here $r_s^0$ is the turn rate in the absence of
chemoattractant, $\alpha_s$ is a chemoattractant sensitivity
parameter, and $r_s^{\high}$ and $r_s^{\low}$ 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,
\eqgone{
\begin{equation}
  r_\torsional(t)=r_\torsional^0.
  \label{rate of cork}
\end{equation}
}

In Sec.\ \ref{Chemotaxis in 3D} we show that assumptions
\eqref{rate of turn} and \eqref{rate of cork} are sufficient to
produce efficient 3D\ chemotaxis.

\subsection{The geometry of simulated chemotaxis assay}

In our 3D\ numerical simulations of chemotaxis of burrowing and
swimming nematodes [see Fig.\ \ref{chemotaxis apparatus}(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
\cite{%
  Bargmann-Horvitz:1991,%
  Nishino-Sato-Ito-Matsuura:2015}.

The chemotaxis efficiency is quantified using the chemotaxis index
\eqgone{
\begin{equation}
  \textit{CI}(t)=\frac{\noworms_\attr(t)-\noworms_\ctrl(t)}
           {\noworms_\attr(t)+\noworms_\ctrl(t)},
           \label{chemotaxis index cumulative}
\end{equation}
}
where $\noworms_\attr(t)$ and $\noworms_\ctrl(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 \cite{Nishino-Sato-Ito-Matsuura:2015}).  A sample trajectory
of a swimming {\it C. elegans}\ reaching the test region is shown in Fig.
\ref{chemotaxis apparatus}(b).

In our simulations we consider a system with the following parameters.
The spatial concentration of chemoattractant  is given by a 3D\
Gaussian distribution
\eqgone{
\begin{equation}
  \label{Gaussian}
  C\sim\exp\left(-r^2/2\sigma^2\right)
\end{equation}
}
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).

\subsection{Simulation results}
\label{Chemotaxis in 3D}

We have performed simulations for two cases: (a) burrowing
\textit{W}-shaped\ nematodes (no transverse slip) and (b) swimming
\textit{C}-shaped\ nematodes.  The nematode gait parameters and parameters of
the chemotaxis model are listed in Table \ref{chemosensing parameters}
(see Sec.\ \ref{Chemotaxis simulation details} for more details).
The gait parameter values are based on observations of crawling and
swimming {\it 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 \eqref{chemotaxis index cumulative}, shown in
Fig.\ \ref{chemotaxis saturation}, is significantly greater than zero,
indicating a clear trajectory bias towards the chemoattractant
concentration peak.  At long times the value of $\textit{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 \eqref{chemo signal} 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.

\ \clearpage

\section{Summary and discussion}
\label{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
{\it 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 {\it C. elegans}\ in 3D\ environments was experimentally
demonstrated
\cite{Pierce-Shimomura-Chen-Mun-Ho-Sarkis-McIntire:2008,%
  Beron-Vidal-Gadea-Cohn-Parikh-Hwang-Pierce-Shimomura:2015},
but, to our knowledge, the formulation presented here is the first
quantitative analysis that takes into account the gait geometry and
swimming mechanics.

\gone{Our numerical simulations of nematode trajectories in which the turn
parameters $s^\curvatural_1$ and $\Deltas^\curvatural$ are randomly
chosen yield many turns that produce only small directional changes.
In contrast, experiments show that most $\Omega$ postures adopted by
{\it C. elegans}\ result in relatively sharp turns.  Our calculations thus
suggest that the nematode is able to choose optimal turn initiation
and termination points.}

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 {\it C. elegans}.

Body movements of {\it 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 \cite{Altun-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
\cite{White-Southgate-Thomson-Brenner:1986,%
Altun-Hall:2009,%
Riddle-Blumenthal-Meyer-Priess: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 \textit{et al.\/}
\cite{Kwon-Pyo-Lee-Je: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.\ \ref{real versus simulated worms 3d turn comparison}, 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${}^\circ$ 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 {\it C. elegans}\ and its
important role in generation of a planar undulation wave have been
demonstrated
\cite{Wen-Po-Hulme-Chen-Liu-Kwok-Gershow-Leifer-Butler-Fang-Yen-Kawano-Schafer-Whitesides-Wyart-Chklovskii-Zhen-Samuel:2012}.
We thus conjecture that {\it 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
\cite{Kwon-Pyo-Lee-Je: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 \textit{et al.}
\cite{Wen-Po-Hulme-Chen-Liu-Kwok-Gershow-Leifer-Butler-Fang-Yen-Kawano-Schafer-Whitesides-Wyart-Chklovskii-Zhen-Samuel:2012},
``[t]he cellular economy of the {\it 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 {\it C. elegans}\ will help determine if such a
control system exists.

\ \clearpage
\section{Methods}
\label{experimental}

\subsection{Imaging}

Wild type {\it C. elegans}\ was used to record the 3D swimming in our
experiments.  Worms were cultured in NGM plate (Nematode Growth
Medium) seeded with bacteria \textit{E.\ coli} OP50 at
20${}^\circ$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${}^\circ$C.

\subsection{Chemotaxis simulation details}
\label{Chemotaxis simulation details}

In our simulations, the planar turns and rolls are initiated randomly
at the rates $r_s$ and $r_\torsional$, 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 \ref{chemosensing parameters}.  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 $\widebar{A}_2/q=1.67$ and
$\sigma_{A_2/q}=0.28$; turning-mode length
$\widebar{\Deltas}^\curvatural/L=0.5$ and
$\sigma_{\Deltas^\curvatural/L}=0.167$; normalized torsion
$\widebar{\tau}_0L=2$ and
$\sigma_{\tau_0L}=1$; and torsional-mode length
$\widebar{\Deltas}^\torsional/L=0.5$ and
$\sigma_{\Deltas^\torsional/L}=0.1$.  The parameter
distributions are the same for burrowing and swimming.

The turning rates listed in Table \ref{chemosensing parameters}
correspond to the directional persistence length of burrowing
trajectories $l_p\approx12L$ (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.

{ \small
\paragraph*{\bf \small {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.
\paragraph*{\bf \small 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.
\paragraph*{\bf \small 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.
\paragraph*{\bf \small Competing interests.}
We have no competing interests.
}

\ \clearpage


\section{Introduction}

Equilibrium and nonequilibrium phase transitions in ensembles of
complex networked systems have been a subject of intense research in
the field of statistical physics and many other disciplines
\cite{PRP06000175,PRP08000093,RMP08001275,PRP2014}. Owing to the
inherent randomness and heterogeneity in the interacting patterns,
phase transitions on complex networks are drastically different from
those on regular lattices in Euclidean space. Examples range from
the anomalous behavior of Ising model
\cite{PHA02000260,PLA02000166,PRE02016104,EPB02000191,PhysRevLett.104.218701}
to a vanishing percolation threshold
\cite{PhysRevLett.85.4626,PRL00005468} and the absence of epidemic
thresholds that separate healthy and endemic phases
\cite{PRL01003200,PhysRevLett.90.028701,PhysRevLett.111.068701} as
well as explosive emergence of phase transitions
\cite{Science323.1453,PhysRevLett.103.255701,PhysRevLett.103.168701,
PhysRevLett.103.135702,PhysRevLett.105.255701,PhysRevLett.106.225701,Science333.322,
Nature2010,PhysRevLett.105.048701,PhysRevLett.107.195701,
PhysRevLett.106.128701,PhysRevLett.108.168702,NatPhys2014}. So far,
unveiling the relationship between the onset of phase transitions
and the topology of the underlying networks is still a topic of
considerable attention.

The majority-voter (MV) model is a simple nonequilibrium Ising-like
system with up-down symmetry that presents an order-disorder phase
transition at a critical value of noise \cite{RMP09000591}. Since
Oliveira pointed out that the MV model on a square lattice belongs
to the universality class of the equilibrium Ising model
\cite{JSP1992}, the model has been extensively studied in the
context of complex networks, including random graphs
\cite{PhysRevE.71.016123,PA2008}, small world networks
\cite{PhysRevE.67.026104,IJMPC2007,PA2015}, scale-free networks
\cite{IJMPC2006(1),IJMPC2006(2)}, and some others
\cite{PhysRevE.75.061110,PhysRevE.77.051122,PhysRevE.81.011133,PA2011,PhysRevE.86.041123,PhysRevE.89.052109}.
These results showed that the critical exponents are generally
dependent on the underlying interacting substrates. However, the
studies of pervious works are mainly based on numerical simulations.
Especially, the analytical determination of the critical noise is
still lacking at present.

For this purpose, in this paper we employ the heterogeneous
mean-field theory to derive the rate equation for governing the MV
model's dynamics on undirected networks. According to linear
stability analysis, we determine the critical point of noise $f_c$,
the onset of an order-disorder phase transition in the limit of
infinite size networks $N \rightarrow \infty$. The analytical result
shows that $f_c$ is related to the ratio of the first moment to the
${3 \mathord{\left/ {\vphantom {3 2}} \right.
\kern-\nulldelimiterspace} 2}$ order moment of degree distribution.
Furthermore, we derive the Langevin equation to study the effect of
stochastic fluctuation on finite size networks. By solving the
potential of the corresponding Fokker-Planck equation, we calculate
the susceptibility as a function of noise and determine the
effective critical noise $f_c(N)$ on finite size networks at which
the susceptibility is maximal. We find that the difference
$f_c-f_c(N)$ decays in power-law ways with $N$. Extensive Monte
Carlo (MC) simulations are performed on diverse network types to
validate the theoretical results.

\section{Model}
We consider the MV model with noise on complex networks defined by a
set of spin variables $\{\sigma_i\}$ $(i = 1, \ldots ,N)$, where
each spin is associated to one node of the underlying network and
can take the values $\pm 1$. The system evolves as follows: for each
spin $i$, we first determine the majority spin of $i'$s
neighborhood. With probability $f$ the node $i$ takes the opposite
sign of the majority spin, otherwise it takes the same spin as the
majority spin. The probability $f$ is called the noise parameter and
plays a similar role of temperature in equilibrium spin systems. In
this way, the single spin flip probability can be written as
\begin{eqnarray}
w({\sigma _i}) = \frac{1}{2}\left[ {1 - \left( {1 - 2f}
\right){\sigma _i}S\left( {\sum\limits_j {{a_{ij}}{\sigma _j}} }
\right)} \right], \label{eq1}
\end{eqnarray}
where $S(x)=sgn(x)$ if $x\neq0$ and $S(0)=0$. In the latter case the
spin $\sigma_i$ is flipped to $\pm 1$ with probability equal to ${1
\mathord{\left/ {\vphantom {1 2}} \right. \kern-\nulldelimiterspace}
2}$. The elements of the adjacency matrix of the underlying network
are defined as $a_{ij}=1$ if nodes $i$ and $j$ are connected and
$a_{ij}=0$ otherwise. In the case $f=0$, the majority-vote model is
identical to the zero temperature Ising model
\cite{PNAS2005,JSM2006}.

\section{Results}
To proceed a mean-field treatment, we first define $q_k$ as the
probability that a node of degree $k$ is in $+1$ state, and $Q$ as
the probability that for any node in the network, a randomly chosen
nearest neighbor node is in $+1$ state. Furthermore, for any node
the probability that a randomly chosen nearest neighbor node has
degree $k$ is ${{kP(k)} \mathord{\left/
 {\vphantom {{kP(k)} {\left\langle k \right\rangle }}} \right.
 \kern-\nulldelimiterspace} {\left\langle k \right\rangle }}$, where $P(k)$ is degree distribution
defined as the probability that a node chosen at random has degree
$k$ and $\left\langle k \right\rangle $ is the average degree
\cite{RMP08001275}. It is supposed to be reasonable only in networks
without degree correlation. The probabilities $q_k$ and $Q$ satisfy
the relation
\begin{eqnarray}
Q = \sum\limits_k {kP(k)} {{{q_k}} \mathord{\left/
 {\vphantom {{{q_k}} {\left\langle k \right\rangle }}} \right.
 \kern-\nulldelimiterspace} {\left\langle k \right\rangle }}. \label{eq2}
\end{eqnarray}
Thus, we can write rate equations for $q_k$ as
\begin{eqnarray}
{{\dot q}_k} &=&  - {q_k}(1 - {\psi _k}) + (1 - {q_k}){\psi _k}
\nonumber \\ &=&  - {q_k} + {\psi _k}, \label{eq3}
\end{eqnarray}
where ${\psi _k}$ is is the probability that a node of degree $k$
takes the $+1$ value, which can be expressed as
\begin{eqnarray}
{\psi _k(Q)} = (1 - f){\varphi _k}(Q) + f({1-\varphi _k}(Q)).
\label{eq4}
\end{eqnarray}
Here, ${\varphi _k}(Q)$ is the probability that a node of degree $k$
with $+1$ state takes the majority rule, which can be written by a
binomial distribution,
\begin{eqnarray}
{\varphi _k}(Q) = \sum\limits_{n = \left\lceil {k/2} \right\rceil
}^k {\left(1 - \frac{1}{2}{\delta _{n,k/2}}\right) C_k^n} {Q^n}{(1 -
Q)^{k - n}}, \label{eq5}
\end{eqnarray}
where $\left\lceil \cdot \right\rceil$ is the ceiling function,
$\delta$ is the Kronecker symbol, and $C_k^n = {{k!} \mathord{\left/
 {\vphantom {{k!} {[n!(k - n)!])}}} \right.
 \kern-\nulldelimiterspace} {[n!(k - n)!]}}$ are the binomial coefficients. By introducing Eq.(\ref{eq3}) into Eq.(\ref{eq2}) we obtain a closed rate equation for
the quantity $Q$,
\begin{eqnarray}
\dot Q =  - Q + \Psi (Q), \label{eq6}
\end{eqnarray}
where
\begin{eqnarray}
\Psi (Q) = \sum\limits_k {kP(k)} {{{\psi _k}(Q)} \mathord{\left/
 {\vphantom {{{\psi _k}(Q)} {\left\langle k \right\rangle }}} \right.
 \kern-\nulldelimiterspace} {\left\langle k \right\rangle }}. \label{eq7}
\end{eqnarray}

\begin{figure}
\centerline{\includegraphics*[width=1.0\columnwidth]{fig1.eps}}
\caption{(color online). Graphic demonstration of steady state
solutions of $Q$. When $f$ is less than a critical value $f_c$,
there are three solutions. One is at $Q =\frac{1}{2}$ corresponding
to a disordered phase and the other two correspond to two symmetric
ordered phases. For $f>f_c$ , there is only one solution at $Q
=\frac{1}{2}$. \label{fig1}}
\end{figure}

In the steady state $\dot Q = 0$, we have $Q_s = \Psi (Q_s)$.
Fig.\ref{fig1} shows that the two typical examples of graphic
solutions of $Q_s$. One can easily find that a trivial stationary
solution, $Q_s=\frac{1}{2}$, always exists irrespective of the value
of $f$ (corresponding to a disordered phase $\left\langle {{\sigma
_i}} \right\rangle  = 0$), as ${\varphi _k}(\frac{1}{2}) =
\frac{1}{2}$ and ${\Psi}(\frac{1}{2}) = \frac{1}{2}$. However, the
other two solutions are possible if $f$ is less than a critical
value $f_c$, and they represent the existence of two ordered phases
with up-down symmetry. Therefore, the critical noise $f_c$ is
determined by the condition that the derivation of $\Psi(Q)$ with
$Q$ equals to one at $f=f_c$, i.e.,
\begin{eqnarray}
{\left. {\frac{{d\Psi (Q)}}{{dQ}}} \right|_{Q = \frac{1}{2}}} = 1.
\label{eq8}
\end{eqnarray}
To do this, we rewrite approximately Eq.(\ref{eq5}) as
\begin{eqnarray}
{\varphi _k}(Q) = {\varphi _k}(\frac{1}{2} + y) = \frac{1}{2} +
\frac{1}{2}erf\left( {y\sqrt {2k} } \right),\label{eq9}
\end{eqnarray}
where $erf(x)$ is the error function. Note that this approximation
is plausible for large values of $k$ as the binomial distribution
can be approximated by a normal distribution and the sum over $n$ in
Eq.(\ref{eq5}) can be substituted by an integral \cite{JSM2006}. The
derivation of $\Psi(Q)$ with $Q$ can be expressed analytically by
\begin{eqnarray}
{\left. {\frac{{d\Psi (Q)}}{{dQ}}} \right|_{Q = \frac{1}{2}}} &=& {\left. {\sum\limits_k {\frac{{kP(k)}}{{\left\langle k \right\rangle }}(1 - 2f)} \frac{{d\varphi (\frac{1}{2} + y)}}{{dy}}} \right|_{y = 0}}  \nonumber \\
& = &{\left. {\sum\limits_k {\frac{{kP(k)}}{{\left\langle k \right\rangle }}(1 - 2f)} \sqrt {\frac{{2k}}{\pi }} {e^{ - 2k{y^2}}}} \right|_{y = 0}} \nonumber \\
 &= &(1 - 2f)\sqrt {\frac{2}{\pi }} \frac{{\left\langle {{k^{{3
\mathord{\left/
 {\vphantom {3 2}} \right.
 \kern-\nulldelimiterspace} 2}}}} \right\rangle }}{{\left\langle k \right\rangle
 }},\label{eq10}
\end{eqnarray}
where $\left\langle {{k^n}} \right\rangle  = \sum\nolimits_k
{{k^n}P(k)}$ is the $n$th moment of degree distribution. Inserting
Eq.(\ref{eq10}) into Eq.(\ref{eq8}), we arrive at the analytical
expression of $f_c$,
\begin{eqnarray}
{f_c} = \frac{1}{2} - \frac{1}{2}\sqrt {\frac{\pi }{2}}
\frac{{\left\langle k \right\rangle }}{{\left\langle {{k^{{3
\mathord{\left/
 {\vphantom {3 2}} \right.
 \kern-\nulldelimiterspace} 2}}}} \right\rangle }}.\label{eq11}
\end{eqnarray}

\begin{figure*}
\centerline{\includegraphics*[width=2.3\columnwidth]{fig2.EPS}}
\caption{(color online). The critical noise $f_c$ obtained from the
theory (lines) and MC simulation (symbols) on (a) RkRN, (b) ERRN,
and (b-c) SFN.\label{fig2}}
\end{figure*}

To validate the theoretical results on $f_c$, we shall consider the
three network types: the random $k$-regular networks (RkRN) and
Erd\"os-R\'enyi random networks (ERRN), as the representations of
degree homogeneous networks, and scale-free networks (SFN) as the
representations of degree heterogeneous networks. For RkRN, each
node has the same degree $k$ and degree distribution follows the
delta-like function $P(k) = \delta (k)$, and thus Eq.(\ref{eq11})
can be reduced to
\begin{eqnarray}
f_c^{RkRN} = \frac{1}{2} - \frac{1}{2}\sqrt {\frac{\pi }{2}}
\frac{1}{{\sqrt { k  } }}.\label{eq12}
\end{eqnarray}
From Eq.(\ref{eq12}), we see that $f_c$ is a decreasing function of
$k$. For $k=N-1\rightarrow \infty$, $f_c=\frac{1}{2}$ and RkRN thus
become the globally connected networks. For ERRN, degree
distribution follows Poisson $P(k) = {\left\langle k \right\rangle
^k}{{{e^{ - \left\langle k \right\rangle }}} \mathord{\left/
{\vphantom {{{e^{ - \left\langle k \right\rangle }}} {k!}}} \right.
\kern-\nulldelimiterspace} {k!}}$ with the average degree
$\left\langle k \right\rangle$, and the theoretical value of $f_c$
for ERRN can be numerically calculated according to Eq.(\ref{eq11}).

For SFN, degree distribution follows a pow-law function $P(k) \sim
{k^{ - \gamma }}$, with degree exponent $\gamma>2$. In the
thermodynamic limit $N \to \infty$, ${\left\langle {{k^{{3
\mathord{\left/ {\vphantom {3 2}} \right.
 \kern-\nulldelimiterspace} 2}}}} \right\rangle }$ diverges for $\gamma \leq {5 \mathord{\left/
{\vphantom {5 2}} \right. \kern-\nulldelimiterspace} 2}$, such that
the critical noise becomes $f_c=\frac{1}{2}{\kern 1pt} {\kern 1pt}$
according to Eq.(\ref{eq11}). For $\gamma > {5 \mathord{\left/
{\vphantom {5 2}} \right. \kern-\nulldelimiterspace} 2}$, both
$\left\langle k \right\rangle$ and ${\left\langle {{k^{{3
\mathord{\left/ {\vphantom {3 2}} \right.
 \kern-\nulldelimiterspace} 2}}}} \right\rangle }$ are finite in the limit of $N \to \infty$, and they
are $\gamma-$dependent given by $\left\langle k \right\rangle  =
{{(\gamma  - 1){k_0}} \mathord{\left/ {\vphantom {{(\gamma  -
1){k_0}} {(\gamma  - 2}}} \right. \kern-\nulldelimiterspace}
{(\gamma  - 2}}{\kern 1pt} )$ and $\left\langle {{k^{{3
\mathord{\left/
 {\vphantom {3 2}} \right.
 \kern-\nulldelimiterspace} 2}}}}
\right\rangle  = {{(\gamma  - 1)k_0^{{3 \mathord{\left/
 {\vphantom {3 2}} \right.
 \kern-\nulldelimiterspace} 2}}} \mathord{\left/
{\vphantom {{(\gamma  - 1)k_0^{1.5}} {(\gamma  - 2.5}}} \right.
\kern-\nulldelimiterspace} {(\gamma  - {5 \mathord{\left/
 {\vphantom {3 2}} \right.
 \kern-\nulldelimiterspace} 2}}}{\kern 1pt} )$ with
$k_0$ being the minimal node degree. By the above analysis, we
immediately obtain the critical noise $f_c$ for SFN
\begin{equation}
f_c^{SFN} = \left\{
\begin{aligned}
&\frac{1}{2} - \frac{1}{2}\sqrt {\frac{\pi }{2}} \frac{{\gamma  - {5
\mathord{\left/ {\vphantom {5 2}} \right. \kern-\nulldelimiterspace}
2}}}{{\gamma  - 2}}\frac{{1{\kern 1pt} {\kern 1pt} {\kern 1pt}
{\kern 1pt} }}{{\sqrt {{k_0}} }},  & \gamma  > {5 \mathord{\left/
{\vphantom {5 2}}
\right. \kern-\nulldelimiterspace} 2}{\kern 1pt} {\kern 1pt} \\
&\frac{1}{2},  & \gamma  \le {5 \mathord{\left/ {\vphantom {5 2}}
\right. \kern-\nulldelimiterspace} 2}
\end{aligned}
\right.\label{eq13}
\end{equation}

We firstly generate the networks according to the Molloy-Reed model
\cite{RSA95000161}: each node is assigned a random number of stubs
$k$ that is drawn from a given degree distribution. Pairs of
unlinked stubs are then randomly joined. We then run the standard MC
simulation: at each MC step, each node is firstly randomly chosen
once on average and then make an attempt to flip spin with the
probability according to Eq.(\ref{eq1}).

In order to numerically obtain the $f_c$, we need to calculate the
Binder's fourth-order cumulant $U$, defined as
\begin{eqnarray}
U = 1 - \frac{{\left[ {\left\langle {{m^4}} \right\rangle }
\right]}}{{3{{\left[ {\left\langle {{m^2}} \right\rangle }
\right]}^2}}},\label{eq14}
\end{eqnarray}
where $m={{\sum\nolimits_{i = 1}^N {{\sigma _i}} } \mathord{\left/
 {\vphantom {{\sum\nolimits_{i = 1}^N {{\sigma _i}} } N}} \right.
 \kern-\nulldelimiterspace} N}$ is the average magnetization per node, $\left\langle  \cdot  \right\rangle$ denotes time averages
taken in the stationary regime, and $\left[  \cdot  \right]$ stands
for the averages over different network configurations. The critical
noise $f_c$ is estimated as the point where the curves $U\sim f$ for
different network sizes $N$ intercept each other. In our
simulations, $f_c$ is determined by five different network sizes:
$N=500$, $1000$, $2000$, $5000$ and $10000$.

For comparison, in Fig.\ref{fig2} we plot the $f_c$ obtained from
the theoretical prediction (lines) and the MC simulation (symbols),
respectively. In Fig.\ref{fig2}(a) and Fig.\ref{fig2}(b), we show
the results on RkRN and on ERRN, respectively. In
Fig.\ref{fig2}(c-d), we show the results on SFN and plot the $f_c$
as a function of $k_0$ for some fixed $\gamma$ in Fig.\ref{fig2}(c)
and of $\gamma$ for some fixed $k_0$ in Fig.\ref{fig2}(d). It is
clearly observed that for large values $k$ there are an excellent
agreements between the theory and simulation. However, for
relatively small $k$ the used approximation in Eq.(\ref{eq9}) is not
very valid, such that the discrepancy between them exists.

So far, we have obtained the analytical expression of $f_c$ and
confirmed its validity by performing MC simulations on different
networks. The expression is only valid for infinite size networks $N
\rightarrow \infty$ where the finite-size fluctuation is ignored.
For finite size networks, the fluctuation is unavoidable and the
actual phase transition never happens. However, one can define an
effective critical noise $f_c(N)$ at which the susceptibility (the
variance of an order parameter) is maximal. Obviously, $f_c(N)$ is
size-dependent and recovers $f_c$ in the limit of $N \rightarrow
\infty$. To get $f_c(N)$, we will derive the fluctuation-driven
Langevin equation for $Q$ \cite{PhysRevE.79.036110,JSM2009}:
\begin{eqnarray}
\dot Q = \mathbb{E}[\Delta Q] + \sqrt {\mathbb{V}[\Delta Q]} \xi (t)
\label{eq15}
\end{eqnarray}
where $\mathbb{E}[\Delta Q]$ and $\mathbb{V}[\Delta Q]$ are the mean
value and the variance of the variation of $Q$, respectively, and
$\xi(t)$ is a Gaussian white noise satisfying $\left\langle {\xi
(t)} \right\rangle  = 0$ and $\left\langle {\xi (t)\xi (t')}
\right\rangle  = \delta (t - t')$. For the present model,
$\mathbb{E}[\Delta Q]$ and $\mathbb{V}[\Delta Q]$ can be computed as
\begin{eqnarray}
  \mathbb{E}[\Delta Q]&=&N\sum\limits_k {P(k){q_k}[1 - {\psi _k}(Q)]\left( { - {k \over {N\left\langle k \right\rangle }}} \right)}  \nonumber\\ &&+ N\sum\limits_k {P(k)(1 - {q_k}){\psi _k}(Q)\left( { {k \over {N\left\langle k \right\rangle }}} \right)} \nonumber \\
   &=&  - Q + \Psi (Q),\label{eq16}
\end{eqnarray}
and
\begin{eqnarray}
\mathbb{V}[\Delta Q] &=& N\sum\limits_k {P(k){q_k}[1 - {\psi
_k}(Q)]} {\left( { - {k \over {N\left\langle k \right\rangle }}}
\right)^2} \nonumber \\&&+ N\sum\limits_k {P(k)(1 - {q_k}){\psi
_k}(Q){{\left( {{k \over {N\left\langle k \right\rangle }}}
\right)}^2}}  \nonumber \\ &=& \sum\limits_k {{{{k^2}P(k)} \over
{N{{\left\langle k \right\rangle }^2}}}} \left[ {q{_k}(1 - {\psi
_k}(Q)) + (1 - {q_k}){\psi _k}(Q)} \right] \nonumber \\ \label{eq17}
\end{eqnarray}
Equation (\ref{eq17}) is not yet a closed equation for $Q$ because
the diffusion term $\mathbb{V}[\Delta Q]$ involve degree-dependent
quantities $q_k$. To close it, we use the quasi-static approximation
obtained from the rate equations (\ref{eq3}) imposing $\dot
q_k\simeq0$, i.e., $q_k\simeq\psi_k(Q)$
\cite{PhysRevE.79.036110,JSM2009}. The approximation assumes that
$Q(t)$ varies much slowly with respect to the dynamics of the
microscopic degrees of freedom $q_k(t)$. By the approximation,
Eq.(17) becomes
\begin{eqnarray}
\mathbb{V}[\Delta Q] = {2 \over {N{{\left\langle k \right\rangle
}^2}}}\sum\limits_k {{k^2}P(k)} {\psi _k}(Q)(1 - {\psi _k}(Q))
\label{eq18}
\end{eqnarray}
Therefore, we obtain the fluctuation-driven Langevin equation with a
closed form,
\begin{eqnarray}
\dot Q =  - Q + \Psi (Q) + \sqrt{2D(Q)} \xi (t) \label{eq19}
\end{eqnarray}
with multiplicative noise $D(Q)={1 \over {N{{\left\langle k
\right\rangle }^2}}}\sum\nolimits_k {{k^2}P(k)} {\psi _k}(Q)(1 -
{\psi _k}(Q))$. Clearly, in the limit of $N\rightarrow\infty$, the
fluctuation term $D(Q)\rightarrow0$, and Eq.(\ref{eq19}) thus
recovers to the mean-field equation derived in Eq.(\ref{eq6}).

Furthermore, let $P(Q,t)$ denote the probability density
distribution of $Q$ at time $t$. Then, the Fokker-Planck equation of
$P(Q,t)$ corresponding to Eq.(\ref{eq19}) can be given by
\begin{eqnarray}
\frac{{\partial P(Q,t)}}{{\partial t}} =  &-& \frac{\partial
}{{\partial Q}}\left[ { - Q + \Psi (Q) + \sqrt {D(Q)D'(Q)} }
\right]P(Q,t) \nonumber \\ &+& \frac{{{\partial ^2}}}{{{\partial
^2}Q}} {D(Q)} P(Q,t) \label{eq20}
\end{eqnarray}
The stationary distribution is $P(Q)=C e^{U_{FP}(Q)}$ where $C$ is
the normalized constant and
\begin{eqnarray}
{U_{FP}}(Q) = \frac{1}{2}\ln \left[ {D(Q)} \right] - \int_{}^Q
{\frac{{ - S + \Psi (S)}}{{D(S)}}} dS \label{eq21}
\end{eqnarray}
is called the potential of the Fokker-Planck equation
\cite{Risken1992}.

The critical noise $f_c(N)$ for finite size networks is determined
using the modified susceptibility $\chi '$ defined as
\begin{eqnarray}
\chi ' = N\left[ {\left\langle {{y^2}} \right\rangle  -
{{\left\langle |y| \right\rangle }^2}} \right], \label{eq22}
\end{eqnarray}
where $y=Q-\frac{1}{2}$, $\left\langle {{|y|}} \right\rangle$ and
$\left\langle {{y^2}} \right\rangle$ are calculated by the integrals
\begin{eqnarray}
\left\langle {\left| y \right|} \right\rangle  = \int_0^1 {\left| {Q
- \frac{1}{2}} \right|P(Q)dQ} \label{eq23}
\end{eqnarray}
and
\begin{eqnarray}
\left\langle {{y^2}} \right\rangle  = \int_0^1 {{{\left( {Q -
\frac{1}{2}} \right)}^2}P(Q)dQ}, \label{eq24}
\end{eqnarray}
respectively. We expect that $\chi'$ have a peak at $f=f_c(N)$ that
diverges and $f_c(N)$ converges to $f_c$ when the network size
increases.

\begin{figure}
\centerline{\includegraphics*[width=1.0\columnwidth]{fig3.EPS}}
\caption{(color online). The susceptibility $\chi '$ as a function
of noise $f$ for some different network sizes: $N=1000, 2000, 5000,
10000$. (a) RkRN: $k=20$, (b) ERRN: $\left\langle k
\right\rangle=20$, (c) SFN: $\gamma=2.4$ and $k_0=20$, and (d)
$\gamma=3$ and $k_0=20$. Lines and symbols indicate the results from
theory and MC simulation, respectively. \label{fig3}}
\end{figure}

\begin{figure}
\centerline{\includegraphics*[width=1.0\columnwidth]{fig4.eps}}
\caption{(color online). Log-log plot of the differences of the
effective critical noise $f_c$ for finite size networks with $f_c$,
$\Delta f(N) = {f_c} - {f_c}(N)$, with $N$ on (a) RkRN: $k=20$,
ERRN: $\left\langle k \right\rangle=20$, and (b) SFN: $k_0=20$.
$\Delta f(N)$ power law decays with $N$: $\Delta f(N) \sim {N^{ -
\nu }}$. Lines and symbols indicate the results from theory and MC
simulation, respectively. The inset of Fig.4(b) shows the exponent
$\nu$ as a function of $\gamma$. \label{fig4}}
\end{figure}

In Fig.\ref{fig3}, we show that $\chi '$ as a function of noise $f$
for some different network sizes $N$ on RkRN, ERRN, and SFN. For
comparison, the values of $\chi '$ obtained from the theoretical
calculations for Eq.(\ref{eq22}) and from MC simulations are
indicated by the lines and symbols, respectively. The theoretical
calculations can give a well prediction for simulation results. As
mentioned above, the point corresponding to the maximal $\chi '$
lies in the effective critical noise $f_c(N)$ for finite size
networks. In addition, we should note that we find that from MC
simulations the commonly used susceptibility $\chi=N\left[
{\left\langle {{m^2}} \right\rangle - {{\left\langle |m|
\right\rangle }^2}} \right]$ and our used susceptibility $\chi'$
defined in Eq.(\ref{eq22}) share the same locations where they are
maximal (results not shown here).

In Fig.\ref{fig4}, we plot the difference $\Delta f(N) = {f_c} -
{f_c}(N)$ as a function of $N$ in double logarithm coordinates. The
lines and symbols also indicate the results of theoretical
calculations from Eq.(\ref{eq22}) and MC simulations, respectively.
As pointed out by many previous studies, $\Delta f(N)$ scale with
$N$ in a power-law way: $\Delta f(N) \sim {N^{ - \nu }}$. With the
increment of $N$, $\Delta f(N)$ decreases and tends to zero in the
limit $N\rightarrow \infty$, recovering the result of
Eq.(\ref{eq11}). For RkRN and ERRN, we find that the exponents $\nu$
are independent of $k$ and $\left\langle k \right\rangle$. For SFN,
$\nu$ is also almost independent of $k_0$ but is an increasing
function of $\gamma$ (see the inset of Fig.\ref{fig4}(b)).

\section{Conclusions}
In conclusion, we have used heterogeneous mean-field theory to
derive the rate equation of an order parameter $Q$ for the MV model
defined on complex networks. By the linear stability analysis, we
have analytically obtained the critical noise $f_c$ at which the
order-disorder phase transition takes place in the limit of infinite
size networks. We find that that $f_c$ is determined by both the
first and ${3 \mathord{\left/ {\vphantom {3 2}} \right.
\kern-\nulldelimiterspace} 2}$ order moments of degree distribution
of the underlying networks. Moreover, we have incorporated the
effect of stochastic fluctuation on finite size networks via the
derivation of the Langevin equation of $Q$. By solving the
corresponding Fokker-Planck equation, we have obtained the effective
critical noise $f_c(N)$ where the susceptibility is maximal. The
results show that $f_c-f_c(N)$ power law decreases with $N$ and
reduces to zero in the limit of $N\rightarrow \infty$. To validate
the theoretical results, we have performed the extensive MC
simulations on RkRN, ERRN, and SFN. There are excellent agreement
between the theory and simulations. However, our theory does not
perform well on very sparse networks. Therefore, in the future it
will be desirable to develop high order theories (such as pair
approximation
\cite{PNAS2002,EPL2009,NJP2014,PhysRevLett.107.068701}) to obtain
more accurate estimation of the critical point of the networked MV
model.

\begin{acknowledgments}
We acknowledge supports from the National Science Foundation of
China (11205002, 61473001, 11475003, 21125313) and ``211 project" of
Anhui University (02303319-33190133).
\end{acknowledgments}