Caenorhabditis elegans Multi-Tracker Based on a Modified Skeleton Algorithm

Automatic tracking of Caenorhabditis elegans (C. egans) in standard Petri dishes is challenging due to high-resolution image requirements when fully monitoring a Petri dish, but mainly due to potential losses of individual worm identity caused by aggregation of worms, overlaps and body contact. To date, trackers only automate tests for individual worm behaviors, canceling data when body contact occurs. However, essays automating contact behaviors still require solutions to this problem. In this work, we propose a solution to this difficulty using computer vision techniques. On the one hand, a skeletonization method is applied to extract skeletons in overlap and contact situations. On the other hand, new optimization methods are proposed to solve the identity problem during these situations. Experiments were performed with 70 tracks and 3779 poses (skeletons) of C. elegans. Several cost functions with different criteria have been evaluated, and the best results gave an accuracy of 99.42% in overlapping with other worms and noise on the plate using the modified skeleton algorithm and 98.73% precision using the classical skeleton algorithm.


Introduction
The nematode Caenorhabditis elegans (C. elegans) is a widely studied animal model, as its diverse age-related behavioral patterns provide valuable information on the function of its nervous system and is, therefore, an attractive model to evaluate the effects of mutations [1]. This facilitates the study and treatment of aging, as well as age-related pathologies and neurodegenerative disorders in humans at advanced ages [2,3]. Many of these studies have shown that automatic tracking applications based on computer vision systems help to reduce the manual cost of data acquisition and research hours, improving potential observation of the effects of drug trials [4] and improvements in lifespan or "shelflife" [5][6][7]. These systems provide quantitative information on alterations in the individual motility and behavior of worms produced by chemical substances in their environment (chemotaxis), providing statistical data that allow further research in the field of health and wellness.
C. elegans demonstrate group behavior [8][9][10], among the best known are courtship, mating, aggression, rearing and foraging. Group behavior assays are currently being performed, for example, research into the effect of 0 2 in food search analysis, and aggregation [11,12]. These assays, like others, are visualized manually, due to the complexity of solving the identification problem during an overlapping or body contact of these worms. Currently, the automatic [13][14][15][16][17][18][19][20][21] or semi-automatic applications discard the data from tracks where there are these particular cases (overlapping and bodies contacts) [22][23][24]. Overlapping can take place among worms or may also be due to plate noise. Plate noise was defined as segmentation errors due to edges, or opaque waste in the plate.
Certain applications use diverse techniques and methods such as length [25], smoothness [26,27], previous segmentation [28], or other complex methods [29][30][31][32][33][34] to solve track-ing problems during aggregation. The aforementioned methods were tested by our team and results indicated that previous segmentation is one of the most significant criteria to help identify worms within an aggregation. This was taken as a starting point to design different trackers, adding and combining numerous criteria, and it was found that the criteria of completeness and color used with the improvement of the form of skeletonizing [35] can help to identity the worms in the next image. Furthermore, we introduced three new criteria: length, smoothness, and noise, revealing that they can be useful in re-identification.
We present different multi-trackers which, unlike any other trackers, are fully automatic (Table 1). Segmentation and identification of the edge, the interior rim of the Petri dish (area of interest), and segmentation tracks made by worms are performed without the intervention of human operators. These new methods show worm-tracking accuracy of above 98% in nematode aggregations and problems related to plate noise. The best result reaches 99.42% accuracy. Table 1. Comparison with other multi-trackers. This table shows the comparison of our multi-tracker with respect to others (tierpsy-tracker [22], WF-NTP.v3 [23]).  [36,37] Fully automatic Improved skeleton and segmentations Individual tracking, overlaps, body contacts, rolled worms and occlusions.

Nematode Strains and Culture
The strains (N2) and CB1370 daf-2 (e1370) used in the study were provided by the University of Minnesota Caenorhabditis Genetics Center. The C. elegans were kept at 20 • C and cultured on 55 mm diameter NGM plates with 1 µg/mL of fungizone to prevent fungal contamination [38]. Escherichia coli (E. coli) strain OP50 was used as standard prey for C. elegans in the laboratory. FUdR (0.2 mM) was used to prevent replication. The cultured plates were 10, 15, 30, 60 and 100 adult-young worms synchronized from worm eggs incubated at 20 • C. This variety in the number of worms per plate allowed to obtain aggregations of two C. elegans or more.

