Multiple image stacks containing portions of the organ of Corti were assembled into a single image stack. Shifts of coordinates between image stacks were calculated based on cross-correlation (MATLAB ‘normxcorr2’ function). After image stitching, a blending algorithm (Rankov et al., 2005) was applied to remove sharp intensity changes in the zone of overlap.

Hair cells in each image stack were detected as local intensity peaks (MATLAB ‘imregionalmax’ function). Single-linkage clustering (maximal distance of connection, 25 μm) was effective for eliminating or reducing the number of false positives. A stretch of local peaks corresponding to the entire row of hair cells were divided into segments of 200–300 μm in length. In each segment, the best-fit plane was calculated (MATLAB ‘pca’ function), together with the best-fit arc along the rows of hair cells. The multiple best-fit arcs were stitched into a continuous curve (Figure 2B). A voxel image containing the entire straightened row of hair cells was reconstructed from the image stacks, based on the stitched-fit curve and the normal vectors of the fit planes.

First, local correlation between the hair cell template and the voxel image of linearized epithelium was calculated by template matching (MATLAB ‘normxcorr2’ function), and the peaks of the correlation were detected. Pixels corresponding to detected peaks were grouped according to the physical size of the IHCs via connected-component labeling. These connected pixel groups (hereinafter called ‘cell candidates’) were used as a first approximation of IHC positions linked to other attributes, including correlation values and local intensity distributions.

The cell candidates were further evaluated to eliminate false positives using two successive machine learning models. The first ensemble learning method created the model for selection with predictor data consisting of areas, barycentric coordinates, correlation values, the intensities of the peaks, and the corresponding values of nearby cell candidates (MATLAB ‘fitensemble’ function with ‘GentleBoost’ method) (Friedman et al., 2000). The model was trained to calculate posterior probability (prediction score), and cell candidates with a high prediction score (~1000 candidates out of initial ~50,000) were selected and further analyzed by the second ensemble learning method (MATLAB ‘fitensemble’ function with ‘Bag’ method) (Breiman, 2001). This method was based on expanded predictors (the prediction score from the first step of the candidate and nearby candidates, and the local intensity distribution centered on the barycentric coordinates of the peaks). The cell candidates after the second selection were connected sequentially, subject to the physical constraint that the IHCs must form a single row with roughly constant intervals of more than 6 μm. The resulting putative positions of IHCs were used for fine readjustment of image linearization and three-dimensional structural analysis.

The image processing applied for IHCs in Step 3 was also applied to OHCs. Detection accuracy was improved by two additional evaluations based on machine learning. First, physical constraints of OHC alignment were introduced into three rows. A multiclass classification model, based on the convolutional neural network method [Neural Network Toolbox of MATLAB (LeCun et al., 1989)], sorted cell candidates into respective rows using input images each containing three rows of four or five OHCs. If the distance between two adjacent cell candidates in the same row exceeded 1.5 times the average distance, the presence of additional cells in the gap was assessed by the fourth model based on the convolutional neural network method. Input images for machine learning were sampled by placing small rectangular areas at equal distances from one another within the gap. If the model predicted the existence of additional cells in the gap, the nearest peaks of the correlation coefficient from the first template matching were recovered.

Mouse husbandry, anesthesia, and euthanasia conformed to related regulations of the government and the institutional guidelines. Protocols related to animal handling were approved by the Animal Care and Use Committee of the Graduate School of Medicine, the University of Tokyo. Male or female wild type C57BL/6J, ICR, and CBA/Ca mice at ages of PND 0 to 360 were used for application of modified Sca/eS with or without antibody labeling. For the detection of fluorescent protein signals, Thy1-GFP M line and a transgenic line expressing a fluorescence resonance energy transfer (FRET)-based indicators (GO-ATeam mouse line) (Imamura et al., 2009; Nakano et al., 2011) were used at ages of PND 0 to 120.

After euthanasia, mice were perfused transcardially with 4% paraformaldehyde (PFA) in phosphate buffered saline (PBS). Osteochondral samples (cochlea embedded in temporal bones and femurs) and other soft tissues (brain, heart, stomach, lung, liver, kidney, intestine, and spleen) were isolated by standard dissection techniques. For mice younger than PND 5, the step of transcardial perfusion was omitted. The isolated tissues were placed in the same fixative (4% PFA in PBS) for 10 to 24 hr. Fixed brain samples were sectioned into 100 μm to 2 mm thick sections using a vibratome.

