A stochastic neuronal model predicts random search behaviors at multiple spatial scales in C. elegans
Abstract
Random search is a behavioral strategy used by organisms from bacteria to humans to locate food that is randomly distributed and undetectable at a distance. We investigated this behavior in the nematode Caenorhabditis elegans, an organism with a small, welldescribed nervous system. Here we formulate a mathematical model of random search abstracted from the C. elegans connectome and fit to a largescale kinematic analysis of C. elegans behavior at submicron resolution. The model predicts behavioral effects of neuronal ablations and genetic perturbations, as well as unexpected aspects of wild type behavior. The predictive success of the model indicates that random search in C. elegans can be understood in terms of a neuronal flipflop circuit involving reciprocal inhibition between two populations of stochastic neurons. Our findings establish a unified theoretical framework for understanding C. elegans locomotion and a testable neuronal model of random search that can be applied to other organisms.
https://doi.org/10.7554/eLife.12572.001eLife digest
An animal’s ability to rapidly and efficiently locate new sources of food in its environment can mean the difference between life and death. As a result, animals have evolved foraging strategies that are adapted to the distribution and detectability of food sources. Organisms ranging from bacteria to humans use one such strategy, called random search, to locate food that cannot be detected at a distance and that is randomly distributed in their surroundings. The biological mechanisms that underpin random search are relatively well understood in singlecell organisms such as bacteria, but this information tells us little about the mechanisms that are used by animals, which use their nervous system to control their foraging behavior.
Roberts et al. have now investigated the biological basis for random search behavior in a tiny roundworm called Caenorhabditis elegans. This worm forages for pockets of bacteria in decaying plant matter and has a simple and wellunderstood nervous system. Roberts et al. used information on how the cells in this worm’s nervous system connect together into socalled “neural circuits” to generate a mathematical model of random searching.
The model revealed that the worm’s neural circuitry for random searching can be understood in terms of two groups of neuronlike components that switch randomly between “ON” and “OFF” states. While one group promotes forward movement, the other promotes backward movement, which is associated with a change in search direction. These two groups inhibit each other so that only one group usually is active at a given time. By adjusting this model to reproduce the behavioral records of real worms searching for food, Roberts et al. could predict the key neuronal connections involved. These predictions were then confirmed by taking electrical recordings from neurons. The model could also account for the unexpected behavioral effects that are seen when a neuron in one of these groups was destroyed or altered by a genetic mutation.
These findings thus reveal a biological mechanism for random search behavior in worms that might operate in other animals as well. The findings might also provide future insight into the neural circuits involved in sleep and wakefulness in mammals, which is organized in a similar way.
https://doi.org/10.7554/eLife.12572.002Introduction
Random search is an evolutionarily ancient set of foraging strategies that evolved as an adaptation to environments in which prey items are undetectable at a distance and occur at unpredictable locations. Rather than attempting to exhaustively search a region of interest, the organism samples the environment at randomly selected points. This is achieved by executing a series of straightline movements, called 'runs,' terminated at random intervals by sampling episodes during which the organism may or may not find prey. Sampling ends in a reorientation event, called a 'turn,' such that the next run is usually in a different direction from the preceding one. In optimal random foraging strategies the probability distribution of run length is matched to the statistical distribution of isolated food patches or prey items (Viswanathan, 2011), with power law distributions predominating when resources are sparsely distributed and exponential distributions predominating when resources are densely distributed (Humphries et al., 2010; Sims et al., 2012; Humphries et al., 2012).
Random search has been documented in a wide range of species including microorganisms, nematodes, insects, mollusks, fish, birds, and mammals including humans (Viswanathan, 2011; Berg and Brown, 1972; PierceShimomura et al., 1999). In humans this strategy is observed in diverse contexts, from traditional huntergatherer societies (Brown et al., 2007; Humphries and Sims, 2014) to technologically enhanced fishing industries (Bertrand et al., 2007). The formal similarities between random search across widely diverse phyla and spatial scales (Viswanathan, 2011) may point to a common mechanism, even in organisms that are highly cognitive. Despite the universality of random search, little is known about its neuronal basis, in part because of the difficulty of recording and manipulating activity in the brain of an unrestrained animal while it explores a large region of space.
The relatively small spatial scale of random search behavior in C. elegans, coupled with the simplicity of its nervous system, provides a unique opportunity to identify the neuronal basis of random search in this species. To the unaided eye, C. elegans search behavior consists of forward runs, each terminated after a variable interval by a briefer period of reverse locomotion, which is also variable in duration (PierceShimomura et al., 1999; Zhao et al., 2003; Wakabayashi et al., 2004), with apparently stochastic switching between these two behavioral states. Reversals are followed by resumption of forward movement that frequently begins with a deep body bend. These bends are highly variable in amplitude and lead to movement in a new direction. Thus, the sequence reverse–forward–deep bend, called a 'pirouette' (PierceShimomura et al., 1999) is the fundamental turning event in C. elegans random search, with functional analogies to tumbles in bacterial chemotaxis (Berg and Brown, 1972). Careful inspection reveals a third state, called “pause,” in which locomotion ceases for a fraction of a second or more (Croll, 1975; Shingai, 2000; Stephens et al., 2008; Rakowski et al., 2013; Salvador et al., 2014). Thus, C. elegans locomotion consists of three main behavioral states – forward, reverse, and pause – together with the transitions between them.
C. elegans subsists on a diet of bacteria that it finds mainly in rotting plant material (Frézal and Félix, 2015). In the laboratory, search behavior is studied in worms foraging on agar plates containing one or more dense bacterial lawns, analogous to food patches in the ethological literature. Like many other organisms, C. elegans can tune the spatial scale of random search to its physiological state, the availability of food (Wakabayashi et al., 2004; Gray et al., 2005), and prior knowledge of its distribution (Calhoun et al., 2014). The lowest values of search scale are observed during “cropping,” (Jander, 1975) the exploitation of a dense food patch. In C. elegans, two substates of cropping have been described: "dwelling," characterized by especially low crawling speed and frequent (presumably short) reversals, and "roaming," characterized by somewhat higher speeds and less frequent reversals. Transitions between dwelling and roaming, like the transitions between forward and reverse locomotion, are stochastic (Ben Arous et al., 2009; Fujiwara et al., 2002; Flavell et al., 2013). Intermediate values of search scale are observed during "local search" (Wakabayashi et al., 2004; Hills et al., 2004) when, for example, the animal is suddenly transferred from a bacterial lawn to a foodless region of the plate. The highest values of search scale are observed during “ranging,” when food is exhausted, starvation sets in, and the need to find a new food patch becomes urgent (Wakabayashi et al., 2004; Gray et al., 2005). Worms sometimes spontaneously leave a food patch well before it is exhausted, with leaving rate inversely related to food quality and food density (Shtonda and Avery, 2006; Harvey, 2009), which may reflect a tradeoff between exploitation and exploration (Bendesky et al., 2011).
At the heart of the C. elegans locomotion circuit are five pairs of premotor 'command' interneurons organized into two functional groups that promote forward and reverse locomotion, respectively (Chalfie et al., 1985; Zheng et al., 1999; Stirman et al., 2010; Schmitt et al., 2012). The two groups are reciprocally connected, and make output synapses onto distinct, nonoverlapping sets of motor neurons that control bodywall muscle. The locomotory state (forward or reverse) is believed to be determined mainly by whichever set of motor neurons is more highly activated by input from the command neurons (Kawano et al., 2011; Xie et al., 2013; Gao et al., 2015; Liu et al., 2014). Command neuron activation depends upon influences that are both external and intrinsic to the command neuron network, and appears to have a strong stochastic component that underlies switching between forward and reverse locomotion. Some command neurons are tightly linked both functionally and synaptically to upstream interneurons that also switch state stochastically in concert and counterpoint to them (Gordus et al., 2015), providing a potential additional source of the stochasticity on which random search depends. At least nine classes of chemosensory neurons and twelve classes of upstream interneurons are required for normal regulation of the duration of forward locomotion (Wakabayashi et al., 2004; Gray et al., 2005; Tsalik and Hobert, 2003; FangYen et al., 2015). Input from these neurons onto the command neuron network modulates the mean run length and, thereby, the spatial scale of random search. Search scale also appears to be modulated by neurons that release biogenic amines (serotonin, dopamine, and tyramine) (Flavell et al., 2013; Hills et al., 2004; Bendesky et al., 2011) or peptides (Ben Arous et al., 2009; Flavell et al., 2013; GloriaSoria and Azevedo, 2008; Styer et al., 2008; Reddy et al., 2009; Bhattacharya et al., 2014). These diverse signaling pathways may provide the means by which the worm optimizes its search strategy in response to feeding history (Gray et al., 2005), the quality, density and spatial distribution of food (Shtonda and Avery, 2006; Calhoun et al., 2015), and other factors that constrain survival and reproduction (GloriaSoria and Azevedo, 2008; Pujol et al., 2001; Pradel et al., 2007; Lipton et al., 2004).
Although the neural circuitry for local search has been described in considerable detail, our understanding of the system remains limited, partly for lack of key physiological data, but also for lack of a model in which to interpret the data. Common sense suggests that the forward and reverse command neurons should inhibit each other to minimize simultaneous occurrences of neuronal states for incompatible behaviors (Zheng et al., 1999). A plausible anatomical substrate for such reciprocal inhibitory connections between command neurons exists in the C. elegans connectome (White et al., 1986), but anatomical data do not specify the signs or strengths of synaptic connections. A quantitative model that incorporates physiological properties of the command neurons and their synaptic connections is needed to interpret experimental results, such as the unexpected observation that silencing some of the reverse command neurons causes a reduction in forward dwell time, and conversely for forward command neurons (Rakowski et al., 2013; Zheng et al., 1999). It is also needed to explain complex patterns of changes in dwell times observed across the three locomotory states caused by introducing or eliminating tonic membrane conductances in the command neurons, and to answer basic mechanistic questions about the control of C. elegans locomotion.
At present, the experimental data are insufficient for creating a neuronbyneuron model of the command network that incorporates biophysical details such as synaptic and membrane conductances without introducing a heavy load of unconstrained parameters (Rakowski et al., 2013). Nor would such a mechanistically detailed model necessarily provide the appropriate level of abstraction in which to intuitively understand C. elegans search behaviors, including their strong stochastic component. Instead, we have kept the level of biological detail to the minimum needed to predict the statistical distributions of dwell times in forward, reverse and pause states, and other fundamental aspects of the behavior. Each of the model's three main assumptions remains within the bounds of widely accepted experimental results; our mathematical analysis simply shows what follows necessarily from these assumptions.
To provide an empirical basis for the model we quantified C. elegans search behavior in terms of tangential velocity, defined as the speed and direction of worm's movement along its sinuous trajectory, which we recorded at higher resolution than previously possible. Behavioral data were then fit to a fourstate hidden Markov model in which each state corresponds to a unique pattern of activation across the command neurons. Importantly, rate constants governing probabilistic transitions between states in the Markov model are expressed in terms of synaptic weights in an analytically tractable version of the model. We were therefore able to validate the model by showing that it correctly predicts phenomena on which it was not fit, such as reciprocal inhibition between forward and reverse command neurons in the biological network and the behavioral effects of perturbations introduced by laser ablations and genetic mutations. Although the model is inherently probabilistic, we found that it also makes accurate predictions concerning deterministic behaviors in C. elegans, indicating a potentially high level of generality. The present findings thus establish a simple theory of C. elegans locomotory control and provide a testable model of random search that can be applied to other organisms.
Results
A neuronal model of random search in C. elegans is a theory of the relationship between activation states of the command neurons and foraging behavior. Methods presently available for observing neuronal activity in freely behaving C. elegans utilize calciumsensitive probes that have insufficient temporal resolution to observe the changes in neuronal activity associated with the rapidly changing behavioral states, especially the frequent brief pauses that are an integral part of the behavior. Therefore, as a proxy for command neuron state, we used the worm's tangential velocity, defined as the speed and direction of worm's movement along its sinuous trajectory. We focused on tangential velocity because in sinusoidal locomotion the net reactive forces produced by bodywall muscle contractions acting against the substrate are tangential to the body surface (Gray et al., 2005). Tangential velocity therefore provides the most direct readout of which group of motor neurons and command neurons (forward or reverse) is more active (Qi et al., 2013). Alternative measures of the rate of translation such as centroid velocity (PierceShimomura et al., 1999) or postural phase velocity (Stephens et al., 2011) have a less direct relationship to command neuron state because these measures either depend in complex ways on the shape of the worm, or rely on a representation of posture that ignores some of the thrustgenerating components of the worm’s shape that come into play unless the worm is moving along a fairly linear trajectory. To monitor tangential velocity as directly as possible, we painted a microscopic black spot on the worm and used a motorized stage controlled by a computer to keep the spot in the field of view (Figure 1A). The most common alternative method for measuring tangential velocity, tracking virtual points obtained by segmenting the worm’s centerline, is subject to segmentation errors introduced by low contrast images of the worm's head and tail (see Cronin et al., 2005 ) which changes the distance between virtual points. This method can also be compromised by dropped frames when the worm's centerline crosses itself during tight turns.
At the start of a 10 min observation period an individual worm was transferred from a foodladen culture plate to a bare agar surface devoid of overt sensory cues, thereby inducing a period of intensive local search behavior (Jander, 1975; Hills et al., 2004). The (x, y)coordinates of the centroid of the spot were recorded with a temporal resolution of 33 ms (i.e., frame rate = 30 Hz) and a spatial resolution of 0.5 μm that was limited mainly by the precision of the stage position encoder; the optical tracking error was much smaller (Figure 1—figure supplement 1). A spatial resolution of approximately 0.5 μm amounts to an approximately 10fold improvement over previously published tracking systems (Kocabas et al., 2012); thus worm speed (Figure 1B) could be extracted with unprecedented accuracy. For statistical analysis, worms were grouped into cohorts having the same genotype or neurons ablated (17–31 worms per cohort), which had been reared together and tested in parallel as young adults within the same 2–3 day period. This approach yielded a comprehensive data set containing a total of 8.3 million position measurements from 501 individuals in 20 cohorts.
Modelindependent identification of locomotory states
Figure 1A– D describes important features of search behavior obtained by regarding the worm as a point moving in an external reference frame (allocentric coordinates) without regard to the orientation of the body axis. The speed distribution was bimodal (Figure 1B) with a broad peak around 180 µm/s that includes both forward and reverse motion, and a narrower peak near zero that corresponds to pauses. The speed autocovariance function had multiple exponential components (Figure 1C), suggesting at least three locomotory states. The average change in heading angle ($\overline{\u2206\phi }$), plotted as a function of the intervening time interval (Figure 1D), showed that worms maintained a nearly constant heading for up to 10 s (Stephens et al., 2010; Peliti et al., 2013), but reoriented randomly within ~30 s, establishing the shortest time scale over which the behavior can be considered a Brownian random walk (Figure 1—figure supplement 2), the simplest form of random search. On shorter time scales the path takes on the character of a truncated Lévy flight (Mantegna and Stanley, 1994).
For more detailed analysis we distinguished forward from reverse movement by visual inspection of the recorded videos, and defined velocity, $v\left(t\right)$, to be a signed scalar value that denotes the speed of movement along the worm’s track in the direction of the head (+) or tail () (Figure 1E; see Materials and methods). The probability distribution of $v\left(t\right)$ (Figure 1F) showed two broad peaks that correspond to forward and reverse movement, and a narrow third peak centered at zero that corresponds to pauses. For the initial analysis we defined pauses using a fixed speed threshold of 0.05 mm/sec (Rakowski et al., 2013). Pauses occurred most frequently as transient interruptions of forward locomotion, causing the worm to stutter as it moves (Figure 1E; Video 1); stuttering also occurred, albeit less frequently, during reverse locomotion (Figure 1E; Video 2). Distinct pauses were also observed during transitions from forward to reverse (Figure 1G; Video 3) and from reverse to forward (Figure 1H; Video 4). Most pauses lasted longer than one video frame, indicating the presence of a locomotory state having a detectable dwell time; thus pauses were not merely zero crossings in plots of velocity versus time. We found that pauses during forward to reverse transitions were on average longer in duration than pauses during reverse to forward transitions (Figure 1I; p<10^{–5} ; MannWhitney Utest). These findings are consistent with the predictions of the model presented below, which uses a probabilistic criterion rather than a fixed velocity threshold to identify pauses.
The stochastic switch model
Based on the results presented in Figure 1 and previous studies noted below, we propose a minimal model for the control of random search behavior that involves two opposing neuronlike “units” that can exist in four distinct states corresponding to forward locomotion, reverse locomotion, and two pause states. This model differs from a previous model that represents the worm as a point in “shape space” (Stephens et al., 2008) in that here velocity is measured directly by observing the motion of a point on the body surface relative to the substrate, rather than indirectly by the temporal progression of shape changes. It also differs from previous models (Rakowski et al., 2013; Zheng et al., 1999; Wicks et al., 1996; Kunert et al., 2014) by representing changes in locomotory state as probabilistic transitions in a Markov process.
Ablation of individual premotor interneurons (Chalfie et al., 1985) has led to the hypothesis that the direction of locomotion is controlled by a network comprising five pairs of premotor command interneurons organized into two functional groups that promote forward and reverse locomotion, respectively. Although the anatomical pattern of synaptic connectivity among these interneurons has been established (White et al., 1986) (Figure 2A), this information does not yield an intuitive explanation of how the direction of locomotion is regulated. Nor, in our view, does the present state of the anatomical connectivity provide the basis for a neuronbyneuron simulation of the network (but see Rakowski et al., 2013), as neither signs nor physiological strengths (weights) of synapses in C. elegans can be inferred reliably from anatomical structure or neurotransmitter type, and almost nothing is known about the intrinsic membrane currents of these neurons or how they shape the inputoutput function of individual command neurons.
To establish a mathematically tractable framework for understanding how the command network functions during search behavior, we created a minimal model based on three simplifying assumptions, each of which was biologically motivated. (i) Command neurons act like binary units (Hopfield, 1982). This assumption was based on voltage recordings from command neurons in which we regularly observed two stable membrane potentials with rapid transitions between them (Figure 2B; also see Kato et al., 2015). It is also supported by the observation of a bimodal distribution of calcium activity in AVA neurons and their upstream partners AIB and RIM (Gordus et al., 2015), and the report of distinct up and down states in voltage recordings from motor neruons (Liu et al., 2014). (ii) Command neurons switch state stochastically. This assumption was based on the observation that C. elegans locomotory behavior has a strong stochastic component, with exponentially distributed dwell times in forward and reverse states (PierceShimomura et al., 1999; Zhao et al., 2003; Flavell et al., 2013; Gordus et al., 2015; Stephens et al., 2011). (iii) Command neurons within the forward pool are coactive, as are command neurons in the reverse pool. This assumption is based on simultaneous calcium imaging data from multiple command neurons in freely moving animals which suggest that the activity of neurons within the reversal pool is tightly correlated (Schrödel et al., 2013; Prevedel et al., 2014). Additionally, neurons in opposing groups are likely to be reciprocally active, as indicated by simultaneous calcium imaging from AVA and AVB (Prevedel et al., 2014; Faumont et al., 2011), as well as AVE and AVB (Kawano et al., 2011). A fourth assumption, concerning the relationship between neuronal states and behavioral states, is introduced below.
The three simplifying assumptions, together with the anatomical data (White et al., 1986), lead to a model that has two binary stochastic elements, $\mathcal{F}$ and $\mathcal{R}$, and six synaptic weights (Figure 2C). Each type of weight has a specific interpretation. The crossconnections (${w}_{\mathcal{F}\mathcal{R}}$, ${w}_{\mathcal{R}\mathcal{F}}$) represent mono and polysynaptic connections between command neurons in different groups. The selfconnections (${w}_{\mathcal{F}\mathcal{F}}$, ${w}_{\mathcal{R}\mathcal{R}}$) represent connections between command neurons in the same group, including recurrent polysynaptic pathways involving neurons outside the command network. Selfconnections also represent possible intrinsic voltage dependent currents within the command neurons, such as C. elegans plateau currents (Mellem et al., 2008). The pair of connections originating from an $\mathcal{F}$or $\mathcal{R}$ unit can have either the same sign or different signs. Allowing a single unit to have opposing effects on different postsynaptic targets is justified by the fact that synaptic weights in the model represent polysynaptic pathways, the effects of which can be excitatory or inhibitory, and by the observation that some C. elegans neurons can monosynpatically excite some postsynaptic neurons while inhibiting others (Chalasani et al., 2007). Two additional weights, ${h}_{\mathcal{F}}$ and ${h}_{\mathcal{R}}$, represent inputs from sensory neurons, interneurons, neural modulators, and any other sources outside the network (Gray et al., 2005; Fry et al., 2014), plus intrinsic membrane conductances that produce sustained effects on membrane potential (Zheng et al., 1999; Gao et al., 2015). The summed synaptic inputs onto $\mathcal{F}$ and $\mathcal{R}$ are, respectively, ${S}_{\mathcal{F}}\left(t\right)={h}_{\mathcal{F}}+{w}_{\mathcal{F}\mathcal{F}}{b}_{\mathcal{F}}\left(t\right)+{w}_{\mathcal{R}\mathcal{F}}{b}_{\mathcal{R}}\left(t\right)$ and ${S}_{\mathcal{R}}\left(t\right)={h}_{\mathcal{R}}+{w}_{\mathcal{R}\mathcal{R}}{b}_{\mathcal{R}}\left(t\right)+{w}_{\mathcal{F}\mathcal{R}}{b}_{\mathcal{F}}\left(t\right)$, where ${b}_{\mathcal{F}}\left(t\right)$ and ${b}_{\mathcal{R}}\left(t\right)$ are the states of $\mathcal{F}$ and $\mathcal{R}$ at time $t$ (1 = ON, 0 = OFF). The quantities ${h}_{\mathcal{F}}$ and ${h}_{\mathcal{R}}$ were assumed to be constant during the 10 min observation period of local search behavior on a bare agar surface.
State transitions of $\mathcal{F}$ and $\mathcal{R}$ were modeled as independent nonhomogeneous Poisson processes in which the transition rates are exponential functions of the summed synaptic input to the units, as shown in Figure 2—figure supplement 1. Changes of the state of $\mathcal{F}$ and $\mathcal{R}$ can be regarded as thermallydriven transitions across energy barriers of height proportional to ${S}_{\mathcal{F}}\left(t\right)$ and ${S}_{\mathcal{R}}\left(t\right)$, respectively. Inhibitory synaptic input increased the height of the barrier for OFF→ON transitions while decreasing the height of the barrier for ON→OFF transitions by the same amount; excitatory synaptic inputs had the opposite effect. The variable $A$ (Materials and methods, Equations 26, 27) represents the fundamental timescale of the system, defined as the rate at which units $\mathcal{F}$ and $\mathcal{R}$ change state when the summed synaptic input is zero. The present model is distinct from deterministic models of the command neuron network (Rakowski et al., 2013; Zheng et al., 1999; Wicks et al., 1996; Kunert et al., 2014) in that it is inherently stochastic, like the behavior it is meant to predict. In particular, the synaptic input to a unit does not immediately determine its state, but instead modifies the transition rates between ON and OFF states.
The two binary units of the model can exist in four states (F, R, X, Y; Figure 2D), and provide the basis for a hidden Markov model having eight transitions in which a single unit changes state. The model was further constrained by the synaptic model, which allows the eight transition rate constants to be specified by only six synaptic weights as shown in Equations 31–35 (Materials and methods). A Markov model was adopted to represent the biological system because dwell times in Markov states, like the observed dwell times in forward and reverse states (Zhao et al., 2003; Wakabayashi et al., 2004), are exponentially distributed. A hidden Markov model was required because, as noted above, states of command neurons cannot be observed directly in freely moving animals, even using optical recording methods.
The fourth assumption is a particular mapping between the states of the two command units and behavioral states of the worm. The command units, $\mathcal{F}$ and $\mathcal{R}$, are intended to represent the two pools of forward and reverse command neurons, respectively, such that the worm moves forward when $\mathcal{F}$ is ON and $\mathcal{R}$ is OFF (state F), backwards when $\mathcal{R}$ is ON and $\mathcal{F}$ is OFF (state R), and pauses when both $\mathcal{F}$ and $\mathcal{R}$ are OFF (state X). These associations between states of the model and activation states of the command neurons are well supported by previous experimental evidence, including studies showing that genetic ablation or silencing of all command interneurons induces prolonged pauses (Zheng et al., 1999; Kawano et al., 2011), but they also assume the major simplification that all command neurons in a given pool act together as a unit.
The model also permits a fourth state, in which $\mathcal{F}$ and $\mathcal{R}$ are simultaneously ON (state Y). Whether the corresponding coactivation state of forward and reverse command neurons normally exists with any significant probability remains to be shown, but it has been observed that their downstream targets, the forward and reverse motor neurons, can be active simultaneously, causing the worm to pause (Kawano et al., 2011). Given the existence of gap junction synapses between the main forward and reverse command neurons and their respective sets of forward and reverse motor neurons, it is reasonable to suppose that forward and reverse command neurons are coactive when their motor neurons are coactive. Thus, there is some evidence to designate state Y as a second pause state, which we consider to be a working hypothesis. Together, states X and Y comprise the phenomenological pause state denoted P. In what follows, we explore the logical consequences of the model’s four assumptions; it remains to be shown experimentally how closely the states of the model correspond to activity states of the command neurons.
We used a maximum likelihood method (Colquhoun and Hawkes, 1995) (Materials and methods) to estimate the set of transition rate constants that had the highest probability of generating the observed time series $v\left(t\right)$. Direct transitions between F and R, and between X and Y, were disallowed because the assumed statistical independence of the two command units implies that the probability of simultaneous transitions in $\mathcal{F}$ and $\mathcal{R}$ is vanishingly small. (Note, however, that the model does allow transitions between any two states during the interval between successive video frames by making two or more nonsimultaneous transitions; see Equation 21). We first fit the velocity distribution for each cohort with three overlapping probability distributions corresponding to forward, reverse and pause states (Figure 2—figure supplement 2), then searched for the set of transition rate constants that maximized the likelihood of the observed $v\left(t\right)$ given the velocity distributions. The resulting rate constants were used to compute the most likely sequence of states via the Viterbi algorithm (Rabiner, 1989; Viterbi, 1967). The agreement between observed velocity data and the sequence of states shown in Figure 2E was typical of the entire data set.
Wild type locomotion
The maximum likelihood rate constants for 5 wildtype cohorts, together with the predicted state probabilities and mean dwell times computed from them, are given in column A of Table 1. The model’s predicted mean dwell time in the reverse state (${d}_{\text{R}}=$ 1.94 ± 0.04 s) agreed with previously reported values (Zheng et al., 1999; Kawano et al., 2011). In contrast, the predicted mean dwell time in the forward state (${d}_{\text{F}}=$ 5.33 ± 0.25 s) was smaller than previously reported when dwell time was measured by eye (13–35 sec) (Zhao et al., 2003; Zheng et al., 1999; Brockie et al., 2001; Ryu and Samuel, 2002) or by velocity threshold crossings (9–16 sec) (Rakowski et al., 2013; Stephens et al., 2011). To determine whether this difference arose because we used a hidden Markov model rather than a fixed velocity threshold, we also identified states based on a fixed velocity threshold of 0.05 mm/s, and calculated the resulting mean dwell times: ${d}_{\text{F},0.05}=$ 1.86 ± 0.03 s; ${d}_{\text{R},0.05}=$ 1.23 ± 0.02 s; ${d}_{P,0.05}=$ 0.14 ± 0.001 s. We attribute the short mean dwell times in state F that we observed using either the hidden Markov model or a fixed velocity threshold to the fact that our tracking system is capable of revealing briefer visits to state P, which interrupt forward runs, than previous methods. Ignoring transient interruptions of forward locomotion (i.e., FPF transitions) and using the fixed velocity threshold of 0.05 mm/s yielded longer a mean forward dwell time of 9.13 ± 0.15 s, which matches the value obtained by others using the same threshold (8.98 ± 0.57 s) (Rakowski et al., 2013). Predicted mean dwell times in the two pause states differed substantially from each other (${d}_{\text{X}}=$0.44 ± 0.03 s, ${d}_{\text{Y}}=$ 0.21 ± 0.02). We assigned the long and short pause states to X and Y, respectively, based on the idea that the energetically expensive state in which both units are on should be relatively shortlived.
In previous work, transitions between locomotory states in C. elegans have been analyzed by choosing a speed threshold to distinguish pause states from the movement states (Rakowski et al., 2013; Salvador et al., 2014; Stephens et al., 2011). The choice of threshold is important because it affects the measured dwell times, yet is necessarily arbitrary because the velocity distributions of the states overlap (Figure 1F). The hidden Markov model used here replaces arbitrary thresholds with empirically determined state transition rates (i.e., the set of rates that maximizes the probability of the observed velocity time series), from which one can determine the sequence of states that is most likely to have generated the data (the Viterbi algorithm). The hidden Markov model thus offers two advantages: (1) it provides a statistical criterion for selecting the best parameter values and (2) it takes into account the uncertainties in identifying the state of the system from velocity data.
Under the assumptions of the hidden Markov model the state of the system cannot be observed directly because the velocity distributions overlap, making it impossible to test directly whether the predicted state probabilities agree with the observed velocity data. Nevertheless, an important test of the model can be obtained using the Viterbi algorithm to identify the most likely sequence of states given the observed velocity data, from which the histogram of dwell times in each state can be computed and compared to the exponential distribution predicted by the Markov model (Figure 2—figure supplement 3). The degree of agreement between the distributions shows that our model provides a good description of the system.
The initial rationale for including two pause states in the hidden Markov model came from our modelindependent analysis of the tracking data (Figure 1I), which showed different dwell time distributions for pauses at FPR and RPF transitions. To test whether having two pause states yielded a statistically significant improvement in the ability of the model to fit the data, we eliminated one of the pause states and asked whether the resulting reduction in likelihood was greater than could be attributed to the reduction in the number of free parameters (see Table 1). For this comparison we constrained the transition rates into state Y to be extremely small (${a}_{FY}={a}_{RY}=$ 10^{10} s^{1}), effectively eliminating state Y and reducing the number of free parameters from six to four. The reduction in likelihood caused by eliminating one of the pause states was highly significant, and cannot be attributed simply to the elimination of two parameters (p<10^{100}; likelihood ratio test). Separately, we considered the most general onepause state model, which allows direct transitions between states F and R and has no constraints on the 6 transition rates other than that they are all ≥0. The fit of this model (Table 1 column C) converged to nearly the same set of transition rates as the onestate model described above (Model B). These comparisons show that our model with two pause states and six free parameters (the six synaptic weights) provides a much better fit to the data than models with only one pause state. We conclude that the tracking data contain a statistically significant signature of two distinct pause states. The model explains the observation that the pause dwell times during FPR transitions are longer than during RPF transitions (Figure 1I) in terms of the different dwell times in states X and Y (${d}_{\text{X}}>{d}_{\text{Y}}$), and the strong tendency to cycle clockwise through state space, exiting from state F to state X and from state R to state Y as shown by the fate diagram (Figure 3).
It has been reported that pauses in C. elegans locomotion occur at specific points in “shape space” (Stephens et al., 2011), suggesting the worm pauses in preferred postures. To investigate this possibility, we analyzed worm tracks before and after pauses, inferring posture from the path of the tracking spot. This inference is justified by the fact that on an agar surface the worm moves without slipping, such that each segment of the body traces the trajectory of the one before it. Thus, the path of the tracking spot leading up to the pause reveals the worm’s posture posterior to the spot during forward locomotion, and anterior to the spot during reverse locomotion (Figure 4).
Plotting mean curvature versus distance along the track (Figure 4A) reveals only a weak tendency to stop in a particular posture in state X (r = 0.14; Figure 4B). Nearly all of the transitions into state X were either stutters during forward locomotion (FXF transitions) or reversals (FXR transitions); when these were analyzed separately, similarly weak postural preferences were found at FXF transitions (r = 0.14) and FXR transitions (r = 0.14). A nearly identical result (r = 0.14 ) was obtained using a fixed velocity threshold of 0.05 mm/s rather than the hidden Markov model to determine state. For the latter case, in which there is only one pause state, we analyzed the posture at all FP transitions, which almost always correspond to FX transitions in the hidden Markov model because FY transitions are extremely rare (see Figure 3 ).To test whether the failure to find a strong postural preference at FX transitions was due to including very short pauses in the analysis, we repeated the analysis after reclassifying all pauses shorter than a minimum duration as a continuation of the previous state, and obtained the same result; we found no strong postural preference at FX transitions for minimum pause durations up to 2 s (r = 0.16, 0.19, 0.23, 0.3 for X dwell times > 0.33 s, 0.67 s, 1 s, and 2 s, respectively); longer dwells in state X were too rare to analyze. Thus, FX transitions can occur at any locomotory phase and do not occur preferentially at a particular posture (Figure 4D); in the case of FXR transitions the worm generally retreats along the same track. In contrast, we found a strong tendency to stop in a particular posture in state Y (Figure 4A,C,E; r = 0.71). Almost all entries into state Y were RYF transitions and these were associated with a ventral bend in the middle of the worm (Figure 4E). These results suggest fundamental differences between the control of forward and reverse locomotion. In our model, forward locomotion terminates when forward command neurons turn OFF, and this can happen at any phase, whereas reverse locomotion terminates when forward neurons turn ON, and this is most likely to happen at a particular phase. The latter could be explained by phasic feedback from the locomotory pattern generator to the forward neurons (Li et al., 2006).
Ablation of command neurons
To determine the contributions of individual command neurons to the overall function of the command network, we separately ablated the pair of neurons that comprises each command neuron class, then tracked ablated and sham operated animals during local search. Mean velocities in F and R, if significantly changed, were reduced (Pokala et al., 2014) (Figure 5A; $\star \star $), as was the frequency of undulations during forward and reverse locomotion (Table 3). In many organisms, the frequency of rhythmic behaviors is regulated by the amplitude of tonic excitatory drive to the associated pattern generator (Weeks, 1978; Satterlie and Norekian, 2001; Böhm and Schildberger, 1992; Deliagina et al., 2000; Dembrow et al., 2003; Hedwig, 2000; Sirota et al., 2000). To explain our results we propose that ablation of the locomotory command neurons reduces tonic drive to the presumptive locomotory pattern generator (Xie et al., 2013; Gao et al., 2015).
A previous study found that ablating a subset of the reverse command neurons (AVAL and AVAR) reduces dwell time in the reverse state but also paradoxically reduces dwell time in the forward state (Zheng et al., 1999). Similarly paradoxical effects have been reported following ablation of the reciprocally connected brain stem nuclei that regulate sleep and wakefulness (Saper et al., 2010). The stochastic switch model predicts and explains such effects. In principle, the ablation of a subset of neurons in a pool of coactive neurons can have widespread effects on the pool’s overall input and output connectivity. Widespread effects can be expected because ablation removes not only the outgoing synaptic connections from the ablated neurons, but also the targets of incoming synaptic connections. In the C. elegans command neuron network, ablating a reverse command neuron such as AVA potentially reduces four of the six weights in the network: ${h}_{\mathcal{R}}$, ${w}_{\mathcal{R}\mathcal{R}}$, ${w}_{\mathcal{R}\mathcal{F}}$, and ${w}_{\mathcal{F}\mathcal{R}}$. Thus, a single ablation can move the system a considerable distance in weight space toward the uncoupled state in which all weights are zero. In the limiting case of a fully uncoupled network, all dwell times approach a value of $1/2A$, where $A$ is the intrinsic switching time of the stochastic units (see Materials and methods, equations 31–34); henceforth we will use ${d}_{0}$ to denote the uncoupled dwell time. Dwell times that in intact animals are greater than ${d}_{0}$ will be reduced by ablation, whereas dwell times that are less than ${d}_{0}$ will be increased. In particular, if ${d}_{\text{F}}$ and ${d}_{\text{R}}$ are both greater than ${d}_{0}$, ablation of a reverse command neuron is expected to reduce both dwell times; the same is true for ablation of a forward command neuron. Thus the observed paradoxical effects of ablations are to be expected if ${d}_{0}$ is below ${d}_{\text{F}}$ and ${d}_{\text{R}}$.
To determine the actual relationship between ${d}_{0}$ and dwell times in the forward and reverse state, we estimated the rate constants in ablated animals and sham operated controls, and computed the corresponding dwell times (Figure 5B; Table 4). Dwell times in F and R, if significantly altered by the ablation ($\star \star $), were reduced, indicating that ${d}_{0}$ is indeed below ${d}_{\text{F}}$ and ${d}_{\text{R}}$. Additionally, dwell times in the pause states ${d}_{\text{X}}$ and ${d}_{\text{Y}}$ were increased, with one exception (${d}_{\text{Y}}$, AVB). Thus, the observed pattern of dwell time changes is consistent, overall, with a value of ${d}_{0}$ that is between the dwell times of the movement states and the dwell times of the pause states. This finding allowed us to place bounds on ${d}_{0}$. Specifically, ${d}_{0}$ must be less than or equal to the lowest postablation value of ${d}_{\text{R}}$, and greater than or equal to the largest postablation value of ${d}_{\text{X}}$; thus, 0.58 ≤ ${d}_{0}$ ≤ 1.24 sec. Furthermore, because $A=1/2{d}_{0}$, we can infer that 0.40 Hz ≤ $A$ ≤ 0.86 Hz. This inequality provides an estimate of the fundamental time scale of stochastic switching in C. elegans locomotion. For subsequent analysis, we defined ${A}_{min}=$ 0.40 Hz and ${A}_{max}=$ 0.86 Hz.
Synaptic weights in the stochastic switch model
Having placed bounds on $A$, we were able to compute synaptic weights in the model (Table 2). This was done by deriving expressions for the weights in terms of the rate constants (Materials and methods, Equations 36–38) and substituting into these equations our estimates of rate constants together with the values of ${A}_{min}$ and ${A}_{max}$. We found that input weights, ${h}_{\mathcal{F}}$ and ${h}_{\mathcal{R}}$, are small and positive, suggesting that these inputs may provide modest but steady excitation to the system (Figure 6A). The selfconnections ${w}_{\mathcal{F}\mathcal{F}}$ and ${w}_{\mathcal{R}\mathcal{R}}$ are also mainly positive, indicating that the ON states may be stabilized by intrinsic or extrinsic positive feedback. The crossconnections ${w}_{\mathcal{F}\mathcal{R}}$ and ${w}_{\mathcal{R}\mathcal{F}}$ are negative, indicating reciprocal inhibition, as expected for neurons that activate opposing behavioral states. Furthermore, the magnitude of ${w}_{\mathcal{F}\mathcal{R}}$ is greater than the magnitude of ${w}_{\mathcal{R}\mathcal{F}}$, suggesting that the animal spends more time in the forward state than the reverse state in part because the forward neurons inhibit the reverse neurons more strongly than the reverse neurons inhibit the forward neurons.
Synaptic weights in an abstract network model such as this one, where neuronal state is activation rather than voltage, are not generally interpretable as synaptic conductances. Rather, they represent the functional effects of one neuron on another, such as the degree of excitation or inhibition produced by a unit change in activation. Thus, synaptic weights in the Stochastic Switch model cannot be said to predict the magnitude of synaptic conductances, but they can be said to predict aspects of functional connectivity in certain cases. For example, as command neurons AVA and AVB are behaviorally much more important than the others (Chalfie et al., 1985) (see also Figure 5A,B), it is reasonable to assume that the signs of their functional synaptic connections match the signs of the net functional connections in the biological network. Thus, the model predicts reciprocal inhibition (Qi et al., 2012) between AVA and AVB under this assumption. We tested this prediction by photoactivating either AVA or AVB with channelRhodopsin2 and recording electrophysiologically from AVB or AVA, respectively (Figure 6B,C). We found that the reversal potential of optically induced synaptic currents in AVA and AVB was more negative than the zerocurrent potential in these neurons (Figure 6B,C,D), indicating synaptic inhibition as predicted by the model. This inhibition is likely to be monosynaptic as C. elegans command neurons are cholinergic and express inhibitory postsynaptic receptors that respond to acetylcholine (Pereira et al., 2015). Additionally, the connection from AVB to AVA appeared to be stronger than the connection from AVA to AVB (Figure 6E), measured in terms of the amplitude of the synaptic current at a holding potential approximately equal to the membrane potential when command neurons are in their depolarized state (Figure 2B). However, we do not exclude the possibility that AVB was more strongly activated than AVA as a result of differential expression of the photoprobe. These findings demonstrate the feasibility of using the worm’s velocity, $v\left(t\right)$, a simple behavioral measure, to predict functional synaptic connections between populations of neurons in a biological neural network, at least under certain assumptions concerning the relationship between model network weights and physiological synaptic strengths.
Genetic effects on command neuron function
Two classes of ion channel mutants that affect membrane conductances in the command neurons are also known to alter locomotory behavior in systematic ways, thus providing key insights into command neuron function (Zheng et al., 1999). The hyperpolarizing class (“HYP”) comprises three genotypes in which release of the excitatory neurotransmitter glutamate, presumed to be tonic, is disrupted by mutations that affect either presynaptic (eat4(ad572), eat4(ky5)) or postsynaptic mechanisms (glr1(n2461)). These mutations are hypothesized to cause chronic hyperpolarization of the command neurons by reducing depolarizing currents. The depolarizing class (“DEP”) comprises two genotypes in which a constitutively activated glutamate receptor is expressed in the command neurons (glr1::glr1(A/T), nmr1::glr1(A/T)). These mutants are hypothesized to chronically depolarize the command neurons.
We found that the frequency of locomotory undulations was decreased in HYP mutants and increased in DEP mutants compared to wildtype controls (Table 5), consistent with the likely effects of respectively increasing and decreasing tonic drive to the presumptive pattern generator for locomotion. Importantly, however, it is possible that both classes of mutation also alter the input resistance of the command neurons. The closure or removal of glutamate receptors in HYP mutants should increase input resistance whereas the introduction of constitutively active glutamate receptors in DEP mutants should decrease it. Thus, the previously observed effects of these mutations on locomotory state transitions (Zheng et al., 1999) could be the result of changes in membrane potential (ΔV), input resistance (Δr), or both.
Changes in membrane potential and input resistance can both be represented in the stochastic switch model by changes in synaptic weights. We modeled the effects of ΔV by adding an increment $\u2206h$ ($1\le \u2206h\le $ 1) to wild type $h$ values, with negative $\u2206h$ for HYP mutations and positive $\u2206h$ for DEP mutations. We modeled the effect of Δr as a change in the magnitude of synaptic weights ($h$ and $w$ quantities). This representation of Δr is appropriate because changes in input resistance alter the magnitude of the voltage change that would be produced by a fixed postsynaptic current. All weights were scaled by a common factor $z$ (1 < $z$ < 2 for HYP mutants; 0 < $z$ < 1 for DEP mutants).
Here we consider the effects of ΔV and Δr on dwell times in the stochastic switch model to enable direct comparison with the original study of HYP and DEP strains (Zheng et al., 1999). Dwell times can be written as functions of weights:
These equations show that the ΔV and Δr hypotheses make qualitatively distinct predictions. The simplest case is dwell ${d}_{\text{X}}$, which depends only on ${h}_{\mathcal{F}}$and ${h}_{\mathcal{R}}$. Equation 1 shows that ${d}_{\text{X}}$ rises and falls as $h$ terms are made more negative or positive, respectively. Thus, under the ΔV hypothesis, ${d}_{\text{X}}$ should rise in HYP mutants and fall in DEP mutants (Figure 7A, row 4). In contrast, under the Δr hypothesis, in which weight magnitudes ($\leftw\right$ and $\lefth\right$) decrease in DEP mutants and increase in HYP mutants, ${d}_{\text{X}}$ should rise in DEP mutants and fall in HYP. To distinguish between these hypotheses, we measured dwell times in mutants and wild type animals during local search. The pattern of observed changes in ${d}_{\text{X}}$ matched the pattern predicted by the ΔV hypothesis but not the Δr hypothesis (Figure 7C, row 4). Thus, the effects of membrane potential appear to dominate the effects of changes in synaptic strength in the case of mutant values of ${d}_{\text{X}}$.
In contrast to ${d}_{\text{X}}$, ${d}_{\text{F}}$ and ${d}_{\text{R}}$ depend on $w$ terms as well as $h$ terms. Under the ΔV hypothesis, the $h$ terms but not the $w$ terms would be affected by the mutations. Positive and negative increments in $h$ have the effects shown in Figure 7A, rows 1 and 2; ${d}_{\text{F}}$ and ${d}_{\text{R}}$ are predicted to shift in opposite directions. Changes in ${d}_{\text{F}}$ are dominated by the effects of ${h}_{\mathcal{F}}$ on the first term in Equation 2 (which represents ${a}_{\text{FX}}$) because the second term in the equation (which represents ${a}_{\text{FY}}$) remains close to zero in the mutants. Analogously, changes in ${d}_{\text{R}}$ are dominated by the effects of ${h}_{\mathcal{F}}$ on the second term in Equation 3 (${a}_{\text{RY}}$) because the first term in the equation (${a}_{\text{RX}}$) remains close to zero in the mutants.
The Δr hypothesis makes a distinctly different prediction. In this version of the model, $w$ terms and $h$ terms would both be affected by the mutations. Now, the predicted pattern of dwell time changes across both ${d}_{\text{F}}$ and ${d}_{\text{R}}$ is such the both dwell times shift in the same direction (Figure 7B, rows 1 and 2); specifically, dwell times in DEP and HYP mutants move toward or away from their uncoupled dwell times, respectively. Taken together, the pattern of observed changes in ${d}_{\text{F}}$ and ${d}_{\text{R}}$ matched the pattern predicted by the Δr hypothesis (Figure 7C, rows 1 and 2) but not the ΔV hypothesis. We conclude that changes in synaptic strength may dominate the effects of changes in membrane potential on mutant values of ${d}_{\text{F}}$ and ${d}_{\text{R}}$.
Neither hypothesis predicts the observed changes in ${d}_{\text{Y}}$ (Figure 7C, row 5) which resembled the pattern of changes in ${d}_{\text{X}}$ (Figure 7C, row 4). However, the ΔV hypothesis correctly predicts observed dwell times in the overall pause state ${d}_{\text{P}}$ (Figure 7C row 3). This is because ${d}_{\text{P}}$ is dominated by ${d}_{\text{X}}$ and changes in ${d}_{\text{X}}$ are wellpredicted by the ΔV model as noted above. Overall, our analysis of the effects of HYP and DEP mutations in terms of the Stochastic Switch Model points to a role for changes in both membrane potential and input resistance in regulating dwell times.
Regulation of search scale
The Stochastic Switch Model immediately suggests a family of models for the regulation of the spatial scale of random search in response to the availability of food and the worm’s physiological state. The scale of random search is determined primarily by ${m}_{\text{F}}$, the mean distance traveled during a forward run. In C. elegans, a run begins with a transition from state R (via P) into state F and continues until the next transition into state R. Any run may include one or more visits to state P, but FPF transitions are not usually associated with changes in heading. In terms of the Stochastic Switch Model, ${m}_{\text{F}}=\overline{{v}_{\text{F}}}{p}_{\text{F}}/{f}_{\text{RPF}}$, where $\overline{{v}_{\text{F}}}$ is the average velocity in state F, ${p}_{\text{F}}$ is the probability of being in state F, and ${f}_{\text{RPF}}$ is the frequency of RPF transitions (Materials and methods, Equation 39), which coincide with random reorientations. Importantly, under the approximation ${a}_{\text{FY}}\cong 0$ (Table 1, column A), ${m}_{\text{F}}$ is can be expressed as a function of just three of the six weights in the network:
We refer to these weights as potential control points in the network. In a minimal model of search scale regulation, ${m}_{\text{F}}$ could be controlled by sensory inputs represented by ${h}_{\mathcal{F}}$ and ${h}_{\mathcal{R}}$ (Figure 8A).
Search scale (${m}_{\text{F}}$) together with the frequency of reversals (FPR transitions), have been used to define the three search modes commonly recognized in C. elegans: cropping, local search, and ranging. To find minimal models for regulation of search mode, we performed exhaustive searches of subregions of network’s sixdimensional weight space. Subspaces, defined by one, two, or three weights, were scanned across a wide range of values (6 ≤ $w$ ≤ 6) while other weights remained fixed at their wild type levels (Figure 7B–H). The performance of each configuration of the network was scored according to whether it matched the range of ${m}_{\text{F}}$ magnitudes and reversal frequencies characteristic of each mode (see Materials and methods). Another consideration was the number of distinct search modes available; accordingly, we also noted the density with which the plane defined by reversal frequency and ${m}_{\text{F}}$ was covered in the scan (Figure. 8BH, gray symbols).
All three search modes were available in the subspace defined by the control points $({h}_{\mathcal{F}},{h}_{\mathcal{R}},{w}_{\mathcal{F}\mathcal{F}})$ (Figure 8B, Figure 8—figure supplement 1). However, only cropping and local search were available in the complementary subspace $({w}_{\mathcal{R}\mathcal{R}},{w}_{\mathcal{F}\mathcal{R}},{w}_{\mathcal{R}\mathcal{F}})$ (Figure 8C); thus, to achieve the full set of search modes, at least one of the weights in equation 5 must be free to change. None of the controlpoint weights was sufficient on its own to produce all three search modes (Figure 8D–F). Scanning the subspaces $({h}_{\mathcal{F}}$, ${w}_{\mathcal{R}\mathcal{R}})$ and $({h}_{\mathcal{R}}$, ${w}_{\mathcal{F}\mathcal{R}})$ showed these pairs of weights to be sufficient for all modes (Figure 8G,H), but a threedimensional subspace containing at least one of the controlpoint weights was a necessary condition for both dense coverage of this plane and the presence of all three search modes (Table 7). We suggest that these threeweight subspaces constitute the most likely minimal models for the regulation of search in C. elegans. They could be tested by chronic manipulation of controlpoint weights utilizing a variety of approaches, such as chemical or optical probes that alter tonic inputs to the command network from sensory neurons and interneurons represented by the parameters ${h}_{\mathcal{F}}$ and ${h}_{\mathcal{R}}$.
Biased random walks
Mean forward run length is also modulated during biased random walks, increasing or decreasing when the animal is moving in a favorable or unfavorable direction, respectively (PierceShimomura et al., 1999; Block et al., 1982; Iino and Yoshida, 2009; Luo et al., 2014). When C. elegans is engaged in chemotaxis toward an attractive substance, the direction of motion relative to the gradient is represented by specialized chemosensory neurons that respond either to increases (ON cells) or decreases in concentration (OFF cells) (Thiele et al., 2009) ; moreover, interventions that activate ON cells or OFF cells promote runs and pirouettes, respectively (Suzuki et al., 2008). Thus, in one simple model of randomwalk chemotaxis, ON cells increase ${h}_{\mathcal{F}}$ and decrease ${h}_{\mathcal{R}}$, whereas OFF cells do the opposite. Simulations show that this model is sufficient to generate realistic chemotaxis in a point model of search behavior in C. elegans (Figure 8—figure supplement 2) when the worm is below the target concentration of attractant. Similar circuitry can explain biased random walks in response to other physical gradients (Lockery, 2011).
The Stochastic Switch Model and deterministic behaviors
In addition to random search, the command neurons in C. elegans are required for a variety of escape responses (Hart et al., 1995) that are deterministic in that ${p}_{\text{R}}$ closely approaches unity for strong stimuli (Wittenburg and Baumeister, 1999; Tobin et al., 2002; Liu et al., 2012; Mohammadi et al., 2013). C. elegans escape responses can be produced by two pathways, one that requires the reverse command neurons (Chalfie et al., 1985) and one that does not (Piggott et al., 2011). Three distinct circuit motifs for the functional connectivity underlying escape responses requiring reverse command neurons are conceivable (Figure 9A). In the Push motif, nociceptive neurons excite reverse command neurons via ${h}_{\mathcal{R}}$ thereby increasing the rate constants for transitions in which $\mathcal{R}$ turns ON (${a}_{\text{XR}}$ and ${a}_{\text{FY}}$), and decreasing the rate constants for transitions in which $\mathcal{R}$ turns OFF (${a}_{\text{RX}}$ and ${a}_{\text{YF}}$). In the limit where ${h}_{\mathcal{R}}$→ ∞, both ${a}_{\text{XR}}$ and ${a}_{\text{FY}}$→ ∞, whereas ${a}_{\text{RX}}$ and ${a}_{\text{YF}}$→ 0 (Figure 9B). The system now inhabits only states R and Y, and ${p}_{\text{R}}={a}_{\text{YR}}/\left({a}_{\text{YR}}+{a}_{\text{RY}}\right)$. In the Pull motif, nociceptive neurons inhibit the forward command neurons via ${h}_{\mathcal{F}}$. In the limit where ${h}_{\mathcal{F}}$→ ∞, the system switches only between states R and X and ${p}_{\text{R}}={a}_{\text{XR}}/\left({a}_{\text{XR}}+{a}_{\text{RX}}\right)$. In the third motif, in which Push and Pull are combined, R becomes an absorbing state (${p}_{\text{R}}=1)$. Using the rate constants shown in column A of Table 1 to compute limiting values of ${p}_{\text{R}}$ in each motif, we found that the Pull and PushPull motifs are sufficient for deterministic escape, whereas the Push motif is not (Figure 9B). Thus, inhibition of forward command neurons is required for deterministic escape, predicting that nociceptive neurons functionally inhibit these neurons.
To test this prediction we examined the ASH neurons, a pair of nociceptive sensory neurons required for the majority of escape responses in C. elegans. ASH neurons have anatomically defined monosynaptic and polysynaptic connections to both the behaviorally dominant command neurons AVB and AVA (Chalfie et al., 1985; White et al., 1986). We have previously shown that the functional connection from ASH to AVA is excitatory (Lindsay et al., 2011). To test whether the functional connection from ASH to AVB is inhibitory, we photoactivated ASH neurons while recording from AVB (Figure 9C,D). The reversal potential of this connection was more negative than the zero current potential, indicating inhibition as predicted by the model. Thus, ASHmediated escape may be controlled by a PushPull motif, further demonstrating the feasibility of using behavioral data to predict populationlevel synaptic connectivity. The source of the AVB inhibition could be the inhibitory connection from AVA, polysynaptic pathways from ASH to AVB, or both.
Notably, the Pull and PushPull motifs are equally effective in driving ${p}_{\text{R}}$ to 1.0 (Figure 9B). Nevertheless, computation of the expected latency to the first reversal event when a forward moving animal suddenly encounters a strong nociceptive stimulus indicates a 2.3fold reduction in latency for the PushPull motif (Figure 9B, parenthetical values). We conclude that the ASH mediated escape circuit in C. elegans may be specialized for short latency escape responses.
Discussion
The Stochastic Switch Model is cast at a level of biological detail that is minimally sufficient to capture the stochastic dynamics of C. elegans locomotion in neuronal terms. Despite its simplicity, the model predicts the unexpected effects of neuronal ablations and genetic manipulations. It also predicts the sign and strengths of key synaptic connections, which were confirmed by combining optogenetics with electrophysiology. The model is immediately extensible to random search at a variety of spatial scales, biased random walks such as chemotaxis, and deterministic escape behaviors. The predictive success of the model indicates that random search in C. elegans can be understood in terms of a neuronal flipflop circuit involving reciprocal inhibition between two populations of stochastic neurons. Two likely sources of stochastic state transitions are quantal synaptic transmission and ion channel gating. Both of these sources derive their randomness from thermal fluctuations at the molecular level, a phenomenon that is common to all nervous systems. The stochasticity underlying search behavior in C. elegans could be intrinsic to the command neurons, their presynaptic neurons (Gordus et al., 2015), or both.
The simplifying assumptions of the model introduce several limitations worth noting. (i) By representing the ten command neurons as only two functional units, the model ignores possible functional differences between individual neurons within each group. (ii) By design, the model predicts exponentially distributed dwell times, but Figure 2—figure supplement 3 shows that this relationship is only approximate. (iii) The model also has no provision to explain the strong correlation between locomotory phase and entry into state Y (Figure 4), although this could be added by modeling feedback from the pattern generator as a timevarying component of ${h}_{\mathcal{F}}$ and ${h}_{\mathcal{R}}$. (iv) The model does not take into account temporal correlations in velocity, but instead uses only the present velocity, along with the present state, to compute transition probabilities. For example, the fact that locomotion gradually slows before the worm enters the pause state (Figure 1G,H) suggests that transition probabilities might be more reliably calculated by including the recent velocity history, rather than just the present velocity. (v) Finally, the model does not attempt to explain the observation that the number of command neurons that are present and the degree of command neuron activation has an effect on velocity and undulation frequency (Figure 5A, Table 3, Table 5). Velocity modulation could be incorporated by relaxing the assumptions that command neurons within pools are coactive and have a single nonzero level of activation.
Although the model correctly predicts several unexpected and even paradoxical observations at the behavioral and electrophysiological levels, it would be premature to conclude that the biological system functions as assumed. This caution extends to all of the model's assumptions, including the mapping relationship between pause states X and Y and their behavioral correlates. We view the pause states as theoretical constructs having an epistemological status akin to theoretical constructs in many widely accepted models, such as the gating particles that were proposed in the HodgkinHuxley model of the squid action potential to explain the voltage sensitivity of ion channels.
An altogether different method for analyzing locomotory states in C. elegans also proposed the existence of two pause states (Stephens et al., 2008). In that work, each pause state was associated with a particular locomotory phase. In contrast, we found that only state Y occurred in association with a particular posture (a ventral bend in the middle of the body), whereas state X occurred with essentially no postural preference. The reason for this discrepancy may be that pauses are identified in different ways in the two studies. Here pauses are identified in terms of tangential velocity. In Stephens et al. (Stephens et al., 2008; 2011; 2010), however, pauses are identified in the phase space defined by the amplitudes of first two principle components of the worm's instantaneous shape. For the two approaches to yield the same result, minima in the magnitudes of tangential and phase velocity would have to be coincident at all times. We believe this outcome is unlikely because the third and fourth principle shape components, which account for approximately 30% of the shape variance (Stephens et al., 2008), meet the necessary and sufficient conditions for generating tangential thrust: a gradient of curvature along the worm’s centerline (Gray, 1946; Gray, 1953; Gray and Lissmann, 1964); this is one way thrust is believed to be generated during omega turns (Stephens et al., 2008). Thus, the worm can be moving with respect to the substrate even when phase velocity is zero. Overall, we speculate that pauses in phase velocity are a subset of pauses in tangential velocity. The extent to which this is true could be determined by performing spot tracking and shape analysis on the same individual worms.
It will be interesting to test several additional predictions of Stochastic Switch Model:
The sign of the input weights, ${h}_{\mathcal{F}}$ and ${h}_{\mathcal{R}}$ predicts tonic excitation of the network. This could be the result of constitutive excitatory synaptic inputs, or depolarizing leakage currents in individual command neurons as has been proposed (Gao et al., 2015).
The sign of the selfconnections ${w}_{\mathcal{F}\mathcal{F}}$ and ${w}_{\mathcal{R}\mathcal{R}}$ predicts one or more mechanisms of selfexcitation within command neuron pools. These might include excitatory connections between command neurons, or intrinsic membrane currents capable of producing plateau potentials (Mellem et al., 2008).
The fate diagram (Figure 3) predicts that forward commands neurons generally lead the changes in direction during spontaneous locomotion. For example, transitions from F to R almost always begin with the $\mathcal{F}$ unit turning off, whereas transitions from R to F almost always begin with the $\mathcal{F}$ unit turning on. This prediction could be tested by calcium imaging in command neurons in freely moving animals (Nguyen et al., 2015); Venkatachalam et al., 2015).
Finally, the prediction that forward command neurons lead the changes in direction, coupled with the observation that transitions from R to Y occur at a particular phase, predicts that the forward command neurons are the predominant site at which phasic feedback from the locomotion pattern generator influences the network. Direct observation of neuronal activity in freely moving animals would be the ideal experiment to confirm the existence of the two pause states proposed in the Stochastic Switch Model (Nguyen et al., 2015; Venkatachalam et al., 2015). In particular, it will be necessary to show that whenever all command neurons are OFF, or all are ON, tangential velocity goes to zero. These experiments will be challenging because they must be done by imaging neuronal activity in freely moving animals at a temporal resolution that exceeds what can be obtained with the current generation of calcium probes. In fact, it may be necessary to use voltage probes rather than calcium indicators because even a very fast calcium probe will be limited by the dynamics of calcium accumulation, which is slow on the time scale of the pause dwell times predicted by the model. Another potential complication is that velocity may not change instantaneously with changes in the state of the command network, but with a delay imposed by time constants in the motor system. A less direct approach, although one with much higher temporal resolution, would be to make whole cell current clamp recordings from command neurons or motor neurons in restrained animals, which cycle through global brain states analogous to forward and reverse locomotion (Kato et al., 2015) even though they cannot move. Instances in which both motor systems are OFF or ON would provide evidence for states X and Y, respectively.
Like the Stochastic Switch Model, a previous model of the command neuron circuit by Rakowski et al. (Rakowski et al., 2013) predicts reciprocal inhibition between command neurons. Although the two models analyze locomotion behavior in terms of the same three behavioral states – forward, reverse, and pause – the models have essentially no points of mathematical contact. In the Rakowski model, neurons are deterministic electrical compartments and only the longterm average state probabilities of the network are computed. In the Stochastic Switch Model, by contrast, neurons are inherently stochastic and instantaneous state is computed. These disparities are significant because only the Stochastic Switch Model can predict temporal phenomena including such fundamental quantities as transition rates and mean dwell times. The fact that the both models predict reciprocal inhibition may reflect that fact that the behavioral signal of reciprocal inhibition is strong enough to transcend large differences between models.
Mammalian sleep, like C. elegans locomotion, is composed of numerous abrupt alternations between opposing behavioral states. Sleep is punctuated frequently by brief periods of wakefulness, and dwell time distributions in sleep and wake states indicate that switching between them is a stochastic process (Lo et al., 2004). Sleep and wakefulness are controlled by mutually inhibitory brainstem nuclei, implying a reciprocal inhibition motif. In a significant parallel to the effects of command neuron ablations on dwell times in C. elegans locomotion (Figure 5B), lesions of sleeprelated nuclei simultaneously reduce the dwell times in both sleep and wake states, as do lesions of wakefulness nuclei (Saper et al., 2010). Thus the relationship between synaptic uncoupling of the circuit and changes in dwell times may be a general principle of reciprocal inhibition in stochastic neuronal networks. Further study of invertebrate models of this circuit motif would be a productive means of identifying the genetic and physiological underpinnings of such circuits.
The debut of the essentially complete wiring diagram of the C. elegans nervous system raised the prospect of the first account of the entire behavioral repertoire of an organism at singleneuron resolution (White et al., 1986; Varshney et al., 2011). To date, the repertoire of behaviors commonly recognized in C. elegans can be divided into three main functional categories, subsuming 23 different elementary actions (Faumont et al., 2012). Because the command neurons considered here are required for almost half of this repertoire, the Stochastic Switch Model is a significant step toward a comprehensive understanding of the neuronal basis of behavior in this animal, bringing us closer to the goal of computing the behavior of an entire organism. Though abstract by design in its representation of individual neurons and synapses, the model accommodates not only random search at multiple spatial scales (Figure 8), but also biased random walks (Figure 8—figure supplement 1) and deterministic escape behaviors (Figure 9). We propose, therefore, that the Stochastic Switch Model could serve as a multipurpose module for computing C. elegans behavior. Combining this mathematically tractable module with others representing sensory inputs, modulatory states, and the presumptive pattern generators for forward and reverse locomotion, could lead to essentially complete models of the C. elegans nervous system that are at once predictive and intuitively comprehensible (Abbott, 2008).
Materials and methods
Strains
All strains were cultivated at 22.5°C on lowdensity NGM (nematode growth medium) agar plates seeded with the E. coli bacteria (OP50) as described by Brenner (Brenner, 1974). Transgenic lines were made using standard protocols (Mello and Fire, 1995).
Experiment  Figure  Strains and genotypes 

