Shared behavioral mechanisms underlie C. elegans aggregation and swarming.

In complex biological systems, simple individual-level behavioral rules can give rise to emergent group-level behavior. While collective behavior has been well studied in cells and larger organisms, the mesoscopic scale is less understood, as it is unclear which sensory inputs and physical processes matter a priori. Here, we investigate collective feeding in the roundworm C. elegans at this intermediate scale, using quantitative phenotyping and agent-based modeling to identify behavioral rules underlying both aggregation and swarming-a dynamic phenotype only observed at longer timescales. Using fluorescence multi-worm tracking, we quantify aggregation in terms of individual dynamics and population-level statistics. Then we use agent-based simulations and approximate Bayesian inference to identify three key behavioral rules for aggregation: cluster-edge reversals, a density-dependent switch between crawling speeds, and taxis towards neighboring worms. Our simulations suggest that swarming is simply driven by local food depletion but otherwise employs the same behavioral mechanisms as the initial aggregation.

Introduction

Collective behavior has been widely studied in living and non-living systems. While very different in their details, shared principles have begun to emerge, such as the importance of alignment for flocking behavior in both theoretical models and birds (Bialek et al., 2012; Pearce et al., 2014; Reynolds, 1987). Until now, the study of collective behavior has mainly focused on cells and active particles at the microscale, controlled by molecule diffusion and direct contact between cells or particles (Köhler et al., 2011; De Palo et al., 2017; Peruani et al., 2012; Starruß et al., 2012), and on animals at the macroscale, aided by long-range visual cues (Bialek et al., 2012; Katz et al., 2011; Pearce et al., 2014). Collective behavior at the intermediate mesoscale is less well-studied, as it is unclear what processes to include a priori. At the mesoscale, sensory cues and motility may still be limited by the physics of diffusion and low Reynolds numbers, respectively, yet the inclusion of nervous systems allows for increased signal processing and a greater behavioral repertoire. Do the rules governing collective behavior at this intermediate scale resemble those at the micro- or the macro-scale, some mixture of both, or are new principles required?

C. elegans collective behavior can contribute to bridging this scale gap. Some strains of this 1 mm-long roundworm are known to aggregate into groups on food (de Bono and Bargmann, 1998); here we also report an additional dynamic swarming phenotype that occurs over longer time periods. C. elegans represents an intermediate scale not only in physical size but also in behavioral complexity—crawling with negligible inertia, limited to touch and chemical sensing, yet possessing a compact nervous system with 302 neurons (White et al., 1986) that supports a complex behavioral repertoire (Hart, 2006; Schwarz et al., 2015). Wild C. elegans form clusters on food at ambient oxygen concentrations, as do loss-of-function neuropeptide receptor 1 (npr-1) mutants. The laboratory reference strain N2, on the other hand, has a gain-of-function mutation in the npr-1 gene that suppresses aggregation (de Bono and Bargmann, 1998), rendering N2 animals solitary feeders. Thus, a small genetic difference (just two base pairs in one gene for the npr-1(ad609lf) mutant) has a big effect on the population-level behavioral phenotype. Previous research on collective feeding has focused primarily on the genetics and neural circuits that govern aggregation (Bretscher et al., 2008; Busch et al., 2012; Chang et al., 2006; Cheung et al., 2005; de Bono et al., 2002; de Bono and Bargmann, 1998; Gray et al., 2004; Jang et al., 2017; Macosko et al., 2009), rather than on a detailed understanding of the behavior itself. Rogers et al. (2006) is a notable exception and includes an investigation of the behavioral motifs that might lead to cluster formation including direction reversals at the edge of clusters. However, we do not know whether these candidate motifs are sufficient to produce aggregation. We also do not know whether aggregation at short times and swarming at longer times are distinct behaviors or different emergent properties of the same underlying phenomenon.

In this paper, we use fluorescence imaging and multi-worm tracking to examine individual behavior inside aggregates. We present new and systematic quantification of the aggregation behavior in hyper-social npr-1(ad609lf) mutants (henceforth referred to as npr-1 mutants) and hypo-social N2 worms. Next, we draw on the concept of motility-induced phase transitions to explain aggregation as an emergent phenomenon by modulating only a few biophysical parameters. Unlike aggregation driven by attractive forces, in motility-induced phase transitions individuals can also aggregate simply due to their active movement and non-attractive interactions, such as volume exclusion (avoidance of direct overlap) (Redner et al., 2013a). For instance, this concept has contributed understanding to the aggregation of rod-shaped Myxococcus xanthus bacteria, which, similar to C. elegans, also exhibit reversals during aggregation (Mercier and Mignot, 2016; Peruani et al., 2012; Starruß et al., 2012). We build an agent-based phenomenological model of simplified worm motility and interactions. By mapping out a phase diagram of behavioral phenotypes, we show that modulating cluster-edge reversals and a density-dependent switch between crawling speeds are sufficient to produce some aggregation, but not the compact clusters observed in experiments. We found that medium-range taxis towards neighboring worms is necessary to tighten clusters and increase persistence. Finally, combining this model with food depletion gives rise to swarming over time, suggesting that the same behavioral rules that lead to the initial formation of aggregates also underlie the dynamic swarming reported here.

Results

Dynamic swarming occurs in social worms at long time scales

Aggregation has most often been characterized as the fraction of worms inside clusters, where individual worms can move in and out of clusters. Here we report an additional dynamic swarming phenotype in aggregating C. elegans that occurs on a timescale of hours. Here, swarming refers to the collective movement of a coherent group of worms across a bacterial lawn (Figure 1A, Video 1). Because of the long timescale, this behavior is not obvious from manual observations of worms on a plate, but becomes clear in time lapse videos (Figure 1B and C, npr-1 panels). Even though N2 worms do not swarm in our experiments (Figure 1B and C, N2 panels), they can swarm under appropriate conditions, such as when a clonal population has depleted almost all food (Hodgkin and Barnes, 1991) or on unpalatable Pseudomonas fluorescens bacterial lawns (personal communication from J. Hodgkin and G.M. Preston). Thus swarming in C. elegans does not require loss of npr-1 function in all environments.

Figure 1.

npr-1 but not N2 worms show swarming behavior over time on thin bacterial lawn.

(A) A few hundred npr-1 mutant worms form dense clusters that move on food over time. Red dashed lines show the food boundary, where area with food is to the right and food-depleted area is to the left; red arrows show the direction of cluster movement. (B) Forty npr-1 mutant worms also cluster and swarm on food. Solid circles encompass the same cluster at different time points; dashed circles show cluster positions prior to the current time point. The same number of N2 worms do not swarm under our experimental conditions, and instead disperse after initial transient aggregation. (C) Visualization of persistent swarming over time. One frame was sampled every 30 s over the duration of the videos and binary segmentation was applied using an intensity threshold to separate worm pixels from the background. Blobs with areas above a threshold value were plotted as clusters to show cluster position over time. The same videos as in (B) were used. Dashed circles show the food boundary. Crosses are cluster centroids at each sample frame. (D) Centroid speed of persistent npr-1 clusters, calculated from centroid positions as indicated in (C) and smoothed over 10 min. Shaded area shows standard deviation across five replicates.

Forty npr-1 animals swarm over a big food patch, reminiscent of a persistent random walk. More than one large moving clusters co-exist towards the end of the video (orange and yellow), and a cluster (orange) disperses and re-forms elsewhere (orange and yellow) when it crosses its previous path (blue), presumably due to local food depletion.

One-hour fluorescence recordings of npr-1 animals under our experimental conditions consist of reproducible temporal dynamics encompassing three phases: transient (animals move about the lawn and start to form clusters), aggregation (clusters largely remain stable with individuals entering and exiting), and swarming (worms move across the lawn in persistent clusters) (see sample Video 2). Percentage of in-cluster worms remain largely consistent throughout the latter two phases, except that clusters remain in place during the aggregation phase and become dynamic during the swarming phase. Error bars represent standard deviation across 13 (npr-1) and 9 (N2) replicates. The average duration of each phase derived from npr-1 experiments are applied to N2 data to maintain temporal consistency, even though N2 does not exhibit aggregation or swarming. Subsequent quantitative analyses for both strains were restricted to using the data from the aggregation (for all but Figure 1D) or swarming (for Figure 1D) phase, in order to reveal the mechanisms necessary for producing aggregation or the dynamics of swarming, respectively.

Video 1.

Sample video showing npr-1 collective feeding dynamics (bright field high-number swarming imaging).

The video plays at 300x the normal speed.

Video 2.

Sample video showing npr-1 collective feeding dynamics (fluorescence 40 worm aggregation imaging).

The video plays at 90x the normal speed.

Dynamic swarming occurs with just 40 npr-1 mutants (Figure 1B, top row), making it experimentally feasible to study. Usually a single npr-1 aggregate forms on the food patch and then moves around the lawn in a persistent but not necessarily directed manner (Figure 1C, left; Figure 1—figure supplement 1), at a steady speed (Figure 1D). The onset of this collective movement appears to coincide with local food depletion, and continues until complete food depletion, at which time the cluster disperses. More than one moving cluster may co-exist, and occasionally a cluster may disperse and form elsewhere when it crosses its previous path (Figure 1—figure supplement 1), presumably due to local food depletion. The observed pattern of npr-1 cluster motion is reminiscent of a self-avoiding, persistent random walk (i.e. not returning to areas that the worms have previously been where there is no food left). By contrast, after initially forming transient clusters on the lawn, N2 worms move radially outwards with no collective movement (Figure 1C, right).

Figure 1—figure supplement 1.

npr-1 swarming on a bigger food patch.

Forty npr-1 animals swarm over a big food patch, reminiscent of a persistent random walk. More than one large moving clusters co-exist towards the end of the video (orange and yellow), and a cluster (orange) disperses and re-forms elsewhere (orange and yellow) when it crosses its previous path (blue), presumably due to local food depletion.

Fluorescence imaging and automated animal tracking allows quantification of dynamics inside and outside of aggregates

Based on our observation that swarming appears to be driven by food depletion, we hypothesize the phenomenon may be a dynamic extension of the initial aggregation that occurs before depletion. To test this idea, we first sought to identify the mechanisms underlying aggregation.

The presence of aggregates is clear in bright field images, but it is difficult to track individual animals in these strongly overlapping groups for quantitative behavioral analysis. We therefore labeled the pharynx of worms with green fluorescent protein (GFP) and used fluorescence imaging in order to minimize overlap between animals (Video 2), making it possible to track most individuals even when they are inside a dense cluster (Figure 2A). We also labeled a small number of worms (1–3 animals out of 40 per experiment) with a red fluorescent protein (RFP)-tagged body wall muscle marker instead of pharynx-GFP. These RFP-labeled worms were recorded on a separate channel during simultaneous two-color imaging (Figure 2B), thus allowing both longer trajectories and the full posture to be obtained in a subset of animals. We wrote a custom module for Tierpsy Tracker (Javer et al., 2018) to segment light objects on a dark background and to identify the anterior end of the marked animals automatically, in order to extract trajectories and skeletons of multiple worms from our data (Figure 2C).

