1.

Introduction

Over the past five years there has been growing interest in biomedical imaging based on the nonlinear optical properties of target tissues and molecules. Such methods include second^{1, 2} and third harmonic generation,³ sum frequency generation,^{4, 5, 6} and coherent anti-Stokes Raman scattering.^{7, 8} Second harmonic generation (SHG) imaging has several features that make it attractive for biological applications: it is nondestructive, noninvasive, and rapid enough for real-time three-dimensional (3D) scanning. Because the signal arises from the intrinsic optical properties of the component macromolecules rather than from dyes or contrast agents, the signal-to-noise ratio is excellent. SHG imaging appears particularly suited for analyzing the structural organization of biological macromolecules. Collagen,^{9, 10} myosin,¹¹ and tubulin¹²—three of the most widely distributed structural biomolecules in higher life forms—are all noncentrosymmetric molecules capable of SHG. Furthermore more, these molecules are organized into complex assemblies whose chiral nature at the supramolecular level most likely contributes to the magnitude of the SHG signal.

Collagen, the most abundant protein in most organisms, produces an exceptionally strong SHG signal; its hyperpolarizability is reported to be within an order of magnitude of that of crystalline quartz.¹³ Although SHG was first observed in collagenous tissue almost $50\text{years}$ ago, significant work on the relationship between collagen structure and SHG was not reported until the seminal studies of Roth and Freund.^{14, 15} Since then, we and others have reported on the use of SHG scanning to obtain structural information about collagenous tissues.^{16, 17} There has been particular interest in determining how SHG can be used to analyze fundamental properties of microfibril structure and organization.

In the present study, we report on a novel signal analysis method in which we use polarization-modulated SHG (PM-SHG) to quantify the level of structural disorder in collagenous tissues. We identify two independent parameters related to the structural organization derived from SHG scan data and define a “disorder index” as a function of these parameters. We validate this method through detailed correlative studies of heat-induced disorder in a rat tail tendon fascicle, a model we have previously described.¹⁸ We measure disorder in this model in two ways: (1) by applying our algorithm to SHG scan data and (2) by morphometric analysis. Correlation between these methods was excellent. We show that our method can be applied to the analysis of pathological changes in structurally complex tissues by assessing disorder in injured intervertebral disks. We use an in vivo model of disk loading developed by one of us (J.L.); this model has been extensively characterized with respect to biomechanical, histochemical, and immunochemical properties.^{7, 8} In the present study, we performed SHG scanning on tissue sections prepared from the original tissue blocks; we compared SHG-derived values for disorder with values obtained by morphometric analysis.⁹ Our data show that PM-SHG analysis represents a robust analytical method that is not affected by potentially confounding conditions, and that measurement of disorder represents a useful approach for discriminating between normal and abnormal tissue. To our knowledge, this is the first report of a method in which the degree of structural disorder, determined by SHG, is used to define both the extent and severity of tissue pathology. In this context, the physical meaning of disorder as a phase state subtends the more general meaning of disorder as a pathological condition.

2.

Experimental Setup

2.1.

Analysis of SHG Scan Data

We were interested in developing an analytical method sensitive to regional variations in levels of disorder. We therefore selected a “nearest neighbor” strategy for analyzing SHG scan data. Data acquired during a scan are stored in two matrices representing signal intensity and polarization angle of maximum signal, respectively. The intensity matrix is used to define a floor value applied to both data sets; the angle data are then used to assess structural order in each $3\times 3\text{pixel}$ “neighborhood.” Each pixel, excluding those in the first and last rows and columns, is analyzed in the context of its own nine-pixel neighborhood, in which it occupies the central position. We do not need to know specific physical correlates—that is, the number of molecules in the area represented by a pixel or the hierarchical level of organization—to analyze the structural relationships inferred from the orientation data. This issue has been explicitly addressed by Roth and Freund.⁹ Order in a type I collagen lattice presupposes parallel alignment of its components: this state is shown schematically in Fig. 1a . Although an ordered state is energetically favorable, its stability depends on the particular confluence of mechanical influences, stochastic influences, and genomic programs occurring locally.^{19, 20} Under disordering conditions, the relationship between the index pixel and its neighbors may change. Mild disorder is shown in Fig. 1b: only two neighbors are not parallel to the index pixel, and their deviation from parallel alignment is relatively small. Severe disorder is shown in Fig. 1c: no neighbor is parallel to the index pixel, and the magnitude of their deviation is large. In some cases, lack of alignment with the index pixel does not represent disorder. In Fig. 1d, several nonparallel neighbors represent the boundary of a region whose components have a different orientation. This condition is identified by an edge detection algorithm, described below.

Fig. 1

Patterns of disorder in nine-element neighborhoods. A neighborhood consists of a $3\times 3$ array of matrix elements; the central element is designated as the index element. Panels (a) to (d) represent patterns of neighborhood organization when the matrix elements consist of SHG angle data. (a) Maximally ordered neighborhood in which all elements are parallel to the index element. (b) Neighborhood with low level of disorder: only two elements are not parallel to the index element and their deviation from parallel orientation is minimal. (c) Neighborhood with high level of disorder: most neighbors are not parallel to the index element; the magnitude of their deviation is high. (d) Mixed neighborhood: three of the elements (elements 1, 4, and 6) that are not parallel to the index element are part of a boundary and therefore do not contribute to disorder. Only element (2) is not parallel to any other element. There are five valid neighbors.

Analysis of data from earlier studies suggested that in most tissues disorder can be characterized as a function of two parameters: (1) the magnitude of the index pixel’s mean deviation from parallel alignment with its neighbors, and (2) the number of nonparallel neighbors. In the present study, we refer to the first parameter as the angle index (AI) and the second parameter as the neighbor index (NI). The magnitude of AI for a given index pixel is expressed as the sum of the differences in the orientation of the index pixel and that of each of its neighbors normalized to $n$ , the number of neighbors used in the computation (i.e., excluding neighbors that represent the boundary of another region). The definition of AI is shown in Eq. 1

Eq. 1

$$AI=\left(\frac{1}{n}\right)\sum _{i=1}^n\left(\mid \alpha -{\theta }_i\mid \right),$$

