Quantifying Intermembrane Distances with Serial Image Dilations

Statistical Methods: Comparisons were made between experimental groups using Student's t-tests. A p-value < 0.05 is considered significant and all values are represented as mean ± standard deviation.

Manual Segmentation. The quantification of the GJ-adjacent perinexus nanodomain width (W_p) is typically accomplished by manual segmentation. This manual segmentation process is demonstrated in Figure 1A and was described previously⁶. The observer identifies the edge of the GJ (Figure 1, red dot), measures 5nm along the center of the perinexus, and measures the distance between the membranes at that point. The process is then repeated at 10, 15, 30, and every 15 nm up to 150 nm. This technique, while effective, has limitations of time and spatial under-sampling along the length of the perinexus.

Mean W_p measurements from previous studies can vary from approximately 10 to 20 nm^2,3,10, and 3 nm appears to be the mean difference needed to detect statistical significance, which is well above the spatial Nyquist frequency of 0.7 nm per measurement based an interpixel resolution of 0.34 nm. Therefore, while manual segmentation is time consuming, the method is sufficient to measure differences in W_p associated with an intervention or disease state.

Serial Image Dilations. In order to measure the perinexus in a faster, reproducible manner with appropriate spatial resolution, we developed a program based on serial image dilations to count the pixels between two manually-traced membranes, which can be seen in Figure 1B.

The serial dilation process is illustrated in Figure 4. As the binary image is dilated (Figure 4A-4D), that dilation is then inverted and added to a working image – the non-binary form of the original outline (Figure 4E-4H). The process is repeated until the outline is entirely filled in (Figure 4D). At this point, the final working image (Figure 4H) is a count of the number of times a particular pixel has remained un-dilated. As such, the values near the outline of the cell membranes are very low, while the values at the center are highest. By counting the number of dilations to fill the area at each point, the distance between the membrane edges can be calculated. The next challenge is to identify and isolate the centerline in order to quantify perinexal width as a function of distance from the GJ, which is done by first applying a spatial derivative to the final working image (Figure 2-last image, and Figure 5A). A second example of a more irregularly-shaped perinexus can be found in supplemental file S3.

Centerline Identification. The gradient of the final working image can be quantified by a spatial derivative, as dilation count values from edge to edge change from high to low to high again (Figure 5A left to right). Considering only the magnitude of the spatial derivative (Figure 5B), the outline and centerline, highlighted with white arrows, are immediately identifiable as areas of discontinuity. At these locations, the gradient direction changes from increasing to decreasing or vice versa. Applying a threshold (Figure 5C) produces a binary image of the centerline and outline, and subtracting the original outline yields the isolated centerline (Figure 5D). While this method of isolating the centerline is computationally efficient, the threshold applied to the spatial derivative creates gaps in the resulting centerline. These gaps (Figure 5D, insert) must be filled to provide an accurate measurement of the distance from the GJ and to ensure the perinexus is measured in its entirety. First, the centerline is dilated to fill in any gaps (Figure 5E), followed by an erosion (Figure 5F) and a "bwmorph" function (operation = 'skel', n = inf) to eliminate as many points as possible while leaving a continuous centerline, thereby increasing the computational efficiency of a subsequent pathfinding algorithm developed by Wasit Limprasert and available on MATLAB Central⁹. This dilation-erosion function produces the completed centerline, which is combined with the final working image (Figure 5G). However, this centerline is often more than one pixel thick and therefore is not a precise isolation of the centerline.

The Wasit Limprasert pathfinding algorithm is used to determine the perinexus centerline. The pathfinding algorithm is able to track the highest values – in this case the values closest to the center which remained un-dilated through the most iterations along the centerline (Figure 5G, insert). The result is an automatic trace of the centerline, as shown in Figure 6. By isolating the centerline, perinexal width can be presented as a function of distance from the end of the GJ, as shown in Figure 6B (top), or as the average width of a specified region of interest.

