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 (${\displaystyle \overline{\left|\mathrm{\Delta }\phi \right|}}$), 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, ${\displaystyle V(t)}$ , 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 ${\displaystyle V(t)}$ (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 ; Mann-Whitney U-test). 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.

Video 1

Download asset

Video 2

Download asset

Video 3

Download asset

Video 4

Download asset

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 neuron-like “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 neuron-by-neuron 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 input-output function of individual command neurons.

Figure 2 with 3 supplements see all

Download asset Open asset

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 (Pierce-Shimomura 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 co-active, 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; Pokala et al., 2014). Additionally, neurons in opposing groups are likely to be reciprocally active, as indicated by simultaneous calcium imaging from AVA and AVB (Pokala 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{ℱ}$ and $\mathcal{ℛ}$, and six synaptic weights (Figure 2C). Each type of weight has a specific interpretation. The cross-connections (${\displaystyle {\mathit{W}}_{\mathcal{ℱ}\mathcal{ℛ}}}$, ${\displaystyle {\mathit{W}}_{\mathcal{ℛ}\mathcal{ℱ}}}$) represent mono- and polysynaptic connections between command neurons in different groups. The self-connections (${\displaystyle {\mathit{W}}_{\mathcal{ℱ}\mathcal{ℱ}}}$, ${\displaystyle {\mathit{W}}_{\mathcal{ℛ}\mathcal{ℛ}}}$) represent connections between command neurons in the same group, including recurrent polysynaptic pathways involving neurons outside the command network. Self-connections 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{ℱ}$or $\mathcal{ℛ}$ 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, ${\displaystyle h_{\mathcal{ℱ}}}$ and ${\displaystyle h_{\mathcal{ℛ}}}$, 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{ℱ}$ and $\mathcal{ℛ}$ are, respectively, ${\displaystyle S_{\mathcal{ℱ}}(t)=h_{\mathcal{ℱ}}+{\mathit{W}}_{\mathcal{ℱ}\mathcal{ℱ}}{b}_{\mathcal{ℱ}}(t)+{\mathit{W}}_{\mathcal{ℛ}\mathcal{ℱ}}{b}_{\mathcal{ℛ}}(t)h_{\mathcal{ℱ}}}$ and ${\displaystyle S_{\mathcal{ℛ}}(t)=h_{\mathcal{ℛ}}+{\mathit{W}}_{\mathcal{ℛ}\mathcal{ℛ}}{b}_{\mathcal{ℛ}}(t)+{\mathit{W}}_{\mathcal{ℱ}\mathcal{ℛ}}{b}_{\mathcal{ℱ}}(t)}$, where ${\displaystyle b_{\mathcal{ℱ}}(t)}$ and ${\displaystyle b_{\mathcal{ℛ}}(t)}$ are the states of $\mathcal{ℱ}$ and $\mathcal{ℛ}$ at time $t$ (1 = ON, 0 = OFF). The quantities ${\displaystyle h_{\mathcal{ℱ}}}$ and ${\displaystyle h_{\mathcal{ℛ}}}$ were assumed to be constant during the 10 min observation period of local search behavior on a bare agar surface.

State transitions of $\mathcal{ℱ}$ and $\mathcal{ℛ}$ were modeled as independent non-homogeneous 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{ℱ}$ and $\mathcal{ℛ}$ can be regarded as thermally-driven transitions across energy barriers of height proportional to ${\displaystyle S_{\mathcal{ℱ}}(t)}$ and ${\displaystyle S_{\mathcal{ℛ}}(t)}$, 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{ℱ}$ and $\mathcal{ℛ}$ 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; Viswanathan, 2011), 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{ℱ}$ and $\mathcal{ℛ}$, are intended to represent the two pools of forward and reverse command neurons, respectively, such that the worm moves forward when $\mathcal{ℱ}$ is ON and $\mathcal{ℛ}$ is OFF (state F), backwards when $\mathcal{ℛ}$ is ON and $\mathcal{ℱ}$ is OFF (state R), and pauses when both $\mathcal{ℱ}$ and $\mathcal{ℛ}$ 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{ℱ}$ and $\mathcal{ℛ}$ are simultaneously ON (state Y). Whether the corresponding co-activation 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 co-active when their motor neurons are co-active. 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{ℱ}$ and $\mathcal{ℛ}$ 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 non-simultaneous 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.