Proposed Tracking Method
At present, automatic or semi-automatic trackers discard the data of tracks where there are overlap or body contacts, due to the difficulty of solving the identity of each individual in these situations. Our method does solve these situations using a fully automated pipeline based on an intelligent active backlight, a modified skeletonization method and new criteria to solve the final optimization problem.
Steps of the proposed tracking method are described in Figure 1a-f, and more detailed in Sections 2.3-2.7, respectively. This method starts with the acquisition of sequences of C. elegans images, Figure 1a. These image sequences are then processed (Figure 1b) to obtain mathematical models of C. elegans Figure 1c, and possible solutions to these in an aggregation. Finally, the result (Figure 1f) is the optimal minimization cost between models and possible solutions. Figure 1. General scheme of image processing. The image shows the different stages that the tracker goes through to obtain the results of the skeletons and to follow all the worms within the plate. (a) Image acquisition [37]. (b) Pre-processing image. (c) Worm models. (d) Improved skeleton using proposed method [35]. (e) Optimization method. (f) Optimization results (predicted poses).

Image Acquisition
The image acquisition process was performed using the capture system [37]. To use this system, first, a laboratory operator removed the plates from the incubator, then the plates were analyzed to find condensation on the covers, if so, they were removed, otherwise, the image sequence was captured. The system [37] is automatic and captures an image every second. Escherichia coli (E. coli) strain OP50 was placed in the center of the plate to capture the worms inside the Petri dish and not scaling the edges or near them. The intelligent active backlighting method was used as the illumination technique [36]. This method is more robust than standard backlighting methods, as it allows us to obtain constant intensity values for the bottom of the Petri dish and the worms (greater than 48 and less than 35, respectively). This facilitated automatic segmentation with fixed thresholds on all images.
The image sequence was acquired at a resolution of 1944 × 1944 pixels and a frequency of 1 Hz (1 image per second) using the system [37] (Table 2). This image acquisition system is open hardware and its assembly procedure, parts and description are described in detail in another work [37]. Worm tracking, using these image acquisition conditions, is a very challenging problem. An image resolution of 1944 × 1944 pixels is the lowest able to detect worms, when a complete Petri dish (55 mm. of diameter) is monitored with a fixed camera. In addition, this problem was solved by using a low frame rate of 1 Hz. The dataset collected was composed by sequences of 30 images where contact events between worms occurred, to perform all the experiments.

Image Processing
Image processing began with the segmentation of the region of interest and C. elegans tracks in the image sequence (white circle and red track in Figure 1b). To obtain the region of interest, first, a segmentation was carried out on all the images of the sequence using a threshold with a fixed intensity value (35 in the gray scale). Then, an AND operation was executed between all the segmented images, the result obtained went through a "Fillhole" operation to fill small holes. The region of interest was selected as the largest connected component of the resulting image. It is important to mention that the use of a fixed threshold for all images is due to the intelligent active backlighting method as a lighting technique [36], as this system allows to conserve background intensity values constantly.
In parallel to the previous step, the C. elegans tracks were segmented with a process using different threshold levels Figure 1b. Threshold levels were below 35 on the gray scale. The resulting segmentations went through two filters in order to eliminate those tracks that did not correspond to worms. In the first stage, those tracks with an area smaller than the minimum area of a worm were filtered. The second filter analyzed the skeleton of each image, if the skeleton did not correspond to a minimum expected length, it was classified as noise. The results of the number of skeletons found in each image were stored in a 30 × N matrix, where N is total number of tracks and 30 is the number of images in the sequence.
Due to low resolution of the worms, a scale factor of 3 was applied to increase resolution. This process was applied by obtaining the model to the end of tracking. The model of each worm was obtained by analyzing 30 × N matrix, finding image "k" of 30 images, where worms were further apart from each other, and their ends, head-tail, were also separate. Tracking of each worm started from image k + 1 to image 30 and from image k − 1 to the first. The skeletonization method proposed in previous work [35] was used in each image. This method used distance transformation [39] to obtain possible worm skeletons, and through of an optimization method using different criteria found the best skeleton prediction. Once the tracking process had finished, the results were reconverted to the original scale to be saved.

Worm Model
At present, there are different skeletonization methods [39,40]. Matlab's bwmorph function was used as a classic skeletonization method to obtain worm skeleton model (color pixels in Figure 2b). The proposed worm model consists of width and color values along the skeleton of each individual, Figure 2c. The width values are obtained using classical skeleton in resulting image after using distance transformation function on segmentation of the image "k". Grayscale image is shown in Figure 2a, while color values of pixels are obtained using the classical skeleton in gray image in Figure 2b. The length value is the total number of pixels in the skeleton. The length model is averaged while the C. elegans are separate and tracking progresses.