Figure 2.

Fluorescence multi-worm tracking.

(A) npr-1 mutant and N2 animals exhibit different social behaviors on food, with the former being hyper-social (top left) and the latter being hypo-social (top right). Using a pharynx-GFP label (bottom row), individual animals may be followed inside a cluster. (B) In two-color experiments, worms are either labeled with pharynx-GFP (left) or body wall muscle-RFP (middle). As the two colors are simultaneously acquired on separate channels, the selected few RFP-labeled individuals are readily segmented and may be tracked for a long time, even inside a dense cluster. (C) Tierpsy Tracker tracks multiple worms simultaneously, generating both centroid trajectories (left, image color inverted for easier visualization; multiple colors show distinct trajectories) and skeletons (middle, pharynx-marked animal; right, body wall muscle-marked animal; red dots denote the head nodes of the skeleton).

Ascarosides and direct adhesion are unlikely to drive different aggregation phenotypes

We first considered long-range chemotaxis driven by food or diffusible ascaroside pheromone signals as a potential behavioral mechanism. Chemotaxis towards food can likely be ignored as our experiments were performed on thin, even bacterial lawns, and worms are mostly on food during the aggregation phase of the experiments (99.7 ± 0.4% for npr-1 and 99.8 ± 0.3% for N2, mean ±S.D.). Although ascarosides are important for processes such as mating and dauer formation in C. elegans (Srinivasan et al., 2008), it is less clear whether long-range signaling via pheromones plays a role in aggregation (de Bono et al., 2002; Macosko et al., 2009). daf-22(m130) mutants do not produce ascarosides, but daf-22;npr-1 double mutants aggregate similarly to npr-1 single mutants (Figure 3—figure supplement 1), consistent with the observation that the hermaphrodite-attractive pheromone icas#3 is attractive to both N2 animals and npr-1 mutants (Srinivasan et al., 2012) and is thus unlikely to explain the difference in their propensity to aggregate. Moreover, attraction between moving objects is known to produce aggregation in active matter systems (Redner et al., 2013a), but it is not known whether this applies to worms. Short-range attraction between worms may exist in the form of adhesion mediated through a liquid film (Gart et al., 2011), but we have no reason to believe this would differ between npr-1 and N2 strains.

Figure 3—figure supplement 1.

Pheromones appear unimportant for aggregation.

npr-1 and N2 animals with pheromones removed by a daf-22 mutation aggregate to similar levels as their pheromone-intact counterparts. Top row: snapshots of 40 worms from each strain behaving on a thin, uniform lawn. Bottom left: quantification of cluster area relative to single worm area for each strain; dashed line shows the cut-off values used to generate the violin plot on the bottom right. Bottom right: probability of having a relative cluster area above the threshold value (dashed line on the bottom left). Blob area were extracted as tracking features. For each recording, a random sample (without replacement) of 500 single worms was used to calculate single-worm mean area, which was used to normalize multi-worm cluster areas from that recording. Relative cluster area values for each strain were pooled across recording replicates, and histograms were created with a bin width of 0.5.

Reversal rates and speed depend on neighbor density more strongly in npr-1 mutants than in N2

Having considered long-range food- or ascaroside-mediated attraction and short-range adhesion, we next focused on behavioral responses to nearby neighbors. While postural changes do not seem to be a main driver of aggregation as principal component analysis of lone versus in-cluster npr-1 worms revealed similar amplitudes in the posture modes (Figure 3—figure supplement 2), we found experimental evidence for density-dependence of both reversal rates and speed and that these differ between the two strains we studied.

Figure 3—figure supplement 2.

Shape analysis for lone and in-cluster npr-1 worms.

Left two panels: first four eigenworms (Stephens et al., 2008) plotted in real space for projections of lone worms and in-cluster worms. Right: variance explained as a function of the number of eigenworms. Eigenworms are based on common reference (Brown et al., 2013) set for both strains and worm categories.

Reversals have been previously suggested as a behavior that may enable npr-1 worms to stay in aggregates (Rogers et al., 2006). To avoid cluster definitions based on thresholding the distance between worms, we quantified individual worm behavior as a function of local density (Figure 3A) instead. Calculating the reversal rates relative to that of worms at low densities, we found that npr-1 mutants reverse more at increased neighbor densities, while N2 animals do not (Figure 3B).

Figure 3.

Individual-level behavioral quantification.

(A) Schematic explaining k-nearest neighbor density estimation. (B) Relative rate of reversals as a function of local density (k-nearest neighbor density estimation with k = 6) for npr-1 (blue) and N2 (orange) strains. Lines show means and shaded area shows the standard error (bootstrap estimate, 100 samples with replacement). (C) Distributions of crawling speeds at different local neighbor densities for both strains. Lines show histograms of speeds for each density bin, and the color of the line indicates the density (blue is high, magenta is low). (D) Midbody absolute speed for manually annotated npr-1 cluster entry (left, n = 28) and exit events (right, n = 29). Each event was manually identified, with time 0 representing the point where the head or tail of a worm starts to enter (left) or exit (right) an existing cluster. Skeleton xy-coordinates were linearly interpolated for missing frames for each event, before being used to calculate midbody speed extending 20 s on both sides of time 0 of the event. Speeds were smoothed over a one-second window. Shading represents standard deviation across events. Each red line shows the midbody absolute speed of a selected event that is shown in Video 3 (left) or Video 4 (right).

npr-1 and N2 animals with pheromones removed by a daf-22 mutation aggregate to similar levels as their pheromone-intact counterparts. Top row: snapshots of 40 worms from each strain behaving on a thin, uniform lawn. Bottom left: quantification of cluster area relative to single worm area for each strain; dashed line shows the cut-off values used to generate the violin plot on the bottom right. Bottom right: probability of having a relative cluster area above the threshold value (dashed line on the bottom left). Blob area were extracted as tracking features. For each recording, a random sample (without replacement) of 500 single worms was used to calculate single-worm mean area, which was used to normalize multi-worm cluster areas from that recording. Relative cluster area values for each strain were pooled across recording replicates, and histograms were created with a bin width of 0.5.

Left two panels: first four eigenworms (Stephens et al., 2008) plotted in real space for projections of lone worms and in-cluster worms. Right: variance explained as a function of the number of eigenworms. Eigenworms are based on common reference (Brown et al., 2013) set for both strains and worm categories.

Video 3.

A single event showing switch from high to low motility state prior to cluster entry (fluorescence 40 worm aggregation imaging).

The red worm at the bottom (arrow) decreases speed before entering a cluster. Inset: midbody absolute speed of that individual with respect to time 0 as the point of the head entering a cluster; open blue circle shows the current speed matched to the video frame.

Video 4.

A single event showing switch from low to high motility state prior to cluster exit (fluorescence 40 worm aggregation imaging).

The red worm increases speed before exiting a cluster. Inset: midbody absolute speed of that individual with respect to time 0 as the point of the head exiting a cluster; open blue circle shows the current speed matched to the video frame.

Next we calculated the speed distributions of individual worms, binned by local neighbor density. We found that both strains slow down when surrounded by many other worms, but the shift is more pronounced for npr-1 animals. npr-1 worms move faster than N2 at low densities, showing a distinct peak at high speeds. As neighbor density increases, this high speed peak gradually becomes replaced by a peak at low speeds, so that the overall speed distribution for npr-1 resembles that of N2 at very high densities. Thus, npr-1 and N2 animals show different density-dependent changes in their respective speed profiles (Figure 3C).

Since the observed transition of the speed profiles could occur due to active behavioral changes as well as restricted movement in clusters, we also considered tracks of individual worms. Using body wall muscle-marked worms allowed us to obtain longer trajectories that could be joined for the duration of an entire video, including cluster entry and exit events. We compared the speed of these tracks with visual assessment of when a worm entered or exited a cluster based on the proximity to pharynx-labeled worms. We found that worms are able to move inside of clusters and observed that speed changes can occur prior to cluster entry and exit events (Figure 3D, Video 3 and Video 4). This change of speed is neither purely mechanical nor a deterministic response to a certain neighbor density, and suggests a mechanism in which worms probabilistically switch between different speeds.

Spatial statistics show group-level differences between npr-1 and N2 animals

The differences in aggregation behavior between npr-1 and N2 are visually striking, but previous quantification has typically been limited to the fraction of animals in clusters. Using the tracked positions of pharynx-labeled worms (Figure 4A), we calculated the pair-correlation function (Figure 4B), commonly used to quantify aggregation in cellular and physical systems (Gurry et al., 2009). We also computed a hierarchical clustering of worm positions (Figure 4C), which is calculated from the same pairwise distances but emphasizes larger scale structure. Using both measures, we found that as a population, npr-1 animals show quantifiably higher levels of aggregation than N2, especially at scales up to 1 mm (pair-correlation ‘S’, Figure 4D) and 2 mm (hierarchical clustering ‘S’, Figure 4E). We also quantified aggregation using scalar spatial statistics, namely the average standard deviation (‘S’) and kurtosis (‘S’) of the distribution of positions. This confirms that the positions of npr-1 worms are less spread-out and more heavy-tailed than those of N2 (Figure 4D).

Figure 4.

Population-level behavioral quantification.

(A) Positions of npr-1 worms in an example frame. (B) Schematic explaining pair correlation function (S), which counts the number of neighbors at a distance r, normalized by the expectation for a uniform distribution. (C) Example dendrogram from which hierarchical clustering branch length distributions (S) can be calculated. (D) Pair correlation function for npr-1 (blue) and N2 (orange). Lines show mean and shaded area shows standard error of the mean. (E) Hierarchical clustering branch length distributions for npr-1 (blue) and N2 (orange). Histograms show relative frequency of inter-cluster distances (single linkage distance in agglomerative hierarchical clustering, equivalent to the branch lengths in the example dendrogram in (C)). (F) Mean standard deviation (S) and kurtosis (S) of the positions of worms, with the mean taken over frames sampled.

Agent-based model captures different aggregation phenotypes