iDISCO; iDISCO was performed according to the published protocol (Renier et al., 2014). Fixed samples were treated with methanol with increasing concentrations, bleached with 5% H₂O₂, and rehydrated by decreasing concentration of methanol in PBS. Anti-MYO7A staining was performed after 0.3 M glycine treatment and blocking with 6% normal goat serum. Tissue clearing was achieved by treatment with increasing concentration of tetrahydrofuran/H₂O up to 80% (v/v), subsequent transfer to dichloromethane, and final incubation with dibenzyl ether.

3DISCO; 3DISCO was performed according to the published protocols (Acar et al., 2015; Ertürk et al., 2012; Yokomizo et al., 2012). Fixed samples were blocked with normal goat serum and then reacted with primary and secondary antibodies. After immunolabeling, samples were treated with increasing concentration of tetrahydrofuran/H₂O up to 80% (v/v), subsequent transfer to dichloromethane, and final incubation with dibenzyl ether.

CLARITY; CLARITY was performed according to the published protocols (Chung et al., 2013; Tomer et al., 2014). Briefly, animals were fixed by perfusion with a solution containing 4% paraformaldehyde, 4% acrylamide, 0.05% bis-acrylamide, 0.25% VA-044 initiator in PBS, followed by induction of hydrogel tissue embedding by raising temperature to 37°C for 3 hr. After electrophoretic tissue clearing in a chamber containing 200 mM boric acid and 4% sodium dodecyl sulfate with 20 V applied across the sample at 37°C for 72 hr, samples were blocked with normal goat serum, reacted with primary and secondary antibodies. After washing, the samples were placed in FocusClear (CelExplorer Labs) medium.

CUBIC; CUBIC was performed according to the published protocols (Susaki et al., 2014; Susaki et al., 2015; Tainaka et al., 2014). Briefly, fixed samples were placed in ScaleCUBIC-1 (25% urea, 25% N,N,N’,N’-tetrakis (2-hydroxypropyl) ethylenediamine, 15% Triton X-100) at 37°C for 24 hr, followed by brief washing and incubation with ScaleCUBIC-2 (50% sucrose, 25% urea, 10% 2,2’,2’’-nitrilotriethanol) for 2 hr. After brief washing, samples were placed in ScaleCUBIC-2.

CB-perfusion; CB-perfusion was performed according to the published protocols (Tainaka et al., 2014). Briefly, animals were first treated with perfusion of an excess amount of 4% paraformaldehyde in PBS, followed by sequential perfusion with PBS and 50% of ScaleCUBIC-1. The target organs were excised and immersed in ScaleCUBIC-1 for 5 days. After this step, samples were transferred to ScaleCUBIC-2 and incubated for 2–3 days. After brief washing, samples were placed in ScaleCUBIC-2.

ScaleS; ScaleS was performed according to the published protocols (Hama et al., 2015). Briefly, inner ear was isolated by a standard dissection procedure. Cochleae were isolated and 4% paraformaldehyde in PBS was perfused from the oval window. ScaleS0 (20% sorbitol, 5% glycerol, 3% dimethylsulfocide, 1% N-acetyl-L-hydroxyproline, 1 nM methyl-β-cyclodextrin, 1 mM γ-cyclodextrin), ScaleA2 (10% glycerol, 4 M urea, 0.1% Triton-X-100), ScaleB4(0) (8 M urea), and ScaleA2 were applied sequentially for permeabilization and tissue clearing. After samples were returned to PBS, immunolabeling was performed. Samples were rinsed by AbScale solution (0.33M urea, 0.1–0.5% Triton X-100 in PBS) and AbScale rinse solution (2.5% BSA, 0.05% Tween-20 in 0.1 × PBS). Finally, samples were placed in ScaleS4 (40% sorbitol, 10% glycerol, 4 M urea, 0.2% Triton X-100, 15–25% dimethylsulfoxide).

Antibodies used in this study were as follows:

MYO7A (rabbit polyclonal, Proteus Bioscience), 200 kDa subunit of neurofilament protein (anti-NF200, mouse monoclonal, clone NE14, Sigma Aldrich), SOX2 (rabbit polyclonal, EMD Millipore), CtBP2 (mouse monoclonal, clone 16/CTBP2, BD Biosciences), VGLUT3 (guinea pig polyclonal, kindly gifted from Dr. Hioki), Alexa Fluor 488 goat anti-rabbit IgG (H + L) (Life Technologies), Alexa Fluor 546 goat anti-mouse IgG (H + L) (Life Technologies), Alexa Fluor 647 goat anti-mouse IgG (H + L) (Life Technologies), and Alexa Fluor 488-conjugated VE cadherin (mouse monoclonal, clone BV13, eBioscience). Concentration of antibodies should be adjusted depending on the sample size. Rhodamine phalloidin (Thermo Fisher) was used for labeling F-actin.