Wild type  18  N2 
AVA → AVB synaptic current  5B  XL238 ntIs[Prig3::ChR2, Punc122::dsRed]; ntIs35[Psra11::tdTomato]; lite1(ce314) 
AVB → AVA synaptic current  5C  XL237 kyEx3801[Psra11::ChR2::GFP, Punc122::dsRed]; ntIs29[Pnmr1::tdTomato]; lite1(ce314) 
AVA ablation  4  N2 
AVD ablation  4  XL59 akIs [lin15(+); Pnmr1::GFP] 
AVE ablation  4  XL59 akIs [lin15(+); Pnmr1::GFP] 
AVB ablation  4  N2 
PVC ablation  4  XL59 akIs [lin15(+); Pnmr1::GFP] 
HYP A^{†}  6  DA572 eat4(ad572) 
HYP B^{†}  6  MT6308 eat4(ky5) 
HYP C^{†}  6  KP4 glr1(n2461) 
DEP A^{†}  6  VM1136 lin15(n765); akIs9 [lin15(+), Pglr1::GLR1(A/T)] 
DEP B^{†}  6  VM188 lin15(n765); akEx52[lin15(+), Pnmr1::GLR1(A/T) ] 
ASH → AVB synaptic current  8  XL194 ntIs27 [ Psra6::ChR2::YFP, Punc122::dsRed]; ntIs35 [Psra11::tdTomato]; lite1(ce314) 

^{†}Internal reference HYP A = HYP16, HYP B = HYP 56, HYP C = HYP20, DEP A = DEP14, DEP B = DEP19.
Physiological solutions
Request a detailed protocolExternal saline for electrophysiology (mM): 5 KCl, 10 HEPES, 8 CaCl2, 143 NaCl, 30 glucose, pH 7.2 (NaOH); internal saline for electrophysiology (mM): 125Kgluconate, 1 CaCl2, 18 KCl, 4 NaCl, 1 MgCl2, 10 HEPES, 10 EGTA, pH 7.2 (KOH). Medium for behavioral assays (mM): NH_{4}Cl 2, CaCl_{2} 1, MgSO_{4} 1, and KPO_{4} 25, pH 6.5; M9 Buffer (grams): 3 KH_{2}PO_{4}, 6 Na_{2}HPO_{4}, 5 NaCl, 1 ml 1 M MgSO_{4}, H_{2}O to 1 liter.
Behavior and tracking system
Request a detailed protocolPrior to each assay, an individual adult hermaphrodite was picked to a bacteriafree agar transfer plate by means of a platinumwire pick. The worm was then washed in M9 to remove excess bacteria, then transferred in a pipette filled with assay medium to a 10 cm petri plate containing 1.7% agarose in assay medium. A black dot approximately 40 microns in diameter was applied to the center of the body as shown in Figure 1A (see Spotting procedure). The worm was allowed to recover from transfer and handling for 2 min., then recorded for 10 min. The assay plate rested on a motorized microscope stage (Applied Scientific Instrumentation MS2000, Eugene, OR USA) fitted with position encoders (Gurely Precision Instruments LE1800, Troy, NY USA) having a resolution of 0.5 μm. Behavior was recorded using an analog video camera (CCD Sony XCST70, 29.97 frames per second) fitted with a 12× zoom lens (Navitar 50486D, Rochester, NY USA). For tracking purposes, video was analyzed in real time by custom software to calculate the eccentricity of the ink spot relative to the center of the field of view, and to compute the stage movements required to recenter the spot. Motion blur was minimized by making stage speed during corrective movements an increasing exponential function of target eccentricity such that small corrections were made more slowly than large corrections. Position encoders were read in synchrony with the video stream and this information was stored for offline analysis. The overall trajectory of the worm was computed by combining the location of the spot in the field of view with stage position in each video frame. The direction of movement (forward or reverse) at the start of each recording was keyed by the observer and subsequent assignments were made automatically by computer. Each recording was spotchecked for correct assignments at four or more points during the recording. In experiments involving neuronal ablations or genetic mutations, recordings of sham operated controls or wild type worms, respectively, were interleaved with worms in each treatment group.
Spotting procedure
Request a detailed protocolThe animal was immobilized by a stream of humidified CO_{2} emitted by a 1.5 mm diameter glass pipette positioned near the worm. The spotting ink was comprised of petroleum jelly (1 ml), mineral oil (1 ml), and black iron oxide (3 g). Ink was applied by means of 1.5 mm diameter glass rod that had been pulled to a fine point, fire polished to produce a bulbous tip, and dipped in the ink. The rod was positioned by means of a micromanipulator. To control for the effects of the spotting procedure, we compared the speed of locomotion of worms that had been immobilized, or immobilized and spotted, to untreated worms. There were no significant differences between these three groups.
Electrophysiology
Request a detailed protocolWorms were glued to an agarose coated coverslip using cyanoacrylate adhesive as previously described (Lindsay et al., 2011). The coverslip formed the bottom of the recording chamber, which was filled with external saline. The cell body of the neuron to be recorded was exposed by making a small slit in the cuticle using a finely drawn glass rod. Recording pipettes had resistances of 10–20 MΩ when filled with internal saline. Voltage and currentclamp recordings were made with a modified Axopatch 200A amplifier (Lockery and Goodman, 1998). In reversal potential measurements, recordings of photostimulationevoked synaptic currents were filtered at 2 kHz and sampled at 10 kHz. Postsynaptic neurons (AVA, AVB) were identified using a combination of fluorescent markers and distinctive voltage clamp currents as described (Lindsay et al., 2011). Presynaptic neurons (AVA, AVB, and ASH) were activated by expression of ChannelRhodopsin2 expressed under the control of neuronspecific promoters as described (see “Strains”). Worms were photostimulated in electrophysiological experiments using the blue channel (470 nm) of a dualwavelength LED module (Rapp OptoElectronic, Wedel, Germany) that was focused by a 63×, 1.4 NA oil immersion objective lens (Zeiss, part number 440762–9904). Irradiance (12.5 mW/mm^{2}) was determined by measuring the power emitted from the objective using an optical power meter placed above the front lens of the objective and dividing by the area of the field of illumination at the focal plane of the preparation.
Ablations
Request a detailed protocolNeurons were ablated using a laser as described previously (Bargmann and Avery, 1995). L1 larvae were mounted on 2.5% agarose pads containing 5–7 mM of the immobilizing agent NaN_{3}. AVA and AVB neurons were ablated in N2 animals and identified by position. AVD, AVE and PVC were ablated in animals expressing nmr1::GFP and identified by a combination of position and GFP expression. To limit potential behavioral differences in the two strains, we outcrossed (4×) the nmr1::GFP strain to the N2 strain used for AVA and AVB ablations. All animals were remounted 1–3 hr after surgery to confirm the ablation; those with collateral damage were discarded. Shamoperated animals were treated in the same manner except that the laser was not fired.
Statistical tests
Request a detailed protocolDescriptive statistics
The worm’s position in video frame $k$ is represented as the row vector:
where $x\left({t}_{k}\right)$ and $y\left({t}_{k}\right)$ are the coordinates of centroid of the tracking spot in the frame of reference of the agar plate, ${t}_{k}=k\u2206t$, $\u2206t=$33 ms, and $N\cong $18000 is the number of video frames analyzed in a continuous recording of one worm. We made the following definitions:
Row vectors
Request a detailed protocolVelocity:
Heading:
Scalar quantities
Request a detailed protocolSpeed:
Mean speed:
Instantaneous turn rate:
Mean heading change:
Speed autocovariance:
Velocity autocorrelation:
Heading autocorrelation:
Mean squared displacement:
Maximum likelihood estimation of state transition rates in a hidden Markov model
Request a detailed protocolTo analyze locomotory states we converted the velocity vector, $V\left(t\right)$, into a signed scalar quantity $v\left(t\right)$ that represents the component of velocity in the direction of the worm’s track, with positive values indicating forward movement. We first smoothed $x\left(t\right)$ and $y\left(t\right)$ using an 11 frame window, assigned a direction to the smoothed track with respect to the head/tail orientation of the worm, and projected $V\left(t\right)$ onto the smoothed track to obtain $v\left(t\right)$. For each cohort of worms we collected all $v\left(t\right)$ values into a single velocity distribution $g\left(v\right)$. The central peak of $g\left(v\right)$ was fit by a Cauchy distribution with median 0 and halfwidth $b$ = 18 µm/s (Figure 2—figure supplement 2), which we used to approximate the pause velocity distribution for states X and Y for all worms:
We used a Cauchy distribution because it has long tails that describe the pause velocity distribution better than a Gaussian distribution (i.e., the worm does not stop instantaneously when it switches from forward or reverse locomotion into one of the pause states). We estimated the forward and reverse velocity distributions ${g}_{\text{F}}\left(v\right)$ and ${g}_{\text{R}}\left(v\right)$ by scaling ${g}_{\text{P}}\left(v\right)$ to fit the peak at $v=0$, subtracting it from the overall distribution and splitting the remaining distribution into ${g}_{\text{F}}\left(v\right)$ for $v>0$ and ${g}_{\text{R}}\left(v\right)$ for $v<0$. Velocity distributions were scaled to be probability densities (area =1) and collected into a row vector:
where ${g}_{i}\left(v\right)$ is the estimated probability density that worms move at velocity $v$ when in state $i$.
The goal of the maximum likelihood fitting procedure is to find the set of state transition rates $\left\{{a}_{\text{XF}},{a}_{\text{FX}},{a}_{\text{XR}},{a}_{\text{RX}},{a}_{\text{FY}},{a}_{\text{YF}},{a}_{\text{RY}},{a}_{\text{YR}}\right\}$ that maximize the probability of the observed velocity time series $v\left(t\right)$ given the velocity distribution $G\left(v\right)$. All transition rates were constrained to be $\ge 0$, and usually were additionally constrained to correspond to valid synaptic weights as described below$.\text{}$ The likelihood is most conveniently calculated using matrix notation as follows; see Colquhoun and Hawkes, 1995 for a more complete explanation of these computations. Let:
Element ${q}_{ij}$ ($i\ne j)$ of matrix $Q$ is the transition rate from state $i$ to state $j$ (i.e., the instantaneous probability per unit time that the system in state $i$ will make a transition to state $j$, and element ${q}_{ii}$ is the negative of the total transition rate out of state $i$, which is related to the mean dwell time in state $i$ by:
Matrix $Q$ is composed of instantaneous transition rates, which can be converted into the matrix of transition probabilities during a brief time interval of duration $\epsilon $ by multiplying $Q$ by $\epsilon $ and adding 1 to each diagonal element (i.e., by calculating $\epsilon \xb7Q+I$, where $I$ is the 4×4 identity matrix). If $\epsilon $ is sufficiently small that multiple state transitions can be ignored, then element $ij$ of matrix $\epsilon \xb7Q+I$ is the probability that the system is in state $j$ at the end of a time interval of duration $\epsilon $ given that it was in state $i$ at the beginning of the interval. For longer time intervals during which multiple state transitions may occur, transition probabilities can be calculated by repeatedly multiplying matrix $\epsilon \xb7Q+I$ by itself. Thus, if
then $M$ is the matrix of transition probabilities during a time interval of duration $K\epsilon $. If $K$ and $\epsilon $ are chosen such that $\u2206t=K\epsilon $, then element $ij$ of matrix $M$ is the transition probability from state $i$ to state $j$ during one video frame of duration $\u2206t$. We chose $K={2}^{30}$ and let $\epsilon =\u2206t/K=30.7$ picoseconds, a time interval during which multiple state transitions can safely be ignored. Since $K$ was chosen to be a power of 2, $M$ could be rapidly and accurately calculated by 30 serial multiplications using 64bit floating point arithmetic.
Let $P\left(t\right)$ be the row vector of historydependent state probabilities:
where ${p}_{i}\left(t\right)$ is the probability of being in state i at time $t$ given $v\left(u\right)$ for all $u$ up to and including the present time ($u\le t$). The matrix product $P\left(t\right)\xb7M$ is the state probability vector at time $t+\u2206t$ prior to accounting for the observed velocity at time $t+\u2206t$. To account for $v(t+\u2206t)$ we used the information contained in $G\left(v(t+\u2206t)\right)$ and applied Bayes theorem:
where $diag\left[G\left(v(t+\u2206t)\right)\right]$ is the 4×4 matrix with the elements of $G\left(v(t+\u2206t)\right)$ along the diagonal, and $l$ is the scalar multiplicative factor required for the sum of the four elements of $P(t+\u2206t)$ to equal 1 (i.e., $P(t+\u2206t)$ is a vector of probabilities). Initially ($t=0$) we set $P\left(0\right)$ equal to the steadystate probability vector ${P}_{\infty}$, which is given by:
where ${U}_{4}$ is the $1\times 4$ row vector of ones and ${Q}_{a}$ is the $4\times 5$ matrix constructed by appending a column of ones to $Q$. To break the symmetry between the behaviorally indistinguishable states X and Y, we identified X as the state with higher steadystate probability.
We then calculated the loglikelihood, summed over all worms in the cohort:
where ${v}_{w}\left(t\right)$ is the velocity and ${P}_{w}\left(t\right)$ is the historydependent state probability vector of worm $w$ at time $t$.
We used a random optimization algorithm to find the set of transition rates that maximized $ln\left(L\right)$. Initial guesses for 6 of the 8 rates were chosen independently from log uniform distribution between 0.01 Hz and 10 Hz. The remaining 2 rates were calculated to satisfy the constraints needed to generate valid synaptic weights (see below). At each iteration, each of the 6 independently chosen rates was altered by adding a random number chosen from a Cauchy distribution with median 0 and width ${b}_{random}$ (initially ${b}_{random}=$ 0.01 Hz), and the remaining 2 rates were recalculated. To avoid getting trapped in local likelihood maxima, the new rates were rejected and another set was calculated if any of the new rates were <0.01 Hz. If the new rates generated an increased likelihood, the new rates were accepted and ${b}_{random}$ was increased by 3%. Otherwise the old rates were retained and ${b}_{random}$ was decreased by 0.5%. The procedure was iterated until ${b}_{random}<$ 0.001 Hz. The random optimization procedure was replicated 10 times for each cohort using different randomly chosen initial guesses. In 71% of the replicates the procedure converged on a set of transition rates in which none of the transition rates differed from the best set by more than 5%. The best set of transition rates was then refined by applying the optimization procedure using a success criterion of ${b}_{random}<$ 10^{5} and constraining transition rates to be ≥ 10^{–4} Hz.
The likelihood calculations described above use only past and present velocity observations to calculate $P\left(t\right)$, but once the optimal transition rates were determined, the ForwardBackward algorithm (Rabiner, 1986) can be used to yield a better estimate of the state probabilities based on past, present and future velocity observations, and the Viterbi Algorithm can be used to find the sequence of states with the highest probability of producing the observed velocities (Figure 2E).
Stochastic model neurons
Request a detailed protocolWe expressed the effect of synaptic inputs to command units $\mathcal{F}$ and $\mathcal{R}$ by equations of the form:
where ${a}_{\text{OFF}}$ is the transition rate from ON to OFF, ${a}_{\text{ON}}$ is the transition rate from OFF to ON, and $S$ is the total synaptic input to the unit. We do not attach any mechanistic significance to these equations, but note that they are analogous to the Arrhenius Equation (Stiller, 1989) an approximation commonly used to describe the rates of chemical reactions in terms of an activation energy, $E$:
where a is the reaction rate constant, $A$ is an empirically determined constant, ${k}_{B}$ is the Boltzmann constant, and $T$ is the absolute temperature. Under this interpretation, $S$ is analogous to activation energy expressed in units of ${k}_{B}T$. Thus, $\mathcal{F}$ and $\mathcal{R}$ are assumed to be symmetrical bistable units that change state at rate $A$ when $S=0$. Deviations from this baseline condition are modelled as external synaptic inputs ${h}_{\mathcal{F}}$ and ${h}_{\mathcal{R}}$.
We represented the total synaptic input onto units $\mathcal{F}$ and $\mathcal{R}$, respectively, by:
where ${b}_{\mathcal{F}}\left(t\right)$ and ${b}_{\mathcal{R}}\left(t\right)$ are the states of $\mathcal{F}$ and $\mathcal{R}$ (1 = ON, 0 = OFF), ${w}_{\mathcal{R}\mathcal{F}}$ and ${w}_{\mathrm{FR}}$ are the synaptic weights from $\mathcal{R}$ onto $\mathcal{F}$ and from $\mathcal{F}$ onto $\mathcal{R}$, respectively, and ${w}_{\mathcal{F}\mathcal{F}}$ and ${w}_{\mathrm{RR}}$ represent synaptic interactions among command neurons of the same class, plus any intrinsic membrane properties that may promote bistability. Applying these definitions to the rate constants in Figure 2C gives:
In these experiments the sensory environment was kept constant (e.g., no chemical or temperature gradients). Therefore ${h}_{\mathcal{F}}$ and ${h}_{\mathcal{R}}$ were assumed to be constant. For simulations of chemotaxis ${h}_{\mathcal{F}}$ and ${h}_{\mathcal{R}}$ varied with position in the chemical gradient.
Equations 31–34 express the 8 transition rates in terms of 6 parameters and yield the following two constraints on the transition rates:
The inverse relations between transition rates and synaptic parameters are:
Derivation of the mean distance traveled during a forward run
Request a detailed protocolThe time series of the worm’s locomotory states can be divided into forward runs, during which the worm is in either the F or P state, and reverse runs, during which the worm is in either the R or P state. Forward runs always begin with an RPF transition and end with the next FPR transition, which marks the beginning of a reverse run. Thus forward runs and reverse runs occur in strict alternation, such that the number of forward runs equals the number of reverse runs.
Let ${m}_{\text{F}}$ denote the mean distance traveled during a single forward run, assuming that forward runs are straight. The value of ${m}_{\text{F}}$ is most easily calculated by dividing time into nonoverlapping epochs, each of which begins with an RPF transition and ends immediately before the next RPF transition. Each epoch thus contains exactly one forward run, which includes all visits to state F during the epoch. Therefore, ${m}_{\text{F}}$ is also equal to the mean distance travelled while in the forward state during a single epoch:
where $\overline{{v}_{\text{F}}}$ is the mean velocity in the forward state and ${f}_{\text{RPF}}$ is the frequency of RPF transitions. Since FPR and RPF transitions occur in strict alternation they must occur in equal numbers: ${f}_{\text{RPF}}={f}_{\text{FPR}}$. Thus, eq. 39 can also be written with ${f}_{\text{FPR}}$ in the denominator, which is more useful for the calculation that follows, although the form shown above is more directly interpreted in terms of the frequency of random reorientations, which occur at the RPF transitions. It is straightforward to calculate ${f}_{\text{FPR}}$ given ${p}_{\text{F}}$, ${a}_{\text{FX}}$, ${a}_{\text{FY}}$, and the probabilities that the transitions out of states X and Y will be into state R:
Combining eqns. 39 and 42 yields:
An approximation to ${m}_{\text{F}}$ in terms of synaptic weights is obtained by noting that transitions from F to Y were extremely rare (${a}_{\text{FY}}=$0.007 s^{1}; Table 1). Setting ${a}_{\text{FY}}\cong 0$ yields:
Simulations of worm behavior
Request a detailed protocolIn Figure 8—figure supplement 1 and 2, the worm was represented as a point that moved forward or backward at speeds of 0.2 and 0.3 mm/sec, respectively, and was stationary during the pause state. Rate constants were calculated according to equations 31–34 based on the weights that pertain under random search or chemotaxis, using either $A={A}_{min}$ or $A={A}_{max}$. Weights were used to compute the state transition matrix $M$. At each time step (∆t = 33 ms), the next state was selected randomly according to the state probabilities given by $M$. When an RPF transition occurred, a new direction of movement (heading) was selected from a uniform distribution. The random component of the heading was modeled as Gaussian noise having a standard deviation of 0.001 degrees. In the case of chemotaxis simulations, the values of ${h}_{\mathcal{F}}$and ${h}_{\mathcal{R}}$ were updated at every time step according to the direction in which the worm was heading, leading to an updated set of weights and a new $M$ matrix.
Definition of modes of random search in C. elegans
Request a detailed protocolTo date, these behaviors have been defined mainly in operational terms. Following the terminology of Jander 1975: (i) cropping is the locomotory behavior exhibited by wellfed worms on plates with densely populated patches of bacteria; (ii) local search (also “area restricted search” (Hills et al., 2004) or “pivoting” (Wakabayashi et al., 2004) is exhibited by wellfed worms within about 10 min after being transferred to a foodless plate; and (iii) ranging (“dispersal” (Gray et al., 2005) or “traveling” (Wakabayashi et al., 2004) is exhibited by wellfed worms tens of minutes after being transferred to a foodless plate. Each mode can be associated with approximate ranges of three parameters: mean forward run length (${m}_{\text{F}}$), mean frequency of reversals (${f}_{\text{FPR}}$), and mean reverse run length (${m}_{\text{R}}$). Local search serves as a useful reference point. During cropping, ${m}_{\text{F}}$ is greatly reduced, ${f}_{\text{FPR}}$ is greatly increased, and ${m}_{\text{R}}$ is also reduced, being limited to “short reversals” (the distance traveled in one or two head sweeps, or about 0.5 mm (Gray et al., 2005); during local search, reverse runs are almost always “long” (the distance traveled in at least three head sweeps). During ranging, ${m}_{\text{F}}$ is greatly increased, ${f}_{\text{RPF}}$ is reduced, and reversals are long. Cutoff values for search modes, inferred from behavioral data (Wakabayashi et al., 2004; Gray et al., 2005; Fujiwara et al., 2002; Hills et al., 2004) were: Dwelling – short forward run length (${m}_{\text{F}}$ < 0.5 mm), high reversal frequency (${f}_{\text{FPR}}$ > 6.0/min), short reversals (${m}_{\text{R}}$ < 0.5 mm); Local search – moderate forward run length (0.5 mm ≤ ${m}_{\text{F}}$ < 5.0 mm), moderate reversal frequency (2.0/min ≤ ${f}_{\text{FPR}}$ < 6.0/min), nonshort reversals (${m}_{\text{R}}$ ≥ 0.5 mm); Ranging – long forward run length (${m}_{\text{F}}$ ≥ 5.0 mm), low reversal frequency (${f}_{\text{FPR}}$< 2/min), nonshort reversals (${m}_{\text{R}}$ ≥ 0.5 mm).
Data archive
Request a detailed protocolAll data and the analysis program are publicly available at doi:10.5061/dryad.35qv6.
Data availability

Data from: A stochastic neuronal model predicts random search behaviors at multiple spatial scales in C. elegansAvailable at Dryad Digital Repository under a CC0 Public Domain Dedication.
References

Laser killing of cells in Caenorhabditis elegansMethods in Cell Biology 48:225–250.https://doi.org/10.1016/S0091679X(08)613904

Scaleinvariant movements of fishermen: the same foraging strategy as natural predatorsEcological Applications : A Publication of the Ecological Society of America 17:331–337.https://doi.org/10.1890/060303

Lvy flights in dobe ju/’hoansi foraging patternssHuman Ecology 35:129–138.https://doi.org/10.1007/s1074500690834

Brain neurones involved in the control of walking in the cricket gryllus bimaculatusJournal of Experimental Biology 166:113–130.

The neural circuit for touch sensitivity in Caenorhabditis elegansThe Journal of Neuroscience 5:956–964.

Activity of reticulospinal neurons during locomotion in the freely behaving lampreyJournal of Neurophysiology 83:853–863.

A newly identified buccal interneuron initiates and modulates feeding motor programs in aplysiaJournal of Neurophysiology 90:2190–2204.https://doi.org/10.1152/jn.00173.2003

Illuminating neural circuits and behaviour in caenorhabditis elegans with optogeneticsPhilosophical Transactions of the Royal Society B: Biological Sciences 370:20140212.https://doi.org/10.1098/rstb.2014.0212

Neuronal microcircuits for decision making in C. elegansCurrent Opinion in Neurobiology 22:580–591.https://doi.org/10.1016/j.conb.2012.05.005

The mechanism of locomotion in snakesThe Journal of Experimental Biology 23:101–120.

A circuit for navigation in Caenorhabditis elegansProceedings of the National Academy of Sciences of the United States of America 102:3184–3191.https://doi.org/10.1073/pnas.0409009101

Nondauer larval dispersal in caenorhabditis elegansJournal of Experimental Zoology Part B: Molecular and Developmental Evolution 312B:224–230.https://doi.org/10.1002/jez.b.21287

Control of cricket stridulation by a command neuron: efficacy depends on the behavioral stateJournal of Neurophysiology 83:712–722.

Dopamine and glutamate control arearestricted search behavior in Caenorhabditis elegansThe Journal of Neuroscience 24:1217–1225.https://doi.org/10.1523/JNEUROSCI.156903.2004

Neural networks and physical systems with emergent collective computational abilitiesProceedings of the National Academy of Sciences of the United States of America 79:2554–2558.https://doi.org/10.1073/pnas.79.8.2554

Optimal foraging strategies: Lévy walks balance searching and patch exploitation under a very broad range of conditionsJournal of Theoretical Biology 358:179–193.https://doi.org/10.1016/j.jtbi.2014.05.032

Foraging success of biological Lévy flights recorded in situProceedings of the National Academy of Sciences of the United States of America 109:7169–7174.https://doi.org/10.1073/pnas.1121201109

Parallel use of two behavioral mechanisms for chemotaxis in Caenorhabditis elegansThe Journal of Neuroscience 29:5370–5380.https://doi.org/10.1523/JNEUROSCI.363308.2009

Ecological aspects of spatial orientationAnnual Review of Ecology and Systematics 6:171–188.https://doi.org/10.1146/annurev.es.06.110175.001131

Lowdimensional functionality of complex network dynamics: neurosensory integration in the Caenorhabditis Elegans connectomePhysical Review. E, Statistical, Nonlinear, and Soft Matter Physics 89:.https://doi.org/10.1103/PhysRevE.89.052805

Mate searching in Caenorhabditis elegans: a genetic model for sex drive in a simple invertebrateThe Journal of Neuroscience 24:7427–7434.https://doi.org/10.1523/JNEUROSCI.174604.2004

Common scaleinvariant patterns of sleepwake transitions across mammalian speciesProceedings of the National Academy of Sciences of the United States of America 101:17545–17548.https://doi.org/10.1073/pnas.0408242101

Tightseal wholecell patch clamping of Caenorhabditis elegans neuronsMethods in Enzymology 293:201–217.https://doi.org/10.1016/S00766879(98)930166

The computational worm: spatial orientation and its neuronal basis in C. elegansCurrent Opinion in Neurobiology 21:782–790.https://doi.org/10.1016/j.conb.2011.06.009

Stochastic process with ultraslow convergence to a Gaussian: The truncated Lévy flightPhysical Review Letters 73:2946–2949.https://doi.org/10.1103/PhysRevLett.73.2946

Action potentials contribute to neuronal signaling in C. elegansNature Neuroscience 11:865–867.https://doi.org/10.1038/nn.2131

DNA transformationMethods in Cell Biology 48:451–482.https://doi.org/10.1016/S0091679X(08)613990

Wholebrain calcium imaging with cellular resolution in freely behaving caenorhabditis elegansProceedings of the National Academy of Sciences of the United States of America 201507110.https://doi.org/10.1073/pnas.1507110112

The fundamental role of pirouettes in Caenorhabditis elegans chemotaxisThe Journal of Neuroscience 19:9557–9569.

Inducible and titratable silencing of Caenorhabditis elegans neurons in vivo with histaminegated chloride channelsProceedings of the National Academy of Sciences of the United States of America 111:2770–2775.https://doi.org/10.1073/pnas.1400615111

Detection and avoidance of a natural product from the pathogenic bacterium Serratia marcescens by Caenorhabditis elegansProceedings of the National Academy of Sciences of the United States of America 104:2295–2300.https://doi.org/10.1073/pnas.0610281104

Photoinducible cell ablation in Caenorhabditis elegans using the genetically encoded singlet oxygen generating protein miniSOGProceedings of the National Academy of Sciences of the United States of America 109:7499–7504.https://doi.org/10.1073/pnas.1204096109

Hyperactivation of Btype motor neurons results in aberrant synchrony of the Caenorhabditis elegans motor circuitThe Journal of Neuroscience 33:5319–5325.https://doi.org/10.1523/JNEUROSCI.401712.2013

An introduction to hidden markov modelsIEEE ASSP Magazine 3:4–16.https://doi.org/10.1109/MASSP.1986.1165342

A tutorial on hidden markov models and selected applications in speech recognitionProceedings of the IEEE 77:257–286.https://doi.org/10.1109/5.18626

Synaptic polarity of the interneuron circuit controlling C. elegans locomotionFrontiers in Computational Neuroscience 7:.https://doi.org/10.3389/fncom.2013.00128

Thermotaxis in Caenorhabditis elegans analyzed by measuring responses to defined Thermal stimuliThe Journal of Neuroscience 22:5727–5733.

Mechanistic analysis of the search behaviour of caenorhabditis elegansJournal of the Royal Society Interface 11:20131092.https://doi.org/10.1098/rsif.2013.1092

Mechanisms of locomotory speed change: the pteropod solutionAmerican Zoologist 41:1001–1008.https://doi.org/10.1093/icb/41.4.1001

Dietary choice behavior in Caenorhabditis elegansThe Journal of Experimental Biology 209:89–102.https://doi.org/10.1242/jeb.01955

Lévy flight and Brownian search patterns of a freeranging predator reflect different prey field characteristicsThe Journal of Animal Ecology 81:432–442.https://doi.org/10.1111/j.13652656.2011.01914.x

Stimulation of the mesencephalic locomotor region elicits controlled swimming in semiintact lampreysThe European Journal of Neuroscience 12:4081–4092.https://doi.org/10.1046/j.14609568.2000.00301.x

Emergence of long timescales and stereotyped behaviors in Caenorhabditis elegansProceedings of the National Academy of Sciences of the United States of America 108:7286–7289.https://doi.org/10.1073/pnas.1007868108

Dimensionality and dynamics in the behavior of C. elegansPLoS Computational Biology 4:e1000028.https://doi.org/10.1371/journal.pcbi.1000028

BookArrhenius Equation and NonEquilibrium Kinetics: 100 Years Arrhenius EquationBSB BG Teubner.

Highthroughput study of synaptic transmission at the neuromuscular junction enabled by optogenetics and microfluidicsJournal of Neuroscience Methods 191:90–93.https://doi.org/10.1016/j.jneumeth.2010.05.019

The neural network for chemotaxis to tastants in Caenorhabditis elegans is specialized for temporal differentiationThe Journal of Neuroscience 29:11904–11911.https://doi.org/10.1523/JNEUROSCI.059409.2009

Functional mapping of neurons that control locomotory behavior in Caenorhabditis elegansJournal of Neurobiology 56:178–197.https://doi.org/10.1002/neu.10245

Structural properties of the Caenorhabditis elegans neuronal networkPLoS Computational Biology 7:e1001066.https://doi.org/10.1371/journal.pcbi.1001066

Panneuronal imaging in roaming Caenorhabditis elegansProceedings of the National Academy of Sciences of the United States of America 113:.https://doi.org/10.1073/pnas.1507109113

The physics of foraging: an introduction to random searches and biological encounters, Cambridge University PressThe physics of foraging: an introduction to random searches and biological encounters, Cambridge University Press.

Error bounds for convolutional codes and an asymptotically optimum decoding algorithmIEEE Transactions on Information Theory 13:260–269.https://doi.org/10.1109/TIT.1967.1054010

Neurons regulating the duration of forward locomotion in Caenorhabditis elegansNeuroscience Research 50:103–111.https://doi.org/10.1016/j.neures.2004.06.005

Initiation maintenance and modulation of swimming in the medicinal leech by the activity of a single neuronethe Journal of Experimental Biology 77:71–88.

The structure of the nervous system of the nematode caenorhabditis elegansPhilosophical Transactions of the Royal Society B: Biological Sciences 314:1–340.https://doi.org/10.1098/rstb.1986.0056

A dynamic network simulation of the nematode tap withdrawal circuit: predictions concerning synaptic function using behavioral criteriaThe Journal of Neuroscience 16:4017–4031.

Thermal avoidance in Caenorhabditis elegans: an approach to the study of nociceptionProceedings of the National Academy of Sciences of the United States of America 96:10477–10482.https://doi.org/10.1073/pnas.96.18.10477

Reversal frequency in Caenorhabditis elegans represents an integrated response to the state of the animal and its environmentThe Journal of Neuroscience 23:5319–5328.
Decision letter

Ronald L CalabreseReviewing Editor; Emory University, United States
In the interests of transparency, eLife includes the editorial decision letter and accompanying author responses. A lightly edited version of the letter sent to the authors after peer review is shown, indicating the most substantive concerns; minor comments are not usually included.
Thank you for submitting your work entitled "A stochastic neuronal model predicts random search behaviors at multiple spatial scales in C. elegans" for consideration by eLife. Your article has been reviewed by favorably evaluated by Eve Marder (Senior Editor) and three reviewers, one of whom, Ronald L. Calabrese, is a member of our Board of Reviewing Editors.
The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.
Summary:
The authors present a kinematic analysis of random local search behavior in C. elegans at high resolution. This analysis leads to a Markov model of underlying neuronal network based on the C. elegans connectome. This inherently stochastic model not only accounts for the data but also predicts the sign and strengths of key connections which are then confirmed by electrophysiology/optogenetics. The model can be adjusted easily to account for other types of random searches by adjusting parameters compatible with sensory input or modulation. Further the model can be expanded to incorporate directed searches as in chemotaxis and can also account for 'deterministic' behavior such as escape movements. The model also accounts for counter intuitive results on dwell times associated with genetic manipulation and laser ablation of key command neurons. These findings conceptually inform experiments in the mammalian sleep network. Altogether this work is an impressive display of the power of modeling when combined with detailed experimental manipulation. The authors make a strong effort to cast their results in the context of other efforts to measure search motion and define pause states. They also explicitly address more neuronalbased determinist models of the network underlying random searches.
Essential revisions:
1) Results, second paragraph: The claim that this study represents "a 10fold improvement over previously published tracking systems" seems to be a stretch. In some sense, centerlinetracking experiments operate at a much finer resolution than tracking a single dot on a worm, with papers such as Brown, et al., PNAS (2013) tracking many more worms with more detail and higher resolution. I agree that the authors appear to have done a fine job at spot tracking, but we don't believe that they have a claim to novelty here.
2) The authors should address the Stephens et al. (2011) paper in a different manner than they do here (particularly in the Discussion section, fifth paragraph). To begin with, the statement that the "mathematical relationship between tangential velocity and phase velocity (so defined) has not been delineated, but is likely to be complex" largely due to slippage between the worm and the substrate seems an overstatement. If the paper is interested in modeling the animal's neural control of motion, then shouldn't the model be more concerned about the dynamics actually being controlled – in this case, the postural dynamics? And while it is indeed theoretically possible to produce thrust without advancing the phase, the collection of papers from Stephens et al. show that the worms' dynamics lie almost entirely on a manifold of remarkably small dimension, showing that these types of potential postural changes occur only rarely. Moreover, the authors themselves admit (in the subsection “Wild type locomotion”, fifth paragraph) that "on an agar surface, the worm moves without slipping."
In addition, the authors state that the Stephens et al. (2010) paper shows that postural modeling does not accurately model the worm when its center of mass trajectory follows an arc. In fact, the cited paper shows the necessity of looking at arcing trajectories between reorientation events and not just using a runandtumble analogy from bacteria, showing that shapespace dynamics form a predictive relationship with foraging trajectories.
All this being said, we are not disputing the authors' modeling choice of only using the midpoint of the worm's body. We have no arguments that there is utility in using a simplified description of a system to gain quantitative insight, but we want to see the authors distinguish themselves from the Stephens (2011) paper differently, focusing instead on the fundamental difference of modeling via a hidden Markov model instead of fitting parameters in a set of deterministic/stochastic differential equations. The choice of measuring Euclidean velocity instead of phase velocity is a modeling choice. Put another way, if someone with comparable amounts of shapespace data were to fit a hidden Markov model of your form to their data using phase velocity instead of velocity, would they likely obtain the same results?
3) In Figure 4, could the result that there is no postural stereotypy entering into state $X$ be the result of an overzealous "pause" caller. Specifically, this model calls many more things pauses than other approaches because of the nature of the subframeresolution hidden Markov fit. What do these plots look like if only the longest visits to the $X$ states are included?
4) Is there any evidence that the time scales of neural activity in C. elegans' motor neurons can be as short as the pause states being predicted? Although the model is supposed to be an abstraction of what's occurring in the worm's nervous system, one should expect to predict reasonable numbers for this nonetheless.
5) In the first paragraph of the subsection “Wild type locomotion” the authors claim that their "tracking system is capable of revealing briefer visits" to the pause state and that this is the reason why their measured dwell time is much shorter than previously measured. Again, we are not yet convinced that this method is much more accurate than the posture tracking used in Stephens et al. (2011) to measure the forward dwell time. The Stephens paper, however, uses a slightly different definition of the dwell time, using cessation of the forward phase velocity instead. Although the paper here does not measure phase velocity, they could measure the dwell time between visits to the reversal state (i.e. ignoring $F$>$P$>$F$ type transitions). What happens if this definition is used? Does the same result emerge?
6) A more general comment is that we would like to understand more what type of neural activity the authors would predict based on their model, connectome data, and polarity results from papers such as Rakowski (2013) despite imaging in a freelybehaving worm being well outside of the scope of this paper (although not nearly as impossible as the authors seem to suggest in the Discussion section). If this experiment were performed, what is the most likely neural instantiation of this model that is consistent with the current literature? The strength of this model is that it makes an attempt at getting at how this circuit may actually function. Accordingly, if possible, a concrete prediction or predictions to this end would greatly increase the value of this paper to researchers and could guide experimentalists performing imaging in freelybehaving animals in their measurements and analysis.
[Editors’ note: a previous version of this study was rejected after peer review, but the authors submitted for reconsideration. The previous decision letter after peer review is shown below.]
Thank you for choosing to send your work entitled "A stochastic neuronal flipflop circuit regulates random search during foraging behavior" for consideration at eLife. Your full submission has been evaluated by Eve Marder (Senior editor) and three peer reviewers, one of whom, Ronald L. Calabrese, is a member of our Board of Reviewing Editors, and the decision was reached after discussions between the reviewers. Based on our discussions and the individual reviews below, we regret to inform you that your submission in its present form is not suitable for publication in eLife. That said, if you feel that you can adequately rebut, answer the reviews with a combination of rewriting or new work, we would be willing to consider a new submission of this material, as it clearly asks an important and interesting set of questions.
Consensus review:
The authors present a kinematic analysis of random local search behavior in C. elegans. The then use a Hidden Markov Model to quantitatively model such searches based their kinematic analyses. They then identify this HMM with known aspects of the worm connectome to understand the neural implementation of switches between different behavioral states. In particular, the model contains two populations of neurons, each controlling either forward or reverse locomotion, with the overall network resulting in four possible behaviors: forward, reverse, and two pause states. They perform a maximum likelihood fit of the model to measured time series, making predictions about the effect of neural ablations, the sign and strength of synaptic connections, effects of perturbations to membrane potentials. Moreover, they phenomenologically model a potential underlying neural mechanism for different search modalities.
The subject matter of the paper is definitely one of broad interest, and if the model is correct it would provide us with a substantive understanding into the neural implementation of the forward/pause/reversal/pause(?) dynamics it would be a major contribution. However, the reviewers are not fully convinced that the model is the simplest possible explanation of the animal's behavior and that the authors have shown that a clear connection can be drawn to neural correlates.
Enthusiasm for the subject matter and the potential neural correlates vs. concerns about the appropriateness of the approach and its rigor were weighed differently by the reviewers initially, but in consultation the concerns predominated. The detailed reviews of the expert reviewers are appended but their concerns are summarized below.
1) The model is not adequately placed within the context of previous work in the field. There have been other papers attempting to understand both exploitation/exploration behavioral generation (e.g. Flavell et al., Cell, 2013) and forward/reversing dynamics (e.g. Stephens, PNAS, 2011).
Moreover, a direct comparison to (Rakowski et al., 2013) is needed. This work, which is cited as Ref. 12, also develops a model of the command neuron network and uses this model to predict behavioral states and the effects of ablations.
2) The existence of two different pause states $X$ and $Y$ is not well supported in the data and the HMM does not have the power to confirm their existence.
3) The mapping of the HMM onto the underlying neuronal network is not convincing. In particular there little data in this paper or in any of the citations to support the contention that pause state $X$ = all command neurons off and pause state $Y$ = all command neurons on. Mapping synaptic weights onto coupling in the HMM is not rigorous.
For the paper to be suitable for eLife the following extensive changes would have to be made.
1) Relax the claims about the model's prediction of synaptic weights.
2) The authors need to establish firmly that tracking a single point provides an unbiased estimate of trajectory in light of Stephens, PLoS One (2010). Given their extensive data set, we presume that converting their point imaging into centerline tracking would be extremely difficult and time consuming, but the authors should at least place their analysis in this light. They should also bolster their arguments that there are two behaviorally distinct pause states and discuss how they are related to the two pause states of Stephens.
3) The authors should differentiate their work with respect to Rakowski (2013). Does their model provide different predictions? Can they be disambiguated? It would be great if they could show this, but minimally, they should suggest what experiments should be performed in the future.
4) Similarly – what measurements should be made to distinguish $X$ and $Y$ pause states map onto the electrical activity of the command neuron network at least at some coarsegrained level?
Reviewer #1:
The authors present a kinematic analysis of random local search behavior in C. elegans at high resolution. This analysis leads to a Markov model of underlying neuronal network based on the C. elegans connectome. This inherently stochastic model not only accounts for the data but also predicts the sign and strengths of key connections which are then confirmed by electrophysiology/optogenetics. The model can be adjusted easily to account for other types of random searches by adjusting parameters compatible with sensory input or modulation. Further the model can be expanded to incorporate directed searches as in chemotaxis and can also account for 'deterministic' behavior such as escape movements. The model also accounts for counter intuitive results on dwell times associated with genetic manipulation and laser ablation of key command neurons. These findings conceptually inform experiments in the mammalian sleep network.
Altogether this work is an impressive display of the power of modeling when combined with detailed experimental manipulation. Moreover the approach complements and expands other efforts using unbiased approaches that make no assumptions about whether and how behavioral categories can be mapped to specific states of the nervous system. This approach allows mapping of states on to the neuronal network. The paper is very clearly and crisply written. The figures are clear and contain useful data well organized. All essential data is presented.
This reviewer is not an expert in Markov models and the associated mathematics, so I would defer to the experts on any potential flaws in this analysis, but the behavioral analyses are very straightforward and carefully done, the electrophysiology is state of the art for C. elegans, and the mutant analysis and cell ablation experiments seem carefully controlled. The statistical analysis seems appropriate but I defer to the other more expert reviewers. Methods are a model of completeness. There are no major concerns noted by this reviewer. I really enjoyed reading this paper and I learned a lot.
Reviewer #2:
Roberts et al. present a study combining behavioral measurements of C. elegans locomotion, modeling of a twostate stochastic system, and measurements of functional connectivity between two command interneurons. The authors claim that using the stochastic model, they can make nontrivial predictions about the state of the underlying neural network, solely based on behavioral measurements. They also claim that using behavioral data and their twostate model they can predict the strengths of synaptic connections between interneurons. If these claims were valid, this would be a significant achievement and worthy of publication in eLife, although the novelty of this work is somewhat compromised by another recent publication (Rakowski et al., 2013).
My major concerns are:
1) The claim that there are two distinct pause states resulting from different levels of activity in the command interneurons is not supported by evidence presented in this paper or cited work.
2) The stochastic switch model is not grounded enough in the biology of the system to justify comparing model fit parameters to measurements made on individual neurons. In particular the comparison of the model fit parameters labeled "synaptic weights" to actual synaptic connections is inappropriate.
Behavioral Measurements and Hidden Markov Model:
In Figure 1, the authors present data obtained by tracking a fiducial marker painted onto the back of a freely moving worm. They show a distribution of "velocities" (1D) that appear to be the sum of 3 distributions, which they label "reverse, pause," and "forward." In Figure 1E, they then show a velocity vs. time graph that has been segmented into these states by thresholding the velocity at +/ 0.05 mm/s. In Figure 2E, they show a similar velocity vs. time graph that has been similarly segmented using a Hidden Markov Model. The claim is that this model does a much better job of describing the behavioral state than the process used in Figure 1E, although I cannot find any quantitative measures in the paper to support this assertion (more on this later).
The HMM differs from simple thresholding by (1) assigning different prior probabilities to transitions between states and staying in the same state and (2) including 2 pause states (called $X$ and $Y$). The states of the HMM are related to behavior and to the stochastic flip flop model as follows.
$F$: worm is moving forward. "$F$" unit is active and "$R$" unit is inactive. What $F$ and $R$ represent is unclear in the paper? They are described as "binary stochastic elements," (subsection “The Stochastic Switch Model”, third paragraph), "single binary neuronlike units" (Figure 2 caption), and as a collection of neurons as in "self connections represent synaptic connections between neurons comprising a given unit." (Figure 2 caption). I will come back to this under Stochastic Switch Model, but for discussion of the hidden Markov model, the distinctions are not important.
$R$: worm is moving backward. "$R$" unit is active and "$F$" unit is inactive.
$X$: worm is not moving/moving slowly. Both units are inactive.
$Y$: worm is not moving/moving slowly. Both units are active.
A great deal of the argument in this paper relies on the existence of two pause states and of the identification of these two behaviorally identical pause states with specific patterns of activity in command interneurons, so I will go over in some detail the arguments the manuscript advances:
1) Using a speed threshold, pauses during forward>reverse transitions were found to be longer than pauses during reverse>forward transitions.
2) "When both $F$ and $R$ are OFF (state $X$) it is assumed that movement ceases, consistent with studies showing that genetic ablation or silencing of all command interneurons induces prolonged pauses (Kawano et al., 2011; Zheng et al., 1999)." Zheng et. al (Zheng et al., 1999) report that glr1::ICE (klys36) "moved significantly slower than wild type worms and also had long pauses during which no movement occurred." Kawano et al. (Kawano et al., 2011) report that silencing all of AVA, AVB, AVE, AVD, and PVC led to prolonged pausing with the body in a straight position and that using a combination of genetic and laser ablation to destroy these same methods resulted in worms that were "kinked" with odd body postures.
3) "In the event that $F$ and $R$ are simultaneously ON (state $Y$), the resulting motor commands are assumed to conflict at the level of the motor neurons of the body wall muscles, resulting in a second motionless state, as has been observed (Kawano et al., 2011)." The way this sentence is structured may leave a mistaken impression of what is reported in Kawano et al., 2011. Kawano et al., 2011 does not show that when $F$&$R$ command interneurons are simultaneously active, worms stop moving. Instead it shows that in innexin mutants (unc7 and unc9), the "kink" state is correlated with the A and B type motorneurons having similar levels of calcium.
4) In the second paragraph of the subsection “Wild type locomotion“. Removing state $Y$ (i.e., setting ${a}_{\mathrm{FY}}={a}_{\mathrm{RY}}=0$) caused a large and highly significant reduction in the summed loglikelihood of the 5 wild type cohorts (𝑝 < 10^{100}) demonstrating that the second pause state greatly improves the fit.
Taken individually or together, these arguments fail to support the existence of multiple pause states or that these states can be reached by coactivation or inactivation of command neurons.
Argument 1 merely shows that the speed of a dot on a worm's midline recovers more quickly in the reverse to forward transition than the forward to reverse. This could be due to measurement artifacts, biomechanical differences, differences in the response latencies of A and B motorneurons, or differences elsewhere in the neural network.
Argument 2 shows that permanent inactivation of all command motor neurons leads to a worm with severe motion defects. It's quite a stretch to go from there to saying that transient decreases in the activities of these neurons will lead to immediate cessation of movement. It does seem at least plausible that if both forward and reverse command neurons are "off," the worm stops moving, although I question whether this cessation would really be immediate, especially since forward motion results from the propagation of a wave down the body due to proprioreceptive feedback (Wen et al., 2012)
Argument 3 does not show that coactivation of command neurons leads to pausing in wild type animals.
Argument 4 suffers from several flaws. First, the hidden Markov model used with the sole observable being the speed of a dot on the worm's midline may not be an appropriate description of the system (discussed in other comments). Second, the fact that one model produces a better fit than another does not show that either model is true. IE: that adding a state to the HMM improves the model's fit does not prove that state exists.
There are other ways we might improve the model fit that are incompatible with the stochastic switch model: Allowing the HMM to have 8 degrees of freedom rather than 6; adding a third pause state or a second forward state with a different velocity distribution; allowing direct transitions between $X$ and $Y$ or $F$ and $R$. If any of these produce a significant improvement in the likelihood of the data given the model, what would that imply for the stochastic switch model? I would also like to see the likelihood ratio test applied to compare a model where the state is determined by simple thresholding (as in Figure 1E) to the 4state HMM.
In summary, the authors present only weak evidence that there even are two distinct pause states. There is no data presented here or in the cited literature that show that simultaneous activation or deactivation of forward and reverse command neurons causes pausing. There is no evidence given for the correspondence of the two pause states in the HMM with two distinct states of the command neuron network. Finally, beyond a very thin justification about metabolic efficiency, there is nothing to support the assignment of state $X$ to the $R$>$F$ transition and state $Y$ to the $F$>$R$.
I think the idea that the same behavioral state (pausing) can be the result of opposite patterns of activity in the command neurons is intriguing and worth pursuing.
The authors write that "The currently available set of optical probes of neuronal activity do not have sufficient temporal resolution to allow direct observation of the underlying states of the command network in intact, freely moving animals." The recovery time of GCaMP3 is ~300 ms and faster dynamics can be inferred by deconvolving the observed signal with the known GCaMP3 filter (Kato et al., 2014). And faster indicators, like GCaMP6f, are now available. The mean duration of inferred state $X$ pauses is 480 ms and 20% of inferred state $Y$ pauses are over 300 ms (Figure 2—figure supplement 3). The authors assert that in state $X$ the activity in all command neurons will be much lower than in state $Y$. Surely in these longer pauses it would be possible to make a measurement that verifies this prediction.
A second approach would be to use optogenetic tools to determine the pause state. In state $X$ (all command neurons off), ChR2 activation of AVB should induce forward locomotion by hyperpolarization of AVA via NpHR/halo should not. In state $Y$ (all command neurons on), hyperpolarization of AVA should lead to forward locomotion.
The "Stochastic Switch Model":
The stochastic switch model uses a nonstandard description of neural state. Instead of parameterizing the state of a neuron by a continuous variable, e.g. membrane potential, the neuron is treated as a binary element, existing in one of two discrete states. Synaptic connections are also represented/modeled in a nonstandard way. Normally, if A is presynaptic to B, activity in A causes a change in the membrane potential of B by an amount set by the weight of the synaptic connection from A to B. In this model, if A is active, the probability that B will spontaneously switch states is modified by an amount that depends exponentially on the "weight" of the connection. As a model of coupled two state systems, there is nothing inherently wrong with this approach, but comparing the resulting model fit parameters to biological parameters is confusing and potentially misleading.
Figure 5, "the stochastic switch model correctly predicts the sign and strength of synaptic connections" shows the confusion generated by the choice of terminology. Figure 5A shows model fit parameters that have no obvious biological interpretation, but the figure caption identifies these fit parameters as "synaptic weights." Then Figure 5B and C show "synaptic currents" (currents that result from functional, but not necessarily direct synaptic, connections) measured by recording the amount of current received by neuron B when neuron A is activated. The figure invites the reader to compare 5A to the 5B5E even though 5A shows fit parameters from a model that does not even have synaptic currents. This presentation is at best confusing and at worst misleading.
In their description of the stochastic switch model, the authors should only use the term "binary stochastic elements" to describe $F$ and $R$, and choose a term other than "synaptic weights" to describe the model fit parameters (hF, wff etc.), which are not synaptic weights as the term is commonly used. The chosen terminology will lead all but the most careful readers to misunderstand the work that has been done.
Measurements of functional connectivity, effect of ablations, and extensions to chemotaxis:
The measurements of functional connectivity appear to be well done; I have no major concerns with these, beyond that already noted in the text: different levels of ChR2 expression in AVA and AVB could complicate the interpretation of these results. Perhaps the current evoked by ChR2 activation could be directly measured and calibrated. In most of the literature about the C. elegans motor circuit, there is assumed to be reciprocal inhibition between AVA and AVB, but I think Figure 5 shows the first direct measurement, a significant contribution.
The rest of the paper (Figures 4, 6–8) relies on the results of the modeling and behavioral experiments.
Reviewer #3:
The submission by Roberts, et al., uses a Hidden Markov Model to quantitatively model the locomotory behavior of C. elegans, with a particular emphasis on understanding the neural implementation of switches between different behavioral states. In particular, the model contains two populations of neurons, each controlling either forward or reverse locomotion, with the overall network resulting in four possible behaviors: forward, reverse, and two pause states. They perform a maximum likelihood fit of the model to measured time series, making predictions about the effect of neural ablations, the sign and strength of synaptic connections, effects of perturbations to membrane potentials. Moreover, the phenomenologically model a potential underlying neural mechanism for different search modalities.
In general, I found the work to be careful done, and understanding the mechanisms of behavioral control is an important question in the field. My major concerns have to do with the novelty of the work.
As far as novelty is concerned, there have been other papers attempting to understand both exploitation/exploration behavioral generation (e.g. Flavell, et al., Cell, 2013) and forward/reversing dynamics (e.g. Stephens, PNAS, 2011), with both of the cited examples leading to the generation of long behavioral time scales. It would be useful for the authors to place their results in the context of these papers. In particular, the Stephens (2011) paper uses an analogous (but nonidentical approach) of fitting a dynamical system to the data and makes (highly accurate) predictions about forward/reversal switching rates using Arrheniuslike stochastic dynamics. In addition, an earlier paper from the same group (Stephens, PLoS Comp Bio, 2008) shows that a 4state system emerges naturally from a data set. Here, this setup is assumed a priori.
Specific comments:
1) In Figure 4B, I find the excuse for omitting error bars to be relatively weak. Minimally, bootstrapped confidence intervals could be found, which would be significantly better than the current lack of any measure of variability.
2) In the fifth paragraph of the subsection “Genetic effects on command neuron function” and Figure 6. The claim is made that "changes in synaptic strength dominate the effects of changes in membrane potential on Δ$F$ and Δ$R$." This certainly seems true for Δ$F$, but both models appear to give moreorless identical predictions here for Δ$R$.
3) In the first paragraph of the subsection “The Stochastic Switch Model and other random search behaviors“. It is claimed that "changes are required in at least three weights to produce the full range of search behaviors." This is not actually shown anywhere in the paper, just that three weights are sufficient to induce the full range. A supplementary figure to this effect would be useful.
4) The predicted value for the dwell time in the forward state is not comparable to the rest of the literature, as a different definition is used here. What would happen if a standard definition was used? Does the same rate result?
https://doi.org/10.7554/eLife.12572.032Author response
Essential revisions: 1) Results, second paragraph: The claim that this study represents "a 10fold improvement over previously published tracking systems" seems to be a stretch. In some sense, centerlinetracking experiments operate at a much finer resolution than tracking a single dot on a worm, with papers such as Brown, et al.,
PNAS (2013) tracking many more worms with more detail and higher resolution. I agree that the authors appear to have done a fine job at spot tracking, but we don't believe that they have a claim to novelty here.
None of the published centerline tracking papers cite the spatial resolution of the method, so the argument against our claim to novelty seems to be unfounded. However, we are willing to limit this claim to measurement oftangential velocityand have adjusted the manuscript accordingly (Introduction, last paragraph). Note that our demonstration of higher velocity resolution than other methods is significant because it provides with an explanation for why we obtained a lower value of mean dwell time in the forward state than previous studies (subsection “Wild type locomotion”, first paragraph).
In preparing this manuscript, we reviewed all C. elegans centerline tracking papers, including Brown et al., PNAS (2013) cited above, and the complete set of the Stephens/Bialek papers. As far as we can tell, none of them presents an explicit statement of spatial resolution or the information needed to compute spatial resolution (pixel size on the camera chip, and overall optical magnification), nor do any of them show data from motionless worms (Figure 1—figure supplement 1), which would at least yield an estimate of resolution in the limiting case of no worm movement. Thus, we find no basis in the literature for the referee’s statement that centerline methods have higher spatial resolution.
Whereas we agree that centerline tracking is superior when the worm’s posture is the key datum, here our focus is on tangential velocity, for which centerline tracking is not the most direct measure possible. As we noted in the main text of the paper (Results, end of first paragraph), the centerline method computes tangential velocity by following virtual points rather than physical marks on the body. This is done by fitting a spline to the worm’s boundary, breaking the spline into a fixed number of segments, and computing the velocity of segment endpoints. The spline method is inaccurate because the tail of the worm, being faint, is differentially resolved from frame to frame. As a result, the length of the spline, and thus the distance between segment endpoints varies significantly (e.g. Figure 1B of Karbowski et al., 2008; Figure 2 of Cronin et al., 2005 where centerline head and tail are chopped off and misaligned in comparison to real worm; Figure 1 panel B of Karbowski et al. 2006). Thus the spline method adds segmentation noise to the velocity measurements. In contrast, the spot tracking method follows a sharply defined, highcontrast mark that is essentially unaffected by image segmentation noise. Our literature review revealed only one other C. elegans tracking paper that cites the accuracy of spot of tracking. Here the tracking object is a fluorescent neuron and the resolution is given as 5 µm (Guo, 2009). This paper is the basis of our statement that our method affords a 10fold improvement, as the stage encoders we used had a resolution of 0.5 µm. Our dead worm experiments verify that encoder precision is the limiting factor in the overall spatial resolution.
2) The authors should address the Stephens et al. (2011) paper in a different manner than they do here (particularly in the Discussion section, fifth paragraph).
We agree with the reviewers’ main points (see details below): (i) Worm's do not slip under our conditions, nor under those of Stephens et al. (ii) Thrust can be produced when phase velocity is zero. (iii) Midpoint tracking is a useful simplified description of the system. The main force of this critique is thus the request to distinguish our work from Stephens' in a different manner – by using the thought experiment in which phase velocity is substituted for tangential velocity when fitting and analyzing our Markov model. Would the two definitions of velocity yield the same results? We have now adopted this excellent suggestion in the Discussion (fourth paragraph).
To begin with, the statement that the "mathematical relationship between tangential velocity and phase velocity (so defined) has not been delineated, but is likely to be complex" largely due to slippage between the worm and the substrate seems an overstatement.
We have removed the quoted statement, and the discussion of slippage. (As point of fact, however, we did not say that the complexity of the relationship between tangential and phase velocities is due to the possibility of slippage. Rather, this complexity is due the fact that tangential velocity in the worm is fully determined by thrust, and the higher order components of the worm’s shape, which are capable of generating thrust, are ignored by definition in computing phase velocity.)
If the paper is interested in modeling the animal's neural control of motion, then shouldn't the model be more concerned about the dynamics actually being controlled – in this case, the postural dynamics?
We could just as easily argue that from the perspective of the command neurons, what is being controlled is the presence or absence of thrust, and its direction.
And while it is indeed theoretically possible to produce thrust without advancing the phase, the collection of papers from Stephens et al.show that the worms' dynamics lie almost entirely on a manifold of remarkably small dimension, showing that these types of potential postural changes occur only rarely.
Shape components 3 and 4 account for approximately 30% of the shape variance (Stephens et al., 2008, Figure 2B), so it is hard to see how they can be characterized as “rare.” Moreover, these two components meet necessary and sufficient conditions for generating thrust: a gradient of curvature along the worm’s centerline (Gray 1946; Gray 1953; Gray 1964); Thus there are nonrare components that are generating thrust.
Moreover, the authors themselves admit (in the subsection “
Wild type locomotion”, fifth paragraph) that "on an agar surface, the worm moves without slipping."
Yes, we do believe that slip is generally minimal. Now, we no longer use slippage as a means of illustrating the differences between the two methods of analyzing behavior.
In addition, the authors state that the Stephens et al.(2010) paper shows that postural modeling does not accurately model the worm when its center of mass trajectory follows an arc.
We are not saying that postural modeling – defined in terms of all shape components – is inaccurate; we are actually making a much more restricted assertion. We are merely saying that phase velocity, because it is based on components 1 and 2 alone, is not necessarily an accurate measure of tangential velocity. In particular, this mismatch is accentuated in arcs, including omega turns, where other thrust producing shape components are strongly engaged.
In fact, the cited paper shows the necessity of looking at arcing trajectories between reorientation events and not just using a runandtumble analogy from bacteria, showing that shapespace dynamics form a predictive relationship with foraging trajectories.
We agree with this point. The work by Iino and colleagues (Iino, 2009) has shown that spatial orientation behavior is the product of two mechanisms: runandtumble (klinokinesis) and weathervaning (klinotaxis). However, klinotaxis, being deterministic rather than random, is outside the scope of the paper. We now make it clear that the type of chemotaxis modeled here is “randomwalk chemotaxis” (Discussion, first paragraph).
All this being said, we are not disputing the authors' modeling choice of only using the midpoint of the worm's body. We have no arguments that there is utility in using a simplified description of a system to gain quantitative insight, but we want to see the authors distinguish themselves from the Stephens (2011) paper differently, focusing instead on the fundamental difference of modeling via a hidden Markov model instead of fitting parameters in a set of deterministic/stochastic differential equations. The choice of measuring Euclidean velocity instead of phase velocity is a modeling choice. Put another way, if someone with comparable amounts of shapespace data were to fit a hidden Markov model of your form to their data using phase velocity instead of velocity, would they likely obtain the same results?
We would not get the same results (rate constants, weights, etc.) because thrust (hence movement) can be generated even when phase velocity is zero (see point (ii) of agreement above). Thus some of the pauses identified phase velocity, would not be identified as pauses by tangential velocity because the worm is still moving. Would all pauses in tangential velocity be detected by the phase velocity method? We speculate that this is mainly true, because when tangential velocity is zero, net thrust is zero, so none of the shape components are generating thrust, including components 1 and 2.
3) In Figure 4, could the result that there is no postural stereotypy entering into state $X$ be the result of an overzealous "pause" caller. Specifically, this model calls many more things pauses than other approaches because of the nature of the subframeresolution hidden Markov fit. What do these plots look like if only the longest visits to the $X$ states are included?
We tested for this possibility by reanalyzing the average posture at $\mathrm{FX}$ transitions after eliminating all but long pauses, and found a similar lack of postural preference at $\mathrm{FX}$ transitions. This result held for any minimum pause duration up to 2 s; longer dwells in state $X$ were too rare to analyze. We have added this result to the text (subsection “Wild type locomotion”, last paragraph).
4) Is there any evidence that the time scales of neural activity in C. elegans' motor neurons can be as short as the pause states being predicted? Although the model is supposed to be an abstraction of what's occurring in the worm's nervous system, one should expect to predict reasonable numbers for this nonetheless.
To our knowledge, there is only one report of voltage recordings from body wall muscle motor neurons in C. elegans (Liu, 2014). Interestingly, these motor neurons, which happen to be in the pool of reverse motor neurons, exhibit distinct up and down states. Theoretically, the down states could be pauses of the state$X$ type in our model when they coincide with down states in forward motor neurons, whereas the up states could be pauses of the state$Y$ type in our model when they coincide with up states in forward motor neurons. Unfortunately, no simultaneous recordings of reverse and forward neurons were made so it is impossible to determine the duration of possible pause states. The Liu paper is now referenced in the third paragraph of the subsection “The Stochastic Switch Model”.
5) In the first paragraph of the subsection “Wild type locomotion” the authors claim that their "tracking system is capable of revealing briefer visits" to the pause state and that this is the reason why their measured dwell time is much shorter than previously measured. Again, we are not yet convinced that this method is much more accurate than the posture tracking used in Stephens et al.
(2011) to measure the forward dwell time. The Stephens paper, however, uses a slightly different definition of the dwell time, using cessation of the forward phase velocity instead. Although the paper here does not measure phase velocity, they could measure the dwell time between visits to the reversal state (i.e. ignoring $F$>$P$>$F$ type transitions). What happens if this definition is used? Does the same result emerge?
We have added this calculation to the text. The resulting mean dwell time in the forward state using the standard fixed velocity threshold of 0.5 mm/s and ignoring $\mathrm{FPF}$ transitions was 9.13 ± 0.15 s. which matches the value obtained by others using the same threshold (8.98 ± 0.57 s) (Rakowski, 2013). This result has been added to the text (subsection “Wild type locomotion”, first paragraph).
6) A more general comment is that we would like to understand more what type of neural activity the authors would predict based on their model, connectome data, and polarity results from papers such as Rakowski (2013) despite imaging in a freelybehaving worm being well outside of the scope of this paper (although not nearly as impossible as the authors seem to suggest in the Discussion section). If this experiment were performed, what is the most likely neural instantiation of this model that is consistent with the current literature? The strength of this model is that it makes an attempt at getting at how this circuit may actually function. Accordingly, if possible, a concrete prediction or predictions to this end would greatly increase the value of this paper to researchers and could guide experimentalists performing imaging in freelybehaving animals in their measurements and analysis.
This is a great idea! We have now included a paragraph of predictions in the Discussion (fifth paragraph).
[Editors’ note: the author responses to the previous round of peer review follow.]
Consensus review:
The authors present a kinematic analysis of random local search behavior in C. elegans. The then use a Hidden Markov Model to quantitatively model such searches based their kinematic analyses. They then identify this HMM with known aspects of the worm connectome to understand the neural implementation of switches between different behavioral states. In particular, the model contains two populations of neurons, each controlling either forward or reverse locomotion, with the overall network resulting in four possible behaviors: forward, reverse, and two pause states. They perform a maximum likelihood fit of the model to measured time series, making predictions about the effect of neural ablations, the sign and strength of synaptic connections, effects of perturbations to membrane potentials. Moreover, they phenomenologically model a potential underlying neural mechanism for different search modalities. The subject matter of the paper is definitely one of broad interest, and if the model is correct it would provide us with a substantive understanding into the neural implementation of the forward/pause/reversal/pause(?) dynamics it would be a major contribution. However, the reviewers are not fully convinced that the model is the simplest possible explanation of the animal's behavior and that the
authors have shown that a clear connection can be drawn to neural correlates. Enthusiasm for the subject matter and the potential neural correlates vs. concerns about the appropriateness of the approach and its rigor were weighed differently by the reviewers initially, but in consultation the concerns predominated. The detailed reviews of the expert reviewers are appended but their concerns are summarized below. 1) The model is not adequately placed within the context of previous work in the field. There have been other papers attempting to understand both exploitation/exploration behavioral generation (e.g. Flavell et al., Cell, 2013) and forward/reversing dynamics (e.g. Stephens, PNAS, 2011).
The requested context is now provided in the Introduction.
Moreover, a direct comparison to (Rakowski et al., 2013) is needed. This work, which is cited as Ref. 12, also develops a model of the command neuron network and uses this model to predict behavioral states and the effects of ablations.
The Rakowski paper is difficult to decipher at several key respects, but as we understand the text, reviewers’ points (i) and (ii) are almost certainly overinterpretations of what the Rakowski model (RM) achieves.
(i) Prediction of behavioral states. We agree that the RM enables one to predict/compute from a given weight matrix the quantity $R$ = Tf/(Tf + Tb), where Tf and Tb are, respectively, the total time the animal spends in the forward and backward states during the observation period. Note, however, that $R$ is not a state. Rather, it is the steadystate probability of forward locomotion (conditional on the animal not being in the pause state). The individual behavioral states themselves – forward, backward, and stopped, in their terminology – are not represented in the mathematics. It would therefore be impossible to derive mathematical expressions for mean dwell times in these states, nor could one simulate the network in time to produce a predicted sequence of behavioral states. Thus, it is simply not true that the RM predicts behavioral states.
(ii) Prediction of ablation effects. We agree that the RM computes the effects of ablations on the quantity R, but we don’t agree it predicts them. This is for the simple reason that the model is fitted to the ablation data, so it can’t be said to predict them. It should have been tested on data to which is was not fitted.
With respect to the question of ablation effects, we must call the reviewers' attention to the fact that the ablation data in Rakowski were not analyzed for statistical significance, and the Sternberg lab has confirmed this in personal communications. The absence of statistics means that the question of which purported "effects" are real, and which are false positives of false negatives simply cannot be addressed. Rakowski’s claims related to behavioral effects of ablations should therefore be put on hold until they have been subjected to the field’s common standard of validity.
In revising our manuscript, we nevertheless made the requested direct comparison be between the two models. We added a paragraph in Discussion that is dedicated to the comparison (sixth paragraph). We chose not to mention the missing statistics.
2) The existence of two different pause states $X$ and $Y$ is not well supported in the data and the HMM does not have the power to confirm their existence.
This comment raises two separate questions: (a) Does the addition of a second pause state significantly improve the ability of the model to fit the data? (b) Do the two pause states correspond to “real’” physiological states? These two questions raise fundamentally different issues, which we consider separately below:
a) Does the addition of a second pause state significantly improve the ability of the model to fit the data? Yes, the second pause state improves the fit well beyond what can be attributed to the addition of two free parameters in the model (p<10^{100} by likelihood ratio test (see subsection “Wild type locomotion“, fourth paragraph and Table 1). For this test we constrained the rates into state $Y$ to be ${a}_{\mathrm{FY}}$ = ${a}_{\mathrm{RY}}$ = 10^{10} s^{1}, which effectively eliminates state $Y$ (${p}_{Y}$ < 10^{18}); model B in Table 1). Because the model with one pause state (model B in Table 1) is a special case of the model with two pause states (model A in Table 1), we could apply the likelihood ratio test. The likelihood ratio test is designed specifically for the purpose of comparing two models, one of which is a constrained version of the other, by comparing the reduction in likelihood in the constrained model with the reduction that can be attributed to having fewer free parameters. We further showed that there is almost no benefit in allowing direct transitions between $F$ and $R$ states (model C in Table 1; equivalent to the most general model with one pause state and all 6 transition rates left unconstrained). We conclude that our synaptic model with two pause states, which has 6 free parameters, fits the data far better than any similar model with only one pause state.
b) Do the two pause states correspond to “real’” physiological states? This is subtle question because “real” physiological state is not easily defined. Therefore, for the purpose of this rebuttal we will consider the strongest version of this question, which we infer to be the question that the reviewer is asking: Does the model have the power to confirm the hypothesis that state $X$ corresponds to both pools of command neurons being “OFF”, whereas state $Y$ corresponds to both pools being “ON”? Our answer is no, and we do not make this claim. Instead, we propose this as a working hypothesis that provides a useful framework in which to interpret results but which will need to be tested directly by physiological experiments when and if such experiments become technically feasible. In support of the working hypothesis, it is important to note that it is consistent with experimental observations beyond the wellknown result that the forward command neurons promote forward locomotion and reverse command neurons promote reverse locomotion. It has also been shown that genetic ablation or silencing of both pools of command neurons induces long pauses (subsection “The Stochastic Switch Model, seventh paragraph), as expected for state $X$, and that worms pause when forward and reverse motor neurons simultaneously enter their activated state (in the eighth paragraph of the aforementioned subsection). Although the latter results are for the motor neurons, not the command neurons that drive them, the gap junction synapses known to exist between the command neurons and their motor neuron targets makes it reasonable to hypothesize that both pools of command neurons are also coactivated during these pauses, as expected for state $Y$. Thus, there is sufficient experimental evidence to justify proposing our working hypothesis. Additional evidence that states $X$ and $Y$ are “real” by physiological criteria comes from our results showing that the two states have different properties. State $Y$ is intimately associated with the termination of reverse locomotion, whereas state $X$ is associated with the end of forward locomotion (Figure 3). Because worms typically end reverse locomotion with a ventral bend, state $Y$ is associated with this posture (Figure 4), whereas state $X$ has no strong postural preference. The two states also have different mean dwell times (Figure 2—figure supplement 3), and state $Y$ almost always leads into state $F$, whereas state $X$ can exit to either state $F$ or $R$ (Figure 3). Thus, pauses that follow reverse locomotion are different from pauses that follow forward locomotion, and can therefore be regarded as two different states of the system. It remains to be shown experimentally how closely the states of the model correspond to activity states of the command neurons
We have made three changes to the manuscript which we hope will limit the risk that other readers will overinterpret what we are claiming with respect to the possible association of states of the model with activity states of the command neurons:
A) We have rewritten the paragraphs that explain the model assumptions and the possible relationships between the four states of the model and the activity states of the command neurons (subsection “The Stochastic Switch Model”, third, seventh and eighth paragraphs).
B) The title has been changed from “A stochastic neuronal flipflop circuit regulates random search during foraging behavior,” to “A stochastic neuronal model predicts random search behaviors at multiple spatial scales in C. elegans.”
C) We have restated the primary conclusion of the paper (Discussion, fourth paragraph): “Nevertheless, the predictive and explanatory successes of the model demonstrate the utility of present conceptualization of the network.”
3) The mapping of the HMM onto the underlying neuronal network is not convincing. In particular there little data in this paper or in any of the citations to support the contention that pause state $X$ = all command neurons off and pause state $Y$ = all command neurons on.
As noted in response to critique 2, this is an assumption, not a contention, and we have taken steps to clarify this point in the Discussion.
Mapping synaptic weights onto coupling in the HMM is not rigorous.
We understand this statement to mean that we have not provided an equation for computing synaptic strength, as would be measured electrophysiologically in pairwise recordings quantifying synaptic currents, in terms of synaptic weights (w) as defined in the Stochastic Switch Model. We agree. As requested below, we have relaxed claims to have demonstrated prediction of synaptic weights in the electrophysiological sense of the term.
For the paper to be suitable for eLife the following extensive changes would have to be made.
1) Relax the claims about the model's prediction of synaptic weights.
We have relaxed these claims by pointing out that the term “synaptic weight” has different meanings in the neuralnetwork versus the electrophysiological literature (neither field, of course, can be said to own the term). We now relax our claims by explicitly stating that the Stochastic Switch Model does not predict the magnitude of synaptic conductances (subsection “Synaptic weights in the stochastic switch model”, second paragraph). But it does predict aspects of functional connectivity, such as excitation or inhibition, because weights in the model are signed quantities. We also relax our claims of correspondence between the model and the electrophysiology by noting that our method of testing the model carries with it the assumption that functional connectivity between populations of neurons reflects the connectivity between the neurons in those pools that have the greatest effects on behavior (e.g. in ablations).
2) The authors need to establish firmly that tracking a single point provides an unbiased estimate of trajectory in light of Stephens, PLoS One (2010). Given their extensive data set, we presume that converting their point imaging into centerline tracking would be extremely difficult and time consuming, but the authors should at least place their analysis in this light.
The field recognizes two definitions of trajectory. One is path of the centroid (Pierce Shimomura et al., 1999), defined as the midpoint of the smallest rectangle that can enclose the worm. The other is the path taken by one or more points on the worm’s centerline (Cronin et al., 2005); spot tracking is an instance of the latter.
Centroid trajectory. Stephens et al. are mainly concerned with the time series of postural measurements, not trajectory. The 2010 paper is an exception, but it considers only centroid trajectory. The main point of that paper is to show that the centroid trajectories of real worms can be recovered from the time series of their posture. The authors report, however, that their “eigenworm model” can be expected to do this only when the trajectory of the real worm is relatively straight. The degree of fit between real and computed trajectories is shown in the supplemental material; the fits may or may not be convincing, depending upon one’s point of view. No attempt is made by Stephens et al. to compute centerline trajectories.
Centerline trajectory. According to the theory of thrust generation in sinusoidal locomotion, established by Gray (1946, 1953) and extended to nematodes by (Gray and Lissman, 1964), the spot tracking method should provide an accurate estimate of centerline trajectory, an estimate that is probably more accurate than previous centerline methods (Cronin et al., 2005). In the case of ideal sinusoidal locomotion (no lateral slip), each segment of the body recapitulates the path of the segment ahead of it, like train cars on a sinuous track; thus the trajectory is just the path taken by the segments. On agar plates, C. elegans plows a furrow as it crawls; surface tension provides the necessary downward pressure. The walls of this furrow, which provide the lateral resistance for thrust generation, are what causes each segment to follow the one ahead of it. It is general practice in the field to dry agar plates slightly before use to reduce slip, ensuring that the body stays within the furrow, and locomotion approaches the ideal case. Under these conditions, a radially symmetric tracking perfectly centered on the body axis introduces no bias into trajectory. Bias occurs only when the center of mass of the tracking spot is eccentric to the centerline. This offset shifts the trajectory by the extent of the eccentricity and can increase/decrease the curvature maxima/minima. But in practice eccentricities were small, probably less than 1/10 of the worm’s diameter, because spots with greater eccentricities were quickly rubbed off by friction with the substrate; such recordings were discarded. Another potential problem is deformation of the spot as the worm’s body undulates, but this was almost certainly not an issue in our experiments because the long axis of the spot was small with respect to the radius of curvature of the worm. Note: The first step in Stephens’ image analysis was to compute the worm’s centerline (Stephens et al., 2008), but centerlines were used only to compute posture. They were not used to compute the centerline velocity, although this would have been possible as previously shown (Cronin et al., 2005).
Conversion of spottracking to centerline tracking. Centerlines can be recovered easily from spot tracking data. We have now done this for the entire wild type data set and included it in the paper (Figure 4). Our goal was to ask whether the worm adopts a particular posture as it enters the pause state, as Stephens et al. report. In this analysis we were somewhat limited by the fact that in the spot tracking method, pause posture can be delineated only up to the anteriorposterior level of the spot (approximately the midpoint). This limitation occurs because the spot’s future trajectory after the pause has not yet occurred, but it is mitigated by the fact that postures in which curvature is discontinuous do not occur in live, wild type worms. We find that pauses following reversals occur at a particular posture whereas pauses following forward locomotion do not. In contrast, Figure 5D of Stephens et al. (2008) indicates pauses occur a particular postures regardless of the direction of movement. One obvious reason for the discrepancy between the two studies is the difference in definitions of velocity used to identify the pause state. This topic is addressed in detail in the Discussion (fifth paragraph).
In summary, the eigenworm model is based on postures, not trajectories. Centroid trajectory can be recovered in the eigenworm model, but it is inaccurate unless the worm is moving comparatively straight. In contrast, the spot tracking method provides an accurate measure of centerline trajectory without significant bias.
They should also bolster their arguments that there are two behaviorally distinct pause states and discuss how they are related to the two pause states of Stephens.
(i) On the existence of two pause states. As stated above, we do not claim that the model has the power to confirm the existence of the two pause states. However, the model does provide reasons to believe such states might exist. As requested, we have reinforced these arguments, now devoting a whole new paragraph to the issue (subsection “Wild type locomotion“, fourth paragraph).
(ii) Discuss relationship to pause states of Stephens. Done. See fifth paragraph of the Discussion.
3) The authors should differentiate their work with respect to Rakowski, 2013. Does their model provide different predictions? Can they be disambiguated? It would be great if they could show this, but minimally, they should suggest what experiments should be performed in the future.
See Discussion, sixth paragraph.
4) Similarly –
what measurements should be made to distinguish $X$ and $Y$ pause states map onto the electrical activity of the command neuron network at least at some coarsegrained level?
See Discussion, third and fourth paragraphs.
Reviewer #3: 1) In Figure 4B, I find the excuse for omitting error bars to be relatively weak. Minimally, bootstrapped confidence intervals could be found, which would be significantly better than the current lack of any measure of variability.
Note: Figure 4 is now Figure 5.
Calculating confidence intervals would be nice for graphical representation, but was unnecessary to compute statistical significance of the differences shown in the figure, which was done using the likelihood ratio test, results of which are shown in Table 4. Computing the results shown in Table 4 was extremely timeconsuming; >95% of the computation time for the entire data analysis presented in the paper was for computing the pvalues in this table. Computing confidence intervals for Figure 5B (formerly Figure 4B) would require another ~10x increase in computation time, thousands of hours. Given that the confidence intervals would not be used to for the key statistical tests, which we already have done using the appropriate test, it would be an unnecessary burden to compute them.
2) In the fifth paragraph of the subsection “Genetic effects on command neuron function “and Figure 6. The claim is made that "changes in synaptic strength dominate the effects of changes in membrane potential on Δ$F$ and
Δ$R$." This certainly seems true for Δ$F$, but both models appear to give moreorless identical predictions here for Δ$R$.
We understand the concern here to be that the ΔV and Δ$R$ models make the same prediction for changes in Δ$R$, meaning that the Δ$R$ data cannot be used to distinguish between the two models. However, taking a step back and looking across the Δ$R$ and Δ$F$ data simultaneously, it can be seen that the Δ$R$ model gets both patterns right, whereas the ΔV model gets only one of them right. It is this fact that we take as evidence in support of the Δ$R$ model for Δ$F$ and Δ$R$. The text has been modified to make it more clear that we are looking across Δ$F$and Δ$R$ when comparing the models (subsection “Genetic effects on command neuron function”, sixth paragraph).
3) In the first paragraph of the subsection “The Stochastic Switch Model and other random search behaviors“. It is claimed that "changes are required in at least three weights to produce the full range of search behaviors." This is not actually shown anywhere in the paper, just that three weights are sufficient to induce the full range. A supplementary figure to this effect would be useful.
The revised text and figures now carefully document weight changes that are minimally required for the full range of search modes (Figure 8 including supplementals, Table 7 (new), and accompanying text (subsection “Regulation of search scale, last two paragraphs and subsection “Biased random walks”).
4) The predicted value for the dwell time in the forward state is not comparable to the rest of the literature, as a different definition is used here. What would happen if a standard definition was used? Does the same rate result?
We now report the mean dwell times for $F$, $R$ and $X$ based on fixed velocity thresholds of ± 0.05 mm/s. = 1.86 ± 0.03; = 1.23 ± 0.02; ${d}_{P,0.05}$= 0.14 ± 0.001. Thus, the short dwell times in the forward state that we observed were not due to the model. We attribute them instead to the higher time resolution of our tracking hardware (subsection “Wild type locomotion”, first paragraph).
https://doi.org/10.7554/eLife.12572.033Article and author information
Author details
Funding
National Institutes of Health (RO1 MH051383)
 Steven B Augustine
 Kristy J Lawton
 Theodore H Lindsay
 Tod R Thiele
 Serge Faumont
 Rebecca A Lindsay
 Matthew Cale Britton
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
This research was supported by a grant MH051383 from the National Institute of Mental Health to SRL.
Reviewing Editor
 Ronald L Calabrese, Emory University, United States