To test whether the individual behavioral differences measured between npr-1 and N2 worms are sufficient to give rise to the observed differences in aggregation, we constructed a phenomenological model of worm movement and interactions. The model is made up of self-propelled agents (Figure 5A), and includes density-dependent interactions motivated by the experimental data, namely reversals at the edge of a cluster (Figure 5B) and a switch between movement at different speeds (Figure 5C). As a model of collective behavior this differs from those commonly considered in the literature, such as the Vicsek model (Vicsek et al., 1995) and its many related variants (Vicsek and Zafeiris, 2012; Yates et al., 2011). Such models typically feature attractive forces or align the direction of motion at ranges much longer than the size of the moving objects, and result in flocking or clustering with global alignment (Figure 5D), which we do not observe in our experimental data. In contrast, our model needs to produce dynamic, disordered aggregates (Figure 1B, Figure 2A and Video 2), and should primarily rely on short-range interactions that are motivated by behaviors measured in our data.

Figure 5.

Agent-based modeling of emergent behavior.

(A) Schematic of individual worm in the agent-based model. Each worm is made up of M nodes (here M = 18), connected by springs to enforce non-extensibility. Each node undergoes self-propelled movement, with the head node (red dot) undergoing a persistent random walk, and the rest of the nodes follow in the direction of the body. (B) Schematic of simulated reversals upon exiting a cluster. Each worm registers contact at the first and last 10% of its nodes within a short interaction radius. If contact is registered at one end but not the other, the worm is leaving a cluster and thus reverses with a Poisson rate dependent on the local density. (C) Schematic of density-dependent switching between movement speeds. Worms stochastically switch between slow and fast movement with Poisson rates kslow and kfast, which increase linearly and decrease exponentially with neighbor density, respectively. (D) Snapshots of simulations with commonly considered aggregation mechanisms, which produce unrealistic behavior for worm simulations, with flocking and highly aligned clustering. Arrows indicate the direction of movement of large clusters. (E) Phase portrait of model simulations, showing snapshots from the last 10% of each simulation, for different values of the two free parameters: density-dependence of the reversal rate and density-dependence of speed-switching (here kslow = kfast). Blue and orange panels highlight best fit for npr-1 and N2 data, respectively. (F) Summary statistics S (pair correlation, top) and S (hierarchical clustering, bottom) for the simulation which most closely matches the experimental data for the npr-1 and N2 strains (blue and orange panels in (E), respectively).

The density-dependence of the reversal rate and speed switching is implemented as follows: The rate of reversals increases linearly with density with slope r’, which is a free parameter, and is thus given by rrev = r’ ρ. The reversal rate at zero density is zero as we ignored spontaneous reversals outside of clusters as these were only rarely observed under our experimental conditions (see Appendix 1 for further discussion of the model construction). This parameterization of the reversal rate may be unbounded, but we can prevent unrealistically high reversal rates for a given maximum worm number by choosing our prior distribution of the parameter r’. The rate of slowing down is similarly approximated as a linear function of density, with free parameter ks’, and is given by kslow = ks0+ks’ ρ, where ks0 is the slowing rate at zero-density. The rate of speeding up is given by kfast = kf0 exp[-kf’ ρ], where the exponential decay is chosen to ensure positivity of the rate, and kf0 is the rate at zero density. The rates of slowing down and speeding up at zero density (ks0, kf0) were obtained from published single-worm experimental data (Javer et al., 2018; Yemini et al., 2013).

We initially ran a coarse parameter sweep, sampling uniformly in the two-dimensional parameter space associated with the density-dependence of reversals and speed switching. As a simplifying assumption, the density-dependence of the speeding-up and slowing-down rates was set equal (k’s = k’f = k’). The remaining parameters, r’ and k’, were varied to explore the global model behavior. This demonstrates that our model can capture different aggregation phenotypes from solitary movement to aggregation (Figure 5E) by varying just two free parameters, and provides important general insights. Inspection of the model simulations shows that each behavior alone (just reversals or slowing) does not give the same level of aggregation as when both parameters are modulated (Figure 5E), so that using both behavioral components proves important. Quantifying the aggregation and comparing it to the npr-1 experiment, however, highlights incomplete quantitative agreement with both the pair correlation function and hierarchical clustering distribution (Figure 5F). Thus, we reasoned additional interactions may be required to match the experimentally observed behaviors.

Adding a medium-range taxis interaction promotes stronger aggregation

To explore improvements in clustering, we extended the model by an attractive taxis interaction. Attraction should intuitively improve clustering, but we knew from our model exploration that an attractive potential between bodies produces undesirable cluster shapes (Figure 5D) and reasoned that a long-range interaction may be unrealistic (Figure 3—figure supplement 1). Thus, we include taxis towards neighboring worms and model worm movement as an attractive persistent random walk. The taxis contribution to a worm’s motile force has an overall strength controlled by parameter ft, with multiple nearby neighbors contributing cumulatively, weighted by 1/r, where r is the distance to a neighboring worm. Neighboring worms beyond a cut-off distance equal to the length of a worm have no contribution. Thus, this taxis interaction is acting at a natural intermediate length scale of our system (see Appendix 1 for details).

The resulting extended model has four free parameters: density-dependent reversals (r′), speed-switching rates (ks′, kf′) and taxis (ft). To find the parameter combinations that best describe each strain, as well as the uncertainty in the parameter values, we used an approximate Bayesian inference approach (see Appendix 1). To increase the computational efficiency of our inference pipeline, we excluded infeasible regions of parameter space to reduce the prior distribution of parameters that we need to sample from (Figure 6—figure supplement 1) (see Appendix 1). We then selected the closest matching simulations from about 27,000 simulations for npr-1 and about 13,000 simulations for N2, equally weighting all four summary statistics. Results from our extended model (Figure 6A, Video 5 and Video 6) show markedly improved quantitative agreement with the experiments (Figure 6B). The approximated posterior distributions of the parameters (Figure 6C–D) show the most likely values of the parameters for each strain, as well as the uncertainty associated with the individual and joint marginal parameter distributions. In particular, to achieve npr-1-like aggregation, the reversal (r’) and taxis (ft) parameters need to be higher than for N2, albeit not too high. The density-dependence of the slowing rate (k’s) is only subtly different between the two strains, while the dependence of the speeding up rate (k’f) is greater in npr-1, but with broader uncertainty.

Figure 6—figure supplement 1.

Reduced prior distribution used for approximate Bayesian inference of extended model.

Marginal and joint prior parameter distributions for npr-1 (A) and N2 (B), that have been constructed from a set of pilot runs excluding any parameter combinations that lead to stable pairs for either strain, unstable clusters for npr-1, and stable clusters for N2. Remaining parameter values were used to construct the prior distributions via kernel-density estimation. See Appendix 1 for details.

Figure 6.

Model with taxis captures quantitative aggregation phenotypes.

(A) Sample snapshot of the closest matching simulations for npr-1 (top) and N2 (bottom). (B) Summary statistics for npr-1 (orange) and N2 (blue): S: pair correlation function; S: hierarchical clustering distribution; S: standard deviation of positions; S: kurtosis of positions. Solid lines show the closest matching simulations; dashed lines show sample mean over the posterior distribution; and dotted lines show experimental means, with error bars showing standard deviation of 13 (npr-1) and 9 (N2) replicates. (C–D) Approximate posterior distribution of parameters for npr-1 (C) and N2 (D). Diagonal plots show marginal distribution of each parameter, off-diagonals show pairwise joint distributions. Parameters are: increase in reversal rate with density, r'; increase in rate to slow down, k's; decrease in rate to speed up, k'f; and contribution of taxis to motile force, ft.

Marginal and joint prior parameter distributions for npr-1 (A) and N2 (B), that have been constructed from a set of pilot runs excluding any parameter combinations that lead to stable pairs for either strain, unstable clusters for npr-1, and stable clusters for N2. Remaining parameter values were used to construct the prior distributions via kernel-density estimation. See Appendix 1 for details.

(A) Simulation for parameters equal to the mean of the posterior distribution for the npr-1 strain (Figure 6C). (A1), Sample snapshot of simulation; (A2), Pair correlation statistic, averaged over ten simulations (solid line), and standard error of the mean (error bars), with experimental reference (dotted line, as in Figure 6B); (A3), Hierarchical clustering distribution, averaged over ten simulations (solid line), and standard error of the mean (error bars), with experimental reference (dotted line, as in Figure 6B); (A4), Combined score for model agreement with experiment (lower score is better) for summary statistics in (A1) and (A2), calculated as the difference in the logarithm of the summary statistics between experiment and simulation (see Appendix 1 for details). (B) As in (A) but with r' = 0 (no reversals). (C) As in (A) but with worms always moving at the faster speed. (D) As in (A) but with ft = 0 (no taxis towards other worms). (E) As in (A) but with η = 0 (no directional noise in movement). (F) As in (A) but with η = 0.005. This η represents directional noise 10 times lower than in (A). (G) As in (A) but with η = 0.08. This η represents higher noise than in (A), and roughly corresponds to the velocity autocorrelation measured for interacting worms in our experiments (Figure 6 – figure supplement 4A2-3). (H) As in (A) but with sinusoidal undulations in the direction of movement, with a frequency similar to that of npr-1 worms (see Appendix 1 for details).