For samples labeled with multiple antibodies shown in Figure 5B, images were obtained by a FV1000 laser scanning microscope (Olympus) with a 25 × objective lens (NA = 1.05). Fluorophores were excited by 488, 564, and 635 nm lines of diode lasers. The sizes of single horizontal images were set to 512 × 512, with pixel sizes of 0.43 × 0.43 μm and z-spacing of 0.75 or 2 μm.

In double-labeling (anti-MYO7A antibody and rhodamine phalloidin) samples, fluorescence intensity of rhodamine phalloidin was measured as follows. First, MYO7A positive cells were automatically detected by the custom-made program. The space without MYO7A staining was simultaneously detected as a putative position of cell loss. Second, rhodamine phalloidin fluorescence intensity was measured using ImageJ software. OHCs surrounded by intact OHCs and nearby OHCs next to the cell loss positions were selected. Finally, the rhodamine phalloidin intensities between different positions of the organ of Corti were normalized by the intensity of OHCs surrounded by intact OHCs (Figure 5A).

Tissue transparency was measured according to the methods described by Hama et al (Hama et al., 2015). Light transmitted through 100 μm-thick mouse brain sections before and after clearing was captured with a CCD camera and quantitated. Transparency was normalized to the samples without clearing treatments. Increase in tissue size was determined by measuring the area of brain sections before and after clearing. The extent of size increase was expressed as a ratio against the pre-cleared area (Figure 1—figure supplement 1).

To assess the imaging depth of bony tissue, 20 μg of Alexa Fluor 488-conjugated anti-mouse VE cadherin antibody was administered into tail vein followed by tissue fixation. The metaphysis of the tubular bone was analyzed. Imaging depth was measured from the surface and the position at which specific fluorescence signal of the vasculature over the background could no longer be detected was taken as the limit of the imaging depth (Figure 1—figure supplement 2).

Mice were placed within steel wire cages (20 cm ×12 cm × 7 cm) in a ventilated sound exposure chamber. The acoustic stimulus was an octave band noise with a center of 4 kHz at 121 dB sound pressure level (SPL) for 4 hr. The sound-delivery speakers (HFD-261–8, TOA) were driven by a power amplifier (IP-600D, TOA) attached to a noise generator (AA-61B, RION) through a programmable filter (3624, NF Corporation). Sound levels were measured (UC-31 and UN14, RION) at multiple locations within the sound chamber and calibrated (NC-74, RION) to ensure uniformity and stability of the stimulus.

ABRs were recorded before and 7 days after noise exposure. Mice were anesthetized with ketamine hydrochloride (50 mg/kg, i.p.) and xylazine hydrochloride (10 mg/kg, i.p.). Before the ABR test, the tympanic membranes were confirmed to be normal. ABRs were evoked with tone-bursts (5 ms duration at 4, 8, 16, and 31.25 kHz) and measured by a recording system (Neuropack ∑ MEB2208, Nihon Kohden). Needle electrodes were placed subcutaneously at the vertex of the skull, under the ear exposed to noise, and under the opposite ear for ground. A speaker was placed 10 cm from the tragus of the stimulated ear. ABRs from 500 trials were averaged at each sound intensity. ABR thresholds were estimated by changing the intensity in 5 bB steps and finding the lowest sound intensity where reliable response peaks were detected

C57BL/6J (PND 5: four samples, PND 30: 10 samples, PND 60: 14 samples, PND 120: nine samples, PND 360: four samples, and PND 60 with noise: 10 samples) mice and CBA/Ca (PND 60: five samples) mice were analyzed. Statistical analysis was performed by paired t-test (Figure 2G), Kruskal-Wallis test with Steel-Dwass multiple comparison test (Figure 3B), paired t-test followed by Bonferroni’s correction (Figure 3D), Welch’s t-test (Figure 4A), paired t-test with Bonferroni’s correction and Welch’s t-test with Bonferroni’s correction (Figure 5), one-way ANOVA followed by Bonferroni’s post hoc test (Figure 1—figure supplement 1 and Figure 2—figure supplement 1A), or a two-tailed unpaired t-test (Figure 1—figure supplement 2A and Figure 2—figure supplement 1B,C). Data are presented as means ± standard errors of the mean (SEM) except standard deviation (SD) in Figure 4E. The p-values are as follows: *p < 0.05, **p < 0.01, ***p < 0.001.

The outline of the protocol, including stitching of multiple image stacks into a single stack, reconstruction of linearized image, automated detection of IHCs, and automated detection of OHCs, was described in the main text. Full description of the algorithm is in Appendix 2.

The outline of the framework was described in the main text. Full description of the models is in Appendix 2.

To evaluate the performance of automated OHC loss counting program, three skilled human operators (A, B, and C) initially detected empty spaces in the same set of the images (n = 4). A human consensus (HC) was established for the presence of empty spaces when more than two operators agreed. The match rates among human operators (A-B, B-C, and A-C) and between the program and HC were comparable (A-B: 91.6, B-C: 89.8, A-C: 89.7, Auto-HC: 93.1, Table 3). The numbers of lost cells in individual empty spaces were then estimated by both the human operators and the program using 161 empty spaces (Figure 2G). The differences between the human operators and the program were defined as absolute errors, and the errors were averaged (mean absolute error) and compared between the program and the human operator (Auto-M1, Auto-M2, and Auto-M3) or between the human operators (M1-M2, M2-M3, and M1-M3).

The longitudinal pattern of cell loss was analyzed among five segments (‘a’, ‘b’, ‘c’, ‘d’, and ‘e’) along the longitudinal axis (x-axis). The two segments from the basal end with their lengths of 800 μm ('a’ and ‘b’), two segments from the apical end with their lengths of 800 μm (‘e’ and ‘d’) were defined, and the middle remaining segment was defined as ‘c’ (Figure 3B). The radial axis (y-axis) was divided into 15 successive zones from the proximal to the distal, covering the width of three rows of OHCs (Figure 3D).

For automated detection and analyses of hair cells, we examined 43 sets of images from 43 samples (PND 30: 10 samples, PND 60: 14 samples, PND 120: nine samples, noise exposure at PND 60: 10 samples). The samples immunostained with MYO7A were imaged with a two-photon microscope. To cover the entire longitudinal length of the organ of Corti, multiple image stacks with overlapping volumes should be obtained. Single image stack consists of 2D images with their pixel sizes of 512 × 512, with image numbers in the range of 36–459 along the z-axis. Large differences in the stack sizes were due to changes in the height of the fluorescent objects along the spiral of the cochlea. The image data set for a single cochlea sample contains 11–16 image stacks, which follow the spiral organization of the organ of Corti from an apical end toward the hook portion. There were 10–40% of overlap in length between adjacent image stacks. Physical longitudinal lengths of the organ of Corti were relatively uniform across samples (Figure 2C in the main text). The voxel sizes in x, y, and z directions were 0.99, 0.99, and 1.0 μm, respectively.

The scripts written for computational detection and analysis of hair cells consist of the following steps:

Step 1: Stitching of multiple image stacks into a single stack

Step 2: Reconstruction of linearized image

Step 3; Automated detection of IHCs

Step 4; Automated detection of OHCs

Multiple image stacks containing parts of the organ of Corti were assembled into a single image stack (Rankov et al., 2005). The cross-correlation method was used to image stitching. For rough estimation of the relative positions of two adjacent image stacks, we used recorded positions of the microscope stage and the overall distribution of fluorescent signals derived from the organ of Corti after appropriate image thresholding. Relative positions of two adjacent images were determined by searching best cross-correlation of images at the putative region of overlap (MATLAB ‘normxcorr2’).

Overlapped images were blended with gradient mixing of intensities. To simplify the explanation, the following protocol is for two-dimensional blending. Blending of image stacks could be performed with similar principles by setting blending gradient in three dimensions. The mixing ratio of the two images was calculated using a sigmoid function as:

where $R_1$ and $R_2$ denote mixing ratios of the first and the second images, respectively. The value $x$ was set to reflect the relative distance of each pixel from the line running through the points with equal distance from two image centers (the center line). The normalized distance was measured from the center line, where both $R_1$ and $R_2$ were 0.5. $R_1$ reaches ~1 as $x$ moves toward the center of the first image, while $R_2$ reaches ~0. Because the blending process ends at the pixel position within the image overlap that has the largest distance from the center line, value $x$ was normalized to be 1 at this largest distance (Appendix 2—figure 1).

Appendix 2—figure 1

Download asset Open asset

After Image processing by a median filter (MATLAB ‘medfilt3’, neighborhood size: 3-by-3-by-3), local intensity peaks were automatically detected (MATLAB 'imregionalmax, 26-connected neighborhoods). Candidate anti-MYO7A-derived signals were selected from the local intensity peaks by eliminating those with their intensities below the threshold value α determined by the following formula.

where $\mathrm{\mu }$ and $\mathrm{\sigma }$ denote the mean and the standard deviation of image intensities in the background. Otsu’s method was used to determine the background thresholds (MATLAB 'multithresh', number of threshold values: 2, lower value was used) (Otsu, 1979). To further eliminate the intensity peaks derived from non-specific fluorescence, single-linkage clustering was performed (maximal distance of connection, 25 μm). This procedure was effective in identifying clustered intensity peaks corresponding to the organ of Corti as the largest assembly of linked intensity peaks.

The stretch of these intensity peaks corresponded to the entire row of hair cells. In the next step of the analysis, more precise determination of the structural parameters for the sensory epithelium is required. For this purpose, the spiral of the hair cell rows should be divided into segments of 200–300 μm in length. Because the sensory epithelium is curved in each original image stack and the spatial distribution of the intensity peaks roughly follows the shape of the sensory epithelium, principal component analysis [PCA (MATLAB ‘pca’)] was effective in extracting the direction and distance in both longitudinal and radial axes of the organ of Corti (Appendix 2—figure 2). Intensity peaks within a single image stack were divided into two equal halves along the first principal component axis, which were generally 200–300 μm in length. Both were projected onto PCA plane defined by the first and second principal component axes, and the best-fit arcs were calculated. The center and the radius of the arc were determined to minimize the mean square error (MATLAB ‘fminsearch’). Repetition of this procedure for all the original image stacks gave rise to a collection of the sensory epithelium segments, that encompassed the entire spiral of the organ of Corti.

Appendix 2—Figure 2

Download asset Open asset

The multiple best-fit arcs were stitched into a continuous curve. The nearest-neighbor points of two adjacent arcs were detected (Appendix 2—figure 3). Starting from these nearest-neighbor points, both of the arcs were converted to the connection of points with regular intervals of 50 μm. For the calculation of a spline curve that fits two adjacent arcs, the points along each arc were divided into two groups at the position of the nearest-neighbor point and the groups of points that belong to the main portion of the sensory epithelium were selected and combined into a single group. The nearest-neighbor points themselves were excluded from the combined group. The combined points were interpolated by using the spline interpolation method (MATLAB ‘interp1’, method: ‘spline’) and principal normal vectors were calculated (hereinafter just called ‘normal vectors’).

Appendix 2—Figure 3

Download asset Open asset

A voxel image containing the entire straightened rows of hair cells was reconstructed based on the stitched fit curve and the normal vectors. The horizontal center line (x axis) of the reconstructed image was set to match the fit curve. The y axis of the reconstructed image corresponds to the normal vector of the fit curve. The sizes of a voxel along the x, y, and z axes in the reconstructed image were all set to 1.0 μm. As each voxel in the linearized image does not show one-to-one correspondence to the original voxel image, interpolation was necessary to estimate the voxel intensity. In general, when a single voxel in the first image was projected to the second image, it was possible to define 2 × 2 × 2 voxels in the second image enclosing the volume projected from the first image. Therefore, the intensity of a voxel in the linearized image was estimated by the trilinear interpolation method using the corresponding 2 × 2 × 2 voxel in the original voxel image. More precise linearization was performed using a fit curve of IHCs detected by the procedure described in the next section. The voxel values are set to zero when the corresponding voxels are outside the original images.

Template matching was performed with a template image of an IHC (MATLAB ‘normxcorr2’, template image size: 11 × 11 pixels). This process puts a correlation coefficient to each element of a three dimensional matrix, which reflects the degree of matching to the template. By conducting two-dimensional peak detection on this matrix (MATLAB ‘imregionalmax’), candidate voxels for IHC positions were selected. The connected-component labeling was then carried out on the binary matrix of candidate IHCs (26-connectivity based). These connected pixel groups (hereinafter called ‘cell candidates’, Appendix 2—figure 4) were used as the first approximation of IHC positions linked to other attributes, including correlation values and local intensity distributions. The cell candidates were further evaluated to eliminate false positives using two successive machine learning models (the overview of our method related to machine learning is described in ‘Principles of auto-detection with machine learning’ section below).

Appendix 2—figure 4

Download asset Open asset

In the first step, the cell candidates were narrowed down to the number approximately twice larger than the expected number of IHCs. Specifically, the typical initial number of cell candidates was ~50,000 per one linearized image and they were narrowed down to ~ 1000 in this step. An ensemble learning method was employed to create the model for this selection (MATLAB ‘fitensemble’, method: ‘GentleBoost’) (Friedman et al., 2000). The predictor data was composed of the area, barycentric coordinates, correlation values, peak intensities, and those of nearby cell candidates. The model was trained to calculate posterior probability (prediction score) that the cell candidates corresponds to IHCs. The cell candidates with high prediction scores were selected for the next step.

In the second step, further selection of the cell candidates was performed by another machine learning model with expanded predictor data. The model was built based on an ensemble learning method (MATLAB ‘fitensemble’ function with ‘Bag’ method) (Breiman, 2001). The expanded predictor data was composed of the prediction score from the first step of the cell candidates and nearby cell candidates, and a small intensity image centered on the barycentric coordinates of the cell candidate.

Finally, the cell candidates were connected sequentially using physical constraints that the IHCs form a single row with roughly constant intervals of more than 6 μm along the x-axis. The obtained coordinates of estimated cells were used for fine readjustment of image linearization and the three dimensional structural analysis of the whole cochlea.

The first several steps of cell detection protocol were similar to those for IHCs described above. In short, template matching was performed with a small template image of OHCs (image size: 7 × 7 pixels) to create the matrix of correlation coefficients. The connected-component labeling was carried out on the binary matrix for the positions of the correlation peaks. The groups of peaks were taken as the cell candidates and narrowed down through two steps of prediction by two machine learning models. Ensemble learning methods were employed to build these models with similar predictors as described above.

For accuracy improvement, several steps were added after the second prediction step. As OHCs are typically distributed in three rows, the selected cell candidates were classified into three rows by a multiclass classification model (Appendix 2—figure 5). The model was built based on the convolutional neural network method (Neural Network Toolbox of MATLAB) (LeCun et al., 1989). The input of the model was a small image centered on the estimated cells. To reduce the chance of missing OHCs, regularity of distances between adjacent cell candidates in row was investigated. If the distance between two adjacent cell candidates in the same row exceeds 1.5 times of average distance, presence of additional cells in the gap was judged by the fourth model.

Appendix 2—figure 5

Download asset Open asset

To evaluate the missing cells in the space, the forth model was used. This model is also built based on the convolutional neural network method and uses small images as an input. Input small images for machine learning were sampled by placing small rectangular areas as region of interest (ROI) with equal distance from each other within the gap between preexisting cell candidates. The number of ROI per space was set to be equal to the integer portion of the quotient with the gap length as a numerator and the average distance of adjacent cells as a denominator. In the case that the model predicted the existence of a cell at the space, the nearest peak of correlation coefficient by the first template matching was added to the list of estimated cells. The obtained coordinates of estimated cells were used for the analysis of cell loss.

Our machine learning based auto-detection system utilizes two principles. The first principle is sorting hair cells from background noise. The models were trained with labeled dataset containing multiple feature values of cell candidates, such as volume, signal intensity, similarity with template image, coordinates, and relative position from neighboring candidates (Table 1). As morphology of hair cells and imaging conditions vary from place to place within the image, it was almost impossible to find the single or a few feature values that could efficiently distinguish cells from background noise. The approach with machine learning models, however, dramatically improved the accuracy of cell detection (Appendix 2—figure 6). Incorporation of excessive feature values may introduce a risk of overfitting. It is important to select machine learning algorithms resistant to the overfitting problem. After evaluation of available supervised learning algorithms in MATLAB (‘Statistics and Machine Learning Toolbox’ in MATLAB R2017b; Decision trees, Discriminant analysis, Nearest neighbors, Naïve Bayes, Support vector machines, and Classification ensembles), we chose classification ensembles (Gentle AdaBoost and Random forest algorithms) because of their high generalization abilities and processing speeds. Training dataset from ten cochlear samples and test dataset from other four samples were used.

The second principle is recovery of false negatives based on auto-recognition of three-row arrangement of outer hair cells (Appendix 2—figure 5). Although the detection based on the first principle shows high specificity of hair cell detection, it sometimes misses cells with atypical morphology or features. The first step of the false negative recovery is grouping of each cell candidate into one of three rows. Then, spatial gaps are detected by checking irregularity of intervals of cell candidates within each row. Finally, the gaps were evaluated whether there were overlooking cells or not by another model. This step of false negative recovery was realized by a convolutional neural network algorithm. The advantage is that it only requires an image as an input without any additional feature extraction. We expected the complementary effect by combining the first and second principles. Indeed, the recovery rate of the detection was further improved by adding the second process (Table 2, Paired t-test, p < 0.005, t = 3.94, df = 9).

We tested another auto-detection method for linearized image implemented in ‘Classic Watershed’ plugin for Fiji (Soille and Vincent, 1990). Several settings should be configured for this approach. First, a threshold value for voxel intensity should be provided. With lower threshold, more cells will be detected, but false positives will also increase. A threshold value for volume of the segmented region is also required for reducing false positives. In addition, inner and outer hair cells were distinguished based on their center coordinates along the radial axis in the linearized image. We optimized these parameters by grid search on the test samples for maximizing the detection efficiency (F score). Fiji macro and MATLAB source code of this method is available on GitHub (Iida, 2018b; copy archived at https://github.com/elifesciences-publications/Watershed).

The models were trained with the dataset obtained from ten cochleae with manual labeling (PND60: two samples, ACL: five samples, NCL: three samples). The numbers of observations of the training datasets are shown in Table 4. Details of the predictors of datasets are shown in Table 1. The datasets for the models IHC1 and OHC1 include all the peak groups obtained from the ten linearized images with template matching (Appendix 2—figure 4). These models are based on a GentleBoost algorithm. To deal with the imbalance between positives and negatives, we set the cost matrix when training IHC1 and OHC1 (100 times and 50 times higher costs for false negatives than false positives, respectively). The hyperparameter was optimized by five-fold cross validation with a Bayesian optimization method on the training dataset (number of the trees, learn rate; ‘OptimizeHyperparameters’ option for ‘fitcensemble’ function of MATLAB).

A subgroup of cell candidates was selected for further processing with the models ICH2 and OHC2, based on their higher posterior probability computed by IHC1 and OHC1. The fraction of this subgroup over the total sample population was fixed. The models ICH2 and OHC2 were based on a Random Forest algorithm. We set the cost matrix to deal with the imbalance of the dataset when training IHC2 and OHC2 (2 times and 0.4 times costs for false negatives, respectively). The hyperparameter was optimized as with IHC1 and OHC1.

The training dataset for the model OHC3 includes small images centered on each outer hair cell extracted from ten linearized images. This model is based on a CNN algorithm. We used a shallow network with typical settings for training. Details of the layer configuration (filter size and number of convolutional layer) were optimized by five-fold cross validation with a Bayesian optimization method on the training dataset. We trained the network with stochastic gradient descent with momentum algorithm, with the maximum number of epochs of 15, and the initial learning rate of 0.001 (Murphy, 2012). We used the default settings of ‘trainingOptions’ function of MATLAB for the other training options.

The dataset for the model OHC4 includes small images of gaps between cell candidates within each row of outer hair cells (Appendix 2—figure 5). The model was based on a shallow CNN with typical settings for training. We optimized the details of the layer configuration of the model and trained it in the same way as OHC3.

The models were tested by ten cochleae that were not used in the training (PND30: two sample, PND60: three sample, ACL: two sample, NCL: three sample). The way of making the test datasets was the same as the training dataset. The numbers of the test dataset of each model are shown in Table 4. The indices of the performance were obtained as follows:

where TP, FP and FN mean True Positive, False Positive and False Negative, respectively. For multiclass classification, the micro-averages of recall and precision were calculated (Sokolova and Lapalme, 2009).

In the previous section, the longitudinal axis (modiolus axis) was set for each sample and utilized for the presentation of the IHC spiral in the cylindrical coordinate system. Registration of multiple samples were performed using these parameters. To reduce sampling bias in the alignment process of multiple samples, we followed the strategy of first creating a generic template of the IHC spiral by averaging all the available samples. After this step, alignment of individual samples to the generic template was performed. Nevertheless, in both cases we utilized the common protocol for alignment. To simplify the following explanation, we refer two spirals to be aligned as spirals A and B (Appendix 2—figure 10).

Appendix 2—figure 10

Download asset Open asset

We first performed alignment between spirals A and B by rotating objects around the longitudinal axis. This can be mathematically achieved by shifting the value $\phi$ of one spiral against the other. To estimate the amount of shift in $\phi$, a two-dimensional plot of $p$against $\phi$ was created and overlaid for spirals A and B. By evaluating the extent of correlation between two curves in the plane of $\phi$-$p$, the optimized shift in $\phi$ can be estimated. To calculate the extent of correlation, two curves were converted to series of points with their angular intervals of 0.1 rad and their cross-correlation was calculated (MATLAB “normxcorr2”). The shift in $\phi$ that gives the highest cross-correlation was selected and the rotational shift between two curves A and B were corrected.

The last step of alignment between spiral A and B was translation along the longitudinal axis. To this end, the spiral A was fixed and the origin of the second spiral B was shifted along the longitudinal axis. Comparison of the two spirals was based on the search for the minimum sum of the squared distances along the longitudinal axis between points from two spirals. These points were selected within the region of overlapping azimuth for the two spirals. Two corresponding points on different spirals share the same $\phi$. The origin of the second spiral was determined by finding the arguments of the minimum for the function of squared distance between corresponding points for two spirals projected to the longitudinal axis as follows:

where $P=\left\{p_1,p_2,\dots \right\}$ is a set of positions on the longitudinal axis $l$, and $\Phi$ is a set of azimuth $\phi$ shared by the points of comparison between spirals A and B. The function $z_o\left(\phi \right)$ denotes an axial coordinate of the point on the spiral A. The function $z_i\left(\phi \right)$ denotes an axial coordinate of the point on the spiral B when the origin of the spiral was set at the position $p_i$. The values 'Distance from Base (μm)' (Figure 2 and Figure 3 in the main text) and 'Cell Number' (Figure 2 in the main text) were also plotted after alignment at the midpoint of the spiral segments.

The outline of the procedure was described in the main text. Loss of OHCs results in the formation of empty spaces in the epithelium. Conventionally, the lost cell number was estimated by the size of the empty space. This method is not reliable when the empty space is large, the rows of OHCs are difficult to define, or cell-to-cell distances vary in different epithelial positions (Appendix 2—figure 11).

Appendix 2—figure 11

Download asset Open asset

In this study, we evaluated the space where cells were lost without assuming OHC rows. First, the x and y coordinates of the centers of detected cells (hereinafter called ‘cell centers’) were measured in the linearized image and created a scatter plot. Second, we obtained parameters necessary for radial alignment (along y-axis) of cell centers. For this purpose, we searched the longitudinal positions (along x-axis) along the organ of Corti, where a zone with the width of 8 μm (a gray zone in Appendix 2—figure 12) contained more than two cell centers (three black points within the gray zone in Appendix 2—figure 12). When the zone was found, two values, the averaged x coordinates (red circle in Appendix 2—figure 12) and the largest and smallest y coordinates (red line in Appendix 2—figure 12) of cell centers, were calculated for this zone. The former value was used as a reference of the center position of the epithelium and the latter value as a reference of its width (Wy). The vertical widths (subtraction of the largest and smallest y coordinates) were interpolated linearly along the x-axis and used for normalizing the y coordinates of cell centers. Third, we obtained parameters necessary for longitudinal alignment (along x-axis) of cell centers. For this purpose, we searched all the cell centers and created rectangles with their edge lengths of 24 μm along the x-axis and 8 μm along the y-axis with their centers aligned with the cell centers (a red point and a gray rectangle in Appendix 2—figure 13). Distances between the original cell center and the center of the nearest cell within the rectangle along x-axis were calculated (Wx). The averages were calculated with 100 μm intervals along the x-axis and interpolated linearly. Fourth, the coordinates of the cell centers were normalized using averaged Wx and Wy values to equalize the horizontal and vertical distances and to keep the distances uniform along the organ of Corti (Figure 2F in the main text).

Appendix 2—figure 12

Download asset Open asset

Appendix 2—figure 13

Download asset Open asset

A binary image of the normalized epithelium was created based on the equalized coordinates of cell centers (Figure 2G in the main text). The coordinates projected onto an image were adjusted to have the average distances between neighbors in x and y as five pixels. The horizontal center line of the image was set to be on the line y = 0. The height of image was set to 15 pixels and the width was adjusted to the range of x coordinates. Then squares of 5 × 5 pixels centered on each cell point were drawn on the image. Small holes were removed by a morphological closing operation. The empty spaces in the image were considered to be the putative cell loss sites.

The estimated amount of cell loss in the entire organ of Corti or in specific areas was shown as either the number of pixels in the empty spaces (Figure 3D and Figure 3—figure supplement 1B) or the cell number obtained by the formula $\mathrm{n}=\raisebox{1ex}{$v$}\!\left/ \!\raisebox{-1ex}{$25$}\right.$, where $\mathrm{n}$ is the estimated number of cell loss and $v$ is the area of each connected region. The estimated number of cell loss was rounded off to the nearest integer (Figure 3A,B and Figure 3—figure supplement 1A,C). For the simulation analysis of cell loss models, the same formula was used for the estimated number of cell loss (Figure 4 and Figure 4—figure supplement 1).

The MATLAB scripts can be viewed on the website of GitHub without installation of MATLAB environment. The scripts have modular architecture, and the inputs and outputs of each module are annotated in the scripts. It would be relatively easy to transfer to the language which can preserve the architecture, such as C/C ++ or Python, by transferring the modules one by one. As there are several built-in functions specific to MATLAB (ex. ‘findpeaks’), however, alternative means would be needed in some cases. Please refer to the website of MathWorks for the syntax and built-in functions of MATLAB in such cases (https://www.mathworks.com). In our scripts, the names of modules we made are written in upper camel case. The variables are written in lower camel case, and the constants in all capital letters. The names of built-in functions are written in lower case. Regarding how to train machine learning models, please refer to the ‘Training of machine learning models’ section above. Feel free to contact the corresponding author in case you have questions regarding the contents of the code. We would consider the transfer to other programming languages as needed.

• Chisato Fujimoto

• Chisato Fujimoto

• Tatsuya Yamasoba

• Tatsuya Yamasoba

• Shigeo Okabe

• Shigeo Okabe

• Shigeo Okabe

• Shigeo Okabe

• Shigeo Okabe