where $\alpha$ is the angle of the index pixel expressed in radians and $\theta$ is the angle of the $i$ ’the nearest neighbor pixel expressed in radians. A neighbor is defined as nonparallel if the difference between its orientation angle and that of the index pixel is greater than an empirically determined threshold value. We determine the threshold value for each type of tissue to be analyzed from SHG scans of samples defined as “controls:” we determine the value for the mean difference in the orientation angles of all pairs of adjacent pixels in each row and each column and set the threshold value as two standard deviations of the mean difference value.

The neighbor index (NI) is defined as the total number of nonparallel neighbors, excluding neighbors that are part of a boundary. Boundaries need to be identified, because many tissues contain well-defined components whose orientation angles differ from each other; for example, lamellar rings in intervertebral disks contain fibrils whose orientation differs from that of fibrils in neighboring rings by $120\mathrm{deg}$ . Our analysis program therefore includes the ability to recognize both straight and curved boundaries. Edges are determined by a regional variant of Canny edge detection. Two conditions must be met for a nearest neighbor to be defined as part of an edge: (1) the angle difference between the nearest neighbor and each of two adjacent pixels on either the horizontal or vertical axis is less than the threshold value; (2) the angle difference between the index pixel and the nearest neighbor is greater than the threshold. Curved boundaries are detected by examining the orientation of pixels on a diagonal axis in relation to the index pixel. Figure 1d illustrates a neighborhood in which the index pixel lies next to a boundary; the neighborhood contains only one nonexcluded neighbor (pixel 2) that is not parallel.

Although we define the degree of “order” for a neighborhood as a function of two independent parameters (AI and NI), it is useful to express disorder as a single value for convenience in graphic representation. To define such a term, we first need to determine the relationship between NI and AI and whether or not it remains consistent in different tissues as well as under different conditions of disorder. Analyses of random matrices showed that the ratio of NI to AI approaches the value of 3.50 with increasingly smaller standard deviations as the matrices increase in size (Table 1 ). Similar values were found in tendon and disk samples. We also found that the NI/AI ratio remained consistent in rat tail tendon samples across a wide range of disordered states (see Sec. 3). Assuming that the angle index and the neighbor index contribute equally to the overall state of disorder, we define the disorder index, referred to as ANI, as the sum of NI and 3.5 times the value of AI

Eq. 2

$$\mathit{ANI}=\mathit{NI}+3.5\left(\mathit{AI}\right).$$

Table 1

Effects of matrix size on variance of mean values of disorder indices in random number models.1

Size (pixels)	Index	Mean±SD
$10\times 10$	$\overline{ANI}$	$6.003\pm 0.148$
	$\overline{AI}$	$0.955\pm 0.161$
	$\overline{NI}$	$2.797\pm 0.054$
$25\times 25$	$\overline{ANI}$	$7.697\pm 0.068$
	$\overline{AI}$	$0.931\pm 0.055$
	$\overline{NI}$	$3.331\pm 0.028$
$50\times 50$	$\overline{ANI}$	$7.842\pm 0.043$
	$\overline{AI}$	$0.991\pm 0.048$
	$\overline{NI}$	$3.390\pm 0.006$
$100\times 100$	$\overline{ANI}$	$7.874\pm 0.007$
	$\overline{AI}$	$1.041\pm 0.048$
	$\overline{NI}$	$3.48\pm 0.005$
$500\times 500$	$\overline{ANI}$	$7.877\pm 0.002$
	$\overline{AI}$	$1.046\pm 0.008$
	$\overline{NI}$	$3.50\pm 0.005$

¹Random number arrays were generated in MATLAB as described in the text. Array size ranged from 10×10 elements to 500×500 elements. The indices of disorder were determined in five matrices of each size; mean values and standard deviation (SD) are shown.

A disorder map, obtained by plotting the value of ANI for each pixel (excluding the outer scan pixels), provides information about regional distribution of disorder. If we want to compare samples with respect to disorder, we can easily compute the mean disorder for the whole scan (the sum of the disorder values for each pixel divided by the number of pixels analyzed). However, if the pattern of disorder in a tissue is heterogeneous or the tissue is structurally complex, we may want to determine mean disorder values for smaller regions of the scan in order to assess intrasample as well as intersample variability.

When we compare groups of scans with respect to their level of disorder, we are primarily interested in determining whether or not differences between the groups are statistically significant. Do the values for the indices of disorder have any intrinsic meaning? To impute physical meaning to the values for disorder, we need to express them in relation to a common frame of reference. The disorder index for a maximally disordered system is a logical choice. We cannot determine this value in a biological sample, because no SHG signal would be detectable.²¹ We can, however, construct a mathematical model that simulates the signal of a completely disordered collagenous array. To do this, we generate a random number matrix whose elements extend over $\pi$ radians (the range of angles in the SHG data set) and compute its disorder indices. Values for disorder indices obtained in biological samples are normalized to the maximum values obtained in the random number matrix. If we designate the term ${\overline{\mathit{ANI}}}_{\mathit{max}}$ as the mean value for the disorder index in a maximally disordered system, the computational procedure for determining the mean disorder index, $\overline{\mathit{ANI}}$ , for a sample whose scan region is $j\times k$ pixels in size is summarized in Eq. 3

Eq. 3

$$\overline{\mathit{ANI}}=\frac{1}{{\overline{\mathit{ANI}}}_{\mathit{max}}}\frac{1}{jk}\left\{\sum _{p=1}^j\sum _{q=1}^k{\mathit{ANI}}_{p,q}\right\}.$$

These equations are implemented in a program originally written in C for this purpose and subsequently incorporated into a MATLAB m-file.

Since completing the experiments reported in this study, we have considered an alternative method of determining ANI that eliminates the potential problems inherent in the use of a “scaling factor,” which cannot be assumed to be valid for all tissues or under all conditions. We propose that AI and NI be expressed relative to the values of these parameters in a maximally disordered pixel. We can determine these values by applying our analytical methods to large random number matrices and determining maximum values for AI and NI in the respective data sets. The value for ANI is defined as the average of AI and NI. The mean value for ANI over a given scan region (Eq. 3) is also simplified as it no longer requires normalization to ${\mathit{ANI}}_{\mathit{max}}$ .