(A) Orientational correlation quantifies the alignment of the pharynxes (experiments, (A1), or first three nodes (simulations, A2-4) between pairs of worms a given distance apart. A value of 1 corresponds to parallel alignment, and −1 to anti-parallel alignment. Solid lines show the average directional correlation and shaded area shows the 95% confidence interval. (A1), Experimental measurements; (A2), Simulation for parameters equal to mean of posterior distribution for npr-1 strain (Figure 6C); (A3), As in (A2) but with ft = 0 (no taxis towards other worms); (A4), As in (A2) but with worms always moving at the faster speed. (B) Velocity correlation quantifies the alignment of movement directions between pairs of worms a given distance apart. A value of 1 corresponds to worms moving in the same direction, and −1 to worms moving in opposite directions. Solid lines show the average directional correlation and shaded area shows the 95% confidence interval. (B1), Experimental measurements; (B2), Simulation for parameters equal to the mean of the posterior distribution for the npr-1 strain (Figure 6C); (B3), As in (B2) but with ft = 0 (no taxis towards other worms). (C) Correlation between velocity of a worm and the direction to each neighbor was calculated to quantify the degree of taxis towards other worms. A value of 1 corresponds to worms moving directly towards a neighbor, and −1 to directly moving away from a neighbor. Solid lines show the average directional correlation and shaded area shows the 95% confidence interval. (C1), Experimental measurements; (C2), Simulation for parameters equal to the mean of the posterior distribution for npr-1 strain (Figure 6C); (C3), As in (C2) but with r' = 0 (no reversals).

(A) Velocity autocorrelation. (A1), From experiments with body wall muscle-tracked single worms on circular food patches; (A2), From experiments with 40 worms, of which a few were body wall muscle-tracked to allow acquisition of longer trajectories; (A3), From simulated, non-interacting worms undergoing a persistent random walk for different parameter values of η, the strength of the angular noise. The dashed line shows a value of 0.23, corresponding approximately to the expected correlation for choosing angles at random, uniformly distributed between −3/4π and 3/4π, thus representing an almost complete reorientation with respect to the original direction of motion. Note that this level is reached after about 15 s for η = 0.05 and for single worms (A1), and after about 8 s for η = 0.08 and interacting worms (A2). (B) Relative reversal rates at various local densities from experiments (solid lines and shaded 95% confidence interval, same data as in Figure 3B) and from model equations for reversal rates, parameterized with the mean of posterior distribution (dotted lines). (C) Speed switching rates at various local densities. (C1), Ratio of worms moving at fast (up to 350 μm/s) versus slow (<100 μm/s for npr-1,<50 μm/s for N2) speeds as measured in experiments; (C2), Ratio of worms moving at fast versus slow (<100 μm/s for npr-1,<50 μm/s for N2) speeds as measured in simulations with posterior mean parameters, showing average over ten (npr-1) and eight (N2) simulations, error bars showing error in the mean. The disagreement may indicate that the exponential form of kf(ρ) (see Figure 5C and main text) is only a rough estimate. For (B) and (CC), inferred model parameters were converted to units of worms/mm2.

(A) Simulations with decreasing length of agents. (A1), Snapshot of simulation with posterior mean parameter values for npr-1, as in Figure 6—figure supplement 2A. Worms have M = 18 nodes and a total length of L = 1.2 mm; (A2), Modified model with M = 9 nodes and shorter worms (but same width) still produces aggregation, without readjusting other parameters; (A3), As in (A2) but with M = 6 nodes per worm and shorter total length; (A4), As in (A2) but with M = 5 nodes per worm and shorter total length; (A5), at M = 4 nodes per worm and corresponding length of L = 0.3211 mm, stable aggregates comprising all worms fail to form. At this worm length, the interaction radii of head and tail nodes start to overlap, and worms require a difference in contact between head and tail to initiate reversals in our simulations. (B) Simulations with volume exclusion. (B1), Snapshot of simulation where volume exclusion is enforced, such that worms cannot overlap (apart from themselves), without adjusting any other parameters. The number of nodes per worm has been increased to M = 45 to ensure sufficient overlap between nodes within a worm. Pair correlation function (B2) and hierarchical clustering distribution (B3) show that aggregate is spread out and less dense compared to experiments (dotted line). Solid lines show mean over three simulations and error bars show standard deviation.

Video 5.

Sample model (with taxis) simulation describing npr-1 mutants.

The video plays at 30x the normal speed.

Video 6.

Sample model (with taxis) simulation describing N2.

The video plays at 30x the normal speed.

Figure 6—figure supplement 4.

Additional comparison of model parameters with experimental measurements.

(A) Velocity autocorrelation. (A1), From experiments with body wall muscle-tracked single worms on circular food patches; (A2), From experiments with 40 worms, of which a few were body wall muscle-tracked to allow acquisition of longer trajectories; (A3), From simulated, non-interacting worms undergoing a persistent random walk for different parameter values of η, the strength of the angular noise. The dashed line shows a value of 0.23, corresponding approximately to the expected correlation for choosing angles at random, uniformly distributed between −3/4π and 3/4π, thus representing an almost complete reorientation with respect to the original direction of motion. Note that this level is reached after about 15 s for η = 0.05 and for single worms (A1), and after about 8 s for η = 0.08 and interacting worms (A2). (B) Relative reversal rates at various local densities from experiments (solid lines and shaded 95% confidence interval, same data as in Figure 3B) and from model equations for reversal rates, parameterized with the mean of posterior distribution (dotted lines). (C) Speed switching rates at various local densities. (C1), Ratio of worms moving at fast (up to 350 μm/s) versus slow (<100 μm/s for npr-1,<50 μm/s for N2) speeds as measured in experiments; (C2), Ratio of worms moving at fast versus slow (<100 μm/s for npr-1,<50 μm/s for N2) speeds as measured in simulations with posterior mean parameters, showing average over ten (npr-1) and eight (N2) simulations, error bars showing error in the mean. The disagreement may indicate that the exponential form of kf(ρ) (see Figure 5C and main text) is only a rough estimate. For (B) and (CC), inferred model parameters were converted to units of worms/mm2.

Figure 6—figure supplement 2.

Core model components, but not noise and undulations in movement, are necessary for quantitative agreement with aggregation summary statistics.

(A) Simulation for parameters equal to the mean of the posterior distribution for the npr-1 strain (Figure 6C). (A1), Sample snapshot of simulation; (A2), Pair correlation statistic, averaged over ten simulations (solid line), and standard error of the mean (error bars), with experimental reference (dotted line, as in Figure 6B); (A3), Hierarchical clustering distribution, averaged over ten simulations (solid line), and standard error of the mean (error bars), with experimental reference (dotted line, as in Figure 6B); (A4), Combined score for model agreement with experiment (lower score is better) for summary statistics in (A1) and (A2), calculated as the difference in the logarithm of the summary statistics between experiment and simulation (see Appendix 1 for details). (B) As in (A) but with r' = 0 (no reversals). (C) As in (A) but with worms always moving at the faster speed. (D) As in (A) but with ft = 0 (no taxis towards other worms). (E) As in (A) but with η = 0 (no directional noise in movement). (F) As in (A) but with η = 0.005. This η represents directional noise 10 times lower than in (A). (G) As in (A) but with η = 0.08. This η represents higher noise than in (A), and roughly corresponds to the velocity autocorrelation measured for interacting worms in our experiments (Figure 6 – figure supplement 4A2-3). (H) As in (A) but with sinusoidal undulations in the direction of movement, with a frequency similar to that of npr-1 worms (see Appendix 1 for details).

To address whether all three behaviors (reversals, speed changes, and taxis) were necessary for aggregation we ran additional simulations: starting from the mean of the posterior distribution for npr-1 (Figure 6C) as a reference, we removed individual model components by setting the corresponding parameters to zero. These perturbed simulations show that removing speed switching or taxis from the model disrupts aggregation, while removing reversals reduces the overall quantitative agreement with experimental data (Figure 6—figure supplement 2, B–D). In some cases, removing individual model behaviors also produced correlations of velocity and orientation between neighbors that are different from what we measure in experiments (Figure 6—figure supplement 3). Thus, we conclude that we have identified sufficient behavioral components for aggregation, and that these are also necessary to quantitatively match aggregation in npr-1 mutants.

Figure 6—figure supplement 3.

Analysis of orientational and velocity correlations in experiments and simulations.

(A) Orientational correlation quantifies the alignment of the pharynxes (experiments, (A1), or first three nodes (simulations, A2-4) between pairs of worms a given distance apart. A value of 1 corresponds to parallel alignment, and −1 to anti-parallel alignment. Solid lines show the average directional correlation and shaded area shows the 95% confidence interval. (A1), Experimental measurements; (A2), Simulation for parameters equal to mean of posterior distribution for npr-1 strain (Figure 6C); (A3), As in (A2) but with ft = 0 (no taxis towards other worms); (A4), As in (A2) but with worms always moving at the faster speed. (B) Velocity correlation quantifies the alignment of movement directions between pairs of worms a given distance apart. A value of 1 corresponds to worms moving in the same direction, and −1 to worms moving in opposite directions. Solid lines show the average directional correlation and shaded area shows the 95% confidence interval. (B1), Experimental measurements; (B2), Simulation for parameters equal to the mean of the posterior distribution for the npr-1 strain (Figure 6C); (B3), As in (B2) but with ft = 0 (no taxis towards other worms). (C) Correlation between velocity of a worm and the direction to each neighbor was calculated to quantify the degree of taxis towards other worms. A value of 1 corresponds to worms moving directly towards a neighbor, and −1 to directly moving away from a neighbor. Solid lines show the average directional correlation and shaded area shows the 95% confidence interval. (C1), Experimental measurements; (C2), Simulation for parameters equal to the mean of the posterior distribution for npr-1 strain (Figure 6C); (C3), As in (C2) but with r' = 0 (no reversals).

Searching for evidence of taxis in the experimental tracking data, we calculated the correlation between worm velocity and the vector towards nearby worms, and found this correlation to be nearly zero in both experiments and simulations for all distances up to 2 mm (Figure 6—figure supplement 3B1–2), which is larger than the size of a typical worm cluster. This may not be intuitive, and we suspect the reason is twofold: (a) the taxis effect is only a small influence on the instantaneous direction of the movement of a worm, compared to persistence and noise; and (b) we only tracked the pharynx in our experiments, and reproduced this restriction in our analysis of simulations, but the whole body of the worm is likely giving relevant cues to any chemical or mechanical taxis. Our methodology that enables us to track inside worm clusters therefore brings with it the caveat that there is unseen worm density that affects any potential taxis behavior, but which remains undetectable in our tracking. Thus, our analysis shows that a taxis behavior similar to our simulations may be present in experiments, even if it is difficult to detect with correlation analysis. We compared the other inferred parameters with experimental measurements: The reversal rate shows a similar increase with density that is greater for npr-1 than N2 (Figure 6—figure supplement 4B). The speed switching rates could only be compared indirectly by calculating the ratio of fraction of worms in fast vs. slow movement in experiments (Figure 6—figure supplement 4C1) and model simulations (Figure 6—figure supplement 4C2). The disagreement may indicate that the exponential form of k(ρ) is only a rough approximation. However, aggregation in the model is not sensitive to speed switching rates, as shown by the broad posterior distributions for the inferred parameters (Figure 6C–D).

Extending the model with food-depletion captures dynamic swarming

Since we hypothesize that the swarming we observed at longer time scales may be explained as aggregation under food depletion conditions, we further extended the model to allow the local depletion of food. Food is initially distributed uniformly, and becomes depleted locally by worm feeding (see Appendix 1 for details). Absence of food suppresses the switch to slow speeds, thus causing worms to speed up when food is locally depleted. As a result, we hypothesize that worm clusters begin to disperse but reform on nearby food, leading to sweeping.

Selecting the parameter combination best matching the npr-1 strain (Figure 6) and an appropriate food depletion rate (chosen such that all food was depleted no faster than observed in experiments), the resulting simulation produced long-time dynamics qualitatively representative of the experimentally observed swarming (Figure 7A–B, Video 7). Worm clusters undergo a persistent but not necessarily directed random walk, can disperse and re-form elsewhere, and multiple clusters may co-exist, all of which we observe experimentally. Tracking the centroid of worms in our simulations, we find a comparable cluster speed as the median experimental value of 172 μm/min (Figure 1D) for a range of feeding rates (Figure 7C) (feeding rate is an unknown parameter as our model only accounts for relative food concentration). Thus, the model indicates that dynamic swarming of npr-1 aggregates may be explained as an emergent phenomenon resulting from individual locomotion, and that the same behavioral mechanisms that produce the initial aggregates, when coupled with local food depletion, give rise to the observed swarming behavior.

Figure 7.

Simulations capture dynamic swarming.

(A) Snapshots of aggregation simulation with food depletion. Background color shows relative food concentration with white indicating high food and black indicating no food. (B) Visualization of worm positions in (A) over time, showing cluster displacement. Note the periodic boundary conditions. (C) Cluster speed at various feeding rates relative to lawn thickness (other parameters equal to mean of posterior distribution for npr-1). The upward trend is expected: smaller lawn thickness leads to faster movement as worms run out of food quicker and need to re-form clusters on nearby food. Cluster speed is calculated the same way as in Figure 1D; error bars show median absolute deviation over five simulations. Dashed line indicates experimentally-derived median cluster speed (from Figure 1D) for comparison.

Video 7.

Sample swarming simulation describing npr-1 mutants.

Background color shows relative food concentration with white indicating high food and black indicating no food.The video plays at 30x the normal speed.

Discussion

We have investigated the mechanisms of aggregation and swarming in C. elegans collective feeding using quantitative imaging and computational modeling. We show that while a combination of increased reversals upon leaving aggregates and a neighbor density-dependent increase in speed switching rates is sufficient to produce aggregation, the addition of taxis towards neighbors improves the quantitative agreement between simulations and experiments. Removing any one of the core behavioral mechanisms (reversals, speed changes, taxis) from our model either disrupts aggregation or otherwise reduces the quantitative agreement with experiments (Figure 6—figure supplement 2–3). The proposed taxis might be driven by a shallow O2 or CO2 gradient created by a worm cluster (discussed further below), to additional chemical signals unaffected by daf-22 loss of function, or to another unknown mechanism. By extending the aggregation model to include food depletion, we show that the same behavioral mechanisms also underlie dynamic swarming in the hyper-social C. elegans strain, reminiscent of wild fires and other self-avoiding dynamics.

We focused on identifying phenomenological behavioral components giving rise to aggregation, while remaining agnostic as to the sensory cues causing the behaviors. The density-dependent interactions could arise from local molecular signaling, or be mediated through contact-sensing, and the 1/r dependence of the taxis interaction is compatible with a diffusible, non-degrading factor (such as CO2, or O2 depletion; dependence would likely be different for a pheromone depending on its degradation rate). Given that aggregates break up when ambient O2 concentration is reduced to 7% (Gray et al., 2004), the preferred concentration of npr-1 mutants, the most obvious candidate for the sensory cue guiding aggregation is O2 (Rogers et al., 2006). A simple hypothesis would be that oxygen consumption by worms locally lowers O2 concentration to the 5–12% preferred by npr-1 mutants, promoting their aggregation. To support this, Rogers et al. (2006) report low O2 concentrations inside worm clusters. However, non-aggregating N2 worms also prefer O2 concentrations lower than atmospheric (5–15%) (Gray et al., 2004). Furthermore, a strong reduction of oxygen concentration inside an aggregate to near 7% is unlikely based on reaction-diffusion calculations: the diffusion of oxygen through worm tissue, or their oxygen consumption, would need to be several orders of magnitude different from estimated values to create O2 gradients as steep as reported by Rogers et al. (Appendix 2—figure 1). However, as worms have been reported to respond even to small changes in oxygen concentration (McGrath et al., 2009), aggregation may still be mediated through a shallower local oxygen gradient.

Appendix 2—figure 1.

Oxygen consumption-diffusion calculations predict shallow O2 concentration gradients.

() Plot of feasible oxygen gradients inside worm aggregates. The oxygen concentration decays with length constant , with diffusion constant (in water) and oxygen consumption rate (estimated as an upper bound for 200 pl/min [Shoyama et al., 2009; Suda et al., 2005] at 21% oxygen and 8000 pl worm volume). The thinnest dimension of a cluster is relevant for diffusion, which is its thickness. We can approximate the cluster geometry either as flat, which results in a 1D diffusion gradient (solid line), or as hemispherical, which we approximate by spherically symmetric diffusion in 3D (dashed line). In either case the reaction-diffusion equation was solved at steady state. (B) Gradient of diffusible, non-degrading signal, ∂t for example CO2, outside a point source. Without decay, this problem is equivalent to calculating the λ potential around a point charge, and the concentration would be , in 3D, where λ is the production rate times the volume of a worm, 0.14/min (equal and opposite to the O2 consumption, based on mass conservation). A point source represents the contribution of a short section of a worm, and the contributions of many worms in an aggregate would integrate to give an approximately logarithmic gradient of signal outside the aggregate.

In this scenario, high ambient O2 concentration serves as a permissive signal for aggregation and a shallow oxygen gradient induces worms to stay inside aggregates. Our agent-based simulations are entirely compatible with this picture. Further experiments would be required to test the hypothesis that oxygen is playing such a dual role. One possibility would be to introduce mutations leading to aerobic metabolism deficiencies into npr-1 mutants. Such mutants would still be able to sense ambient oxygen, but are expected to produce an even weaker oxygen gradient in an aggregate. The resulting phenotype could then be compared quantitatively to model predictions, for example with reduced taxis and/or modified rates of density-dependent reversal and speed switching. Additionally, one may seek evidence for the ability of worms to sense a shallow oxygen gradient by repeating the gas-phase aerotaxis experiment described in Gray et al. (2004), but with a much smaller gradient (19–21%) in the light of our new calculations, to see if worms can sense and move towards environments where oxygen levels are only slightly below ambient concentrations. Further work quantifying the behavior of individual worms at different oxygen concentrations, such as during oxygen-shift experiments inside flow chambers where single animals experience acute switches between 21% and 19% oxygen, may also help to distinguish oxygen as a direct cue or part of the ‘sensory triggers that can initiate social behavior by activating chemotaxis or mechanotaxis’ (Gray et al., 2004).

The model of worm movement and interactions presented here was chosen for a balance of simplicity and realism, and is not necessarily unique. Our model comprises a persistent random walk of chain-like worms, which were loosely inspired by work on bacterial systems (Balagam and Igoshin, 2015). We have adopted Bayesian parameter inference to capture the uncertainty in our parameter estimates, and to enable flexible extension to additional experimental data or comparison of different models in future work. An alternative approach is to be entirely data-driven in the construction of the model and compute interactions between worms directly based on their tracked positions at every time step, as has been done in collective behavior of Myxococcus xanthus (Cotter et al., 2017; Zhang et al., 2018). This approach may require higher worm numbers and improved tracking, to ensure comparably large statistical sample sizes with bacterial studies. We have used experimental data to inform our modeling framework where appropriate (size, shape, speed of agents, and reversal and speed change rates at zero density), and verified that the aggregation outcome is robust and quantitatively similar to experimental results regardless of the amount of noise in the persistent random walk (Figure 6—figure supplement 2E–G), or the presence of undulations in agent movement (Figure 6—figure supplement 2H). We have further verified that aggregation still occurs with shorter simulated worms (and fewer nodes per worm), given they are long enough to detect a contact difference between head and tail when exiting a cluster, which is required to initiate reversals (Figure 6—figure supplement 5A). Lastly, in the model presented here, we have allowed for overlap between worms to reflect a degree of overlap in clusters when worms can crawl over each other. With volume exclusion our model still produces aggregation, although the clusters are less dense and more extended (Figure 6—figure supplement 5B).

Figure 6—figure supplement 5.

Aggregation model requires minimum length of simulated worms, and is robust to introducing volume exclusion.

(A) Simulations with decreasing length of agents. (A1), Snapshot of simulation with posterior mean parameter values for npr-1, as in Figure 6—figure supplement 2A. Worms have M = 18 nodes and a total length of L = 1.2 mm; (A2), Modified model with M = 9 nodes and shorter worms (but same width) still produces aggregation, without readjusting other parameters; (A3), As in (A2) but with M = 6 nodes per worm and shorter total length; (A4), As in (A2) but with M = 5 nodes per worm and shorter total length; (A5), at M = 4 nodes per worm and corresponding length of L = 0.3211 mm, stable aggregates comprising all worms fail to form. At this worm length, the interaction radii of head and tail nodes start to overlap, and worms require a difference in contact between head and tail to initiate reversals in our simulations. (B) Simulations with volume exclusion. (B1), Snapshot of simulation where volume exclusion is enforced, such that worms cannot overlap (apart from themselves), without adjusting any other parameters. The number of nodes per worm has been increased to M = 45 to ensure sufficient overlap between nodes within a worm. Pair correlation function (B2) and hierarchical clustering distribution (B3) show that aggregate is spread out and less dense compared to experiments (dotted line). Solid lines show mean over three simulations and error bars show standard deviation.

One advantage of using C. elegans to study animal collective behavior is the opportunity to experimentally control and perturb the system. It should be possible to experimentally modify the key behavioral parameters identified in this paper with mutations or acute stimulus delivery in order to test our model. For example, one can introduce a reversal phenotype with unc-4 mutations, or alter the speed switching rates with mutations that affect the roaming-dwelling transition. Controlled stimulus delivery has already been used in previous oxygen-shift experiments. The resultant experimental outcomes may then be compared to theoretical predictions. Thus, there are ample opportunities for future studies to further integrate experimental and theoretical methods in the study of C. elegans collective behavior.

Despite its extensive study in the lab, it is still uncertain whether aggregation and swarming have a function in the wild. Aggregation may serve to protect C. elegans from desiccation or UV radiation associated with the surface environment (Busch and Olofsson, 2012). C. elegans swarming on unpalatable bacteria may also facilitate predation, perhaps through the collective action of secreted molecules that overcome bacterial defenses (personal communication from J. Hodgkin and G.M. Preston) in a manner similar to the well-described cooperative predation strategy used by Myxobacteria xanthus (Muñoz-Dorado et al., 2016; Pérez et al., 2016). Moreover, social versus solitary foraging strategies may confer selective advantages in different food abundance, food quality, and population density environments (de Bono and Bargmann, 1998). The observation that aggregating strains are less fit in laboratory conditions (Andersen et al., 2014) suggested that social feeding is not an efficient strategy at least in abundant food conditions. However, the observed fitness difference between aggregating and non-aggregating strains is actually dissociable from the feeding strategy in the lab (Zhao et al., 2018), leaving the question unresolved. Furthermore, in other systems, social feeding can increase fitness in natural environments via improved food detection and intake (Cvikel et al., 2015; Li et al., 2014; Snijders et al., 2018). It would be time consuming to experimentally measure the feeding efficiency of different behavioral strategies for a wide range of food patch sizes, distributions, and qualities. The agent-based model used in this study presents an opportunity to use a complementary approach to finding conditions that may favor social feeding.

C. elegans bridges the gap between the commonly studied micro- and macro-scales, and finding the behavioral rules underlying this mesoscale system allows us to consider principles governing collective behavior across scales. Indeed, key behavioral rules identified here for C. elegans aggregation have been observed at other scales. Spontaneous reversals have been implicated in bacterial aggregation at the microscale (Mercier and Mignot, 2016; Starruß et al., 2012; Thutupalli et al., 2015). By contrast, aggregating worms reverse mainly in response to leaving a cluster rather than spontaneously, thus requiring more complex sensory processing and behavioral response than seen in bacterial systems. Furthermore, changes in movement speed are a common feature in motility-induced phase transitions (Großmann et al., 2016; Redner et al., 2013b; Abaurrea Velasco et al., 2018). The emergent phenomena observed in models of interacting particles generally range from diffusion-limited aggregation to jamming at high volume fractions to flocking of self-propelled rods through volume exclusion (in two-dimensions). In contrast, aggregation in C. elegans occurs at much lower numbers of objects (tens of worms) and lower densities (area fraction of 4–6%) than typically studied in this field (thousands of objects at area fractions of 20–80%), and the density dependence of motility changes again emphasizes the role of more complex sensing and behavioral modulations common in macroscale animal groups such as fish shoals (Ward et al., 2011). Thus, collective behavior of C. elegans at the mesoscale indeed draws from both ends of the size scale and complexity spectrum, linking the physical mechanisms familiar from microscopic cellular and active matter systems with the behavioral repertoire of larger multicellular organisms.

Our approach of decomposing aggregation into component behaviors through modeling may also have applications in quantitative genetics beyond the scope of our current study. While hyper-social npr-1 mutants and hypo-social N2 worms show phenotypic extremes, wild isolates of C. elegans aggregate to different degrees (de Bono and Bargmann, 1998). Previous work has shown that even a very small increase in the phenotypic dimensionality (from one to two) can reveal independent behavior-modifying loci (Bendesky et al., 2012). Thus inferring model parameters for data from multiple wild C. elegans strains would produce behavioral parameterizations that might serve as a powerful set of traits for finding further behavior-modifying loci.

Materials and methods

Animal maintenance and synchronization

C. elegans strains used in this study are listed in Key Resources Table above. All animals were grown on E. coli OP50 at 20°C as mixed-stage cultures and maintained as described (Brenner, 1974). All animals used in imaging experiments were synchronized young adults obtained by bleaching gravid hermaphrodites grown on E. coli OP50 under uncrowded and unstarved conditions, allowing isolated eggs to hatch and enter L1 diapause on unseeded plates overnight, and re-feeding starved L1’s for 65–72 hr on OP50.

Bright field high-number swarming imaging

The strain used here (Figure 1A and Video 1) is DA609. On imaging day, synchronized adults were collected and washed in M9 buffer twice before several hundred animals were transferred to a seeded 90 mm NGM plate using a glass pipette. After M9 is absorbed into the media, ten-hour time-lapse recordings were taken with a Dino-Lite camera (AM-7013MT) at room temperature (20°C) using the DinoCapture 2.0 software (v1.5.3.c) for maximal field of view. Two independent replicates were performed.

Bright field standard swarming imaging

Step-by-step protocol is available at dx.doi.org/10.17504/protocols.io.vybe7sn. All recordings from this dataset are listed in Supplementary file 2.

The strains used here (Figure 1B) are DA609 and N2. Prior to collecting the full dataset, a single batch of OP50 was grown overnight, diluted to OD600 = 0.75, aliquoted for use on each imaging day, and stored at 4°C until use. Imaging plates were 35 mm Petri dishes containing 3.5 mL low peptone (0.013% Difco Bacto) NGM agar (2% Bio/Agar, BioGene) to limit bacteria growth. A separate batch of plates was poured exactly seven days before each imaging day, stored at 4°C, and dried at 37°C overnight with the agar side down before imaging. The center of an imaging plate was seeded with a single 20 μL spot of cold diluted OP50 one to three hours before imaging. The overnight plate drying step allowed the bacteria to quickly dry atop the media in order to achieve a more uniform lawn by minimizing the ‘coffee ring’ effect that would thicken the circular edge of the bacterial lawn. For each imaging day, synchronized young adults were collected and washed in M9 buffer twice before 40 animals were transferred to a seeded imaging plate using a glass pipette.

Imaging commenced immediately following animal transfer in a liquid drop, on a custom-built six-camera rig equipped with Dalsa Genie cameras (G2-GM10-T2041). Seven-hour recordings with red illumination (630 nm LED illumination, CCS Inc) were taken at 25 Hz using Gecko software (v2.0.3.1), whilst the rig maintained the imaging plates at 20°C throughout the recording durations. Images were segmented in real time by the Gecko software. The recordings were manually truncated post-acquisition to retain aggregation and swarming dynamics only. The start time was defined as the moment when the liquid dried and the all the worms crawled out from the initial location of the drop, and the end time was when the food was depleted and worms dispersed with increased crawling speed. Twelve independent replicates were performed for each strain.

Bright field big patch swarming imaging

Step-by-step protocol is available at dx.doi.org/10.17504/protocols.io.vyhe7t6. All recordings from this dataset are listed in Supplementary file 2.

The experiments here (Figure 1—figure supplement 1) are identical to those in the bright field standard swarming imaging, except for two differences. First, the imaging plates were seeded with a 75 μL spot of diluted OP50 (OD600 = 0.38) and allowed to inoculate overnight at room temperature before being used for imaging the next day. Second, recordings were taken over 20 hr instead of seven. Eight independent replicates were performed for each strain.

Bright field pheromone imaging

Step-by-step protocol is available at dx.doi.org/10.17504/protocols.io.vyie7ue. All recordings from this dataset are listed in Supplementary file 2.

The strains used here (Figure 3—figure supplement 1) are DA609, N2, DR476, and AX994. Bacteria aliquots and imaging plates were prepared as in the bright field standard swarming imaging assay. For each imaging day, synchronized young adults were collected and washed in M9 buffer twice before 40 animals were transferred to a seeded imaging plate using a glass pipette. After M9 was absorbed into the media following worm transfer in liquid, imaging plates containing the animals were subjected to a gentle vibration at 600 rpm for 10 s on a Vortex Genie two shaker (Scientific Industries) to disperse animals and synchronize aggregation start across replicates. Imaging commenced 20 s after the vibration finish, using the same rig set-up as swarming imaging above, except one-hour recordings were taken. Images were segmented in real time by the Gecko software. At least eight independent replicates were performed for each strain. Automated animal tracking was performed post-acquisition using Tierpsy Tracker software (http://ver228.github.io/tierpsy-tracker/, v1.3), which we developed in-house (Javer et al., 2018). Images with were tracked with customized parameters to create centroid trajectories, 49-point worm skeletons, and a battery of features.

Fluorescence aggregation imaging

Step-by-step protocol is available at dx.doi.org/10.17504/protocols.io.vzje74n. All recordings from this dataset are listed in Supplementary file 2.

The strains used here (Figure 2, Videos 2–4) are OMG2, OMG10, OMG19, and OMG24. One-color imaging consisted of pharynx-GFP labeled worms only, whereas two-color imaging also included a small number of body wall muscle-RFP labeled worms that were recorded simultaneously on a separate channel (thus readily segmented from the rest of the worms). The latter was necessary to follow individuals over a long period of time, particularly while inside a cluster, as frequent pharynx collisions inside clusters lead to lost individual identities and broken trajectories. For two-color imaging, animals with different fluorescent markers were mixed in desired proportion (1–3 red animals out of 40 per experiment) during the washing stage before being transferred together for imaging.

The data collection paradigm was identical to the bright field pheromone imaging assay in terms of bacteria aliquots, imaging plate preparation, and vibration implementation following animal transfer. The difference is that image acquisition was performed on a DMI6000 inverted microscope (Leica) equipped with a 1.25x PL Fluotar objective (Leica), a TwinCam LS image splitter (Cairn) with a dichroic cube (Cairn), and two Zyla 5.5 cameras (Andor) to enable simultaneous green-red imaging with maximal field of view. One-hour recordings were taken with constant blue (470 nm, 0.8A) and green (cool white, 1.4A) OptoLED illumination (Cairn), and images were acquired with 100 ms exposure at 9 Hz using Andor Solis software (v4.29.30005.0). The microscopy room was maintained at 21°C throughout the recording durations. Ten or more independent replicates were performed for each strain. We were able to reproduce stereotyped aggregation dynamics across replicates under our experimental paradigm (Figure 1—figure supplement 2). Image segmentation and automated animal tracking was performed post-acquisition using Tierpsy Tracker software (v1.3) with customized parameters, to create centroid trajectories, obtain two-point skeleton from pharynx-labeled individuals and 49-point midline skeletons from body wall muscle-marked ones, and extract various features. For body wall muscle-marked animals, trajectories were manually joined where broken due to tracking errors.

Figure 1—figure supplement 2.

Stereotyped temporal dynamics.

One-hour fluorescence recordings of npr-1 animals under our experimental conditions consist of reproducible temporal dynamics encompassing three phases: transient (animals move about the lawn and start to form clusters), aggregation (clusters largely remain stable with individuals entering and exiting), and swarming (worms move across the lawn in persistent clusters) (see sample Video 2). Percentage of in-cluster worms remain largely consistent throughout the latter two phases, except that clusters remain in place during the aggregation phase and become dynamic during the swarming phase. Error bars represent standard deviation across 13 (npr-1) and 9 (N2) replicates. The average duration of each phase derived from npr-1 experiments are applied to N2 data to maintain temporal consistency, even though N2 does not exhibit aggregation or swarming. Subsequent quantitative analyses for both strains were restricted to using the data from the aggregation (for all but Figure 1D) or swarming (for Figure 1D) phase, in order to reveal the mechanisms necessary for producing aggregation or the dynamics of swarming, respectively.

Fluorescence aggregation tracking data analysis

The code for tracking data analysis is available at https://github.com/ljschumacher/wormTrackingAnalysis (Schumacher et al., 2019; copy archived at https://github.com/elifesciences-publications/wormTrackingAnalysis).

Tracked blobs were filtered for minimum fluorescence intensity and maximum area, to exclude any larvae and tracking artifacts, respectively, which appeared on the occasional plate. Local worm densities around each individual were calculated using k-nearest neighbor density estimation, where the density is k divided by the area of a circle encompassing the k-th nearest neighbor. We chose  and verified based on visual assessment that the overall distribution of local densities changes very little with increasing k.

Reversals were detected based on a change of sign of speed from positive to negative, which was calculated from the dot-product of the skeleton vector (of the pharynx) and the velocity vector, and smoothed with a moving average over half a second. We only counted reversals that were at least 50 µm in length, and that moved at least half a pixel per frame before and after the reversal. Reversal events thus detected where binned by their local density. For each density bin, reversal rate was estimated as the number of events divided by the time spent in forward motion for that bin. The variability was estimated using a subsampling bootstrap: the reversal rate was estimated 100 times, sampling worm-frames with replacement, and estimating mean and standard deviation.

Speed profiles were generated by binning the measured speed values for local density, and then creating a histogram of speed values for each density bin.

Summary statistics of aggregation, such as pair-correlation and hierarchical clustering, where calculated as described in Appendix 1.

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 article "Shared behavioral mechanisms underlie C. elegans aggregation and swarming" for consideration by eLife. Your article has been reviewed by two peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Naama Barkai as the Senior Editor. The following individuals involved in review of your submission have agreed to reveal their identity: Jonathan Hodgkin and Oleg A Igoshin.

The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.

The manuscript of Ding et al. investigates collective feeding in the roundworm C. elegans with combination of time-lapse imaging, image analysis and mathematical modeling. It's a very well-written paper that presents broadly interesting and intriguing set of observation. Combination of experimental, quantitative analysis and modeling is the major strength of the paper. The authors study a prey-induced worm collective motility behavior called "swarming" and propose a mechanistic basis for this phenomenon. Both individual and population-level behaviors are quantified. Subsequent modeling is used to identify three key behaviors that drive aggregation: edge reversals, a density-dependent switches in speeds, and taxis towards neighboring worms. The same models can account for swarming when local food depletion is taken into account.

Major points:

1) Importantly, the modeling approach could be improved given the wealth of the experimental data. Specifically, it is surprising that the authors chose a fully phenomenological approach to construct their agent-based models that uses very little of the quantified data. Furthermore, authors do little to quantitatively compare resulting behaviors of agents in the model to those of tracked worms. Therefore, the reviewers are not fully convinced that identified behaviors in the model are either necessary or sufficient (given considerable noise in individual behaviors) for the observed population behaviors.

To circumvent this issue, measured individual worm behaviors should be used in constructing and/or constraining their model (see for example References PMID: 30514635, PMID: 28533367). Importantly, the noise and variability of the movement of agent match those of worms. It is not clear if this is the case in the current model. Furthermore, comparison of the resulting behaviors of agents and worms would allow testing the model. For example, is there correlation (e.g. between the agent velocity and the vector to nearby worm) that can quantify the degree of taxis? Can this correlation be compared between the model and experiments? Similar quantitative characteristics of other proposed 'key behaviors" need to be compared. Can changes in worm behaviors in nutrient-rich vs depleted regions be quantified and fit into the model? Either of the two approaches or some combination of them is needed to show convincingly that the postulated behaviors are indeed observed and responsible for the observed population self-organization.

2) The authors discuss possible neuronal mechanisms underlying their behavioural rules. Ideally, it would be good to apply or at least suggest possible experimental tests of these mechanisms, in particular of the proposed important medium-range taxis between neighbouring worms. It is suggested that this might be O2 tension mediated: how might this be tested? They do provide convincing experiments that exclude any significant contribution of ascaroside pheromones to swarming.

Other comments:

- The work also mainly compares behaviour in the standard lab strain N2, which is 'hypo-social', to behaviour in a 'hyper-social' npr-1 mutant of N2. Neither of these exactly corresponds to the natural state of almost all wild C. elegans, so the analysis is somewhat artificial; this might be more explicitly admitted.

- The authors state that to their knowledge "this swarming phenotype has not been reported in C. elegans previously", which betrays ignorance. For example, the following observation in Hodgkin and Barnes (2002, PMID: 1684664): "When a population has consumed almost all the bacteria, the worms… form a swarming mass that moves as a wave across the remaining bacterial lawn, and then disperses when all the bacteria have been consumed." This indicates that N2 swarms under appropriate conditions, perhaps of higher bacterial food density. The authors' experiments were performed only on "thin, even bacterial lawns".

- They also state that "it is still unclear whether aggregation and swarming have a function in the wild". Swarming is almost certainly naturally advantageous in predation on pathogenic or repellent bacteria, as a result of the 'wolfpack' effect that has been studied in myxobacterial swarms, to which some reference should be made. This reviewer, and others working on interactions between C. elegans and soil bacteria, have often observed N2 swarming at the edge of unpalatable bacterial lawns. These swarms appear to allow the worms to overcome bacterial defences by collective action of digestive enzymes or other secretions.

- Do worms in the experiments move in purely 2D, i.e. never crawl on top of one another and never overlap? Is this true for the agents in the simulations?

- Better justification of the level of complexity chosen to model worms with agent-based approach is necessary. Why would simple "vector particle" models fail? How important is that each worm is multi-segmented? How important are the elastic properties of the worm? Does the assumption of the propulsive force generation at the head vs throughout the body affects the motion?

- By comparison with the bacterial literature, the term swarming used by the authors could be somewhat confusing (e.g. by comparison with the widespread phenomenon referred to as bacterial swarming, e.g. PMID: 20694026). It is also probably acceptable to use the term "swarming" in the worm field, but a clear definition should be added in the Introduction.

Major points:

1) Importantly, the modeling approach could be improved given the wealth of the experimental data. Specifically, it is surprising that the authors chose a fully phenomenological approach to construct their agent-based models that uses very little of the quantified data.