Extraction of Possible Solutions
The skeletonization method proposed in the previous work [35], unlike classical methods, enables the separation of aggregated worms (Figure 3a), creating new paths in the skeleton, and some possible solutions for each worm (Figure 3b-e). Maximum and minimum values of the width vector are used in the distance to transform images of each segmentation to find this new skeleton as mentioned in [35]. The possible solution skeletons are obtained from a recursive function, which runs through ends and branch points of the new skeleton that overlap the previous segmentation of the worm's body (red circle in Figure 3b-e).

Optimization Method
The prediction "S" of the following postures will be the possible solution with the minimum value within all the possible "P" combinations of skeletons in one segmentation (1). The value of each possible combination "C P " of skeletons is obtained from the sum of criteria "m", for the number of worms "n" in an aggregation (2). The criterion "C p " is evaluated for each possible worm "i" and for each criterion "j". The criteria analyzed were the length of skeleton, overlap with the previous body, completeness, smoothness of the skeleton, noise in segmentation and the colors of each worm.
The length and color criteria prevent the prediction differing from the model in length and color. The smoothness criterion prevents sudden changes in the direction of the skeletons. The completeness criterion prevents the current segmentation from being incomplete. The overlap criterion with the previous body prevents the identity change during aggregation. Furthermore, the noise criterion prevents the skeleton prediction falling on the plate noise.
The reconstruction of the body of each worm was used for the evaluation of the different criteria. This was performed by using the skeleton pixels in each possible prediction with the width and color values obtained in the model (prediction start), respectively, for each individual.

Length Criterion
The length criterion "C L " is obtained from the sum of the multiplication of average squared width [W 2 i ] with the difference in length (∆ i ), (3). This difference is obtained from subtraction between the model length of each worm (L i ) and length of the skeleton obtained (W L i ). The average squared of the width was used so that the resulting length criterion is as significant as the rest of the criteria where the error in pixels was evaluated.

Overlap Criterion
The overlap criterion "C O " is obtained from the sum of the absolute difference of the reconstruction of the worm's body in the previous state, B Ppx , (white dashed line in Figure 5b,c,e,f) and the current state, B Cpx , (green and blue segmentation in Figure 5b,c,e,f) for each possible worm "i" in the aggregation (4) and (5). This is done for all "m" pixels in each reconstruction. Figure 5a shows in gray scale the aggregation of two worms, Figure 5d,g shows the next postures predictions (blue and green segmentation), parallel to these in dashed lines showing the previous state (B Ppx ) and the rest show the current state (B Cpx ). In Figure 5d it can be seen that the overlap criterion is low because it belongs to the same worms, while in Figure 5g it is higher due to the identity change.

Completeness Criterion
The completeness criterion "C Cp " is obtained from the sum of the absolute difference of the current segmentation of the image (A S ) and the reconstruction of each possible worm body "i" in the current state (B Cpx ) (6) and (7). This is done for all "m" pixels in each reconstruction or current segmentation. Figure 6a shows the aggregation of two worms in gray scale, Figure 6b-d shows completeness values for each image. Figure 6b is the correct prediction, while Figure 6c shows incorrect prediction, due to the identity change and final extremes. Figure 6d shows low completeness criterion and an incorrect prediction too, due to the change of identity in both worms.

Smoothness Criterion
The smoothness criterion "C S " is obtained from the average of the absolute values of the angles obtained for the "nk" pixels of the skeleton and for each worm "i" in the segmentation (9). For each pixel of the skeleton, there is an angle (θ px ), which is obtained by an average of the sum of "nA" angles before and "nA" angles after divided by the total pixels found "c" (8). The value of "c" will be 2*nA if there are "nA" pixels before and after the pixel to be evaluated (8). Figure 7a shows the aggregation of two worms in gray scale, Figure 7b shows a possible prediction of skeletons with low softness criterion, while Figure 7c,d show possible predictions of skeletons with high softness criterion.

Noise Criterion
The noise criterion "C N " is obtained from the intersection of the noise segmentation (N S ) and the reconstruction of each body of a possible worm "i" in the current state (B Cpx ) (10) and (11). This is done for all "m" pixels in each reconstruction or noise segmentation. Figure 8a shows the aggregation of one worm with noise in a gray scale, while Figure 8b shows the correct prediction (segmentation in blue) with a low noise criterion (segmentation in magenta). Figure 8c shows an incorrect prediction (blue-magenta segmentation) with a high noise index (magenta segmentation).

