Computer vision hardware
A Dalsa Falcon 4M30 camera (8 bits; 2352 × 1728 pixels, 31 Hz) was used with a Rodenstock 60 mm f/4.0 Rodagon lens to image a 5 cm plate of worms (resulting in an image with 24.3 μm per pixel) under brightfield illumination (backlit with 4″ × 4.9″ light plate, Schott A08925 with ACE I illuminator). All tracking was performed on computers with 3 GHz Intel Core 2 Duo processors and 4 GB of RAM. Images were captured using National Instruments’ PCIe-1427 CameraLink capture card.
Online image processing
Custom software written in LabView (National Instruments) presents an interface to set up image processing and experimental parameters; once tracking is started, the software captures frames from the camera via CameraLink. (The software also supports GigE, Firewire, and USB input, and can read saved video from AVIs.) Captured frames are passed, on-line, to a custom image analysis library written in C++. Collectively, this software is called the MWT (Multi-Worm Tracker) and is available under an open-source license. The version used in this study is available as Supplementary Software 1
and the most recent version can be found at http://sourceforge.net/projects/mwt
. LabView is required to modify the MWT code but not to run it. However, a runtime license for the LabView Vision package is required in all cases.
To identify moving objects, the MWT searches for portions of the image that are darker (optionally, lighter) than background by some threshold T. Objects are segmented using flood-fill from dark pixels. To avoid segmenting single pixels due to camera noise, T is set stringently but objects are filled out to a less stringent threshold βTT (typical βT ≈ 0.8).
Objects that fall within user-specified size thresholds (with more stringent thresholds S−, S+ for initial capture and less stringent ones βSS− and βT−1S+ thereafter) are presumed to be worms, and are followed by rectangular subregions 10 pixels (by default) larger on each side than the bounding box of the object. In the next frame, the object is sought only in this subregion. Thus, only a small portion of each frame needs to be analyzed. In addition, new objects are sought in a moving band covering 1/16 (by default) of the full image.
The background estimate is updated with a decaying average: B(t + 1) = (1 − α)B(t) + αI(t) where B is the estimate of the background, I is the measured intensity of a certain pixel at a given time, and by default α = 2−5. This estimate is only used to update pixels within the band used to search for new objects, and it only updates those pixels that are not in any animal’s subregion. Permanently motionless animals are treated as part of the background.
Objects are considered to persist from one frame to the next if four criteria are satisfied: (1) the new object is within the (less stringent) size threshold; (2) the new object overlaps an old object with at least 50% of its pixels; (3) only one old object overlaps the new object; and (4) only one new object overlaps that old object. Objects failing the first two criteria are discarded. Objects failing the third criteria (plus several others) are considered “collided” and are tracked as a new joint object; objects failing (4) are considered “split” and are tracked as two new objects. Due to rules (3) and (4), collisions can be robustly detected and isolated from analysis; when animals separate they can be immediately tracked again (original identities are lost, as robust maintenance of object identity is too computationally intensive to solve in real time).
The centroid position, area, best-fit line (least squares), and exterior contour of each object is computed from the binary segmentation. A spine (branchless skeleton) is computed using a fast algorithm that assumes that the objects have two pointed ends. One end isthe sharpest angle:
is the number of points in the contour, c
is a vector between contour points, ci,k
, and xi
is the ith
contour point. The other end is the sharpest weighted by distance from the first:
is the index of the selected end and j
is the index being tested. Eleven spine points are then placed at the midpoint of equal divisions along the outline from the endpoints.
All quantities are computed at each timepoint for each object and accumulated in memory. When an object is lost, all of its data is written to disk. On-line estimates of speed are computed and can be displayed during tracking as a diagnostic tool.
The image processing and analysis are single-threaded (although largely parallelizable in principle), and are not strictly real-time. In particular, if a heavy workload or operating system call causes enough delay so that several images arrive while the previous one is being processed, only the latest image is used while older ones are skipped.
By default, only the first image of an experiment is saved to disk as a reference, since disk I/O would be a limiting factor unless specialized hardware was used. Timelapse images can optionally be saved to disk; on our system, aframe every two or three seconds is tolerated well.
To refine the data saved by the MWT and to quantify specific behavioral parameters, we created a second software package, Choreography. Choreography is written in Java and is available under an open source license as part of the MWT package at http://sourceforge.net/projects/mwt
A variety of conditions can be set to select which animals are suitable for analysis; the primary selection criteria are a minimum duration of tracking (typically 20 seconds) and a minimum movement (here, three times the animal’s length).
Choreography computes a variety of metrics and can provide population-averaged or animal-by-animal output, as appropriate. Metrics include area covered by the animal, bearing of the animal as defined by a least-squares line fit to the segmented pixels, speed and angular speed computed using boxcar averaging with user-defined width, length measured along the spine, average curvature along the spine, changes in direction, and—assuming that the animal’s dominant movement direction is forward—direction of motion. Note that the resolution used does not permit direction of motion to be inferred from shape analysis; head and tail are largely indistinguishable.
For noise estimates, we assumed that high-order derivatives of behavioral parameters should tend towards zero and therefore any residual is caused by measurement noise. Specifically, we assumed Gaussian noise and thus
→ ∞. In most cases, the estimate converges by n
= 2 (second derivative). The estimate is conservative for stochastic camera noise alone since behavior probably contributes somewhat to high-order derivatives. However, it is imperfect: it does not take into account rare large changes due to aliasing along pixel boundaries.
Metrics that rely upon the spine, such as length and curvature, are typically computed by first low-pass filtering the outline using an exponentially decaying kernel (2−1/2 decrease per pixel) before applying a spine-finding algorithm similar to that used in the MWT. If the plate is particularly messy and the outlines are therefore not smooth, a heuristic is applied to the outline to detect T-shapes and remove the stem of the T (presumed to be a spurious projection).
Metrics that rely upon the geometry of the animal’s path are not computed from raw centroid positions. Instead, we perform a greedy geometric decomposition of the path defined by the centroid positions. First, we classify those portions statistically indistinguishable from random jitter. Since our noise-measurement routines do not test for normality, we view the tests not as a guarantee of statistical relevance, but rather as a principled heuristic to obtain a decomposition of the animal’s path that is sufficient for robust behavioral analysis. Remaining points are fit by straight lines that are extended until new points are statistically distinguishable. Straight lines are converted to circular arcs where doing so provides statistically significant improvement (in mean squared displacement of points from the line). To speed computation, circular arcs are fit using the “HyperFit” method with zero essential bias25
. This method requires finding the eigenvalues of a 4 × 4 fitting matrix, which we compute byusing an analytic solution of the quartic equation.
A change of direction is considered to be either two adjacent line segments that meet with a negative dot product, or motion along a single line that backtracks a statistically significant distance (P < 0.05 over the entire track, assuming Gaussian noise). Movement direction is scored by first splitting the path into segments where ends can be unambiguously followed---by rejecting segments whose perimeter to area ratio is statistically consistent with an animal folded back on itself along at least 1/3 of the body---and then scoring direction changes and assuming that any sufficiently long path (by default, three times body length) is predominantly forward.
Averages of computed metrics can be saved to disk or viewed as graphs within Choreography. In addition, metrics can be saved on an animal-by-animal basis or can be viewed as a color-coding of the tracks of the animal on a 2D map. Along with the tracks, the animals’ body postures can be replayed at various zoom levels and speeds to allow manual validation of the results. The replay feature also allows the creation of multiple minimaps that provide an overview of different regions on the plate; by switching between these, one can rapidly move the main viewback and forth between different locations and scales.
Organism- and experiment-specific functionality is implemented as plugins for Choreography. Plugins are tightly integrated in that they have access to all internal data used by the software, and can intervene at various stages of analysis. Although Choreography is open source and therefore can be modified directly, users are encouraged to write plugins when extending the software to reduce the difficulty of moving to new versions.
All assays were conducted with the wild-type Bristol isolate of Caenorhabditis elegans (N2) unless otherwise stated. The mutant strains CB1033 (che-2), PR811 (osm-6), RB2464 (tax-2), and CB1611 (mec-4), as well as the strains of wild-type Caenorhabditis species N2 (C. elegans, Bristol isolate) HK104 (C. briggsae), PB4641 (C. remanei), and PB2801 (C. brenneri), were obtained from the Caenorhabditis Genomic Center (CGC). These strains were maintained unfrozen for less than 30 generations before testing. The strain XJ1 is descended from N2 but was maintained unfrozen for at least several years (hundreds of generations). The 33 strains tested for tap habituation defects () were also obtained from the CGC.
Baseline recordings on food were performed on 5 cm NGM plates seeded 24 hours earlier with 50 μl of OP50 bacteria. Worms were transferred to the plate at 80–90 hours of age and allowed to acclimate for 4 hours before being tracked for 60 minutes. Unless otherwise specified, all analysis was performed with the Choreography options --shadowless -M 3 -t 20 -S --plugin Reoutline::exp --plugin Respine (see MTW User Guide in Supplementary Software 1
or at http://sourceforge.net/projects/mwt
To assess whether plate-to-plate differences in movement rates were a result of the statistics of behavior of individual worms, we pooled all worms from eight plates and randomly sampled (with replacement) to generate data for 1,000 virtual plates. The observed differences between real plates and their mean was 60% larger than the difference between virtual plates and their mean (20.7 μm s−1 vs. 12.9 μm s−1); five of eight real plates had difference scores that were statistical outliers (P < 0.01) of the virtual plates’ distribution. Thus, the observed variability between plates was not solely due to statistical sampling. Other likely causes include different environmental conditions on the plates, subtle perturbations during recording, or behavioral states that last longer than identity is maintained (205 s on average here). Fortunately, plate-to-plate variability does not obscure the advantage of multi-worm tracking: monte carlo generation of virtual plates with 10 tracked worms produced a data set with slightly larger plate-to-plate variability than was observed for real plates.
Omega turns were detected with a combination of three factors. A spine-length-to-perimeter metric flagged regions where the worm had folded back to touch itself. To detect tight turns without self-contact, we used both the third principal component of the spine shape (which correlates with turns), and a folding score
is the vector from spine point a
. Intuitively, this last metric is small when the worm’s spine is folded exactly in half (all dot products will be close to −1); constants are chosen so that one gets a minimum of 0.5 for each pair of spine points folded onto each other. Thresholds for these metrics were set by hand; a wide range of parameters gave results that agreed with human judgment in obvious cases. This method detected 74 of 76 hand-annotated omega turns in our reference data set (3% false negatives), and 98 of 100 automatically-scored omega turns were confirmed by hand in wild-type movement data on relatively clean plates (2% false positives, though the error rate is higher on a sub-optimal plate such as our reference data set). Thus, the overall error rate is ~5%.
Food chemotaxis plates were prepared by placing a 50 μl drop of OP50 bacteria at x,y coordinates (r/2,0) where (0,0) represents the center of a plate of radius r. In practice, on the 50 mm inner diameter plates we used, the actual coordinates (measured by hand) were (12.4,0.2) mm with standard deviation (±0.5,±0.9); the diameter of the spot of food was 3.7 ± 0.3 mm. A control spot of LB growth medium was placed at (−r/2,0). The OP50 was grown for 48 hours. For testing, 6–12 worms were picked from food to an empty plate to reduce the amount of stray food, and then picked again from the empty plate to coordinate (0,0). Tracking was started within 60 s of placement and left to run for 20 minutes.
To quantify food chemotaxis, we counted crossings of the animals’ centroids along a ring 1 mm outside of the spot of food or the control spot. We computed Bayesian preference estimates P
) and the associated standard error σP
) for each spot size, based on a binomial distribution with uniform prior26
) is the cumulative number of (signed) crossings into the food spot at time t
) is the same for the control spot.
To quantify pirouettes, we monitored the animal at 0.5 s intervals and calculated the angle between the direction of movement and the direction to the food. We then pooled the quadrant heading towards the food (angles from −π/4 to π/4) and the quadrant heading away, and counted reversal frequency (probability per unit time of executing a reversal). We compared the behavior of worms in a 1 cm radius circle about the center of the food spot to those in a 1 cm radius circle about the center of the control spot. Within-strain statistical comparisons were performed using a χ2 test; P-values for hypotheses about the similarity of strains were computed using monte carlo sampling.
Salt chemotaxis plates were prepared following previous work9
. In brief, we placed a 5 μl drop of 500 mM NaCl at the center of a low salt plate (1 mM CaCl2
, 1 mM MgSO4
, 5 mM potassium phosphate, pH 6, 2% agar) three hours before tracking. Worms were washed 30 minutes before tracking in low salt buffer (as plate, less agar) and rinsed in low salt buffer plus 0.05% glycerol. 8–12 worms were transferred from the rinse bath to a < 0.5 μl droplet held in a glass loop, which was then dotted onto the plate 12.5 mm from the center. The droplet dried within 60 s; tracking was started within 120 s. To quantify weathervane turning, we segmented the path into one-movement-cycle chunks (two body bends), calculated the angle between the animal’s movement and the direction to the food (“bearing”), and how that angle changed over one cycle (“curve”). We considered only animals between 2.5 mm and 12.5 mm from the location of the NaCl drop, as this was predicted to be where the gradient was steepest (using the method of 8
) and provided some robustness to inaccurate placement of the drop.
Tap habituation experiments
Tap habituation assays were performed on food after a minimum of 6 hours’ recovery from transfer; animals were tracked for 10 minutes to allow them to approach steady-state behavior, and then were tapped 20 or 30 times ( respectively) at a 10 s inter-stimulus interval with a custom-built solenoid tapper that drives a small plunger into the side of the plate. A reversal was scored when a worm was still or moving forward at the time of the tap and moved backwards within one second of the impact; the reversal was considered to be complete when the worm began a pause or forward motion lasting for more than 0.2 s. Cases where the animal was already reversing were removed from analysis. For reversal probability computations, reversal distances of less than 30 μm were considered “no reversal”, and we converted counts of reversals and failures into probabilities using a Bayesian estimate of probability with uniform prior.
Validation was performed by comparing with manually annotated reversals. Out of 148 responses scored by hand as either yes, no, or already-reversing, the tracker disagreed with eight (5% error rate). All but one discrepancy was due to small disagreements over the timing or size of a response.
Putative tap habituation mutants were prepared using a higher-throughput method: five gravid adults were placed on a plate seeded with 50 μl E. coli 24 hours before. The adults were left to lay eggs and then removed from the plate 3 hours later, leaving approximately 60–80 eggs. The plate was tracked once the animals reached adulthood (80 ± 3 hours later, at 20 °C). The recording protocol consisted of 100 s of baseline followed by 30 taps applied to the side of the plate at a 10 s inter-stimulus interval. The data was analyzed with Choreography’s reversal detection plugin: --plugin MeasureReversal::tap::dt=1::collect=0.5. For each mutant strain and wild-type replicate, three or four plates were tested (except CB1220: n = 2); on average, 34 worms per plate were suitable for analysis. A wild-type distribution for all wild-type plates tested on the same day was created for the initial and habituated (28th–30th) responses. All strains and wild-type replicates were standardized to the distribution for the day they were tested. Monte Carlo sampling was used to generate the null hypothesis distribution given this scheme, and the resulting P-values were re-interpreted as an effective Z-score for plotting. Values outside ±3.13 were considered to be significantly different than the wild-type distribution (two-tailed, P < 0.05, with Bonferroni correction for multiple comparisons). CX20 and NM1815 were confirmed using an unpaired two-tailed t-test between the raw reversal probabilities of the last stimulus versus the wild-type replicate that was tested most closely chronologically.
The reference data set used for most hand-annotation and for speed and accuracy tests was a tap habituation screen assay with two differences. First, the worms were given only 5 minutes to recover before tapping instead of 10. Second, the worms were allowed to grow up for 100 hours. These conditions provided an extra challenge for the tracking algorithms, as there were a number of L1 progeny on the plate by that time. To improve the outlines in this more-difficult case, we used the::despike option of the Reoutline plugin.
Since the MWT cannot save video at the full rate (7.2 GB/minute), we used a custom LabView program to stream video to a solid state drive so we had a complete record of the experiment. We then ranthe MWT on the saved video.