We thank the reviewers for their suggestions, however we believe that a phenomenological approach is more appropriate in our case. We lack a single-worm model that is accurate and simple enough for inference. We have recently shown that we can simulate accurate locomotion using a neural network (https://www.biorxiv.org/content/10.1101/222208v1), but achieving this accuracy requires hundreds of parameters. We feel a more appropriate level of description for modeling collective worm behavior is the persistent random walk, which has previously been shown to work well for single-worm foraging (Salvador et al., 2014 – now cited in Appendix 1 under “self-propelled movement”). At the level of this more tractable effective model, we thoroughly considered previously available single worm data to parametrize the size of our agents, rates of slowing, and the duration of slow movement periods at zero neighbor density. Other data on the frequency of undulations and worm shapes were ultimately not considered, as we did not find them relevant to the aggregation mechanism (see for example Figure 6—figure supplement 2H, where we add undulations).

Furthermore, authors do little to quantitatively compare resulting behaviors of agents in the model to those of tracked worms. Therefore, the reviewers are not fully convinced that identified behaviors in the model are either necessary or sufficient (given considerable noise in individual behaviors) for the observed population behaviors.

To address the question of necessity of the observed behaviors, we have run additional simulations and added Figure 6—figure supplements 2 and 3 in response to the reviewers’ suggestion. Starting from parameter values obtained from the mean of the posterior distribution as a reference, we removed or perturbed individual model behaviors. These simulations show that removing speed switching or taxis from the model disrupts aggregation, while removing reversals reduces the overall quantitative agreement with experimental data (Figure 6—figure supplement 2B-D). In some cases, removing individual model behaviors also produces a different correlation of velocity and orientation between neighbors from what we measure in experiments (Figure 6—figure supplement 3).

To circumvent this issue, measured individual worm behaviors should be used in constructing and/or constraining their model (see for example References PMID: 30514635, PMID: 28533367).

We appreciate the reviewers’ suggestions and have considered this approach. However, we favor a Bayesian inference approach instead of feeding the model at every time-step with observed tracking data (which in our case is imperfect), as other work on bacterial systems has done, including the suggested references. The inference approach not only gives us the most likely model parameters and their uncertainty, but also has several additional advantages, such as providing information on model identifiability or “sloppiness”, automatically adjusting for model complexity when sampling in higher dimensional parameter spaces, and being easily extensible to incorporate multiple sources of experimental data. In response to your suggestions, we now discuss the data-driven approach as an alternative (Discussion, fourth paragraph).

Importantly, the noise and variability of the movement of agent match those of worms. It is not clear if this is the case in the current model.

Noise and variability of movement have indeed already been considered when parameterizing the angular noise of movement relative to the size of food patch. In our experiments, long-term persistence of movement is disrupted when worms reach the edge of the food lawn. Thus the noise parameter was chosen to ensure an almost complete reorientation (velocity autocorrelation of 0.23, corresponding to choosing a random angle between ± 3/4π) after a time it takes a single worm to move across a food lawn (about 25 seconds for npr-1 mutants moving at 330 μm/s). This is described in Appendix 1 under “self-propelled movement” and the velocity autocorrelation for our simulation is in Figure 6—figure supplement 4A3. To further compare whether the noise and variability of movement matches that of worms in our experiments, we analyzed the velocity autocorrelation in our worm tracking data. These show similar decay in velocity autocorrelation for individually tracked npr-1 and N2 worms (Figure 6—figure supplement 4A1), and a stronger decay for npr-1 worms in 40-worm experiments (Figure 6—figure supplement 4A2). The latter are thus interacting and expected to have less persistent movement as a result. We have further verified that the model behavior is robust to changing the noise parameter to as low as 0 and as high as 0.08 (which roughly corresponds to the autocorrelation measured for interacting worms in our experiments (Figure 6—figure supplement 4A2); higher noise values not tested), as can be seen in Figure 6—figure supplement 2E-G. These results show that aggregation in our model is robust to noise.

Furthermore, comparison of the resulting behaviors of agents and worms would allow testing the model. For example, is there correlation (e.g. between the agent velocity and the vector to nearby worm) that can quantify the degree of taxis? Can this correlation be compared between the model and experiments?

We have calculated the correlation of agent velocity and the vector to nearby worms to quantify the degree of taxis (Figure 6—figure supplement 3B). In our experiments, this correlation is nearly zero for all distances up to 2 mm, which is larger than the size of a typical worm cluster. This may not be intuitive, but the correlation was comparably low when we quantified the same statistic in our simulations, which have a low but non-zero taxis parameter. We suspect the reason is twofold: a) the taxis effect is only a small influence on the instantaneous direction of movement of a worm, compared to persistence and noise; and b) we only tracked the pharynx in our experiments, and reproduced this restriction in our analysis of simulations, but the whole body of the worm is likely giving relevant cues to any chemical or mechanical taxis. This brings with it the caveat that there is unseen worm density that affects any potential taxis behavior, but which remains undetectable in our tracking. Thus, our analysis shows that the attractive taxis is important for aggregation and may be non-zero (as it is in the parameterized simulations), even if it is difficult to detect with correlation analysis. Our modeling approach thus clearly predicts taxis, or a phenomenologically equivalent behavior, to be an important factor in aggregation, while our analysis also shows that this would have likely been missed by a purely data-driven approach.