Publication history
 Received: October 26, 2015
 Accepted: January 19, 2016
 Accepted Manuscript published: January 29, 2016 (version 1)
 Version of Record published: March 8, 2016 (version 2)
 Version of Record updated: October 11, 2018 (version 3)
Copyright
© 2016, Roberts et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 5,786
 Page views

 1,285
 Downloads

 60
 Citations
Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading

 Computational and Systems Biology
Biological age, distinct from an individual's chronological age, has been studied extensively through predictive aging clocks. However, these clocks have limited accuracy in short timescales. Here we trained deep learning models on fundus images from the EyePACS dataset to predict individuals' chronological age. Our retinal aging clocking, 'eyeAge', predicted chronological age more accurately than other aging clocks (mean absolute error of 2.86 and 3.30 years on qualityfiltered data from EyePACS and UK Biobank, respectively). Additionally, eyeAge was independent of blood markerbased measures of biological age, maintaining an allcause mortality hazard ratio of 1.026 even when adjusted for phenotypic age. The individualspecific nature of eyeAge was reinforced via multiple GWAS hits in the UK Biobank cohort. The top GWAS locus was further validated via knockdown of the fly homolog, Alk, which slowed agerelated decline in vision in flies. This study demonstrates the potential utility of a retinal aging clock for studying aging and agerelated diseases and quantitatively measuring aging on very short timescales, opening avenues for quick and actionable evaluation of geroprotective therapeutics.

 Computational and Systems Biology
Computational models starting from large ensembles of evolutionarily related protein sequences capture a representation of protein families and learn constraints associated to protein structure and function. They thus open the possibility for generating novel sequences belonging to protein families. Protein language models trained on multiple sequence alignments, such as MSA Transformer, are highly attractive candidates to this end. We propose and test an iterative method that directly employs the masked language modeling objective to generate sequences using MSA Transformer. We demonstrate that the resulting sequences score as well as natural sequences, for homology, coevolution, and structurebased measures. For large protein families, our synthetic sequences have similar or better properties compared to sequences generated by Potts models, including experimentally validated ones. Moreover, for small protein families, our generation method based on MSA Transformer outperforms Potts models. Our method also more accurately reproduces the higherorder statistics and the distribution of sequences in sequence space of natural data than Potts models. MSA Transformer is thus a strong candidate for protein sequence generation and protein design.