Kernel Analysis. It is important to note that digitized images are based on square arrays and dilation kernels are likewise based on square matrices. This means that dilation distance across a diagonal is greater than orthogonal. Therefore, we next sought to determine whether the kernel affected the results of the algorithm. In order to quantify kernel-specific variability, five different kernel shapes were analyzed: "Plus" (the shape used in the above analyses), "X", "Box", and "Line", as described in Figure 7A. The kernel is applied at each non-zero point of a binary image. The star in each kernel of Figure 7A represents the center, where white is a value of 1, and black is a value of 0 for the dilation kernel.

Each kernel's influence on the mean W_p measurement of a single approximately horizontal perinexal image (Figure 7B, top), quantified by an experienced user, was determined by rotating the image with Matlab's "imrotate" command and computing W_p in steps of 10°. The W_p measurement values (Figure 7B, bottom) fluctuate with image orientation in a rectified sinusoidal fashion with a Plus shaped kernel. The lowest values occur when a relatively straight perinexus is oriented vertically or horizontally. Neither the X, Box, nor Line kernels provided any advantage over the plus-shaped kernel. The X and Box kernels produced identical results, but the values of mean W_p were out of phase with th_e Plus kernel by 45°. The Line kernel failed to fully dilate the image at certain angles as can be seen by the absence of data in the green trace for images rotated less than 30 or more than 145°. Thus, the orthogonal Plus dilation kernel overestimated membrane separation when dilating a perinexus with an axis oriented diagonally for example at approximately 45°, and the X and Box kernels underestimated mean W_p when the long axis of the perinexus was also at 45°. Based on this analysis, we developed a correction factor applied to the values generated from dilating with the plus-shaped kernel. To account for the overestimation of membrane separation associated with image orientation, this correction factor multiplied by the measured width value depending on the orientation of the image (Equation 1).

If θ <45°
     W_{p corrected} = cosd(θ) * W_{p measured}
If θ ≥45°
     W_{p corrected} = cosd(θ) * W_{p measured}     (Equation 1)

In this equation, W_p measured is the original W_p value generated by the above algorithm and θ is the calculated angle from horizontal, in degrees. θ is calculated by taking the inverse tangent of the total change in the horizontal direction divided by the total change in the vertical direction of the perinexal centerline. The above correction approximates the average angle, relative to horizontal, of the perinexus (Figure 8A, left-top) and results in a measurement as if obtained from a horizontal perinexus (Figure 8A, left-bottom). The rationale behind this equation comes from the fact that the plus-shaped kernel (Figure 7A) is essentially two line-shaped kernels arranged orthogonally to each other. As such, below 45° (closer to horizontal), the dilations occur vertically and therefore multiplying by the cosine of the angle gives the correct measurement. Conversely, for angles above 45° (closer to vertical), dilations occur horizontally and the sine is used to determine the correct measurement. At precisely 45°, the sine and cosine are equal. Supplemental file S4 provides a depiction of this concept. Note that this correction is based on the average angle and caution should be used when analyzing substantially non-linear shapes. This process was repeated on 20 randomly-selected perinexi and the corrected measurements correlated strongly with measurements obtained by manually rotating and re-analyzing the images (Figure 8A, right). To confirm the accurate correction for image orientation, two sets of phantom edges were generated (Figure 8B, left) and rotated 180°. With the trigonometric correction, the algorithm accurately returned the correct value at each orientation, regardless of spatial resolution or image size (Figure 8B, right).

Analytical Application and Reproducibility with Orientation Correction. Recalling that previous studies using manual segmentation report statistically significant mean W_p differences greater than or equal to 3 nm, it was important to determine whether to the algorithm could be used to recapitulate previous findings using a complete dataset. Using the new algorithm, two observers – one experienced and one unexperienced with perinexal analysis (Obs. 1 and Obs. 2, respectively) – analyzed the same images from a previous study⁶ that included 12 patients who were diagnosed with atrial fibrillation (AF) prior to tissue collection and 29 patients that did not have pre-existing AF (No-AF). The experienced user found that W_p was significantly wider in patients with AF than without AF (21.9±2.5 and 18.4±2.0 nm, respectively, Figure 9A). These values with the correction factor applied are similar to those reported previously (24.4±2.2 nm and 20.7±2.4 nm, respectively)⁶. Importantly, the inexperienced user found the same significant difference (22.1±2.8 nm and 20.1±2.6 nm, respectively) between disease states with the automated program. Additionally, the standard deviation of the W_p values did not change with the correction factor, indicating the standard deviation of 2-3 nm is not an artifact of the algorithm but of the structure itself and tissue processing. These results demonstrate that the proposed automated method is capable of recapitulating the results of previous studies.