Color Criterion
The color criterion "C Cl " is obtained by adding all the pixels with an absolute difference greater than the threshold value (U = 1) of the color model (C M ) with respect to the reconstruction of the current body (C A ) (12) for each worm in the aggregation (13). This means that the error only increases the value by 1 if the absolute difference between a pixel of the reconstruction of the body of the possible worm "i" using the color model (C Mpx ) and the reconstruction of the same worm "i" in the current segmentation (C Apx ) is greater than the threshold value. This is done for all "m" pixels in the reconstruction. For both worm body reconstructions (model and current prediction), the same width and length values are used, in order to compare pixels 1 to 1. Figure 9a shows the aggregation of two worms in gray scale, while Figure 9b,c shows the comparison of each prediction with its respective model.

Evaluation Method
Manually labeled skeletons were used as a reference to compare all the results. These skeletons were obtained using an application designed to select each pixel belonging to the skeleton of each worm one by one in the image sequence. This operation was performed for all 3779 postures of the 70 plates used. The shape of these nematodes was recovered using a disk-shaped dilation operation of radius equal to half the width (approx. 2 pixels) on the skeletons obtained.
The Jaccard coefficient, or intersection over the union (IoU), was used to measure the degree of precision in locating worms (14). As its name indicates, it is obtained by dividing the total area of the intersection by the union of the elements [41]. For the evaluation, the area of the reconstructed bodies of the manually labeled skeletons was used, skeletons using the skeletonization method proposed in [35] and the classical skeletonization method using the Matlab bwmorph command.
The IoU index was expected to be higher because a predicted pose (Figure 10b) is compared to an annotated ground-true pose (Figure 10a), which must overlap (Figure 10c-e). The results for the example below are IoU = 0.9784, 0.5667, and 0.2649, respectively.
Matlab 2018b Machine Learning Toolbox was used to obtain the comparison statistics between prediction models (Appendix A Figures A1a,b-A4a,b) and the two skeletonization methods (Appendix A Figures A5a,b and A6a,b). The Kolmogorov-Smirnov normality test was used for large samples (n > 50) and the Wilcoxon Signed Ranks test.

Results
Experiments were performed with 70 plates. Of these, 54 corresponded to plates with 10 and 15 worms, one plate with 30 worms, four plates with 60 worms, and 11 plates with 100 worms, totaling 65,400 worm poses. All these data were analyzed to obtain contact between worms and noise. As demonstrated in [5,14], a higher number of worms per plate will increase the likelihood of contact between them. Nematodes studied were young-adult wild-type (N2) and CB1370, daf-2 (e1370), as mentioned above. 53 tracks with 3240 poses were used to evaluate aggregation between worms, and 17 tracks with 509 poses for aggregation between worms and noise. The IoU index was used to evaluate the percentage of success in tracking the worms and also to compare both skeletonization methods. The area of worm bodies reconstructed from skeletons obtained manually and skeletons obtained with the two skeletonization methods (new and classical) was used to evaluate the IoU index.
Different prediction models were implemented in order to find the most significant criteria. The name of each model has been coded using letters from the criteria names. "O" for overlap, "L" for length, "Cp" for completeness, "N" for noise, "S" for smoothness, and "Cl" for color. The model with the best results was model 7 (OCpCl) with a 99.42% percentage accuracy in aggregated worm tracks and an IoU value of 0.70 in average.  To measure the percentage accuracy in tracking the worms, non-zero IoU values were used, from the beginning to the end of the tracks. The accuracy of the results of the different prediction models using the 2 skeletonization methods are shown in Table 3. In addition, the average IoU value was obtained for all the prediction models (see Table 4), from the beginning of the aggregation to the end. Then, 790 poses were used in aggregations between worms and 509 poses for worms aggregated with noise, 1299 poses in total.

Comparison with Other Trackers
At present, automatic or semi-automatic trackers discard the data of tracks where there are overlap or body contacts, due to the difficulty of solving the identity of each individual in these situations. Our method does solve these situations, so a direct comparison cannot be made with the results of other trackers.
To compare our method graphically with other trackers (tierpsy-tracker [22], WF-NTP.v3 [23]), labeled data was first shaded in different grays and overlapped predicted data in colors (Figure 12a-c). Errors for each comparison are shown in grayscale. The comparison was made with bodies reconstructed from skeletons (skeleton dilation with disk equal to 2), except for WF-NTP.v3 [23], whose results are skeleton centroid points. The tierpsy-tracker [22] multi-tracker, Figure 12a, did not resolve path 2 where aggregation occurs, while the other tracks (1,3,4) were partially resolved, due to the occlusions that the worms make on themselves. The multi-tacker WF-NTP.v3 [23], Figure 12b, did not solve the identity problem in track 2. Furthermore, like the previous multi-tracker it presented problems in the tracks with occlusions (1, 3, 4). Model 7, on the other hand, had almost zero errors, tracked all worms, resolved the identity in the aggregation of track 2 and the occlusions in the remaining tracks, Figure 12c. Results obtained using tierpsy-tracker [22]. (b) Labeled data (grays) compared with colored lines that connect centroids obtained using WF-NTP.v3 [23]. (c) Results obtained using model 7.