3.

Results

3.1.

Mathematical Models

We generated random number matrices to test the behavior of our computational method when used to assess samples in a state of maximum disorder. The SHG signal becomes undetectable before maximum disorder is reached in biological samples. We analyzed random number matrices of different dimensions using our disorder algorithm; in these matrices, the mean values for the three parameters of disorder represent the values characteristic of a maximally disordered system. We therefore append the subscript max to their designation. Values for ${\overline{AI}}_{\mathit{max}}$ , ${\overline{NI}}_{\mathit{max}}$ , and ${\overline{ANI}}_{\mathit{max}}$ are shown in Table 1. As described in Sec. 2.1, we normalized mean values for disorder obtained in biological samples to these values. Thus, the value for $\overline{ANI}$ in a biological sample tells us what percent of maximum possible disorder is present in the scanned region of the sample.

Our algorithm is designed to recognize edges: the boundaries that occur in normal tissues between regions whose components differ in orientation. For example, intervertebral disks comprise 10 to 12 lamellar rings: the orientation of fibrils within any given ring are offset by $120\mathrm{deg}$ from the orientation of fibrils in the two adjacent rings. We have incorporated a simple edge detection algorithm that recognizes organized orientation shifts that about the index pixel. We tested our algorithm in models containing straight boundaries as well as curved boundaries simulating the geometry of lamellar rings and confirmed its ability to recognize these features (data not shown).

Models were constructed to simulate different patterns of disorder. Figure 2a shows the orientation map of a model of heterogeneous disorder: a large, ordered region contains a small, central disordered region. A homogeneous pattern is shown in Fig. 2c. Every pixel has one or more deviant neighbors; the magnitude of the angle deviation for all pixels is relatively low. Disorder maps for the two models are shown in Figs. 2b and 2d. The mean value for disorder $\left(\overline{ANI}\right)$ in the model of heterogeneous disorder is 14.6: this value represents the mean value for the pixels in both the ordered outer region and the disordered central region. The value for $\overline{ANI}$ in the homogeneous model of disorder, shown in Fig. 2d, is 16.1. Although the values for $\overline{ANI}$ are similar in the two models, values for $\overline{AI}$ and $\overline{NI}$ are not. The model of homogeneous disorder has a low value for $\overline{AI}$ coupled with a relatively high $\overline{NI}$ , but the model of heterogeneous disorder has the reverse. These differences are reflected in the values for the ratio of these two indices: the ratio is 1.02 in heterogeneous disorder and 0.21 in homogeneous disorder. Thus, the ratio of $\overline{AI}$ to $\overline{NI}$ may be useful in discriminating between different patterns of disorder.

Fig. 2

Comparison of two models of disorder. Two models representing different patterns of disorder were constructed in MATLAB. (a) Heterogeneous disorder was simulated by constructing a small, highly disordered region within a large, ordered region. The small disordered region in panel (a) was constructed with a random number generator set to include the range of possible angles present in a SHG data set. The graphic representation in panel (a) was generated with the MATLAB quiver function. (b) Disorder map of the model in panel (a), obtained by plotting the value for ANI for each pixel. The mean value for disorder $\left(\overline{ANI}\right)$ for the whole panel (ordered plus disordered region) is 14.6. (c) Orientation map of a model of homogeneous disorder, constructed with a random number generator encompassing only $20\text{-}\mathrm{deg}$ angle variation. All pixels are surrounded by nonparallel neighbors, but the magnitude of deviation is small. (d) Disorder map of same region. $\left(\overline{ANI}\right)$ is similar in the two models—however, the ratio of the independent parameters AI and NI is very different (1.02 versus 0.21).

3.2.

Heat-Induced Disorder in Fascicles

Rat tail tendon fascicles undergo a predictable sequence of changes during heat-induced structural disorganization, with initial changes detectable at the tips of the fascicles and progressing axially. The effects of heating on the SHG signal are shown in Fig. 3 : after $10\text{-}\min$ exposure to $58°\mathrm{C}$ , the outermost fibrils on the fascicle tip show a marked decrease in signal intensity, but have not yet separated into divergent bundles [Fig. 3a]; after $20\min$ , the region of decreased SHG signal has extended axially, and discrete fibril bundles can be seen fanning out [Fig. 3c]. Disorder maps are shown in Figs. 3b and 3d. After $10\text{-}\min$ heating, note how the regions of detectable disorder are confined to the outer rims of the fascicle tips and the outer edges. In contrast, disorder after $20\text{-}\min$ heating has extended across the entire tip.

Fig. 3

Heat-induced disorder in rat tail tendon fascicle. Progressive changes in SHG data in rat tail tendon fascicles exposed to $58°\mathrm{C}$ heat. (a) SHG intensity, $10\text{-}\min$ heating: note SHG intensity is decreased only at the very tip of fascicle even though fibril bundles are beginning to fan out. (b) SHG disorder map (ANI values), $10\text{-}\min$ heating: disorder is detectable only at the fascicle tips and sides. (c) Orientation map: note that disruptions in parallel orientation of fibrils are associated with regions of detectable disorder. (d) SHG intensity, $20\text{-}\min$ heating: note how SHG intensity is now decreased across the entire tip, not just at the edges, as fibril bundles begin to separate. (e) SHG disorder map, $20\text{-}\min$ heating: disorder is detectable across the tip, consistent with loss of SHG signal. (f) Orientation map: note that loss of parallel alignment occurs in some fibril bundles faster than in others and is associated with increased levels of disorder.

Photomicrographs were obtained of fasicles heated for 20 or $40\min$ (data not shown) that corresponded to the regions used for SHG scans. Scans were selected in which there were progressive changes in SHG signal and in disorder along the length of the scan. Morphometric analysis was performed on circular regions, $100\mu \mathrm{m}$ in diameter, at intervals of $150\mu \mathrm{m}$ , with the initial region representing the least damaged area. A total of four separate fascicles was used. The images were binarized and disorder was analyzed using the method of MIL determination. Data are expressed as the aspect ratio of the ellipse best fitted to the data. Typical results are shown in Fig. 4a ; changes in SHG signal intensity over the same region of the fascicle are shown in Fig. 4b. Note that intensity progressively decreases as the aspect ratio increases. Changes in the disorder index for the same scan are shown in Fig. 4(b-2). Disorder increases as signal decreases. To examine correlations between aspect ratio data and the SHG parameters, mean SHG values were obtained in the regions of the scan corresponding to the regions used for MIL detection. The mean values for the SHG indices of disorder and for the corresponding aspect ratios are summarized in Table 2 . Control values are also shown in the table. Correlations between the aspect ratio data and SHG indices of disorder were analyzed using regression analysis; results are shown in Fig. 4c. Significant correlations were found between the mean aspect ratio and each disorder index $(p<0.001)$ . Additional statistical analyses were done to analyze the relationships among the indices of disorder and SHG intensity; data are shown in Fig. 5 . Regression analysis showed that SHG intensity was inversely correlated with all of the indices of disorder [Figs. 4a, 4b, 4c]; regression ANOVA showed that all correlations were highly significant $(p<0.0001)$ . Statistical analysis of residuals was conducted for all correlations; no additional data transforms were indicated. Analysis of the correlation between the two independent indices of disorder, AI and NI, is shown in Fig. 4d. The relationship between the two parameters was linear over much of the disorder spectrum; however, the best-fit curve was a second-order polynomial when data for high levels of disorder were included.

Fig. 4

Correlations between morphometry and SHG indices of disorder. (a) Mean aspect ratio values of a fascicle subjected to heating for $20\min$ at $58°\mathrm{C}$ . Photomicrographs of regions of fascicle $100\mu \mathrm{m}$ in diameter were binarized and analyzed as described in the text. Four regions were analyzed: the value for distance refers to the distance from the top of the SHG scan. The top of the scan is furthest from the fascicle tip. Mean values ± SD for the aspect ratio are shown; each group was significantly different from the others by ANOVA $(p<0.001)$ . (b) Selected data from a SHG scan of the same region of fascicle subjected to morphometry. [(b)-1] SHG intensity values are highest at the start of the scan and decrease nearer the tip. [(b)-2] ANI values from the same scan: disorder increases as the scan approaches the tip. (c) These three panels show correlations between values for the aspect ratio and each of the disorder indices for the length of the scan. Correlations were examined using linear regression analysis. For each region of the fascicle correlations between the aspect ratio and the corresponding SHG index of disorder were highly significant $(p<0.001)$ .

Fig. 5

Correlations among indices of disorder in heated rat tail fascicle. Values for AI, NI, ANI, and SHG intensity were obtained from analyses of fascicles heated for varying amounts of time at temperatures from $56\text{to}58°\mathrm{C}$ . Correlations among the indices of disorder and between disorder and intensity were investigated using linear and nonlinear regression. Best-fit curves and regression coefficients for representative analyses are shown in panels (a) to (d). (a) AI versus SHG intensity. (b) NI versus SHG intensity. (c) ANI versus SHG intensity. (d) AI versus NI: correlation between these independent variables remains linear over the range of disorder detectable in biological samples.

3.3.

Intervertebral Disk Injury

We have previously reported that we could observe differences between SHG scans of control disks and scans of disks subjected for one week to in vivo loading.²⁷ In the present study, we sought to determine if we could quantify these differences and, further, if such differences could be correlated with other measures of disk structure and function. We obtained disk sections prepared from tissue blocks available from studies of in vivo disk loading performed by one of us previously,²⁵ in which compressive loading was compared with tensile loading. Five samples were obtained from each group. The samples analyzed were midsagittal sections, a cutting plane that allows visualization of all of the main structural features in one section. The geometry of the disk components in a midsagittal section is illustrated in Fig. 6 . The complex comprising the vertebral body, end plate, rings, and nucleus pulposus is sometimes referred to as the “spinal motion segment” to denote the structural or functional interplay between these diverse components. All of the components in this region, with the exception of the nucleus pulposus, generate strong SHG signal. The ability of SHG scanning to provide highly detailed images of a structurally complex region—with no exogenous dyes or contrast material—is demonstrated in the SHG intensity scan in Fig. 7 . This scan was extended beyond the “motion segment” to include peripheral supporting tissues, such as muscles, ligaments, fascia, and skin. Although in this study, we confine our analysis to structural changes in the lamellar rings (the region indicated by the white rectangle), other regions may play an equally important role in the mechanisms of injury and repair. One such region is the interface between the cartilage end plate, the lamellar rings, and the vertebral body; we have successfully analyzed structural organization and levels of disorder in this region using the methods described (data not shown).

Fig. 6

Geometry of the intervertebral disk at the midsagittal plane. The midsagittal plane is often used as the cutting plane for analyzing intervertebral disks, as it provides visualization of the major structures: inner and outer lamellar rings and the nucleus pulposus. (a) Schematic view of how the section plane cuts through the middle of the disk. (b) Cutaway view. (c) The section on the disk and the incident laser light present a complex geometrical interface. Legend: NP—nucleus pulposus, OR—outer rings, IR—inner rings.

Fig. 7

Polarization-modulated SHG imaging of intervertebral disk and associated structures. SHG intensity map of the disk section showing all of the structural components in the “spinal motion segment”—the intervertebral disk and associated structures—except the nucleus pulposus. Additional structural components including muscles, ligaments, and skin are also shown. Note that the structures can be clearly distinguished from each other; the signal has not been processed and shows excellent signal-to-noise ratio. Legend: NP—nucleus pulposus, IR—inner rings, OR—outer rings, EP—cartilage endplate, L—ligament, M—muscle, VB—vertebral body, E—epithelial layer.

As described by us earlier, disk loading at $1\mathrm{MPa}$ for one week resulted in deformation of normal ring architecture and qualitative changes in SHG image maps.²⁷ In the present study, we also see marked changes in SHG scans of the annulus fibrosus from disks subjected to compressive or tensile loading as compared with control disks, as well as differences between tensile and compressive loading. As shown in Fig. 8 , compressive loading results in marked deformation of the rings and decreased disk height. Values for disorder $\left(\overline{ANI}\right)$ were obtained in regions measuring $50\times 50\text{pixels}$ ; selected values are shown on the disorder maps in Figs. 8b and 8d. Regional disorder values for control disks never exceeded 4.5%. In contrast, values for $\overline{ANI}$ in loaded disks reached values of 20 to 25%. Disorder in loaded disks had a strikingly heterogeneous pattern. Figure 8d shows a disk subjected to compressive loading characterized by high disorder values in the upper half of the disk; regions in the lower half are much lower. Compressive loading resulted in significantly higher disorder values than tensile loading. Both load and location significantly affected disorder values in the disks. The highest values for disorder were seen in the upper half of disks subjected to compressive loading, with a mean value of 17%. Tensile loading produced slightly greater disorder in the lower half of the disk, but the difference between the upper and lower halves was not significant. Values for $\overline{ANI}$ for all of the groups, as well as statistical analysis by ANOVA, are summarized in Fig. 9 . Although signal intensity was, in general, inversely correlated with disorder, correlations may be absent for several reasons: (a) values for disorder and intensity that are close to median values may fail to show a clear inverse correlation; (b) color scaling in architecturally complex tissues may fail to show subtle differences in some regions; (c) components in a specific region may affect signal intensity but not orientation.

Fig. 8

Comparison of SHG imaging in control and loaded disks. Effects of compressive loading on structural organization in the annulus fibrosus are clearly shown. (a) Control disk, SHG intensity. (b) Control disk, disorder map. Representative values for disorder in $25\times 25\text{pixel}$ regions are shown. (c) Disk subjected to compressive loading, SHG intensity; arrow shows region of ring deformation (arrow) associated with loss of ring height. (d) Compressed disk, disorder map. Note that disorder is predominantly located in the upper half of the disk; high disorder in general corresponds to low intensity. Typical values for $\left(\overline{ANI}\right)$ in $25\text{-}\times \text{-}25\text{pixel}$ regions in this area are 17 to 23%.

Fig. 9

Comparison of disorder in different regions of control and loaded disk. Values for disorder indices obtained by SHG scanning of control and loaded disks ( $n=5$ for each group) were analyzed by ANOVA. Fischer’s PLSD was used for post hoc analysis. Values for $\left(\overline{ANI}\right)$ were obtained from regions of $25\times 25\text{pixels}$ . Disorder is expressed as a percent age of maximum disorder. There were three groups (control, flexion, tension) and four regions analyzed (upper inner, upper outer, lower inner, lower outer). Legend: Bars with differing superscripts are significantly different from each other $(p<0.05)$ .

We performed histomorphometric analyses of the sections subjected to SHG scanning. Typical results from our procedure for generating binary images from grayscale, $40\times$ photomicrographs of disk samples are shown in Fig. 10 . Procedures for determining the mean intercept length in these binary images were described in Sec. 2. To enhance contrast, we used a polarizing lens. We selected 4 to 5 circular regions of the disk from the inner, middle, and outer rings, each approximately $50\mu \mathrm{m}$ in diameter, for analysis (Fig. 11 ). Although methods for image preparation had by far the greatest effect on aspect ratio determination, we also considered whether or not other factors might have an effect. As described in Sec. 2, we analyzed the effects of the number of rotations and the number of grid lines on MIL determination. The number of rotations had a small effect—slightly different results were obtained if only 8 rotations were used or if more than 16 were used—but these differences were not enough to affect discrimination between groups. These results are consistent with the report by Hipp ²⁹

Fig. 10

Morphometric analysis of disorder in disk. Photomicrograph of unstained intervertebral disk used for histomorphometric analysis of disorder using MIL analysis, as described in the text. (a) Photomicrograph, $40\times$ , of lamellar rings of control disk. A $50\text{-}\times \text{-}50\text{-}\mu \mathrm{m}$ circular region (arrow) is used for analysis. (b) The circular region is converted to a binary image.

Fig. 11

Correlation between aspect ratio and disorder index. Values for aspect ratios were obtained in the same sections used for disorder analysis based on SHG scanning. Four to five circular regions, $50\mu \mathrm{m}$ in diameter, were obtained for each computation. Mean values for the groups were used for determining correlations. Correlations between these two measures of disorder were assessed by regression analysis. Linear regression analysis showed there was a highly significant correlation between the parameters $(p<0.001)$ .

We performed regression analysis to determine if correlations existed between SHG indices of disorder and the aspect ratios obtained from morphometric analysis. We found a highly significant correlation $(p<0.001)$ between values for $\overline{ANI}$ and for the aspect ratio obtained in the same sample region (Fig. 10). Our morphometric data are consistent with previously reported data from these models.²⁵

4.

Discussion

In this paper, we describe a novel method for quantifying structural disorder in collagenous tissues using data obtained from PM-SHG scanning. Our algorithm allows us to quantify the magnitude of the disorder and to classify the pattern of disorder. We have validated our method by conducting detailed correlative studies in rat tail tendon fascicles subjected to controlled heating: we found that SHG-derived indices of disorder are significantly correlated with morphometric analyses of disorder. We have also shown that our method can be used to study pathological changes in complex tissues, such as intervertebral disks subjected to loading. Our study raises several important questions. First, what is the physical meaning of disorder, as determined with our method, in a given sample? Second, does the ability to quantify disorder in collagenous matrices have broader implications—beyond identifying pathophysiological changes—for understanding biological systems? Third, does our method have advantages over techniques currently being used to measure disorder? Fourth, what distinguishes our method from techniques described by other groups working with PM-SHG imaging? We address these issues below.