Importantly, the perinexus is a recently-defined structure and no consensus has been reached on the range of absolute values of membrane separation adjacent to the GJ^2,3. Since outer-membrane-to-outer-membrane GJ width has been previously estimated at 20 nm¹³, we sought to determine the algorithm's efficacy by also measuring GJ width. Both observers found no significant difference between gap junction widths (GJW) of patients with or without pre-existing AF (Figure 9B). Absolute GJW values for AF and non-AF patients were 20.5 ± 2.5 nm and 20.3 ± 1.9 nm, respectively, for the experienced observer and 21.0 ± 3.1 nm and 20.0 ± 2.2 nm for the inexperienced observer, similar to what has been reported previously.

To determine whether the automated algorithm required less time to analyze data than manual segmentation, both the experienced and inexperienced user recorded the time needed to quantify a 10-image training set (Supplemental file S5). Table 1 demonstrates that the experienced and inexperienced user decreased analysis time by 4.7- and 8.3-fold respectively using the automated algorithm relative to the manual segmentation approach, with an approximately 43-fold increase in spatial resolution along the perinexus.

Algorithm Troubleshooting. The most common error when running the algorithm occurs when the final centerline does not end at the edge of the image. In such cases, not enough points were selected from the spatial derivative map, causing the program to fail and produce an error message advising the user to select a larger crop area or increase the spatial derivative threshold. Drawing a larger crop box will improve the reliability of the program in some cases as the spatial derivative changes drastically near the edges of the figure, which can disrupt the pathfinding or edge detection algorithms.

It is also possible for the pathfinding algorithm to fail to properly identify the centerline, even if the centerline reaches the edge of the image, particularly if the gradient threshold is too low (Figure 2A). If the gradient threshold is set too high, there will be more unnecessary points incorporated into the pathfinding algorithm (Figure 2B), decreasing computational efficiency. If the user is unable to determine an appropriate threshold, the image array "GMag" (Figure 2C), which is generated by the program and can be found in the Workspace, can help the user determine the threshold. Find points along the centerline and set the threshold slightly above their Index value to ensure these points are selected. In the given example, an appropriate threshold would be above ~5.1 (Figure 2C, insert).

The start point may also fail to reach the beginning of the perinexus (Figure 2D). In this case, re-run the program and set the Manual Start value to 1. After the centerline has been isolated, the user selects a point outside the perinexus and the centerline point closest to the selected pixel (Figure 2E, red square) will be set as the start point. The result is the full centerline (Figure 2F).

Figure 1: TEM images with quantification processes. The manual segmentation process (A) requires the user to perform 12 individual intermembrane measurements while estimating the centerline. The automatic process (B) requires a manual, continuous trace of the outline of the perinexus. The red dot in each image represents the user-identified end of the GJ and the beginning of the perinexus. Please click here to view a larger version of this figure.

Figure 2: Centerline troubleshooting. Two primary modes of failure for the centerline identification and their solutions: Each image is labelled with the array name in Matlab. If the gradient threshold is too low (A, threshold 0.2) the centerline algorithm can fail. Setting the threshold too high (B, threshold 70) can reduce the computational efficiency of the pathfinding algorithm. An appropriate gradient threshold can be determined from the GMag array (C, insert). If the centerline fails to reach the beginning edge of the perinexus (D), the user can opt to manually select the start point. After the "Start Point" option is set to 1 in the opening GUI, the user then selects a point outside the perinexus of interest (E). The end result should be a centerline that accurately depicts the entirety of the perinexus (F). All labels in quotations (A-F) correspond to the variable names in Matlab. Please click here to view a larger version of this figure.