Discussion
Caenorhabditis elegans are aggregated in different ways, such as aggregation of end parts (head or tail), partial aggregation of bodies, and aggregation of parallel bodies, among others. The experiments were conducted with the above mentioned methods with a single criterion and found that the overlapping criterion with the previous prediction is the most significant. This criterion allows part of the previous state to be preserved, helping to solve the next state. A result of 98.49% was obtained using this criterion individually. However, when aggregation gives rise to an overlap in most worms, it is difficult to identify them, even for human observers. To solve this problem, different tests were designed, combining the criteria mentioned in the optimization methods section.
The completeness in cases of aggregation, where there are overlaps between aggregated bodies, partial and total identity changes could be obtained, as observed in the optimization method Figure 6c,d. When using the overlap criterion with completeness criterion, a percentage accuracy of 99.13% was achieved. As shown in Appendix A ( Figures A1a,b and A2a,b) this criterion is statistically significant, helping to solve some cases where the overlap criterion presented problems, Figure 5b.
The color criterion also helped to improve worm prediction, although not so significant (Appendix A Figures A3a,b and A4a,b). The end portion of C. elegan (tail) has a greater number of light pixels (higher gray levels) than the head, where there are fewer. These features are important because before, during and after an aggregation one of the final parts remains visible, and with this color singularity a more accurate prediction can be obtained, helping to solve the identity problem after an aggregation between worms or noise. The use of these three criteria (overlap, completeness, and color) allowed us to obtain a percentage accuracy of 99.42%.
The length, smoothness, and noise criteria were discarded because when they were used individually or together with the overlap criterion, their percentage accuracy decreased. On the one hand, the problem with the smoothness criterion was due to the low resolution of the worms in our images, which provided a poor estimate of this criterion. On the other hand, the problem with the length criterion was due to errors in the length model and the increase or decrease in the length of worms owing to overlaps. The problem with the noise criterion was due to, when the worm is visualized against background noise, many possible solutions were generated.
Finally, it is worth mentioning that the best combination of criteria depends on image quality. In this work, we demonstrated that the combination of the three criteria mentioned above (overlap, completeness, and color) was the best option for automatic tracking of interaction behaviors among C. elegans (contacts or overlapping) with our low-resolution dataset.

Conclusions
This paper presents a method for tracking multiple C. elegans in standard Petri dishes where some worms can come into contact or overlap. This method was evaluated in a difficult scenario using a low-image resolution and a low frame rate. Using an optimizer with the appropriate criteria (overlap, completeness, and color) was shown to solve many worm overlap and contact situations. The accuracy obtained under these conditions was 99.42% and 98.73% using the modified skeleton algorithm and the classical skeleton algorithm, respectively. In addition, the proposed method employs an improved active backlight system and an improved skeletonization algorithm. The active backlight system provides fixed gray levels in all captured images, which allows automatic segmentation using a fixed threshold. The improved skeletonization algorithm uses width information from each worm model to extract skeletons, enabling the tracking of worms moving in parallel (side by side). Our proposal, unlike other trackers that discard worm overlaps and contacts, solves many of these situations, increasing the number of worms tracked in a Petri dish and therefore paving the way to automate new assays where interaction between worms occurs. . The p-value obtained was 8.24E-96 less than the significance value of 0.05, so the null hypothesis was rejected and the alternative hypothesis H1 was accepted (data did not come from normal distribution). Once the alternative hypothesis was accepted, Wilcoxon signed ranks test was used to evaluate both methods. . The p-value obtained was 1.48E-151 less than the significance value of 0.05, so the null hypothesis was rejected and the alternative hypothesis H1 was accepted (data did not come from normal distribution). Once the alternative hypothesis was accepted, Wilcoxon signed ranks test was used to evaluate both methods. Figure A4. Wilcoxon signed rank test. (a) The Wilcoxon signed rank test table shows the difference that exists in two related samples through positive, negative and tie ranges. (b) p-value obtained with Wilcoxon rank test was 0.1054 less than the significance value of 0.05, so it was concluded that there was not a statistically significant difference between both models. Normality test on the difference of methods (New-classical). The p-value obtained was 1.26E-152 less than the significance value of 0.05, so the null hypothesis was rejected and the alternative hypothesis H1 was accepted (data did not come from normal distribution). Once the alternative hypothesis was accepted, Wilcoxon signed ranks test was used to evaluate both methods.