Morphological characteristics of muscle fibers, such as fiber size, are critical factors that determine the health and function of the muscle. However, at this time, quantification of muscle fiber cross-sectional area is still a manual or, at best, a semiautomated process. This process is labor intensive, time consuming, and prone to errors, leading to high interobserver variability. We have developed and validated an automatic image segmentation algorithm and compared it directly with commercially available semiautomatic software currently considered state of the art. The proposed automatic segmentation algorithm was evaluated against a semiautomatic method with manual annotation using 35 randomly selected cross-sectional muscle histochemical images. The proposed algorithm begins with ridge detection to enhance the muscle fiber boundaries, followed by robust seed detection based on concave area identification to find initial seeds for muscle fibers. The final muscle fiber boundaries are automatically delineated using a gradient vector flow deformable model. Our automatic approach is accurate and represents a significant advancement in efficiency; quantification of fiber area in muscle cross sections was reduced from 25–40 min/image to 15 s/image, while accommodating common quantification obstacles including morphological variation (e.g., heterogeneity in fiber size and fibrosis) and technical artifacts (e.g., processing defects and poor staining quality). Automatic quantification of muscle fiber cross-sectional area using the proposed method is a powerful tool that will increase sensitivity, objectivity, and efficiency in measuring muscle adaptation.
- automated software
- cross-sectional area
- image analysis
skeletal muscle is an extremely adaptive tissue that is able to change size depending upon external stimuli. Increases in muscle mass in response to resistance training and losses in mass caused by disuse or associated with chronic diseases such as cancer and HIV are primarily due to hypertrophy or atrophy of individual muscle fibers, respectively, rather than addition or loss of fibers (6). Thus the ability to accurately and efficiently quantify muscle fiber size is vital for assessing muscle function, because muscle mass is the primary determinant of muscle strength (2, 9).
Currently, muscle fiber size is most commonly assessed by quantifying fiber cross-sectional area (CSA) in fixed or frozen samples. To delineate individual muscle fibers, immunohistochemical approaches targeting different antigens in the sarcolemma or basal lamina surrounding the fibers are commonly used. For example, dystrophin immunoreactivity clearly defines the fiber sarcolemma (7); caveolin-3 staining is more patchy, but also localizes to the sarcolemma (4); laminin staining clearly defines the basal lamina surrounding each fiber (3); and wheat germ agglutinin staining reacts with proteoglycans in the basal lamina and extracellular matrix (5, 11). These protocols are favored over more basic histochemical approaches because they offer precisely distinguishable fiber units. However, other simple staining procedures, such as hematoxylin and eosin staining also delineate fiber perimeter and morphology. Even with these techniques that enable visualization of fiber boundaries, quantification of fiber CSA in microscopic images is still most often accomplished by manually tracing the delineated fibers followed by calculation of fiber CSA with software programs such as ImageJ (8, 10). Even with software programs that can semiautomatically aid in identification of fiber units, such as Zeiss AxioVision, manual correction of errors is required that is both time consuming and labor intensive. The purpose of this study was to develop and validate a high throughput, automatic image segmentation algorithm to quantify fiber CSA to increase objectivity, accuracy, and efficiency in muscle image analysis.
MATERIALS AND METHODS
Collecting and processing muscles for immunohistochemistry.
All animal procedures were conducted in accordance with institutional guidelines for the care and use of laboratory animals as approved by the Animal Care and Use Committee at the University of Kentucky. Adult (4 mo of age) male C57B1/6 mice were injected intraperitoneally with ketamine (100 mg/kg) and xylazine (10 mg/kg) and euthanized by cervical dislocation. Plantaris muscles were excised, pinned to a cork block at resting length, covered with a thin layer of Tissue Tek OCT compound (Sakura Finetek), and then quickly frozen in liquid nitrogen-cooled isopentane and stored at −80°C until processing.
To alter muscle morphology, mice were subjected to synergist ablation surgery in which removal of the gastrocnemius and soleus muscles of the lower limb mechanically overloads the remaining plantaris muscle (7). After 2 wk of overload, plantaris muscles were harvested as described above.
Antibodies and reagents used were anti-dystrophin (cat. no. VP D505); mouse IgG blocking reagent (cat. no. MKB-2213); Vectashield mounting medium (cat. no. H-1000), all from Vector Laboratories (Burlingame, CA); Texas Red-conjugated goat anti-mouse (cat. no. 610-109-121, Rockland Immunochemicals, Gilbertsville, PA); and DAPI (cat. no. D3571, Invitrogen, Carlsbad, CA).
Frozen muscles were sectioned (7 μm), air dried, and stored at −20°C. Two different methods were used to delineate fibers. Fresh frozen sections were immunoreacted with the dystrophin antibody, followed by Texas Red-conjugated secondary antibody. Sections were postfixed in 4% paraformaldehyde and then stained with DAPI to visualize nuclei. Alternatively, frozen sections were fixed in 4% paraformaldehyde followed by incubation with Texas Red conjugated wheat germ agglutinin (WGA) (cat. no. W21405, Invitrogen). Three to five washes with phosphate-buffered saline were performed between each step. All images of cross sections were captured with a Zeiss AxioImager M1 upright fluorescent microscope at ×10 or ×20 magnifications.
Analysis of fiber CSA using manually corrected AxioVision software.
CSA was quantified on muscle cross sections following dystrophin immunohistochemistry (Fig. 1) with AxioVision Rel software (v4.8, Zeiss) using the sequence described below, with default settings for brightness, gamma, contrast, sigma filter, and edge enhancement. First, individual fibers were identified based on the intensity and continuity of the dystrophin-stained sarcolemma surrounding each fiber using the segmentation function in the program. A single representative fiber was selected manually and the program automatically defined and filled all remaining fibers in the field of view (Fig. 2A; high magnification image is shown in Fig. 2B). Second, additional parameters within the program were selected manually, in the following order: automatic object separation, which identified separate fibers; deletion of artifacts, which identified connective tissue and other interfibrillar spaces for exclusion; and filling of holes, which filled spaces within fibers left over from segmentation. The resulting image is shown in Fig. 2C, and a high magnification is shown in 2D. Third, interactive processing of the measurement mask was performed to manually correct errors in segmenting the image, such as delineating missed fibers, separating fibers identified as single units, and smoothing wavy or jagged edges of fibers (Fig. 2E). This step is the most labor- and time-intensive procedure depending on the number of errors. Fourth, the set measurement properties function was used to complete the analysis. The measure features function calculated area and number of regions, and the calculated CSA values were converted to a text file format. The draw features functions area and ID region display CSA values for each fiber in the image. Finally, the frame function deleted incomplete fibers around the edge of the image to generate final semiautomatic segmentation results (Fig. 3A; high magnification is shown in Fig. 3B). It is obvious that the whole procedure requires a tremendous amount of human intervention. Even for a very experienced researcher, the whole procedure requires 25–40 min to process a single image containing approximately 100 muscle fibers.
Automated CSA quantification.
Observing the labor- and time-intensive procedure described above, we proposed a completely automatic algorithm (Fig. 4) for muscle fiber segmentation. The automatic image segmentation algorithm and CSA calculation contains three steps: 1) ridge detection to extract the fiber boundaries; 2) mathematical morphology operations to postprocess the detected seeds; and 3) a deformable model to drive the contour to converge to the muscle fiber boundaries.
Muscle fiber boundaries, visualized by dystrophin staining (Fig. 5A), were modeled as intensity ridges. Because both the fibers and edges vary in size, multiscale ridge detection was used to detect these edges. Given image I(x,y), we first computed Hessian matrix of the image, where Iσ(x,y)= I(x,y)*G(x,y,σ), * is used to represent convolution, and G(x,y,σ) denotes the Gaussian kernel with scale σ. We then applied an eigenvalue decomposition operation to the Hessian matrix, with λ1, λ2 representing two eigenvalues of H and choose |λ1| > |λ2|. Then we computed the likelihood measurement for the ridges: where RB = |λ2/λ1|, S = λ12 + λ22. As pointed out in Ref. 2a, the negative sign of λ1 indicates bright ridge structure in a dark environment and vice versa. S is used to suppress those spurious responses with small λ1 and λ2. Considering different scales, the final likelihood of a ridge was determined as: where i represents the i-th scale. A threshold μ was computed with an automatic threshold algorithm (9a) from the likelihood map generated. The pixels with a likelihood larger than the threshold are treated as ridges between different muscle fibers (Fig. 5B).
A morphological closing operation was used to connect nearby ridges. Muscle fiber seeds were generated by using the inverse of the ridge detection map (Fig. 5C). To separate touching fibers, an iterative erosion method was used. First, fiber seeds calculated less than a given threshold (for example, 5,000 pixels) were filtered out; second, the remaining seeds were eroded with a disk-shaped morphological structuring element. These two steps were repeated until convergence was reached. In this way, seeds smaller than a given threshold were kept intact, and seeds larger than the threshold were appropriately separated.
During the postprocessing step, we defined roundness as follows:
Those seeds with small roundness were removed in post processing.
Deformable model to calculate the final boundaries of the fibers.
The gradient vector field (GVF) deformable model algorithm (12) was used to delineate the final muscle fiber boundaries, as illustrated in Fig. 6. A traditional deformable model (snake) is defined as X(s) = [x(s), y(s)], where x(s) and y(s) are x and y coordinates along the contour, and s represents in arc length with value from 0 to 1 to minimize an energy function
The first term Eint[X(s)] represents the internal energy, defined as where Xs(s) is the first derivative of X(s), and Xss(s) is the second derivative of X(s) with respect to s. The external energy term is defined in different format for different variants of deformable model; the specific definition for GVF snake will be mentioned later.
The goal of the computation is to find the local minima of Esnake. On the basis of the Euler-Lagrange equation, minimized when where Xss(s) and Xsss(s) represent the second and fourth derivatives of curve with respect to parameter s, and ∇Eext[X(s)] denotes the gradient of external energy term
To effectively lead the snakes into the concave region, Xu (8) proposed a gradient vector flow based external energy term where Θ represents the GVF field defined as Θ(x,y) = [u(x,y), v(x,y)]. The ∇[Gσ(x,y) * f(x,y)] represents the gradient of the input image after Gaussian smoothing with variance σ and zero mean.
The snake is made dynamic by defining X as the function t and s
The solution is achieved by solving the discrete equations iteratively.
The final segmentation results are shown in Fig. 5D. The number of iterations was chosen as 40 throughout the analyses. Given the segmentation results from GVF snake, we calculated the area of each fiber using region filling. The actual physical size of each individual muscle fiber (μm2) was calculated from the automatic segmented area (in pixels) multiplied by the physical resolution of the microscope (μm/pixel), displayed on the final quantified image (Fig. 3C).
Average muscle fiber CSA in 30 randomly selected images from normal plantaris muscle was determined using both the proposed automatic quantification algorithm and the semiautomatic method using AxioVision software with careful manual annotations, as illustrated for one image in Fig. 3, A and C. A direct comparison of CSA values obtained from this testing dataset demonstrated consistency in CSA values between the automatic algorithm and the AxioVision analysis with manual corrections (Table 1). The percent difference ranged from −4.0 to +2.3%, with an average difference of ±1.3%. The percent coefficient of variations ranged from 0.04 to 2.88%, with an average number <1% (0.93%). The automated algorithm is equally effective in quantifying fiber CSA in sections treated with wheat germ agglutinin (Fig. 7, A–C), that binds proteoglycans in the extracellular matrix surrounding fibers.
Fiber CSA determination is often required following an intervention that may alter the highly organized structure of muscle. The automated algorithm is effective in images of muscle undergoing adaptation in response to synergistic ablation that disrupts normal muscle morphology (Fig. 8). As shown in Fig. 8A, the overloaded plantaris muscle is less organized, demonstrating cellular infiltration, increased noncontractile area, and fiber size heterogeneity due to fiber hypertrophy, fiber splitting, and damage-induced regeneration. The results from the automated approach (Fig. 8C) correlated well with the semiautomated process (Fig. 8B). The latter required more than 40 min of manual correction, whereas the proposed automatic segmentation algorithm required ∼15 s. A direct comparison of the two methods in five randomly selected images demonstrated a percent difference range from −1.1 to +3.8%, with an average difference of ±1.6%. The percent coefficient of variations ranged from 0.4 to 2.7%, with an average of 1.3%.
Finally, sections that demonstrate poor staining quality, illustrated in Fig. 9A, often cannot be quantified, so that the experiment must be repeated, wasting valuable tissue. Poor staining is one of the most common technical problems encountered, requiring extensive manual correction for the software to identify fibers with weak boundaries. The semiautomatic method, followed by extensive manual annotation, required 40 min to obtain the result shown in Fig. 9B. However, it required only 15 s using the automatic segmentation algorithm to process this challenging image (Fig. 9C). The final CSA values differed by only 1.24% (1,611 compared with 1,631 μm2).
Skeletal muscle is uniquely sensitive to mechanical load, with individual muscle fibers able to increase and decrease in size. Due to this robust adaptive property, muscle biologists have long understood the importance of accurately quantifying muscle fiber size (1). Current state-of-the-art image analysis methods (manual and semiautomated) are extremely inefficient and labor intensive, especially for challenging images containing biological inconsistencies and/or technical errors in processing, including differential or incomplete staining, poorly defined fiber boundaries, freeze artifacts, and sample defects. These lead to errors that must be manually corrected. The novel automatic image analysis algorithm proposed here attempts to mitigate such inconsistencies through the use of advanced machine learning and biomedical image analysis technology innovations. Results presented here show that the proposed method provides excellent automatic segmentation results, even in images that show heterogeneity in fiber CSA induced by synergist ablation surgery (Fig. 8). This suggests that fiber CSA can also be accurately and efficiently quantified with this new method in other conditions such as denervation, spinal cord injury, or genetic manipulation (i.e., mdx or myostatin knockout mice).
We acknowledge that other analytical software programs, such as Cyteseer, are currently on the market. However, compared with Cyteseer, our system is designed for processing digitized skeletal muscle images specifically. We have downloaded the Cyteseer software and compared the Cyteseer segmentation results with the automated method described in this paper. In our test, Cyteseer failed to segment many testing images, especially those images exhibiting poor quality with weak boundaries and high background, leading us to conclude that in our hands, Cyteseer does not perform as well as our algorithm (data not shown).
The novel segmentation and quantification method described here is superior based on a data-driven algorithm, providing a general solution to automate muscle fiber image analysis that eliminates subjective variability while dramatically improving efficiency. The algorithm was extensively tested using dystrophin immunohistochemistry and WGA staining. It is currently being adapted to other staining methods, such and hematoxylin and eosin, which would be more widely applicable, particularly in diagnosis of muscle pathology. Moreover, algorithms are currently being developed that will enable quantification of additional morphological properties of muscle, such as fiber shape, fiber type, the number and position of myonuclei and satellite cells, capillary density, cellular infiltration, and fibrosis. Together, these measures will enable a more comprehensive, sensitive, and quantitative evaluation of muscle adaptation.
In summary, we have developed an automatic image segmentation approach that can accurately and efficiently quantify muscle fiber size in cross sections. It is very common for a biology lab to analyze hundreds of histochemical images and stains to appropriately represent tissue-specific adaptations in a single research project. We argue that the greatest advantage of this software is its potential for an extremely high-throughput analysis that will transform days of labor-intensive work into minutes, saving valuable research time. It is even more exciting that a similar approach can be extended to other tissues, such as adipose, where cell size is of interest. Currently, we are working on developing a user friendly interface that will make this specific software widely available to the muscle research community for a broader application in evaluating muscle health and function.
This project was supported by National Institutes of Health (NIH) Grants AG34453 and AR60701 to C.A.P. The project described was also partially supported by the National Center for Research Resources and the National Center for Advancing Translational Sciences, NIH, through Grant UL1TR000117. The content is solely the responsibility of the authors and does not necessarily represent the official views of the NIH.
No conflicts of interest, financial or otherwise, are declared by the authors.
Author contributions: J.M., J.D.L., F.L., L.Y., and C.A.P. conception and design of research; J.M., J.D.L., and F.L. performed experiments; J.M., J.D.L., F.L., L.Y., and C.A.P. analyzed data; J.M., J.D.L., F.L., L.Y., and C.A.P. interpreted results of experiments; J.M., J.D.L., F.L., and L.Y. prepared figures; J.M., J.D.L., F.L., L.Y., and C.A.P. drafted manuscript; J.M., J.D.L., L.Y., and C.A.P. edited and revised manuscript; J.M., J.D.L., F.L., L.Y., and C.A.P. approved final version of manuscript.
The authors acknowledge Cynthia Byars for expertise with the figures.
- Copyright © 2013 the American Physiological Society