Figure 3: Perinexus Selection. To crop the perinexus, click-and-hold to drag a box around it (the crop tool is selected automatically) as shown by the blue arrow. This box can be adjusted by using the squares on the sides and corners to make it larger or smaller. The green arrow represents the end of the perinexus, which the user should ensure remains "open." Please click here to view a larger version of this figure.

Figure 4: Serial image dilations. The binary outline is repeatedly dilated in one-pixel increments (A-D) and added to the working image (the non-binary form of the image, E-H) after each dilation. Please click here to view a larger version of this figure.

Figure 5: Centerline isolation and pathfinding. A spatial derivative is calculated from the final working image (A) and the magnitude of that spatial derivative (B) is used to isolate the outline and centerline (white arrows). A user-defined threshold identifies the outline and centerline and subtracting the original outline yields the centerline (D). However, gaps appear in the centerline as a result of the threshold (D – insert). In order to produce a continuous centerline, a secondary dilation is performed on the isolated centerline (E), followed by a secondary erosion to increase the computational efficiency of a subsequent pathfinding algorithm. This eroded image (F) is then combined with the final working image, allowing for the identification of a continuous, one-pixel-thick centerline (G – insert). Please click here to view a larger version of this figure.

Figure 6: Final data presentation. The program outputs the final outline on top of the original TEM image (A). The line is color-coded green for 0-150 nm, blue for the user-defined region of interest and red for the area beyond 150 nm. Additionally, the program outputs a similarly color-coded chart representing W_p as a function of distance from the start of the perinexus and for the region of interest (B), as well as the average W_p (inset in respective graph). Please click here to view a larger version of this figure.

Figure 7: Dilation kernel shape analysis. Shapes for dilation kernels (A): the star in the middle represents the dilated pixels. White boxes are pixels affected by the dilation, in the shape of a plus, X, box or line. An approximately horizontal perinexus (B, top, with a red line indicating 0°) was rotated clockwise from 0 to 180° in steps of 10° and repeatedly dilated using the different kernel shapes (B, bottom). Plus- and line-shaped kernels produce similar results, although the line kernel fails at certain orientations, while the box- and X-shaped kernels are out of phase by 45°. Please click here to view a larger version of this figure.

Figure 8: Image Orientation Correction. To correct for the orientation of the image, the average angle of the perinexus is calculated from the positions of the start- and end-points (A, left, beginning and end of black line). The inverse tangent of the change in the y-direction (A, left, green line) divided by the change in the x-direction (A, left, red line) yields the correction angle, θ (A, left, yellow). The goal is to then correct the average W_p value to give the minimum distance between the edges, as if the perinexus had been imaged approximately horizontally (A, left below). The application of the correction factor described by Equation 1 was compared to manually rotating each of 15 randomly-selected images by the calculated θ before analysis. The corrected values strongly correlated with the rotated-image values (R² = 0.991, A, Right), indicating Equation 1 is a valid correction factor for image orientation. To confirm the correction factor is appropriate, two phantoms were generated of perfectly-parallel edges with a known distance between them (B, Left). Phantom 1 and Phantom 2 have spatial resolutions of 2.833 pixels/mm. and 71.6 pixels/in., respectively. As shown by the blue diamonds and red squares in B, Right, the algorithm accurately calculates their widths across 180 degrees of image rotation. Please click here to view a larger version of this figure.

Figure 9: Algorithm reproducibility. Using the automatic process with image orientation correction, both an experienced and inexperienced observer found significant differences between AF and No-AF groups (A), consistent with a previous study, detecting a minimum difference of 2.6 nm. Additionally, neither observer found a significant difference in GJW (B). Please click here to view a larger version of this figure.

Table 1: Comparison of manual and automatic processes. Both observers required less time per image to trace the outline than to perform the manual segmentation process for a 10-image training set. Additionally, the automatic process has a higher sampling frequency, recording 3.45 measurements per nm, compared to an average of 1 measurement every 12.5 nm for the manual process. The training set images can be found in Supplemental File S5, along with outlines and measurements as performed by an experienced user.