To draw inferences concerning the physical meaning of disorder in a given sample, such as the length scale, we need to know how the tissue is organized, particularly with respect to the size of its component structures. The fundamental unit of organization in all fibrillar collagen matrices is considered to be the microfibril, a structure consisting of five adjacent collagen molecules offset from each other by one-quarter of their length. Only recently have definitive crystallographic data been published detailing the packing structure of collagen molecules within the microfibril as well as the relationship between adjacent microfibrils.³⁰ The microfibril, approximately $300\text{-}\mathrm{nm}$ long and $7\mathrm{nm}$ in diameter, consists of five staggered, pleated collagen molecules with a right-handed supertwist. Neighboring microfibrils are physically connected to each other through interdigitation of their nonhelical C- and N-terminal extensions. Microfibrils tend to form lateral associations during fibrillogenesis, resulting in the formation of discrete bundles up to $70\mathrm{nm}$ in diameter that are often, but not invariably, referred to as fibrils. In each tissue, fibrils tend to have a characteristic size; fibril diameter is closely regulated by environmental conditions as well as by other components present in the extracellular matrix, particuliarly decorin. Tissues vary considerably in structural complexity. Some tissues are characterized by two or three hierarchical levels of organization, at progressively larger scales, beyond the fibril, as exemplified by rat tail tendon. In tendon, the lateral aggregation of microfibrils produces discrete bundles, approximately $30\mathrm{nm}$ in diameter, termed subfibrils. The subfibrils are organized into larger bundles called fibrils, approximately $500\mathrm{nm}$ in diameter. Fibrils, separated from each other by fascial layers are organized into larger structures called fascicles $\left(50\text{to}500\mu \mathrm{m}\right)$ ; several fascicles, encased within fascia, form a single tendon.

We can interpret disorder data obtained from heated tail tendon in relation to these organizational structures. Figure 12 summarizes the pertinent data obtained from control and disordered tail tendon: note that only small sections of the scans are shown so that individual pixels can be clearly seen. In each scan, the tendon sample is placed on the slide such that its axial orientation is parallel to the $y$ axis. Each pixel represents an area of $1\times 1\mu \mathrm{m}$ .

Fig. 12

Structural correlates of disorder in heated rat tail tendon. Small regions, $20\times 20\text{pixels}$ , were selected from PM-SHG scans of control and heated fascicles. Each pixel represents a $1\times 1\text{-}\mu \mathrm{m}$ region of tissue. The first row shows data from a control sample; the second row from a moderately disordered sample (heated at $58°\mathrm{C}$ for $10\min$ ), and the third row from a severely disordered sample (heated at $58°\mathrm{C}$ for $40\min$ ). In each row, the first panel shows SHG intensity, the second panel shows fibril orientation, the third panel shows disorder, and the fourth panel shows the tissue structures present in a typical pixel within the scanned region. (a) SHG intensity in the control tendon shows a strong signal. (b) Fibril orientation in the control tendon is parallel. (c) Disorder in the control tendon is less than 1%. (d) A typical pixel contains two fibrils that are parallel and undamaged. (e) SHG intensity in moderate disorder is decreased but detectable in all pixels. (f) Fibril orientation in moderate disorder shows regions of nonparallel alignment. (g) Disorder index shows mean value is 12%. (h) A typical pixel in moderate disorder contains two fibrils that are no longer parallel but whose components are not damaged. (i) SHG in severe disorder is markedly diminished and some pixels have no detectable signal. (j) Fibril orientation shows loss of parallel orientation and loss of signal. (k) Disorder map shows that the mean value (for those pixels with signal) is 21%. (l) In a typical pixel, two fibrils are distorted and torn and may include damaged fibril components at each level of organization. Legend: Col—collagen molecules, MF—microfibrils, SF—subfibrils.

Table 2

Aspect ratio values and disorder indices in heated rat tail tendon.1

The panels in the first row of Fig. 12 show data from a normal tendon. As we would expect, SHG signal is high [Fig. 12a], mean orientation angles for each pixel are parallel to the axial orientation of the tendon, and also to each other [Fig. 12b], and the disorder index in the region is extremely low [Fig. 12c]. Figure 12d shows a schematic representation of the structural integrity of the components likely to be found within the area of a single pixel. In a $1\times 1\text{-}\mu \mathrm{m}\text{pixel}$ , we might scan across 1 to 3 fibrils; the length of the fibrils, however, most likely extends through the entire sample. Figure 12d shows that for each pixel the fibrils lie parallel to each other and there is no evidence of disruption in the subfibrillar components.

Figure 12’s second row presents analogous scan data obtained from a tendon subjected to brief heating. The SHG signal [Fig. 12e] is detectable in all pixels, but values are lower than for the control tendon. The mean orientation angle [Fig. 12f] in many pixels deviates from the axial direction of the fascicle; angles in some adjacent pixels are not parallel. The mean disorder index for the $20\times 20\text{-}\mu \mathrm{m}$ scan region is 12% [Fig. 12g]. At what level of organization has the disorder occurred in this sample? If disorder had occurred at the level of the microfibril, then by definition it would no longer be energetically favorable for the collagen molecules to remain aligned in their quarter-stagger conformation. They will become randomly oriented with respect to each other, resulting in destructive interference and loss of SHG signal. This outcome is not consistent with the data shown in Fig. 12d: all of the pixels in the region have a detectable signal. Could partial disassembly of the microfibrils occur? This seems unlikely—the interdigitation of the microfibrils would tend to prevent isolated loss of microfibril integrity. More significantly, although the decreased SHG signal could result from partial loss of microfibril structures, the change in the orientation angle of individual pixels cannot result from their loss. Any collagen molecules freed from the lattice would become randomly oriented with respect to each other and would not change the pixel orientation angle. Similar arguments apply to postulating that disorder occurs at the level of subfibril organization. Significant loss of subfibril structures would likely result in loss of the SHG signal, which we do not observe in this sample. The presence of some disaggregated subfibrils could result in decreased SHG signal, consistent with our observations. However, this phenomenon could not account for the large deviations in the mean orientation angle observed in many pixels. To produce a change in the detected orientation angle for a given pixel, a relatively large number of the disaggregated microfibrils would have to simultaneously alter their orientation (relative to the axial orientation of the fascicle) in the same direction. This is not a credible outcome of thermal denaturation. It is more likely that any disaggregated microfibrils would become randomly oriented with respect to each other; the net effect of these microfibrils on the mean orientation angle of the microfibrils that remained aggregated would be negligible, and we would see no effect on the pixel’s mean orientation angle. The only structural change consistent with our scan data would be disorder at the level of the fibrils such that they no longer lie straight—and/or parallel to each other—within the fascicle. This conformational change could result in significant changes in the mean orientation angle for a given pixel (relative to the fascicle orientation) as well as different orientation angles in adjacent pixels, because only one or two fibrils are likely to pass through any given pixel. The associated decrease in SHG signal could be a result of the altered fibril orientation (as there is less contribution from adjacent structures to signal magnitude) and/or partial disruption of microfibrils. This conclusion is summarized in the schematic shown in Fig. 12d: the most likely disruption occurs primarily at the level of fibril organization.