Similar quantitative characteristics of other proposed 'key behaviors" need to be compared.

We now compare our inferred parameters for reversal and speed switching rates with our measured data (Figure 6—figure supplement 4B-C). The agreement of the reversal rate is imperfect, which may be in part due to the sensitivity of experimental analysis. When quantifying reversals, we had to make a choice as to when a reversal is substantial enough, for which we imposed a minimum absolute speed and a total path length of the reversal of 50 μm (roughly one worm width). Shorter reversals may more likely be due to noisy tracking or small head-movements. However, many reversals are not tracked from the beginning to the end, especially at higher densities; therefore any threshold on path length may falsely exclude partially tracked reversals.

To compare speed switching rates kf and ks, we unfortunately cannot reliably measure these experimentally. However, we have quantified the ratio of worms moving at fast (up to 350 μm/s) versus slow (<100 μm/s for npr-1, <50 μm/s for N2) speeds in experiments and simulations at various local densities. The disagreement (Figure 6—figure supplement 4C) may indicate that the exponential form of k) (see Figure 5C and main text) is not the best approximation. However, aggregation in the model is not sensitive to speed switching rates, as shown by the broad posterior distributions for the inferred parameters (Figure 6C-D), thus we leave it for future work to consider different functional forms of kf(ρ).

Can changes in worm behaviors in nutrient-rich vs depleted regions be quantified and fit into the model?

We have used the fact that worms speed up when off food to demonstrate that the aggregation mechanism also leads to swarming, and the fact that movement on food is persistent in order to parameterize noise for the model. Otherwise, our model mostly concerns on-food behavior, as 99.7 ± 0.4% (npr-1, mean ± S.D.) and 99.8 ± 0.3% (N2, mean ± S.D.) of the worms are on food during the aggregation phase of the experiments.

Either of the two approaches or some combination of them is needed to show convincingly that the postulated behaviors are indeed observed and responsible for the observed population self-organization.

As described above, we have kept the phenomenological model construction and Bayesian parameter inference approaches, but have now added substantial quantitative comparisons between experimental and theoretical outcomes to better illustrate the utility and limitations of our model.

2) The authors discuss possible neuronal mechanisms underlying their behavioural rules. Ideally, it would be good to apply or at least suggest possible experimental tests of these mechanisms, in particular of the proposed important medium-range taxis between neighbouring worms. It is suggested that this might be O2 tension mediated: how might this be tested? They do provide convincing experiments that exclude any significant contribution of ascaroside pheromones to swarming.

We have now suggested several experiments to test the neuronal mechanisms underlying our core behavioral components (Discussion, third paragraph). To elucidate the role of a shallow oxygen gradient in taxis, we suggest introducing mutations leading to aerobic metabolism deficiencies into npr-1 mutants. Such mutants would still be able to sense ambient oxygen, but are expected to produce an even weaker oxygen gradient in an aggregate. A quantitative comparison of the resulting phenotype to model predictions (e.g. with reduced taxis and/or modified rates of density-dependent reversal and speed switching) would be informative. Additionally, one may perform aerotaxis and oxygen-shift experiments (using a much smaller oxygen gradient than previous works have used) to seek evidence for the ability of worms to sense and move towards environments where oxygen levels are only slightly below ambient concentrations. We also discuss the possibility of perturbing reversals and speed-switching rates using mutations.

