An annotation dataset facilitates automatic annotation of whole-brain activity imaging of C. elegans

Annotation of cell identity is an essential process in neuroscience that allows for comparing neural activities across different animals. In C. elegans, although unique identities have been assigned to all neurons, the number of annotatable neurons in an intact animal is limited in practice and comprehensive methods for cell annotation are required. Here we propose an efficient annotation method that can be integrated with the whole-brain imaging technique. We systematically identified neurons in the head region of 311 adult worms using 35 cell-specific promoters and created a dataset of the expression patterns and the positions of the neurons. The large positional variations illustrated the difficulty of the annotation task. We investigated multiple combinations of cell-specific promoters to tackle this problem. We also developed an automatic annotation method with human interaction functionality that facilitates annotation for whole-brain imaging.


Introduction
Identification of the cell is an essential process in the broad fields of biology including neuroscience and developmental biology. For example, identification of cells where a gene is expressed can often be the first step in analyzing functions and interactions of the gene. Also, the identity information is required for comparing cellular activities across different animals. In order to annotate cell identities in microscopic images, features of the cells such as positions and morphologies are often compared between the samples and a reference atlas.
The nematode C. elegans has a unique property that all cells and their lineages have been identified in this animal (Sulston and Horvitz 1977;Sulston et al. 1983). Additionally the morphology and the connections between all 302 neurons in adult hermaphrodite were also identified by electron microscopy reconstruction (White et al. 1986). Such detailed knowledge opens up unique opportunities in neuroscience at both single-cell and network levels. Recent advances in microscopy techniques also enable whole-brain activity imaging of the worm (Schrödel et al. 2013;Prevedel et al. 2014;Kato et al. 2015;Nichols et al. 2017;Kotera et al. 2016;Tokunaga et al. 2014;Toyoshima et al. 2016;Hirose et al. 2017), even for the free-moving worms (Nguyen et al. 2016(Nguyen et al. , 2017Venkatachalam et al. 2016). The neural activities were obtained at single-cell resolution and the identities of limited numbers of neurons were annotated manually in some of the studies (Schrödel et al. 2013;Prevedel et al. 2014;Kato et al. 2015;Nichols et al. 2017;Kotera et al. 2016;Venkatachalam et al. 2016;Toyoshima et al. 2016). However, there is no systematic and comprehensive method to annotate the neurons in whole-brain activity data (Nguyen et al. 2016).
Annotation of neuronal identities in C. elegans is often performed based on positions of the neurons, especially for larval animals in which the neurons are located at stereotyped positions (Bargmann and Horvitz 1991;Bargmann and Avery 1995). However, for adult animals, the positions of the neurons are highly variable between animals (Nguyen et al. 2016). Several additional pieces of information can be used such as superimposed cell identity markers and morphological information of the neurons. Currently, superimposing cell identity markers, such as fluorescent proteins expressed by well-characterized cell-specific promoters, is the most popular and reliable method for neural identification. For example, Serrano-Saiz et al. 2013 showed that such methods are effective when the number of the target neurons is limited (Serrano-Saiz et al. 2013). However, integrating this approach with the whole-brain activity imaging seems difficult because it requires different markers and different fluorescent channels for every neuron in principle. Morphological information is also useful when the number of target neurons is very limited but is not readily utilized for whole-brain imaging because the neurons are distributed densely in the head region of the worms and the morphological information cannot be obtained accurately.
Several efforts for developing automatic annotation methods were reported. In order to annotate the neurons based on their positions, the information of the positions and their variations will be required. Long et al (Long et al. 2008(Long et al. , 2009 produced 3D digital atlas for 357 out of 558 cells from several tens of L1 animals, and related works also used the atlas (Qu et al. 2011;Kainmueller et al. 2014). The atlas consists of positions and their deviations of the cell nuclei of body wall muscles, intestine, pharyngeal neurons, and neurons posterior to the retrovesicular ganglion, as well as some other cell types. However the neurons anterior to the retrovesicular ganglion are omitted because of their dense distribution (Long et al. 2009), and the atlas cannot be applicable to the neurons in head region important for neural information processing. Aerni et al (Aerni et al. 2013) reported positions of 154 out of 959 cells from 25 adult hermaphrodites, including intestinal, muscle, and hypodermal cells, and introduced a method that integrates useful features including fluorescent landmarks and morphological information with the cell positions. Nevertheless, the positions of neurons were not reported. As far as we know, the information of the positions of the neurons in adult worms can be obtained only from the atlas produced by the EM reconstruction work (White et al. 1986).
Unfortunately, the White atlas does not have the information about the variety of the positions between individual animals. Additionally, the atlas may be deformed because of inherent characteristics of the sample preparation methods for electron microscopy. Thus, experimental data of positions of neurons in adult animals were very limited. Also automatic annotation methods for neurons in the head regions of the worm have not been reported.
Here we measure the positions of the neurons in adult animals by using multiple cell-specific promoters and create a dataset. We evaluate the variations of the positions and obtain an optimal combination of the cell-specific promoters for annotation tasks based on accumulated information of cell positions. We also develop and validate an efficient annotation tool that includes both automated annotation and human interaction functionalities.

An annotation dataset of head neurons
In this study, we focused on the head neurons of an adult animal of the soil nematode C. elegans, which constitute the major neuronal ensemble of this animal (White et al. 1986). The expression patterns of cell-specific promoters were used as landmarks for cell identification ( Figure 1A). The fluorescent calcium indicator Yellow-Cameleon 2.60 was expressed in a cell-specific manner by using one of the cell-specific promoters and used as a fluorescent landmark. All the neuronal nuclei in these strains were visualized by the red fluorescent protein mCherry. Additionally, the animals were stained by a fluorescent dye, DiR, to label specified 12 sensory neurons following a standard method (Shaham 2006). The worms were anesthetized by sodium azide and mounted on the agar pad. The volumetric images of head region of the worm were obtained with a laser scanning confocal microscope. All the nuclei in the images were detected by our image analysis pipeline roiedit3D (Toyoshima et al. 2016) and corrected manually. The nuclei were annotated based on the expression patterns of fluorescent landmarks.
Finally, we obtained volumetric images of 311 animals with 35 cell-specific promoters in total ( Figure 1B). On average, 203.7 nuclei were found and 44.2 nuclei were identified ( Figure 1C).
These positions and promoter expression information are hereafter called the annotation dataset. Figure 1D shows names and counts of the identified cells. In most animals, 12 dye-stained cells and 25 pharyngeal cells were identified. Finally, we identified a total of 175 cells in the head region anterior to the retrovesicular ganglion. We didn't identify URA class (4 cells), RIS cell, SIA class (2 of 4 cells) AVK class (2 cells), and RMG class (2 cells) because of the lack of suitable cell-specific promoter.
Please note that we use H20 promoter as a pan-neuronal promoter (Shioi et al. 2001). We confirmed H20 promoter was expressed in the GLR cells and XXX cells by expressing the cell-specific promoters nep-2sp and sdf-9p, respectively. We estimated that H20 promoter was expressed in pharyngeal gland cells and HMC cell, based on their positions. Also we estimated that H20 promoter was expressed weakly in the hypodermal cells, based on their positions and shape of the nuclei, but we remove these cells from our annotation dataset. We confirmed the promoter was not expressed in the socket cells nor the sheath cells by expressing the cell-specific promoter ptr-10p.     After removing contribution of (1) and (2)

Optimal combination of the cell specific promoters
In order to reduce the error rate of the annotation method, one may want to use the information of fluorescent landmarks (Kotera et al. 2016;Nguyen et al. 2016). Using multiple landmarks will reduce the error rate. One or two fluorescent channels are often available for the landmarks in addition to the channels required for the whole-brain activity imaging. We therefore sought for the optimal combination of cell-specific promoters for two-channel landmark observation using the annotation dataset.
Several properties of the promoters were evaluated in order to choose the optimal combination; how many number of cells are labelled ( Figure 3A), stability of expression ( Figure 3A), sparseness of the expression pattern ( Figure 3B, see Methods for definition), and overlap of expression patterns in the case of combinations (See Supplementary Dataset 1). Among the 35 tested promoters, eat-4p was selected because it was expressed in the most numerous cells in the head region ( Figure   3A). The promoters dyf-11p and glr-1p were expressed in numerous cells, and glr-1p was selected as the second promoter because the sparseness of the expression patterns of glr-1p is higher than that of dyf-11p ( Figure 3B) and because the expression patterns of dyf-11p highly overlapped with that of eat-4p. Additionally, ser-2p2 was selected based on the stability of the expression and low overlaps with eat-4p and glr-1p. Thus the combination of eat-4p, glr-1p and ser-2p2 was selected ( Figure 3C).
The latter two promoters were used with the same fluorescent protein assuming only two fluorescent channels can be used for the landmarks as is the case for our experimental setup for whole-brain imaging. In the annotation dataset, eat-4p was expressed in 69 cells and glr-1p + ser-2p2 were expressed in 50 cells out of 196 cells in the head region of adult worms.
All combinations of the promoters could be evaluated by an algorithm that considers the number of expression, sparseness and overlap of expression patterns (see Methods and Supplementary   Table 1). In brief, the algorithm highly evaluated a combination when two neighboring cells were in different colors. In the case of three promoters and two fluorescent channels, the combination consisting of eat-4p, glr-1p and ser-2p2 was placed in the 18th rank out of the possible 20825 combinations.
Here we produced a strain JN3039 as follows. The far-red fluorescent protein tagRFP675 was expressed using eat-4p, and the blue fluorescent protein tagBFP was expressed using glr-1p and ser-2p2. The red fluorescent protein mCherry was expressed using the pan-neuronal promoter H20p.
This strain does not use fluorescent channels of CFP, GFP, and YFP and is useful for the cell identification tasks. For example, if there is a strain that express these fluorescent proteins with a promoter whose expression patterns should be identified, one can do it by just crossing the strain with our standard strain JN3039.
Additionally a strain JN3038 was made from the strain JN3039 by expressing fluorescent calcium indicator Yellow-Cameleon 2.60 with the pan-neuronal promoter H20p. This five-colored strain will enable whole-brain activity imaging with annotation.    Figure 2C.

Automatic annotation using bipartite matching and majority voting
Having a set of atlases that capture the positional variations of the cells, we can account for the spatial uncertainty of cell annotation using majority voting. We propose an automatic annotation method that utilizes bipartite matching and majority voting ( Figure 5A). Each part will be described below in brief (see Method for detail).
In the bipartite matching step, the cells in a target animal were assigned to the cells in an atlas. An assignment of a cell in the target animal to a cell in the atlas has a cost based on the similarity (or dissimilarity) of the two cells including Euclidean distance, expressions of landmark promoters, and the feedback from human annotation. The optimal combination of the assignments that minimizes the sum of the costs was obtained by using Hungarian algorithm. The name of the cell in the target animal can be estimated as the name of the assigned cell in the atlas in this step.
To handle the configurational variations of cells to be annotated, we used the majority voting technique. Assuming the generated atlases could capture the positional variations of the cells, we assigned unannotated cells in the target animal to those in atlases, then giving annotation results. Each assignment of a cell is considered as a vote, and the most voted assignment was considered as the top rank estimation of annotation.
In order to validate our automatic annotation method, a 5-fold cross validation test was performed. All the animals in the annotation data set were randomly divided into five subsets. We perform a total of 5 tests. For each test, we exclude one of the subset from training of atlases, and use it to estimate the annotation performance based on the trained atlases. The error rate of bipartite matching was relatively high, and the majority voting could deliver significant improvements of the annotation accuracy ( Figure 5B). On average, 78.3 nuclei were annotated and 46.2 nuclei were successfully estimated as the top rank, and the error rates of the top rank estimation was 41.4% ( Figure   5B and C). As a control, two methods are introduced; one method only considers the mean and covariance of the cell positions of raw data (without using the atlases and voting, see Figure 2D). The other method considers the mean and covariance of the cell positions in the atlases (without using majority voting). The error rate of the two methods were higher than the proposed method, indicating that the majority voting step in the proposed method contribute to the correct estimation. If we consider the accuracy for the top 5 voted estimations (shown as rank 5), the error rate decreased to 7.3%.
The automatic annotation method was applied to the animals with fluorescent landmarks (strain JN3039, see Figure 3C-E). With the help of the optimized expression of landmark fluorescent proteins, the number of identified cells in an animal will increase compared to the strains used to make the annotation dataset. On average, 202 nuclei were found and 156.3 nuclei were identified from 15 adult animals. The error rates of the top rank estimation with and without fluorescent landmark were 37.7% and 51.5%, respectively, indicating that utilizing the fluorescent landmark also contribute the correct estimation ( Figure 5D). If we consider the accuracy for the top 5 voted estimations, the error rate decreased to 8.1%. These error rates were comparable to the cross-validation results for the annotation dataset, suggesting that our annotation framework will work correctly for the whole-brain activity imaging.
The automatic annotation method was also applied to the animals in a microfluidic chip for whole-brain activity imaging ( Figure 5figure supplement 1). The error rates of the top rank estimation with and without fluorescent landmark were 52.1% and 72.8%, respectively, and that of the top 5 voted estimations was 12.2%. The worms were compressed and distorted to be held in the microfluidic chips, and the distortion of the worm may increase the error rates. During whole-brain imaging for free-moving animals (Nguyen et al. 2016(Nguyen et al. , 2017Venkatachalam et al. 2016), the worms will be less compressed and less distorted, and our algorithm may works better.
Additionally, our algorithm is implemented in the GUI roiedit3d (Toyoshima et al. 2016), and it can handle feedback information from the human annotations. Once annotations are corrected manually, our method can accept corrections and uses them to improve the results. For example, one can identify neurons manually by using other information including the neural activity or morphology, and the automatic estimation for the other neurons will be improved. The final results can be added to the annotation dataset and the annotation algorithm will work more accurately. Thus the feedback system incorporates tacit knowledge into the automatic annotation method. Through the interactive process our algorithm will make human annotation tasks more efficient. (C) Error rates of the automatic annotation method for the animals in the promoter dataset. The error rates were evaluated by cross-validation.
(D) Error rates of the automatic annotation method for the strain JN3039 that expresses the fluorescent landmarks.
(E) The automatic annotation method was integrated in the graphical user interface roiedit3d that enables feedback between automatic and manual annotation.

Figure 5 -Figure Supplement 1
Error rates of the automatic annotation method for the animals in a microfluidic chip for whole-brain activity imaging (JN3038 strain, n=12).

Discussion
In this study, we obtained volumetric fluorescent image of 311 animals using 35 promoters, and created an annotation dataset that contains the positions of the identified cells and expression patterns of promoters in respective animals. Utilizing the annotation dataset we evaluate the variation of the positions of the cells and choose the combination of the promoters optimal for our annotation tasks. We proposed an automatic annotation method and validated its performance on head neurons of adult worms for whole-brain imaging. Thus, we successfully integrate the annotation techniques with the whole-brain activity imaging.
The cell positions of real animals and its variation will be the most important information for the cell identification. As far as we know, this might be the first report about the large-scale information of the positions of the cells in the head region of adult C. elegans, which lead to systematic and comprehensive method to annotation of the head neurons. The error rate of the automatic annotation might be slightly high for fully automatic annotation. The integration of the automatic annotation method to the GUI enables machine-assisted annotation and enhances the process of whole-brain image annotation. Increasing the number of animals and promoters will improve the accuracy and objectivity of the automatic annotation method.
Increasing the number of fluorescent channels and landmarks will also improve the accuracy.
Long-stokes shift fluorescent proteins might be good candidates because they use irregular fluorescent channels that will not be used in standard application. In our case, however, these proteins disrupted the neighbor fluorescent channels by leaking-out. Employing color deconvolution techniques will increase the number of substantial fluorescence channels and may improve the accuracy.
The images of the animals we recorded will have useful information for annotation including size of the nuclei and intensities of the fluorescence. In the manual annotation process we utilized these pieces of information for improving accuracy. On the other hand our automatic annotation algorithm does not utilizes these pieces of information and it may be one of the causes of relatively low accuracy of the algorithm. Recent advances in artificial neural networks especially in the field of image analysis will enable to utilize such information for automatic annotation. It is well known that artificial neural networks require large amount of training data composed of images and the corresponding grand truth. Our annotation dataset contains images with identity information and will be ideal for the training data, but the number of data may not be enough. Our method that makes annotation more efficient will play an important role for opening up the path to utilization of artificial neural networks in the future.
There are no dataset of cell positions that can be used as a benchmark of cell identification methods. For example, a new method that solve the cell identification problem as a nonlinear assignment problem was reported recently (Bubnis et al. 2019). The report utilizes synthesized data and does not use real data. To evaluate the real performance of new methods, the method should be tested on real data. Our annotation dataset will be an ideal benchmark of newly developed cell identification methods. Thus our study will facilitate the future studies for automatic annotation methods.
In order to identify the expression patterns of the promoters, the most accurate method is testing whether the fluorescence of promoter overlaps with the fluorescence of the neuronal identity markers (Serrano-Saiz et al. 2013). In such cases our standard strain and automatic annotation method will help the selection of the markers through objective estimation of cell identities.

Microscopy
A set of static 3D multi-channel images of C. elegans strains ranging from JN3000 to JN3036 were obtained as follows. Day 1 adult animals were stained by the fluorescent dye DiR (D12731, Thermo Fisher Scientific) with the standard method (Shaham 2006). The stained animals were mounted on a 2% agar pad and paralyzed by sodium azide. The fluorescence of the fluorescence proteins and the dye was observed sequentially using laser scanning confocal microscopy (Leica SP5 with 63× water immersion lens and 2× zoom). The sizes of the images along the x1 and x2 axes were 512 and 256 voxels, respectively, and the size along the x3 axis varied depending on the diameter of the animal.
A set of 3D multi-channel images of strain JN3039 was obtained as described above without using the fluorescent dye DiR.
A set of 3D multi-channel images of strain JN3038 was obtained as follows. Day 1 adult animals were conditioned on NGM plate with OP50 (Kunitomo et al. 2013). The conditioned animals were introduced and held in a microfluidic device called olfactory chip (Chronis, Zimmer, and Bargmann 2007). The depth and width of the fluid channel in the chip were modified in order to reduce the distortion of the worms. The animals and their head neurons moved to some extent in the device because the animals were not paralyzed. The fluorescence of the tagBFP, tagRFP675, and mCherry channels was observed simultaneously using customized spinning disk confocal microscopy. The sizes of the image along the x1 and x2, and x3 axes were 512, 256, and 50 voxels, respectively. The sizes of a voxel along the x1, x2, and x3 axes were 0.28, 0.28, and about 0.77 μm, respectively.

Image analysis for the annotation dataset
All the nuclei in the images were detected by our image analysis pipeline roiedit3D (Toyoshima et al. 2016) and corrected manually. The cells stained by the chemical dye were identified as reported (Shaham 2006  For the pairs of 4 ≤ ≤ 10 , all the permutations ( equal to or less than 10! ~ 3.6 × 10 6 combinations) were calculated. For the pairs of > 10, 1 × 10 7 permutations were randomly selected and calculated.

The algorithm for searching optimal combination of cell-specific promoters and the definition of the sparseness
The most important factor for selecting promoters in order to improve annotation accuracy is to achieve a checkerboard-like coloring pattern for the ease of separating neighboring cells. A simple metric to account for this factor is to sum the number of neighboring cell pairs that exhibit a different color based on cell-specific promoters, where each pair is inversely weighted by the distance between the two neurons. Such a metric can be considered as a modification to an Ising model in physics. We choose a Gaussian probability model for the weighting function with an empirically chosen value of the standard deviation to be 9.6 μm. The metric can be written as where is a set of all cells in an animal. and are positions of cell X and cell Y, respectively.
is label of cell X and = (1,0) means that landmark protein of color 1 is expressed in the cell X but that of color 2 is not expressed. Because the experimental setup has a limited amount of channels, we are able to perform an exhaustive search for all possible combination of the available promoters, and compare the final values of the metric as a reference for choosing the combination of cell-specific promoters used in our experiment. We evaluated all the combinations for 3 promoters and 2 colors (20825 combinations). The scores of the single promoter for single color were used as the index of sparseness.

Generating atlases
To obtain an atlas with fully annotated cells, we need to combine positional information of cells from multiple partially annotated images while maintaining the relative position between the cells as much as possible. We achieve this goal by maximizing the consistency (or smoothness) of a displacement flow when combining different images, for which the displacement flow is defined as follows.
Suppose that in two images, denoted by 0 and 1 , there coexist annotated cells. The displacement of cell is denoted by = 1 − 0 where 0 and 1 denote the positions of the cell in 0 and 1 , respectively. Then, we define a displacement flow field 0→1 ( ) from 0 to 1 on the entire space ∈ ℝ 3 : (1) Here, ( | , ) denotes the density function of the normal distribution with mean and covariance Σ (Please note that is 3 × 3 covariance matrix and determine effective range of the displacement of a cell). This represents a flow field function interpolated by the given displacements of the cells in the two images. When taking the weighted average in the calculation of 0→1 ( ), larger weights are assigned to the displacements of more neighboring cells with respect to in 0 .
To generate an atlas, we conducted the following steps (see Figure 4A): 1. Set a randomly ordered sequence { 1 , … , 311 } of the 311 partially annotated animals.
We discard the sequence if the 1 has less than 60 annotated cells.
2. For ∈ {2 … ,311}, cells in were sequentially aligned to those in 1 as follows: A) The positions of all annotated cells in 1 were unchanged.
B) All annotated cells that coexisted in both 1 and were used to calculate the displacement field →1 ( ) with a pre-determined (Eq 1).
C) All cells annotated in but not in 1 with their positions denoted by were shifted and aligned to 1 according to 1 ← + →1 ( ). Add them to the annotated cells in 1 .

D) Terminate the iteration if all annotated cells have been aligned in the synthesized
reference image.
In this scheme, a spatial pattern of produced cells was largely affected by the interpolated flow fields.
In general, the performance will be poor if the number of observed source displacements was small.
To reduce such instability, we skipped and used it later when the shared less than half cells annotated in common with 1 . Repeating this procedure, we generated 3,000 reference samples.
The generated reference samples serve as a set of virtual atlases that imitate observed topological variations of cellular positions across different worm samples. To obtain more realistic atlases, we optimized = diag( 1 , 2 , 3 ) in Eq 1, which is the parameter to control the smoothness of displacements in the sequential alignments. We defined an objective function to reflect the similarity of the topological variations between our raw data set and the generated atlas.
By optimizing such objective function and taking the optimal values of the parameters as a reference, we selected an empirical value of = diag(9.6 μm, 9.6 μm, 9.6 μm). Details of the objective function and optimization is in Supplementary Note 1.

Bipartite graph matching
Detected cells in a target animal and an atlas were matched using the Hungarian algorithm to solve the bipartite graph matching problem. The matching was achieved by comparing one or more selected features between cells. Here, features refer to some quantitative properties for the cells that can be used to distinguish the identity of a cell from another. The most fundamental feature is the positions of cells. Other typical features include cell volume, fluorescence intensities, and so on. We use expression of landmark proteins (i.e. binarized fluorescent intensities) and feedback from human annotation. With such features, the dissimilarity of cells was represented by a matrix , where the { , } entry is the distance of the feature values between the th cell in the target and the th cell in the atlas. When there are features chosen, we can assemble them into a single matrix BGM through a weighted sum: where is the weight for each feature. We use = 1, = 20. For feedback from human annotation, the assignments incompatible with the human annotation have infinity dissimilarity.
With a given assignment, we can calculate the sum of the dissimilarity values in BGM that correspond to the selected matching. A modified Hungarian algorithm (Jonker and Volgenant 1987) was used to minimize the total distance with respect to all possible assignments under the constraint of one-to-one matching.

Majority voting
Multiple name assignments of a cell in the subjective animal were obtained by repeating the bipartite graph matching using 500 different atlases. Each assignment was considered as one vote, and the estimated names for a target cell was ranked by vote counts. The estimation for a cell was independent of each other and multiple cells may have the same estimated names. If non-overlapping result is required, one can assemble cost matrix based on vote counts and apply the Hungarian algorithm.

Calculation of error rate of automatic annotation
All the detected cells in a target animal other than hypodermal cells were used as a target. The names of the cells were estimated by our automatic annotation method based on their positions. Expression of landmark promoters were also used for Figure