The third row in Fig. 12 presents data from a region of heated tendon characterized by high disorder. The SHG signal is significantly lower throughout the scan than it was in the moderately disordered sample, in addition, almost half of the pixels are below detectable levels, as is clearly seen in Fig. 12f. The mean orientation angles in pixels with detectable signal show marked deviation from the axial orientation of the fascicle. These data suggest that disorder is occurring at two different levels of structural organization. In pixels with detectable SHG signal and deviation of the orientation angle, disorder is likely occurring at the fibril level, for the reasons described above. Only disruption of the fibril structure as a whole could produce such marked changes in the mean orientation angle seen in many of the pixels. However, in pixels with no SHG signal (but which correspond to regions where tissue is present, as confirmed by light microscopy), disorder is likely occurring at the microfibril or subfibril level of organization. As described above, this type of disruption would result in destructive interference and complete loss of signal. Figure 12l shows these two types of disorder schematically: disruption of fibril relationships and disaggregation of fibril subcomponents in approximately half of the scan region. This interpretation of the physical meaning of disorder in rat tail tendon subjected to heat is consistent with a very recent study by Chen ³¹ on the effects of thermal denaturation on rat tail tendon structure. They observed that the earliest changes occurred at the level of fibril organization. The changes consisted of transverse breaks or tears, causing regions of the fibril to become twisted. This change was associated with some decrease in SHG signal magnitude. They attributed this phenomenon to the presence of regions within tendon that are particularly vulnerable to thermal effects. Later changes consisted of global loss of SHG signal, which they attributed to thermal denaturation. In summary, the length scale of disorder can be inferred from the data by considering the size of the organizational structures within the sample in relation to pixel area.

The ability to quantify disorder in collagenous tissues may also contribute to a deeper understanding of the role of collagenous lattices in biological systems. Order-disorder transitions in physical systems represent a fundamental problem of interest in statistical mechanics. In a crystalline lattice, order can be assumed to represent the net effects of the interactions between nearest neighbors, as expressed by the Ising model and its more generalized form, the Potts model. Although these models were developed to analyze spin interactions between particles, their use has been extended to analyze states of complex biological systems.^{32, 33, 34} Collagenous matrices are particularly suited to such analytical methods, as they possess many of the attributes of a crystalline lattice.³⁵ The remarkable strength of collagenous tissues depends on properties resulting from its supramolecular organization, such as interdigitation of microfibrils, supertwist conformation, secondary twisting produced by covalent crosslinking, and others, over the region of required structural functioning, not the properties of the individual molecules. Thus, a small amount of disorder may produce a disproportionately large decrement in function if lattice continuity or strength is sufficiently disrupted in key locations. Very recent studies suggest that these “crystalline” attributes of collagen may be more dynamic than was originally postulated. Gutsmann ³⁶ have presented data suggesting that collagen fibrils are best modeled as hollow tubes rather than as solid rods and that the molecules within the tube may flow from one region to another in response to physical perturbations. The ability of our analytical method to assess small, local changes in lattice structure may be well-suited to studying interactions between the matrix and its environment from the perspective of this new structural model of collagen. The ability of the collagen lattice to transmit information within a biological organism has been a topic of recurring interest, particularly in relation to inflammatory response, would healing, bone repair, fibrosis, tumorigenesis, and other pathological conditions. Of particular interest are studies suggesting that conformational changes within the lattice structure itself may allow signalling via piezoelectricity.³⁷ Disruptions in lattice continuity represent a common event in many pathological processes, including mechanical trauma (as in our disk model), stochastic perturbations (as in our heat-treated fascicles), tissue repair and regeneration, and neoplastic invasion. We are currently conducting studies in which we are applying our method to some of these other pathological processes.

There are surprisingly few methods designed to measure disorder in complex biological systems such as tissues. The most commonly used approach is the method of directed secants, originally described by Saltikov in 1958.³⁸ A significant modification, developed by Whitehouse,²⁸ is known as the the mean intercept length (MIL) method. An object is analyzed by obtaining an image of it and converting the image to a binary form. It is a form of topological analysis designed to determine if a predominant orientation of the image components exists. This method continues to be used widely in many different applications, as described by Stroeven, ³⁹ including biomedical image analysis.^{40, 41, 42, 43, 44, 45} Pattern isotropy is not, of course, synonymous with disorder, and there are many cases where this approach will yield misleading results. For example, an intervertebral disk sectioned along the posterior sagittal plane will transect several adjacent lamellar rings. Morphometric analysis of such a section using MIL determination shows a high level of pattern isotropy, with an aspect ratio over 0.9 (unpublished data); the tissue is, in fact, highly ordered. Image analysis in the MIL method requires replicable methods for thresholding and binarization of images, a process that has inherent problems in standardization.⁴⁶ Disorder in tissues is also assessed using indirect methods that rely on the presence of correlative phenomena known to be associated with disorganization. Examples include the Thompson scale for measuring degenerative changes in intervertebral disks and the Mankin scale for measuring changes in cartilage.⁴⁷ Such “semiquantitative” scoring systems are limited in several respects: most of the markers are unrelated to organization of the matrix, observer bias is difficult to avoid, and the scoring system provides only three or four scoring levels. All of these methods require destruction of the tissue and require long preparation steps. All have significant confounding variables. SHG-derived analysis of disorder offers significant advantages over these methods.