The maximum likelihood rate constants for 5 wild-type 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 (${\displaystyle d_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 (${\displaystyle d_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: ${\displaystyle d_{F,0.005}=}$ 1.86 ± 0.03 s; ${\displaystyle d_{R,0.005}=}$ 1.23 ± 0.02 s; ${\displaystyle d_{P,0.005}=}$ 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 (${\displaystyle d_X=}$0.44 ± 0.03 s, ${\displaystyle d_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 short-lived.

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 model-independent 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 (${\displaystyle 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 one-pause 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 one-state 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 (${\displaystyle d_X>d_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).

Figure 3

Download asset Open asset

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).

Figure 4

Download asset Open asset

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).

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 and Kristan, 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).

Figure 5

Download asset Open asset

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 co-active 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: ${\displaystyle h_{\mathcal{ℛ}}}$, ${\displaystyle {\mathit{W}}_{\mathcal{ℛ}\mathcal{ℛ}}}$, ${\displaystyle {\mathit{W}}_{\mathcal{ℛ}\mathcal{ℱ}}}$, and ${\displaystyle {\mathit{W}}_{\mathcal{ℱ}\mathcal{ℛ}}}$. 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 ${\displaystyle d_0}$ to denote the uncoupled dwell time. Dwell times that in intact animals are greater than ${\displaystyle d_0}$ will be reduced by ablation, whereas dwell times that are less than ${\displaystyle d_0}$ will be increased. In particular, if ${\displaystyle d_F}$ and ${\displaystyle d_R}$ are both greater than ${\displaystyle 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 ${\displaystyle d_0}$ is below ${\displaystyle d_F}$ and ${\displaystyle d_R}$.

To determine the actual relationship between ${\displaystyle 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 ${\displaystyle d_0}$ is indeed below ${\displaystyle d_F}$ and ${\displaystyle d_R}$. Additionally, dwell times in the pause states ${\displaystyle d_X}$ and ${\displaystyle d_Y}$ were increased, with one exception (${\displaystyle d_Y}$, AVB). Thus, the observed pattern of dwell time changes is consistent, overall, with a value of ${\displaystyle 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 ${\displaystyle d_0}$. Specifically, ${\displaystyle d_0}$ must be less than or equal to the lowest post-ablation value of ${\displaystyle d_R}$, and greater than or equal to the largest post-ablation value of ${\displaystyle d_X}$; thus, 0.58 ≤ ${\displaystyle d_0}$ ≤ 1.24 sec. Furthermore, because ${\displaystyle A=1/2d_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 ${\displaystyle A_{min}=}$ 0.40 Hz and ${\displaystyle A_{max}=}$ 0.86 Hz.

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 ${\displaystyle A_{min}}$ and ${\displaystyle A_{max}}$. We found that input weights, ${\displaystyle h_{\mathcal{ℱ}}}$ and ${\displaystyle h_{\mathcal{ℛ}}}$, are small and positive, suggesting that these inputs may provide modest but steady excitation to the system (Figure 6A). The self-connections ${\displaystyle {\mathit{W}}_{\mathcal{ℱ}\mathcal{ℱ}}}$ and ${\displaystyle {\mathit{W}}_{\mathcal{ℛ}\mathcal{ℛ}}}$ are also mainly positive, indicating that the ON states may be stabilized by intrinsic or extrinsic positive feedback. The cross-connections ${\displaystyle {\mathit{W}}_{\mathcal{ℱ}\mathcal{ℛ}}}$ and ${\displaystyle {\mathit{W}}_{\mathcal{ℛ}\mathcal{ℱ}}}$ are negative, indicating reciprocal inhibition, as expected for neurons that activate opposing behavioral states. Furthermore, the magnitude of ${\displaystyle {\mathit{W}}_{\mathcal{ℱ}\mathcal{ℛ}}}$ is greater than the magnitude of ${\displaystyle {\mathit{W}}_{\mathcal{ℛ}\mathcal{ℱ}}}$, 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.

Figure 6

Download asset Open asset

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 channelRhodopsin-2 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 zero-current 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 thaCroninn 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, ${\displaystyle v(t)}$, 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.

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 (eat-4(ad572), eat-4(ky5)) or postsynaptic mechanisms (glr-1(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 (glr-1::glr-1(A/T), nmr-1::glr-1(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 wild-type 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 ${\displaystyle \mathrm{\Delta }h}$ (${\displaystyle -1\le \mathrm{\Delta }h\le }$ 1) to wild type ${\displaystyle h}$ values, with negative ${\displaystyle \mathrm{\Delta }h}$ for HYP mutations and positive ${\displaystyle \mathrm{\Delta }h}$ for DEP mutations. We modeled the effect of Δr as a change in the magnitude of synaptic weights (${\displaystyle 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 ${\displaystyle d_X}$, which depends only on ${\displaystyle h_{\mathcal{ℱ}}}$ and ${\displaystyle h_{\mathcal{ℛ}}}$. Equation 1 shows that ${\displaystyle d_X}$ rises and falls as ${\displaystyle h}$ terms are made more negative or positive, respectively. Thus, under the ΔV hypothesis, ${\displaystyle d_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 (${\displaystyle \left|\mathit{W}\right|}$ and ${\displaystyle \left|h\right|}$) decrease in DEP mutants and increase in HYP mutants, ${\displaystyle d_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 ${\displaystyle d_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 ${\displaystyle d_X}$.

Figure 7

Download asset Open asset

In contrast to ${\displaystyle d_X}$, ${\displaystyle d_F}$ and ${\displaystyle d_R}$ depend on $w$ terms as well as ${\displaystyle h}$ terms. Under the ΔV hypothesis, the ${\displaystyle h}$ terms but not the $w$ terms would be affected by the mutations. Positive and negative increments in ${\displaystyle h}$ have the effects shown in Figure 7A, rows 1 and 2; ${\displaystyle d_F}$ and ${\displaystyle d_R}$ are predicted to shift in opposite directions. Changes in ${\displaystyle d_F}$ are dominated by the effects of ${\displaystyle h_{\mathcal{ℱ}}}$ on the first term in Equation 2 (which represents ${\displaystyle a_{FX}}$) because the second term in the equation (which represents ${\displaystyle a_{FY}}$) remains close to zero in the mutants. Analogously, changes in ${\displaystyle d_R}$ are dominated by the effects of ${\displaystyle h_{\mathcal{ℱ}}}$ on the second term in Equation 3 (${\displaystyle a_{RY}}$) because the first term in the equation (${\displaystyle a_{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 ${\displaystyle h}$ terms would both be affected by the mutations. Now, the predicted pattern of dwell time changes across both ${\displaystyle d_F}$ and ${\displaystyle d_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 ${\displaystyle d_F}$ and ${\displaystyle d_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 ${\displaystyle d_F}$ and ${\displaystyle d_R}$.

Neither hypothesis predicts the observed changes in ${\displaystyle d_Y}$ (Figure 7C, row 5) which resembled the pattern of changes in ${\displaystyle d_X}$ (Figure 7C, row 4). However, the ΔV hypothesis correctly predicts observed dwell times in the overall pause state ${\displaystyle d_p}$ (Figure 7C row 3). This is because ${\displaystyle d_p}$ is dominated by ${\displaystyle d_X}$ and changes in ${\displaystyle d_X}$ are well-predicted 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.

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 ${\displaystyle m_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, ${\displaystyle m_F={\overline{V}}_F{p}_F/f_{RPF}}$, where ${\displaystyle {\overline{V}}_F}$ is the average velocity in state F, ${\displaystyle p_F}$ is the probability of being in state F, and ${\displaystyle f_{RPF}}$ is the frequency of RPF transitions (Materials and methods, Equation 39), which coincide with random reorientations. Importantly, under the approximation ${\displaystyle a_{FY}\cong 0}$ (Table 1, column A), ${\displaystyle m_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, ${\displaystyle m_F}$ could be controlled by sensory inputs represented by ${\displaystyle h_{\mathcal{ℱ}}}$ and ${\displaystyle h_{\mathcal{ℛ}}}$ (Figure 8A).

Figure 8 with 2 supplements see all

Download asset Open asset

Search scale (${\displaystyle m_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 six-dimensional 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 ${\displaystyle m_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 ${\displaystyle m_F}$ was covered in the scan (Figure. 8B-H, gray symbols).

All three search modes were available in the subspace defined by the control points (${\displaystyle h_{\mathcal{ℱ}}}$, ${\displaystyle h_{\mathcal{ℛ}}}$, ${\displaystyle {\mathit{W}}_{\mathcal{ℱ}\mathcal{ℱ}}}$) (Figure 8B, Figure 8—figure supplement 1). However, only cropping and local search were available in the complementary subspace (${\displaystyle {\mathit{W}}_{\mathcal{ℛ}\mathcal{ℛ}}}$, ${\displaystyle {\mathit{W}}_{\mathcal{ℛ}\mathcal{ℱ}}}$, ${\displaystyle {\mathit{W}}_{\mathcal{ℱ}\mathcal{ℛ}}}$) (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 control-point weights was sufficient on its own to produce all three search modes (Figure 8D–F). Scanning the subspaces (${\displaystyle h_{\mathcal{ℱ}}}$, ${\displaystyle {\mathit{W}}_{\mathcal{ℱ}\mathcal{ℱ}}}$) and (${\displaystyle h_{\mathcal{ℛ}}}$, ${\displaystyle {\mathit{W}}_{\mathcal{ℱ}\mathcal{ℛ}}}$) showed these pairs of weights to be sufficient for all modes (Figure 8G,H), but a three-dimensional subspace containing at least one of the control-point 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 three-weight subspaces constitute the most likely minimal models for the regulation of search in C. elegans. They could be tested by chronic manipulation of control-point 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 ${\displaystyle h_{\mathcal{ℱ}}}$ and ${\displaystyle h_{\mathcal{ℛ}}}$.

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 (Pierce-Shimomura 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 random-walk chemotaxis, ON cells increase ${\displaystyle h_{\mathcal{ℱ}}}$ and decrease ${\displaystyle h_{\mathcal{ℛ}}}$, 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).

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 ${\displaystyle p_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 ${\displaystyle h_{\mathcal{ℛ}}}$ thereby increasing the rate constants for transitions in which ${\displaystyle \mathcal{ℛ}}$ turns ON (${\displaystyle a_{XR}}$ and ${\displaystyle a_{FY}}$), and decreasing the rate constants for transitions in which ${\displaystyle \mathcal{ℛ}}$ turns OFF (${\displaystyle a_{XR}}$ and ${\displaystyle a_{FY}}$ ). In the limit where $h_{\mathcal{ℛ}}$→ ∞, both ${\displaystyle a_{XR}}$ and ${\displaystyle a_{FY}}$→ ∞, whereas ${\displaystyle a_{RX}}$ and ${\displaystyle a_{YF}}$→ 0 (Figure 9B). The system now inhabits only states R and Y, and ${\displaystyle p_R=a_{YR}/(a_{YR}+a_{RY})}$. In the Pull motif, nociceptive neurons inhibit the forward command neurons via $h_{\mathcal{ℱ}}$. In the limit where $h_{\mathcal{ℱ}}$→ -∞, the system switches only between states R and X and ${\displaystyle p_R=a_{XR}/(a_{XR}+a_{RX})}$. In the third motif, in which Push and Pull are combined, R becomes an absorbing state ${\displaystyle (p_R=1)}$. Using the rate constants shown in column A of Table 1 to compute limiting values of ${\displaystyle p_R}$ in each motif, we found that the Pull and Push-Pull 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.

Figure 9

Download asset Open asset

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, ASH-mediated escape may be controlled by a Push-Pull motif, further demonstrating the feasibility of using behavioral data to predict population-level 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 Push-Pull motifs are equally effective in driving ${\displaystyle p_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.3-fold reduction in latency for the Push-Pull motif (Figure 9B, parenthetical values). We conclude that the ASH mediated escape circuit in C. elegans may be specialized for short latency escape responses.

All strains were cultivated at 22.5°C on low-density 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).

External saline for electrophysiology (mM): 5 KCl, 10 HEPES, 8 CaCl2, 143 NaCl, 30 glucose, pH 7.2 (NaOH); internal saline for electrophysiology (mM): 125K-gluconate, 1 CaCl2, 18 KCl, 4 NaCl, 1 MgCl2, 10 HEPES, 10 EGTA, pH 7.2 (KOH). Medium for behavioral assays (mM): NH₄Cl 2, CaCl₂ 1, MgSO₄ 1, and KPO₄ 25, pH 6.5; M9 Buffer (grams): 3 KH₂PO₄, 6 Na₂HPO₄, 5 NaCl, 1 ml 1 M MgSO₄, H₂O to 1 liter.

Prior to each assay, an individual adult hermaphrodite was picked to a bacteria-free agar transfer plate by means of a platinum-wire 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 MS-2000, Eugene, OR USA) fitted with position encoders (Gurely Precision Instruments LE-1800, Troy, NY USA) having a resolution of 0.5 μm. Behavior was recorded using an analog video camera (CCD Sony XC-ST70, 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 re-center 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 off-line 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 spot-checked 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.

The animal was immobilized by a stream of humidified CO₂ 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.

Worms 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 current-clamp recordings were made with a modified Axopatch 200A amplifier (Lockery and Goodman, 1998). In reversal potential measurements, recordings of photostimulation-evoked 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 ChannelRhodopsin-2 expressed under the control of neuron-specific promoters as described (see “Strains”). Worms were photostimulated in electrophysiological experiments using the blue channel (470 nm) of a dual-wavelength 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²) 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.

Neurons 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₃. AVA and AVB neurons were ablated in N2 animals and identified by position. AVD, AVE and PVC were ablated in animals expressing nmr-1::GFP and identified by a combination of position and GFP expression. To limit potential behavioral differences in the two strains, we outcrossed (4×) the nmr-1::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. Sham-operated animals were treated in the same manner except that the laser was not fired.

Statistical significance for the results shown in Figures 4B and 6C, and in Tables 4 and 6 were obtained using the likelihood ratio test (see Table 1 and 4 legends). Otherwise, two-tailed t-tests or 2-tailed Mann-Whitney U tests were used.

The worm’s position in video frame $k$ is represented as the row vector:

where ${\displaystyle X(t_k)}$ and ${\displaystyle Y(t_k)}$ are the coordinates of centroid of the tracking spot in the frame of reference of the agar plate, ${\displaystyle t_k=k\mathrm{\Delta }t}$, $∆t=$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:

To analyze locomotory states we converted the velocity vector, $\mathit{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 $\mathit{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 half-width $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 ${\displaystyle g_F(v)}$ and ${\displaystyle g_R(v)}$ by scaling ${\displaystyle g_P(v)}$ to fit the peak at $v=0$, subtracting it from the overall distribution and splitting the remaining distribution into ${\displaystyle g_F(v)}$ for $v>0$ and ${\displaystyle g_R(v)}$ for $v<0$. Velocity distributions were scaled to be probability densities (area =1) and collected into a row vector:

where ${\displaystyle g_i(v)}$ 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 {${\displaystyle a_{XF},a_{FX},a_{XR},a_{RX},a_{FY},a_{YF},a_{RY},a_{YR}}$}that maximize the probability of the observed velocity time series $v\left(t\right)$ given the velocity distribution $\mathit{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$.$ 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 ${\displaystyle q_{ij}}$ ($i\ne j)$ of matrix $\mathit{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 ${\displaystyle 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 $\mathit{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 ·\mathit{Q}\mathbf{+}\mathit{I}$, where $\mathit{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 ·\mathit{Q}+\mathit{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 ·\mathit{Q}+\mathit{I}$ by itself. Thus, if

then $\mathit{M}$ is the matrix of transition probabilities during a time interval of duration $K\epsilon$. If $K$ and $\epsilon$ are chosen such that $∆t=K\epsilon$, then element $ij$ of matrix $\mathit{M}$ is the transition probability from state $i$ to state $j$ during one video frame of duration $∆t$. We chose $K=2^{30}$ and let $\epsilon =∆t/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, $\mathit{M}$ could be rapidly and accurately calculated by 30 serial multiplications using 64-bit floating point arithmetic.

Let $\mathit{P}\left(t\right)$ be the row vector of history-dependent 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 $\mathit{P}\left(t\right)·\mathit{M}$ is the state probability vector at time $t+∆t$ prior to accounting for the observed velocity at time $t+∆t$. To account for $v(t+∆t)$ we used the information contained in $\mathit{G}\left(v(t+∆t)\right)$ and applied Bayes theorem:

where $diag\left[G\left(v(t+∆t)\right)\right]$ is the 4×4 matrix with the elements of $\mathit{G}\left(v(t+∆t)\right)$ along the diagonal, and $l$ is the scalar multiplicative factor required for the sum of the four elements of $\mathit{P}(t+∆t)$ to equal 1 (i.e., $\mathit{P}(t+∆t)$ is a vector of probabilities). Initially ($t=0$) we set $\mathit{P}\left(0\right)$ equal to the steady-state probability vector ${\displaystyle {\mathit{P}}_{\mathrm{\infty }}}$, which is given by:

where ${\mathit{U}}_4$ is the $1\times 4$ row vector of ones and ${\mathit{Q}}_a$ is the $4\times 5$ matrix constructed by appending a column of ones to $\mathit{Q}$. To break the symmetry between the behaviorally indistinguishable states X and Y, we identified X as the state with higher steady-state probability.

We then calculated the log-likelihood, summed over all worms in the cohort:

where ${\displaystyle V_W(t)}$ is the velocity and ${\displaystyle {\mathit{P}}_W(t)}$ is the history-dependent 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 ${\displaystyle b_{random}}$ (initially ${\displaystyle 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 ${\displaystyle 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 ${\displaystyle 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 ${\displaystyle 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 $\mathit{P}\left(t\right)$, but once the optimal transition rates were determined, the Forward-Backward algorithm (Rabiner and Juang, 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).

We expressed the effect of synaptic inputs to command units $\mathcal{ℱ}$ and $\mathcal{ℛ}$ by equations of the form:

where ${\displaystyle a_{OFF}}$ is the transition rate from ON to OFF, ${\displaystyle a_{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, ${\displaystyle 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 ${\displaystyle K_{B}T}$. Thus, ${\displaystyle \mathcal{ℱ}}$ and ${\displaystyle \mathcal{ℛ}}$ are assumed to be symmetrical bi-stable units that change state at rate $A$ when $S=0$. Deviations from this baseline condition are modelled as external synaptic inputs ${\displaystyle h_{\mathcal{ℱ}}}$ and ${\displaystyle h_{\mathcal{ℛ}}}$.

We represented the total synaptic input onto units ${\displaystyle \mathcal{ℱ}}$ and $\mathcal{ℛ}$, respectively, by:

where ${\displaystyle b_{\mathcal{ℱ}}(t)}$ and ${\displaystyle b_{\mathcal{ℛ}}(t)}$ are the states of ${\displaystyle \mathcal{ℱ}}$ and ${\displaystyle \mathcal{ℛ}}$ (1 = ON, 0 = OFF), ${\displaystyle {\mathit{W}}_{\mathcal{ℛ}\mathcal{ℱ}}}$ and ${\displaystyle {\mathit{W}}_{\mathcal{ℱ}\mathcal{ℛ}}}$ are the synaptic weights from ${\displaystyle \mathcal{ℛ}}$ onto ${\displaystyle \mathcal{ℱ}}$ and from ${\displaystyle \mathcal{ℱ}}$ onto ${\displaystyle \mathcal{ℛ}}$, respectively, and ${\displaystyle {\mathit{W}}_{\mathcal{ℱ}\mathcal{ℱ}}}$ and ${\displaystyle {\mathit{W}}_{\mathcal{ℛ}\mathcal{ℛ}}}$ 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 ${\displaystyle h_{\mathcal{ℱ}}}$ and ${\displaystyle h_{\mathcal{ℛ}}}$ were assumed to be constant. For simulations of chemotaxis ${\displaystyle h_{\mathcal{ℱ}}}$ and ${\displaystyle h_{\mathcal{ℛ}}}$ 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: