|Home | About | Journals | Submit | Contact Us | Français|
Spinal cord injury (SCI) is a debilitating event with multiple mechanisms of degeneration leading to life-long paralysis. Biomaterial strategies, including bridges that span the injury and provide a pathway to reconnect severed regions of the spinal cord, can promote partial restoration of motor function following SCI. Axon growth through the bridge is essential to characterizing regeneration, as recovery can occur via other mechanisms such as plasticity. Quantitative analysis of axons by manual counting of histological sections can be slow, which can limit the number of bridge designs evaluated. In this study, we report a semi-automated process to resolve axon numbers in histological sections, which allows for efficient analysis of large data sets.
Axon numbers were estimated in SCI cross-sections from animals implanted with poly(lactide co-glycolide) (PLG) bridges with multiple channels for guiding axons. Immunofluorescence images of histological sections were filtered using a Hessian-based approach prior to threshold detection to improve the signal-to-noise ratio and filter out background staining associated with PLG polymer.
Semi-automated counting successfully recapitulated average axon densities and myelination in a blinded PLG bridge implantation study.
Axon counts obtained with the semi-automated technique correlated well with manual axon counts from blinded independent observers across sections with a wide range of total axons.
This semi-automated detection of Hessian-filtered axons provides an accurate and significantly faster alternative to manual counting of axons for quantitative analysis of regeneration following SCI.
Spinal cord injury (SCI) is a traumatic event that affects approximately 12,000 new individuals in the United States each year (Devivo, 2012; National Spinal Cord Injury Statistical, 2013). Ascending and descending axon tracts are disrupted during SCI resulting in loss of sensory and motor function (McDonald and Sadowsky, 2002). Many pre-clinical strategies, including peripheral nerve and cell grafts, have been developed to reconnect severed regions of the cord after SCI (Cheng et al., 1996; David and Aguayo, 1981; Lu et al., 2012; Thuret et al., 2006). Biomaterial bridges have been shown to guide axons across the injury site (Hurtado et al., 2006; King et al., 2003; Moore et al., 2006; Novikov et al., 2002; Pawar et al., 2015; Stokols et al., 2006; Thomas et al., 2013; Tuinstra et al., 2014; Wong et al., 2008) and provide modest improvement in locomotion when coupled with cell transplantation or growth factor delivery (Han et al., 2009; Teng et al., 2002; Tsai et al., 2006). The extent of regeneration through the bridge, however, is often not characterized and recovery may be due to indirect mechanisms including modulation of the immune response, spinal plasticity, and altered cell migration. To examine the role of regeneration in recovery, and to compare different studies, quantitative analysis of axon regrowth through the injury site is needed.
Manual counting of axons in serial sections through the injury site can be time-intensive. Assessing myelination of regenerating fibers requires additional counting and coordination of multi-color images to determine co-localization of axonal and myelin staining. The time necessary to obtain these data can hinder larger studies, such as combinatorial strategies where several factors are tested simultaneously. Automated image analysis has the potential to expedite this process and allow for large, multi-variable studies. Few automated techniques, however, have been developed for quantification of axons or myelination of regenerating axons in the injured spinal cord. Axon densities have been estimated using line profiles, though only features that intersect with the line profile region are analyzed (Grider et al., 2006; Sathyanesan et al., 2012). Variation in axon density across two dimensions in a section cannot be well captured. Regeneration following SCI is often inhomogeneous, with regions of high and low axon density within the same tissue section, requiring methods to approximate total axon numbers rather than densities of axons. Another major obstacle to automated analysis of regeneration is inconsistency in fluorescence due to non-specific background staining and autofluorescence from implanted materials, such as bridges, that limits the use of simple threshold techniques for object detection and counting.
In this study, we describe a semi-automated approach to estimate axon numbers and numbers of remyelinated axons that is consistent with counts performed by blinded observers, yet can be performed in significantly less time. This analysis is applied to quantify axon growth through PLG spinal cord bridges that contain multiple channels to promote directed regrowth through a lateral hemisection model of spinal cord injury (Thomas et al., 2013; Tuinstra et al., 2014; Yang et al., 2009). Histological sections were stained for neurofilament, and images were captured and filtered using a Hessian-based technique to overcome low signal-to-noise ratios. Importantly, the eigenvalues of the Hessian matrix are used to detect curvilinear structures within the image. Similar techniques have been developed to improve the contrast of blood vessels in magnetic resonance imaging for two- and three-dimensional morphological analysis (Frangi et al., 1998; Sato et al., 1998). Hessian-based filtering followed by object detection provides a time-saving tool for semi-automated analysis of axon regeneration and myelination following SCI.
Porous multiple channel bridges were fabricated from PLG (75: 25 ratio of D,L-lactide to L-glycolide; inherent viscosity: 0.76 dL g−1; Lakeshore Biomaterials; Birmingham, AL, USA) and 63-106 μm salt as previously described (Thomas et al., 2013). Briefly, sugar fibers were drawn from 220 °C caramelized blend of dextran (MW ~100,000; Sigma Aldrich; St. Louis, MO, USA), glucose (Sigma Aldrich), and sucrose (Fisher Scientific; Pittsburgh, PA, USA) in a 1/2.5/5.3 ratio by mass. Seven fibers were dipped into a 2.5:1 mixture by mass of PLG microspheres to salt, packed into a custom mold, and then equilibrated under 800 psi CO2 for 16 hours. The CO2 was released over 30 minutes, foaming the bridge into its final structure. The bridges were then cut to 2 mm in length, leached for 4 hours in ultrapure water, dried, and stored in a desiccator until use. Prior to implantation, bridges were immersed in 70% ethanol and then ultrapure water for 30 seconds each.
Female C57Bl/6 mice (age 8-10 weeks) were used for spinal cord injuries studies in accordance with the ACUC at Northwestern University. Animals were anaesthetized using 2% isoflurane and the surgery site was scrubbed by repeated betadine and ethanol exposure. A single incision was made in the skin to expose the underlying muscle. The muscle tissue next to the spinal column was cut from T7-T11. A double laminectomy was performed from T9-T10. Two lateral hemisection incisions were created between T9 and T10, approximately 2 mm apart. The tissue between each incision was cleared to create a gap hemisection injury. A PLG spinal cord bridge was then implanted into the injury site and covered with gel foam (Henry Schein; Dublin, OH, USA). The musculature was then sutured closed and the skin closed using surgical clips. Animals were expressed twice daily and given enrofloxacin (Baytril, Henry Schein, 2.5 mg/kg per day) for two weeks following injury. For analgesia, buprenorphine (0.1 mg/kg) was administered every 12 hours subcutaneously for 60 hours following surgical procedures.
At 8 weeks post injury, animals were perfused and the spinal cord tissue from T7-T12 was harvested. Spinal cords were embedded and 12 μm cryostat sections were obtained in the cross-sectional plane. Sections were batch stained for each analysis. Sections were post-fixed for 10 min in 4% paraformaldehyde (Sigma) followed by permeabilization in 0.5% triton-X (Sigma) in Trizma solution (Sigma) for 15 min. Sections were then blocked in 10% normal donkey serum (Jackson ImmunoResearch; West Grove, PA, USA) and 5% bovine serum albumin (Sigma) in 0.1% triton in Trizma solution for 1 hour. Sections were incubated with the following primary antibodies overnight: rabbit anti-neurofilament (NF200, 1:200, Sigma) and goat anti-myelin basic protein (MBP; 1:500; Santa Cruz, Dallas, TX, USA). Sections were washed 3 times, incubated with the donkey anti-rabbit Alexa Fluor 555 and donkey anti-goat Alexa Fluor 647 secondary antibodies (1:1000; Life Technologies, Grand Island, NY, USA) for 1 hour, and then mounted with coverslips using Fluoromount G (Electron Microscopy Sciences, Hatfield, PA, USA). A series of 200× images encompassing the entire cross-section were obtained for NF200 and MBP stains on a widefield fluorescent microscope (DMIRB; Leica, Buffalo Grove, IL, USA) with standard fluorescence filter cube sets, a 0.4 numerical aperture 20× objective, and a 14-bit CCD camera (CoolSnap HQ-2; Photometrics, Tucson, AZ, USA) resulting in an approximate pixel size of 0.67μm2. Images for each stain corresponding to the same 200× region were converted to 8-bit and merged first prior to stitching. Merged images from the same spinal cord cross-section were then stitched in Photoshop (Adobe, San Jose, CA, USA), resulting in a single composite image of NF200 and MBP staining in the entire section. All sections were stained and imaged under the same conditions.
Prior to importing into MATLAB (Mathworks, Natick, MA, USA), stitched composite images were normalized in Photoshop to utilize the full dynamic range of intensity. The maximum intensity value found in all images was determined and set to a value of 230 which is approximately 90% of the maximum potential value in the 8-bit images. Intensity values for all pixels in all images were scaled accordingly. For manual counting, scaled images were counted by two blinded observers in ImageJ (Schneider et al., 2012) using the cell counter plugin. For automated analysis, scaled images were imported into MATLAB (imread). The area of the section corresponding to the PLG bridge was defined using the polygon region of interested tool (roipoly). A Hessian matrix of the NF200 channel was created by convolution filtering using a 7×7 kernel of the second derivative of the Gaussian function in the x-direction (Gxx), y-direction (Gyy) and xy-direction (Gxy,Gyx) with a scaling factor (σ) of 1. Eigenvalues were determined by taking the determinant of the Hessian matrix consisting of the filtered images in each direction (Gxx, Gyy, Gxy, Gyx). Eigenvalues were then sorted based on their absolute value and the original image was then restructured according to the equation published by Frangi and colleagues (Frangi et al., 1998). The parameters β and c were set to 0.5 and 15 respectively for all analyses on NF200 staining, consistent with previously published values (Frangi et al., 1998). Following filtering, images were counted manually in ImageJ or positive NF200 events were identified by intensity threshold using the im2bw function in Matlab (threshold value = 0.1), single pixel events were removed using the ‘clean’ operation of the bwmorph function, and the remaining events were counted using the bwconncomp function, which determines the number of connected objects identified following intensity threshold. To obtain axon densities, total NF200 counts were divided by the area outlined with the polygon region of interest tool in Matlab.
To determine the extent of myelination of regenerating axons within the PLG bridge, overlapping images of MBP and NF200 staining were examined concurrently. Color channels containing MBP and NF200 fluorescence individually were separated and filtered as described in the previous section. The values for β and c were set to 0.5 and 15 respectively for the NF200 channel, and 0.5 and 10 for the MBP channel. Following intensity threshold (threshold value = 0.1, im2bw), NF200 objects were identified (bwconncomp) and the pixel locations for each object were identified (PixelIdxList). Positive MBP staining was determined by intensity threshold (threshold value = 0.05, im2bw). NF200 objects containing pixel locations that overlapped with positive MBP staining were counted and compared to total NF200 counts.
Statistical significance was determined by one way analysis of variance (ANOVA) using Prism software (GraphPad, La Jolla, CA, USA). Correlation between manual and semi-automated counts was determined in MATLAB. Statistical significance was set at p< 0.05.
The ability for Hessian filtering to improve detection of regenerating axons and reduce background fluorescence was investigated within images of fluorescently labelled axons. Spinal cord cross sections were obtained at 8 weeks post-implantation of a PLG bridge into a lateral hemi-section model of spinal cord injury (Figure 1). Fluorescence images demonstrated bright NF200 antibody staining associated with regenerating axons, as well as background fluorescence likely associated with the PLG polymer (Figure 2A-B). Binary filtering of the unfiltered image by threshold intensity results in selection of many broad areas not specific to axons and grouping of several axons into one object leading to false-positives and false-negatives respectively, which could contribute to inaccurate assessment of axon numbers (Figure 2 C, grey arrows). The area of the spinal cord cross-section occupied by the PLG bridge was selected in MATLAB and convolution filtered using the second derivative of the Gaussian function to establish a Hessian matrix. The eigenvalues of the Hessian matrix were then used to filter the original image. Hessian-based filtering significantly improved the resolution (Figure 2D, white arrows) and detection (Figure 2E) of fluorescently labelled axons and appeared to reduce the detection of non-specific labelling (Figure 2E, grey arrows) compared to the corresponding unfiltered images shown in Fig 2B-C. Features selected by intensity threshold overwhelmingly localized with the axonal structures.
Hessian filtering was also performed using ImageJ through the FeatureJ plugin. The best results were obtained with the largest eigenvalue of the Hessian tensor, absolute eigenvalue comparison, and a smoothing factor of 0.5. Improved contrast in the axonal structures was observed, however, many of the structures were not fully resolved and the grey level of background pixels was higher than in images filtered with the Hessian-based approach in Matlab.
Axons regenerate primarily through the channels of the bridge; thus, axons are not homogeneously distributed within the tissue section and are located within specific regions, requiring the semi-automated analysis to be accurate in regions with variable axon densities. Furthermore, axons are present at many orientations with respect to the plane of the section resulting in different observed morphologies. Initial studies thus aimed to develop the technique by analyzing small regions within sections with various numbers, sizes, and orientations of axons (Figures 3A-B). Two blinded observers manually counted the number of axons in 175 regions from unfiltered images, and automated counting was performed in the same regions following Hessian-based filtering. Manual counts performed by the two blinded observers were found to correlate with a Pearson’s coefficient of 0.832. Automated counts correlated well with the average of the manual counts from the blinded observers (slope = 0.9861, Pearson’s r = 0.788, p < 0.001, Figure 3C). The log derivatives of manual and automated counts were used to obtain a normal distribution (Lillier’s test, p = 0.4343), and the Bland-Altman plot indicated that automated counting had good agreement with manual counting (Figure 3D). Automated counts also correlated well with manual counts of Hessian filtered images of the same regions (slope = 0.996, Pearson’s r = 0.823, p <0.001, Figure 3E-F). While manual counting took up to 1 hour to complete, automated counting was significantly faster (approximately 10 seconds) and worked well for approximating the number of regenerating axons compared to manual counting over a wide range of axon numbers, sizes, and orientations.
The ability for the semi-automated counting technique to recapitulate axon numbers across multiple animals in a spinal cord injury study was investigated. Animals with bridges implanted into a lateral hemisection were randomly divided into 4 cohorts (n = 3 each). At 8 weeks, the spinal cords were harvested and serial cross-sections through the bridge were stained with antibodies against NF200 to identify regenerating axons (Figure 4A). The bridge area from a minimum of 3 sections per animal was counted by two blinded manual observers followed by an automated analysis of the same sections (Figure 4B-C). No statistical differences were observed in average axon densities between the manual and automated counts for each of the 4 groups (n=3 each, Figure 4D). Axon densities were in the range previously reported for this system (1000-2000 axons/mm2) (Thomas et al., 2013; Thomas et al., 2014). Variability within and between the simulated experimental groups, despite all animals receiving the same treatment condition (unmodified PLG bridges), likely arises from subject variability in response to SCI, which may result from minor variations in the surgical procedure. The average of the manual axon counts for each animal statistically correlated with the average of the axon counts obtained from the automated analysis (slope=0.998, Pearson’s r=0.853, p<0.001, Figure 4E).
Automated quantification for myelination of the regenerating axons was investigated within the bridges. Spinal cord cross-sections were stained with antibodies against NF200 and MBP to determine co-localization of axons and myelin. Due to the structural similarity of myelin staining to NF200+ axons (Figure 5A-C), a similar Hessian-based filtering technique was applied. Filtering of MBP+ myelin resulted in reduced background and well defined curvilinear structures that co-localized with filtered NF200+ axons (Figure 5D-F). The number of NF200+ axons that co-localized with MBP+ myelin following intensity threshold was determined and normalized to the total number of NF200+ axons. Manual counts for myelination were obtained in the same region of each section, and were consistent with previous reports (Thomas et al., 2013; Thomas et al., 2014). Manual counts correlated well with counts from the automated process (slope = 0.99, Pearson’s r = 0.823, p < 0.001, Figure 5G). No statistical differences were observed between the averages of manual and automated myelinated axon counts for each of the 4 simulated experimental groups (Figure 5H).
Regenerating the injured spinal cord has been a century long endeavor with various strategies including peripheral nerve grafts, cell suspensions, and biomaterial scaffolds (Geller and Fawcett, 2002; Wilson, 1984). A key aspect in characterizing regeneration is the number of regenerating axons. Counting regenerating axons manually can be time intensive as hundreds of axons can be present even with modest regeneration (Auer, 1994; Grider et al., 2006; Stokols et al., 2006; Tuinstra et al., 2012; Tuinstra et al., 2014). Manual counting typically involves multiple blinded observers that introduces a source of variability due to individual interpretation of fluorescence staining despite specific criteria for counting axons. This inherent subjectivity between multiple individuals introduces an observer-dependent bias. Automated methods for quantifying histology can substantially decrease processing times and increase objectivity. Automated or semi-automated quantification of axon counting or myelination is typically performed in high magnification electron microscopy images (Auer, 1994; Isaacs et al., 2014; Kim et al., 2015; McGavern et al., 1999). Programs for neuron tracing have been developed to study morphology and neurite outgrowth, but have not been fully adapted to provide accurate axon counts in regenerating tissues (Gensel et al., 2010; Koppes et al., 2014; Meijering, 2010; Meijering et al., 2004).
In this study, a semi-automated technique to estimate axon counts was developed based on Hessian-based image filtering to isolate axon structures prior to threshold detection. Convolution filtering, like deconvolution microscopy, can improve the resolution of common widefield fluorescent images in substantially less time than is required for higher resolution imaging techniques such as confocal microscopy. Hessian-based filtering has been employed to specifically enhance the contrast of curvilinear features, including axons in brain and spinal cord tissues, but has yet to be utilized in automated methods to determine total axon number (Grider et al., 2006; Sathyanesan et al., 2012). Several programs were tested in the course of this study to improve axon detection for computer-assisted counting. Hessian-based filtering in Image J, available through the FeatureJ plugin, resulted in images with enhanced contrast. The basal level of background fluorescence following filtering, however, was greater in Image J compared to the Hessian-based filtering technique performed in MATLAB, even with previously published parameter sets. This increased background with ImageJ resulted in decreased ability to detect axons using an intensity threshold. The find edges tool in Photoshop, along with the edge function in MATLAB, resulted in outlines of axons rather than enhancing the contrast of the axon structures. One common difficulty with the programs that were tested is the limit in adjustable parameters available to control the image restructuring.
The method published by Frangi and colleagues provides several adjustable variables including a scaling factor (σ) and kernel size for determining the Hessian matrix as well as the parameters β and c for adjusting how the eigenvalues of the Hessian matrix impact the subsequent image filtering (Frangi et al., 1998). An excellent example of this method in MATLAB is available from Dirk-Jan Kroon through the Mathworks website (www.mathworks.com). In this study, the scaling factor and kernel size were determined empirically based on qualitative assessment of the restructured images while the value for β was the same as previously published. The values for c and intensity threshold were obtained by minimizing false-positive and false-negative events in a small subset of images. While the parameters β and c were consistent with previously published studies, this condition may not persist for all cases. Changing the filtering parameters (σ, β, and c) drastically altered the brightness and morphology of filtered axons. Choosing these parameters was the greatest potential source of error we observed in this study. As with any parameter dependent on image intensity, these values will need to be determined for each study and will require batch staining, imaging, and processing of tissues within the study.
Hessian-based filtering of spinal cord cross-sections stained with antibodies against the axonal marker neurofilament (NF200) resulted in improved contrast and detection of axons regenerating through bridges. The autofluorescence and non-specific background staining from the PLG bridge was effectively removed. Using the Hessian-based filtering technique, biomaterial-associated fluorescent noise is minimized allowing for identification and detection of intended curvilinear features. Axons with a variety of orientations were detected following filtering and intensity threshold, providing a significant advantage over other techniques, such as line profile analysis, that predominantly measure features parallel to the plane of the section. While a few overlapping axons were observed in the semi-automated approach and likely counted as a single event, we often count these structures as one axon during manual counting due to the difficulty in differentiating between one long continuous axon, branching points of a single axon, or two overlapping axons. Hessian-based filtering thus allowed for semi-automated counting by intensity threshold and object detection in MATLAB of axon numbers that had good agreement with manual counts from unfiltered and filtered images. Importantly, the semi-automatic process was approximately 2 orders of magnitude faster than manual counting.
The semi-automated process was able to accurately capture axon numbers across a range of densities, which can be useful for characterizing regeneration across the tissue section that may have inhomogeneous distributed axons. Bridges, in particular, can have an inhomogeneous distribution due to providing a controlled architecture and defined surface to guide regenerating axons (Sakiyama-Elbert et al., 2012). Axons may be concentrated in specific regions (e.g., channels) of the bridge, with sparse axons elsewhere. This inhomogeneous distribution of axons would require significant time for manual counting. Stereology is often used as a means to provide unbiased and quantitative data about tissue sections (Kaplan et al., 2010); however, this method utilizes random, systematic sampling that may not be readily applicable to tissues with an inhomogeneous axon distribution. The speed by which the semi-automatic process demonstrated in this study can quantify axon numbers enables analysis of the entire region. The Bland-Altman plot representing the difference between manual and automated counts was employed to identify any potential bias in the semi-automated technique. Positive and negative normalized residual values of comparative amplitude were found across a wide range of total axon numbers, demonstrating no bias and similar performance at various densities.
In summary, we demonstrate that Hessian-based filtering of fluorescence images can facilitate counting of regenerating axons in cross-sectional images following SCI. Axon counts obtained from the semi-automated analysis statistically had good agreement with manual counts within a designated bridge area in spinal cord cross-sections. The semi-automated analysis also recapitulated the results of a multi-variable experimental study simulated by grouping subjects into 4 distinct cohorts. Co-localization of multiple axonal and axon-related structures can also be quantified to provide additional data concerning regenerating axons such as myelination. We used antibodies against neurofilament to identify axon fibers as it is a commonly used marker in SCI studies, however, this technique could be applied to other axonal markers such as Tau (Hurtado et al., 2006; Novikov et al., 2002; Stokols et al., 2006; Tsai et al., 2006). Importantly, this technique drastically reduced the time necessary to analyze large sets of data in a manner. Based on these results, Hessian-based filtering followed by object detection provides a time-saving tool for semi-automated analysis of axon regeneration and myelination following SCI. Additional applications may benefit from the methods developed here, therefore the code used in this study will be freely provided to any interested researcher. This work provides a significant step towards a more rapid characterization of axon growth in the injured spinal cords, and through engineered biomaterial environments.
The authors would like to acknowledge Dr. Bea Penalver for her expertise in data analysis and Dirk-Jan Kroon for online examples of Hessian-based filtering methods. The published work was funded by the National Institutes of Health (R01 EB005678).
Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.