The possibility of using PM-SHG scanning to assess disorder in macromolecular arrays was first suggested by Roth and Freund.⁹ In a detailed study of age-associated changes in the ratio of the hyperpolarizability tensor elements in rat tail tendon collagen, they found that this ratio decreased with age. They noted that SHG polarization measurements are sensitive to “the degree and kind of molecular order.” In several elegant studies, Webb and coworkers have also raised the possibility of using SHG imaging to examine structural features, including the degree of disorganization, in collagenous matrices.^{2, 48} Recently, they reported on a detailed study of the ratio of forward-scattering to backscattering SHG in fascicles from young and old rats; their data provide further evidence that fascicles may exist as “hollow tubes” rather than as “solid rods.”³⁶ Although these workers raised the issue of studying disorganization with SHG imaging, they did not describe specific methods for doing so.

Other investigators have reported on applications of SHG polarimetry. Yasui ^{49, 50} mapped the orientation of collagen fibrils in cadaveric skin, using a backscattered geometry. Polarization angle was rotated by mechanical rotation of a half-wave plate. The authors describe a value they refer to as “polarization anisotropy” (PA) that they describe as a measurement of organization: it represents the ratio of SHG intensity obtained at orthogonal orientations averaged over the entire scan. They conclude that tissues such as reticular dermis contain progressively “entangled” fibers, in contrast to tissues such as Achilles tendon. PA appears to be a measurement of structural organization at a more macroscopic tissue level than the level of fibrillar organization associated with functional properties of the matrix. Applegate ⁵¹ describe a similar approach for assessing organization in biological samples, based on the ratio of SHG intensities at orthogonal polarization angles. This study represents a somewhat different imaging method, in that SHG is combined with optical coherence tomography. Data are presented for collagen orientation in fish scales. In neither of these studies are the anisotropy parameters correlated with other measures of disorder. Further, it is unclear at what scale of organization order or disorder will be represented. Theodossiou ⁵² used SHG intensity measurements to track thermal denaturation in collagen samples; fluorescence was detected concomitantly. The incident laser beam was not polarized. The authors state that this approach has direct clinical implications, in that the ability to detect phase transitions can be used to direct laser ablative procedures, and thus prevent thermal damage to the surrounding tissue. Although this method is certainly identifies the final stage of tissue disorganization, it does not appear to provide insight into initial and intermediate states of disorganization in tissues.

Several additional methods have been described for inferring structural information from SHG scanning. Cox have presented data purporting to use SHG to discriminate between genetically distinct fibrillar collagens; these data are based on properties characteristic of arrays of these collagens, such as fibrillar diameter, rather than on specific differences in primary structure.^{53, 54} In analysis of several pathological conditions, we did not find these criteria sufficient to discriminate reliably among types I, II, and III. More recently, Cox ⁵⁴ have addressed issues such as colocalization of SHG with fluorescence signal and identification of collagen degradation in pathological conditions. Plotnikov ⁵⁵ have reported on a method for analyzing sarcomere periodicity using fast Fourier transformation analyses of $15\times 15\text{-}\mu \mathrm{m}$ regions. Data are expressed as a frequency distribution of sarcomere period length; the authors suggest that these spatial frequencies can be used for quantitative comparisons. This method could potentially be used to analyze disorder. Boulesteix ⁵⁶ describe a method for analyzing the SHG signal from sarcomeres that achieves a computed resolution of $20\mathrm{nm}$ . They used a $0.9\mathrm{NA}$ , $60\times$ objective to achieve a spatial resolution of $0.45\times 2\mu \mathrm{m}$ , which was used to visualize sarcomeres within living myocytes exposed to a paralytic toxin capable of inducing sarcomere contracture. Both of these approaches appear to be specific for analyzing changes in myosin.

Finally, we consider the rationale for our signal analysis method in the context of the many new developments in signal processing. Our goal was to develop a scalable analytical method capable of (1) detecting disorder that corresponded to disruption of matrix stability, (2) identifying small regions of disorder within large ordered areas, (3) discriminating among different patterns of disorder, (4) providing quantitative information about disorder so that real-time changes could be tracked. These features seemed to be particularly desirable in analyzing effects of aging, certain diseases, trauma, and tissue engineering strategies on tissue response. To our knowledge, no available method offered all of these features. In developing our own algorithm, we identified two independent parameters associated with disorder—a feature that allows us to analyze patterns of disorder by looking at changes in the ratios of these parameters. These features are not available in methods that rely on, for example, the comparison of SHG intensity at two orthogonal polarization angles.

Our analytical method does not address changes occurring at higher levels of matrix organization, unless such changes directly affect normal fibril alignment. Disks subjected to loading, for example, had regions where the borders of the lamellar rings appeared undulant. We could detect disorder in these regions if it was associated with derangements in fibrillar alignment; however, such macroscopic changes were also observed to occur in the absence of disorder as we defined it. Although these phenomena are readily observed in our orientation maps, these changes are not quantified by our disorder algorithm. We explored a number of data analysis methods capable of detecting disorder at this higher level of structural organization. Texture analysis seeks to define texture features in image data using co-occurrence matrices or application of mathematical filters.⁵⁷ Moment invariant analysis of images can be accomplished through any of several techniques, such as Gabor filters,⁵⁸ and represents another approach for discriminating between normal and abnormal tissues. A wide range of image morphometry techniques are available such as dilation, erosion, segmentation, watershed, digital skeletonization, and so forth. Sometimes these methods are used for converting grayscale images to binary images for automated processing; however, they may be used as stand-alone methods to enhance certain features of scan data.^{57, 59} A novel Bayesian segmentation technique has been described for segmentation of microscopic nuclear magnetic resonance data for accurate classification.⁶⁰ Ultimately, we found that quantitative analysis of structural changes at higher levels of organization are most effectively handled using wavelet decomposition, a technique that can be applied to the same SHG data sets analyzed by our current algorithm. We have presented a preliminary report of our findings and are currently developing an integrated analytical system that incorporates both methods.

5.

Conclusion

We have shown that we can use PM-SHG imaging to obtain quantitative information about the magnitude and the pattern of structural disorder in biological tissues, even those with complex geometry containing mixtures of collagen types as well as other matrix components. We have validated our analytical method through detailed correlations between SHG-based data and conventional morphometric analysis in two tissue models of disorder: (1) heat-induced disruption of fascicle and (2) load injury in intervertebral disks. Our studies suggest that the ability to detect very mild disorder, including transient changes, may be of particular value in identifying pathological disturbances at a stage when they may still be reversible.