Other comments:

- The work also mainly compares behaviour in the standard lab strain N2, which is 'hypo-social', to behaviour in a 'hyper-social' npr-1 mutant of N2. Neither of these exactly corresponds to the natural state of almost all wild C. elegans, so the analysis is somewhat artificial; this might be more explicitly admitted.

We have now added relevant wording to emphasize that npr-1 and N2 worms show extreme aggregation phenotypes (Discussion, final paragraph).

- The authors state that to their knowledge "this swarming phenotype has not been reported in C. elegans previously", which betrays ignorance. For example, the following observation in Hodgkin and Barnes (2002, PMID: 1684664): "When a population has consumed almost all the bacteria, the worms… form a swarming mass that moves as a wave across the remaining bacterial lawn, and then disperses when all the bacteria have been consumed." This indicates that N2 swarms under appropriate conditions, perhaps of higher bacterial food density. The authors' experiments were performed only on "thin, even bacterial lawns".

We have now included the fact that N2 swarms under appropriate conditions in our text (Results, first paragraph). We have also modified our text in appropriate places to clarify that swarming in this paper refers to that of npr-1 mutants.

- They also state that "it is still unclear whether aggregation and swarming have a function in the wild". Swarming is almost certainly naturally advantageous in predation on pathogenic or repellent bacteria, as a result of the 'wolfpack' effect that has been studied in myxobacterial swarms, to which some reference should be made. This reviewer, and others working on interactions between C. elegans and soil bacteria, have often observed N2 swarming at the edge of unpalatable bacterial lawns. These swarms appear to allow the worms to overcome bacterial defences by collective action of digestive enzymes or other secretions.

We thank the reviewers for sharing these insights and have now included them (Discussion, sixth paragraph). While swarming may be likely to have an ecological role, without a better understanding of the fitness effects and indeed whether it occurs in nature, we would prefer to remain cautious about the function of aggregation and swarming.

- Do worms in the experiments move in purely 2D, i.e. never crawl on top of one another and never overlap? Is this true for the agents in the simulations?

In our experiments, worms move in 2D on agar and rarely burrow. When multiple worms are interacting, they can be observed to crawl over each other, whilst also restricting each other’s movement. This phenomenology does not perfectly correspond to full volume exclusion or free overlap, and thus is difficult to implement in simulations. We chose not to enforce volume exclusion as we suspect the degree of overlap between worms in aggregates to be considerable. Thus, the worms can also overlap in the model, although with taxis direct overlap is discouraged – the direction of taxis reverses (becomes repulsive) if worms overlap, which was chosen to reflect some extent of physical hindrance. We had originally also tried simulations with volume exclusion, but had chosen not to show these as it was not required to achieve aggregation and the aggregates seemed more spread out than in experiments. We now show an example of this in Figure 6—figure supplement 5B.

- Better justification of the level of complexity chosen to model worms with agent-based approach is necessary. Why would simple "vector particle" models fail? How important is that each worm is multi-segmented? How important are the elastic properties of the worm? Does the assumption of the propulsive force generation at the head vs throughout the body affects the motion?

Simple vector-particle models would not have a front and back to detect when particles are exiting a cluster, which seems to be an important behavioral trigger. Though this could be fixed by including a memory, a particle leaving an aggregate ceases to be an interaction partner for particles in the aggregate unless there are long-range interactions. We chose agents with a long shape as we were initially looking for purely local (i.e. short range) interactions that lead to clustering. However, we believe this is still close in simplicity to the “vector particle” model, as our simulated worms really are just “vector particles with a tail”, or equivalently, a chain of vector particles.

The multi-segmented representation of the worms in our model may not be essential, and we suspect that equivalent phenomenology could be captured with rigid rods and different inferred parameters, as long as the rods are allowed to overlap. However, rigid rods have the undesirable property that the tail swings out upon turning, which is not observed in worms. Hence, we chose a persistent random walker with a tail as a more plausible representation in our current approach.

In the model presented here, elastic effects play no importance apart from enforcing non-extensibility of the chain-like worms.

The propulsive force is not generated only at the head; rather, each node of the chain is a self-propelled particle. The head node undergoes a persistent random walk, whereas the other nodes follow the line of the body (and all nodes are subject to non-extensibility constraints). These are described in Figure 5A legend and Appendix 1 under “SPP worm model”.

To investigate whether our model still works with simplified agents, we ran further simulations (starting from the posterior mean parameterization for npr-1) with fewer nodes per worm (and correspondingly shorter total worm length, because the size of the nodes was not changed to keep the width the same). We found that aggregation was robust to shortening worms up to a point (Figure 6—figure supplement 5A). With 4 nodes per worm aggregation only formed transiently and with fewer worms per aggregate, indicating that the worms need to be longer than this threshold length. This can be explained by the fact that at this worm length, the interaction radii (set to 3 node radii) of the head and tail nodes start to overlap, and in our simulations worms require a difference in contact between head and tail to initiate reversals. Thus, as worms become shorter than M=5 nodes, the head and tail nodes increasingly have the same degree of contact within a cluster, and thus more often fail to sense when they are exiting a cluster (while also having a shorter time window in which to initiate a reversal whilst leaving a cluster).

- By comparison with the bacterial literature, the term swarming used by the authors could be somewhat confusing (e.g. by comparison with the widespread phenomenon referred to as bacterial swarming, e.g. PMID: 20694026). It is also probably acceptable to use the term "swarming" in the worm field, but a clear definition should be added in the Introduction.

We have added an operational definition for C. elegans swarming (Results, first paragraph).

Key resources table

Resource	Designation	Source or reference	Identifiers	Additional information
Strain (C. elegans)	N2	Caenorhabditis Genetics Centre	RRID:WB-STRAIN:N2	Laboratory reference strain.
Strain (C. elegans)	DA609	Caenorhabditis Genetics Centre	RRID:WB-STRAIN:DA609	Genotype: npr-1(ad609)X.
Strain (C. elegans)	OMG2	this paper		Genotype: mIs12[myo-2p::GFP]II; npr-1(ad609)X. Originated from CB5584 and DA609.
Strain (C. elegans)	OMG10	this paper		Genotype: mIs12[myo-2p::GFP]II. Originated from CB5584; outcrossed 6x to CGC N2.
Strain (C. elegans)	OMG19	this paper		Genotype: rmIs349[myo3p::RFP]; npr-1(ad609)X. Originated from AM1065 and DA609.
Strain (C. elegans)	OMG24	this paper		Genotype: rmIs349[myo3p::RFP]. Originated from AM1065; outcrossed 6x to CGC N2.
Strain (C. elegans)	DR476	Caenorhabditis Genetics Centre	RRID:WB-STRAIN:DR476	Genotype: daf-22(m130)II.
Strain (C. elegans)	AX994	Mario de Bono (MRC Laboratory of Molecular Biology)		Genotype: daf-22(m130)II; npr-1(ad609)X.
Software	Tierpsy Tracker (v 1.3)	Javer et al., 2018		Software available at ver228. github.io/tierpsy-tracker.
Software	wormTrackingAnalysis	this paper		Software available at github.com/ljschumacher/wormTrackingAnalysis.
Software	sworm-model	this paper		Software available at github.com/ljschumacher/sworm-model.