Open Access
Add to BibSonomyAbstract
Image correlation spectroscopy (ICS) is known to be a useful tool for the evaluation of fiber width in the extracellular matrix. Here we evaluate a more general from of ICS fit parameters for fiber networks and arrive at a means of quantifying the fiber density, pore size and length which facilitates the characterization of the extracellular matrix. A simulation package was made to create images with different structural properties of fiber networks such as fiber orientation, width, fiber density and length. A pore finding algorithm was developed which determines the distribution of circular voids in the image. Collagen I hydrogels were prepared under different polymerization conditions for validation of our pore size algorithm with microscopy data. ICS parameters included amplitude, standard deviation and ellipticity and are shown to predict the structural properties of fiber networks in a quantitative manner. While the fiber width is related to the ICS sigma; the fiber density relates to the pore size distribution which correlates with the ICS amplitude in thresholded images. Fiber length is related to ICS ellipticity if the fibers have a preferred orientation. Findings from ICS and pore distribution algorithms were verified for both simulated and microscopy data. Based on these findings, we conclude that ICS can be used in the assessment of the extracellular matrix and the prediction of fiber orientation, width, density, length and matrix pore size.
©2012 Optical Society of America
1. Introduction
Correlation spectroscopy uses autocorrelation to detect the likelihood of repetitions in a signal at a certain lag. In the case of two dimensional and spatially resolved data, this technique is called image correlation spectroscopy (ICS) [1,2] and has been used in numerous applications such as the estimation of fluorescent particle diffusion [3,4], degree of aggregation of plasma proteins [5] and the characterization of the extracellular matrix by analyzing the collagen fiber network. Previously Raub et al. [6,7] linked parameters from ICS analysis of multi-photon microscopy images with the bulk mechanical properties of collagen hydrogels. Here we evaluate how the properties of virtually computed fiber matrices affect parameters extracted from ICS and whether other parameters from the ICS analysis correlate with properties of fibrous networks.
In ICS, once the autocorrelation function has been calculated, it is fitted to a suitable function to obtain a set of parameters. Here we discuss how the fit parameters of a more generalized Gaussian function correlate with variations in fiber density (fill factor), fiber length, fiber width and average orientation of the fibers. Our fit parameters include the amplitude, orientation, standard deviation (sigma), constant offset and ellipticity. With the help of simulated data, we establish a relationship between the fit parameters and the fiber properties.
Raub et al. have shown how the collagen fiber width, length and pore size are related to polymerization conditions [7]. They found that ICS sigma correlates well with average fiber diameters while modest correlation was found in an attempt to estimate pore size based on ICS sigma [7]. The ICS’s ability to create other parameters that are related to the voids within the network (pore size), fiber length and the average fiber orientation has not yet been shown.
When employing fibrous networks in biologic experiments, fiber diameter and length are important parameters to characterize the network’s mechanical properties. When evaluating the ability of cells or particles to migrate through the network, its pore size plays an important role [8,9]. In addition to the ICS technique, we have developed an algorithm through which we assess the pore distribution of the fibrous networks and we show that with appropriate image pre-processing steps the ICS fit parameters become correlated with the pore distribution.
Our findings are further verified by applying the technique to second harmonic generation (SHG) images of collagen type I hydrogels. By using polymerization pH, we altered the structural properties of collagen gels [6,7,10] and applied our analysis technique to those images.
2. Materials and methods
2.1. ICS theory
ICS entails the calculation of the autocorrelation function (ACF) of an image. The ACF is calculated more efficiently in the Fourier domain where the multiplication of the Fourier transform of the image (i) and its complex conjugate gives the power spectral density (PSD) and the inverse Fourier transform of the PSD the ACF as shown in Eq. (1):
where u, v are pixel locations and x, y lags in pixels. As the ACF is concentrated about the center of the image, for structural analysis higher spatial lags can be removed by cropping a central portion of the ACF. Since the ACF at zero lag is related to noise, it is also removed in our analysis. In our case, we analyzed a central 12*12 portion of the ACF for a 256*256 pixel image. We determined the size of the 12*12 portion visually by comparing the ACF to a fitted function as illustrated in Fig. (1) . The crop size should be minimized with the goal to cover large ACF values and representing the ACF drop-off with least amount of error. The crop size was kept constant for our study.
Fig. 1 Illustration of Image Correlation Spectroscopy and the dependence of the autocorrelation on orientation. The image with fibers oriented horizontally and its autocorrelation is shown. The Gaussian fit of the centrally cropped ACF is to the right with the wireframe representing the fit.
2.2. Generalized Gaussian fit to describe the ACF
Since the ACF is affected by the predominant orientation of the image’s structural element, we modified the existing approach of parameter extraction [11] and included orientation as an additional analysis parameter. For illustration, Fig. (1) shows the dependence of the ACF on the oriented structural elements where those elements are fibers oriented primarily horizontally. In our modified approach, the ACF is fitted to a generalized 2-D Gaussian function shown in Eqs. (2).
In Eqs. (2), x0, y0 are the center of the ACF, σM, σm the standard deviation along the major and minor axis, theta the orientation of the axis with respect to the horizontal, A the peak amplitude and B the average offset. This allows us to summarize average structural properties of an image with five parameters. We term the ratio of the major axis standard deviation to the minor axis standard deviation as ellipticity which is another parameter in our ICS analysis.
We modified an existing implementation of ICS [11] with the above function. The fitting is optimized in a least squares sense and the fit initialization conditions were adjusted to facilitate convergence.
2.3. Simulation of images
To verify the effectiveness of using ICS as a tool to provide a quantitative assessment of a fiber network, we simulated binary images with fibers of various properties. Lines represented fibers and were varied in length, thickness, orientation and density. Two types of random number generators were used, one creating values with a normal distribution (length and width) and the other with a uniform distribution (position, orientation). The starting point of each fiber was randomly selected over a range representing the desired image plus additional padding accommodating for average fiber length and variation. The orientation was uniformly selected in the specified range and an offset of 90 degrees was randomly added allowing the fibers to extend “forward” but also “backward” from the starting point.
Table (1) describes the simulated parameter space. For each simulated condition we generated ten random images and results are reported as the average and standard deviation of those images. For example, we kept the mean fiber length, fiber density constant and varied the mean fiber width with the fibers following a normal distribution. This gave us images to evaluate the influence of fiber width on the fit parameters. This approach was repeated for varying mean fiber content and length. To evaluate fiber length, fibers oriented with different uniform distributions ranging from a dominant to completely random orientation were simulated for varying fiber lengths. This enabled us to show the dependence of fiber orientation and length on ellipticity.
Table 1. Simulation variables
2.4. Method for estimating the pore distribution
Since our fiber network simulation program produced a random image without knowing the void size and distribution we developed a technique to estimate the average void size. The void size is approximated by finding the largest circular object fitting into a void (pore).
The flowchart to calculate the pore distribution is shown in the Fig. (2a) and can be illustrated by the following example: The input image is a binary image which has been simulated with a fill factor of 0.5, mean fiber length of 50 pixels, mean fiber width of 2 pixels and the orientation following a uniform distribution between [-90, 90] as shown in Fig. (2b). A circular mask with a diameter of 10 pixels is convolved with the image (Fig. 2b). Areas with zero intensity are identified in the convolved image, which corresponds to voids where the mask is capable of fitting. Figure (2c) illustrates all regions with intensity of zero (black area) however, independent regions with an area of one pixel are the locations where the circular mask fits exactly (marked with crosses). When increasing the mask size, the independent regions shrink. They often do not shrink to exactly one pixel and the smallest area is several pixels in size. Therefore we choose an arbitrary threshold of 10 pixels and recorded the centers of all independent regions with a size smaller than this threshold. Independent regions were identified with connected component labeling (Image processing toolbox, Matlab, MathWorks, Natick, MA) [12] which provides the centroid and area for each region. We carried out the entire process for masks representing all diameters of interest and maintained a list of possible pore centers. Centroid coordinates were rounded to the nearest integer pixel and if more than one pore size was found at the same location, the one associated with the largest diameter was kept (Fig. 2d). The list of centroids and diameters enables us to create the pore distribution as shown in Fig. 2(e). Our experimental results indicated that a lognormal function fits the pore distribution best.
Fig. 2 Flow chart and illustration of pore finding algorithm. (A) Steps entailed in finding the pore distribution of an image. The algorithm produces a list of pore centers and the associated diameters. (B) The image is convolved with a circular mask. To find the pore distribution, circular masks with different diameters are generated. Zeros of the convolved image are identified which corresponds to the points where the circular mask fits into the pore. (C) Connected component labeling is used to identify independent regions that have an area smaller than a set threshold. For illustrative purpose, centers of those regions are marked with crosses (C left) and the corresponding circular pores are drawn in the subfigure to the right. (D) Simulated image with all pores identified. (E) Pore diameters [pixel] follow a lognormal distribution.
2.5. Microscopy of collagen type I hydrogels
Imaging was performed with a single beam intravital microscope (TrimScopeI, LaVisionBioTec, Bielefeld, Germany).The Titanium Sapphire laser was tuned to 780nm and the field of view was about 350 microns. SHG images of collagen I hydrogels were recorded similar to previous publication of our group [13]. An image stack was formed by varying focus in the sample and five consecutive images were analyzed for each sample. The average and standard deviation of those analyses is reported for each analyzed hydrogel. 5ml collagen mixture was prepared and with 35ul 1MNaOH added vortexed to mix. 250ul mixture was removed and added to the first well of a 48 well plate sitting on ice. 5ul NaOH was added to the remaining mixture and vortexed, while another 250ul was removed to the second well. This was repeated for 13 wells and the volume of NaOH was reduced to 2.5ul for wells 14 through 15 resulting in collagen gel pH varying from 6 to 9. For this study repeatable polymerization outcome under different pH conditions was not necessary as we primarily wanted to study hydrogels with varying structural properties.
2.6. Preprocessing of microscope images for ICS and pore finding algorithm
Our microscopy grayscale data was subjected to Wiener filtering for noise removal (3x3 kernel, Image processing toolbox [14], ) which improves the output of the pore distribution algorithm for our imaging data. The raw grayscale images are not suitable for thresholding directly as the intensity on the optical center is brighter than at its edges. Adaptive histogram equalization (Image processing toolbox) [15] was applied to create uniform grayscale distribution in 8x8 sub areas of the image (Fig. 3 ). Alternatively the images could be corrected for nonuniform illumination by measuring the spatial system response. In order to create a binary image, the images were then thresholded using Otsu’s method (Image processing toolbox) [16]. As shown in the flowchart in Fig. (3), once the binary image is obtained we conduct ICS and Pore distribution analysis.
Fig. 3 Steps involved in ICS and pore distribution analysis of measured data. The pore finding routine and ICS is applied on images after contrast adjustment, adaptive histogram equalization and thresholding. Image (1) illustrates the raw data, (2) the histogram adjusted data, (3) the thresholded image, and (4) illustrates the pore centers and diameters.
It is possible to conduct ICS directly on a gray scale image. We have tested ICS on raw data, adaptive histogram corrected gray scale images and on threshold images and found that for comparing ICS with our pore routine, analysis on thresholded and histogram adjusted images created most reproducible results (data not shown).
3. Results
Representative images simulated for different widths, fiber densities and lengths are shown in the Fig. (4) . ICS sensitivity analysis was conducted on those images as well as the pore finding algorithm was applied to identify if any of the ICS parameters correlated with the pore diameter.
Fig. 4 Images of virtually computed fiber matrices under different simulation conditions. The top row contains images of different widths with a constant mean fiber density and length. The middle row has images of different fiber densities with constant mean width and length. The bottom row illustrates images of different lengths with constant mean fiber density and width.
3.1. ICS parameter sensitivity
When plotting ICS fit parameters versus fiber properties, it can be observed in Fig. (5) that with increase in fill factor; the amplitude, sigma and offset decreases. However the change in amplitude was more strongly affected by fill factor then the other parameters. Similarly an increase in fiber width increases amplitude, sigma and offset; however the increase in sigma is more pronounced than amplitude and offset. It can also be observed that while the fiber length correlates weakly with the amplitude and offset, it is a strong function of ellipticity when fibers show a preferred orientation. We can conclude that the amplitude is influenced the most by the fill factor, sigma by the width and ellipticity by the length of the fibers subject to a dominant orientation. When comparing amplitude, sigma, offset and ellipticity versus preferred orientation, we found that amplitude is slightly increased with preferred orientation (data not shown), but the change is more than an order of magnitude smaller than the change induced by the fill factor. Offset was very little affected by preferred orientation as changes were within the error of our simulations. Ellipticity was strongly affected by preferred orientation because sigma along the major axis was strongly affected by preferred orientation.
Fig. 5 ICS parameter sensitivity analysis. The top left figure illustrates sensitivity of fit amplitude with fill factor, width and length. Symbols represent and average of 10 simulations. Standard deviation is also plotted but often smaller than the symbol labeling the plot. The amplitude is highly correlated with fill factor. The top right figure illustrates the variation of average sigma, with sigma shown to be highly correlated with fiber width. The average sigma is the mean of sigma on the minor and major axis. The bottom left figure illustrates the offset with offset shown to be highly correlated with fill factor. The bottom right figure illustrates ellipticity to be highly correlated with fiber length if the fibers have a preferred orientation. Orientations follow uniform distributions from [-5°, 5°] (Orient = 5) to [-80°, 80°] (Orient = 80). Uniform distribution of [-80°, 80°] indicates almost randomly oriented fibers and the ellipticity (σM/σm) is almost flat around 1.
3.2. Pore size distribution
The simplest approach to simulate a fiber network with smaller pore size is to increase the fill factor in our simulations. We calculated the pore size distribution for different fill factors ranging from 10% to 80%. We found that the fill factor correlates well with average pore size which was obtained from a log normal fit as shown in Fig. (2e). We also observed that the pore size variation within an image did not show a trend with the fill factor (data not shown). Then ICS analysis was conducted on the same images to confirm whether ICS parameters correlate with the pore size. We found that among all the ICS parameters, ICS amplitude varies the most with changing mean pore diameter (Fig. 6 , left). That sensitivity was stronger than the influence of offset and any of the other ICS parameters (data not shown).
Fig. 6 Correlation of ICS amplitude with average pore size of simulated (left) and microscopy data (right). The microscopy data (n = 15) spans a pore size range that is much smaller than the simulated data.
3.3. ICS parameters from microscopy data
In order to verify our findings on microscopy data, we prepared collagen gels under different polymerization conditions. Our raw grayscale images showed that for low and high pH values, the SHG signal level is low and higher signal levels were obtained when polymerization was carried out at room temperature as compared to 37 degrees C (data not shown). The raw grayscale images were preprocessed as described above and subjected to ICS as well as pore distribution analysis. After preprocessing the fill factor varied between 0.5 and 0.8. For that limited range, the ICS amplitude correlated as expected with fill (R2 = 0.92) and ICS offset had little correlation with fill factor or (R2 = 0.23) and sigma (R2 = 0.36) although trends were consistent with the simulated data (data not shown). ICS amplitude was plotted versus average pore size (Fig. 6, right) and it can be observed that the mean pore diameter follows a linear trend with the ICS amplitude. Also, our experimental results confirmed that the pore size had the best correlation with ICS amplitude as previously found in the simulated data. Using a linear model, the goodness of fit (R2 = 0.901) indicated that the pore size correlates reasonably well with the ICS amplitude as opposed to for example the ICS offset (R2 = 0.01).
4. Discussion and conclusion
With the results on simulated images we have demonstrated how our ICS parameters correlate with the structural properties of virtual fiber networks. While the Gaussian fit was dependent on the orientation of the structural elements in the image, the width of the fibers was found to strongly correlate with the sigma. The fiber length is a major structural parameter and we found a strong dependence of fiber length on ellipticity when fibers tend to align in a particular direction. In simulated data amplitude is decreasing the most at small fill factors while ellipticity changes the most for shorter fibers.
When comparing the ICS fit parameters from microscopy data we found the same trend as observed in our simulated data. There was a strong correlation of ICS amplitude with fill factor. We did not find correlation between ICS offset and any of our computed parameters. Our microscopy data covered a fill factor of 0.5-0.8 where ICS amplitude change is smaller, never the less we were able to obtain a high linear correlation coefficient. Ellipticity in our microscopy data varied between 1 and 1.6 indicating we had a small preferred orientation with relatively short fibers. In that range ellipticity changes the most which might explain why we could observe an increase in ellipticity with increased pH above pH 7 (data not shown). Others have reported that with increase in the basic environment, fibers become longer [7] which can lead to an increased ellipticity if there is a preferred orientation in the field of view.. This is interesting as a preferred orientation is not obvious in our images. This suggests that only a small deviation from uniform orientation is necessary to observe this effect. During cell migration in collagen hydrogels, it is likely that cells will create preferentially oriented fibers along their migration track and ICS ellipticity might be useful in analyzing fibers along those tracks [13]. Additionally, collagen fibers in tissue are generally not randomly oriented and ICS ellipticity might be a useful quantitative parameter for tissue imaging.
Our pore finding algorithm is an effective method to describe pore size distribution. Because ICS amplitude decreases with an increase in fiber density we were able to compare the pore size results with ICS. We found a correlation between the mean pore size and the ICS amplitude. This has been validated with our simulations as well as our SHG image data. Pore size of our microscopy data was at the lower range of the simulated data where it is more difficult to observe a change, however we still found a good correlation coefficient. In order for the ICS amplitude to correlate with pore size it was important to correct the images for nonuniform system response and to threshold our data.
All our analysis was conducted on two dimensional images. The ACF and the Gaussian fit can be expanded to three dimensions. Similarly the pore finding algorithm can be expanded with a third dimension. While analysis with a third dimension will not change the ability to assess fiber diameter, it would affect predicted fiber length. In confocal microscopy the resolution along the optical axis is about three times worse than the planar resolution which allows estimating fibers’ lengths even when they are at oblique angles in respect to the measurement plane, however their length is underestimated. When comparing microscopy images in Fig. 3 with the simulated images in Fig. 4 one can observe that fibers can exhibit tapering towards the end as well as splaying. Our line algorithm does not simulate these characteristics but it is likely that they lead to an underestimation of fiber length.
We can conclude that the structural properties of collagen fibers have been evaluated with Image Correlation Spectroscopy and we have shown that ICS parameters correlate with fiber width, length, average orientation and pore size of simulated fiber networks. Pore size calculations were verified on collagen I hydrogels and showed good correlation.
Acknowledgments
This research was supported by NIH R01 HL077683.
References and links
1. A. Nohe and N. O. Petersen, “Image correlation spectroscopy,” Sci. STKE 2007(417), pl7 (2007). [CrossRef] [PubMed]
2. D. L. Kolin and P. W. Wiseman, “Advances in image correlation spectroscopy: measuring number densities, aggregation states, and dynamics of fluorescently labeled macromolecules in cells,” Cell Biochem. Biophys. 49(3), 141–164 (2007). [CrossRef] [PubMed]
3. M. Srivastava and N. O. Petersen, “Diffusion of transferrin receptor clusters,” Biophys. Chem. 75(3), 201–211 (1998). [CrossRef] [PubMed]
4. P. W. Wiseman, J. A. Squier, M. H. Ellisman, and K. R. Wilson, “Two-photon image correlation spectroscopy and image cross-correlation spectroscopy,” J. Microsc. 200(1), 14–25 (2000). [CrossRef] [PubMed]
5. A. Nohe, E. Keating, T. M. Underhill, P. Knaus, and N. O. Petersen, “Dynamics and interaction of caveolin-1 isoforms with BMP-receptors,” J. Cell Sci. 118(3), 643–650 (2005). [CrossRef] [PubMed]
6. C. B. Raub, V. Suresh, T. Krasieva, J. Lyubovitsky, J. D. Mih, A. J. Putnam, B. J. Tromberg, and S. C. George, “Noninvasive assessment of collagen gel microstructure and mechanics using multiphoton microscopy,” Biophys. J. 92(6), 2212–2222 (2007). [CrossRef] [PubMed]
7. C. B. Raub, J. Unruh, V. Suresh, T. Krasieva, T. Lindmo, E. Gratton, B. J. Tromberg, and S. C. George, “Image correlation spectroscopy of multiphoton images correlates with collagen mechanical properties,” Biophys. J. 94(6), 2361–2373 (2008). [CrossRef] [PubMed]
8. P. Friedl and K. Wolf, “Plasticity of cell migration: a multiscale tuning model,” J. Cell Biol. 188(1), 11–19 (2010). [CrossRef] [PubMed]
9. F. Sabeh, R. Shimizu-Hirota, and S. J. Weiss, “Protease-dependent versus -independent cancer cell invasion programs: three-dimensional amoeboid movement revisited,” J. Cell Biol. 185(1), 11–19 (2009). [CrossRef] [PubMed]
10. B. A. Roeder, K. Kokini, J. E. Sturgis, J. P. Robinson, and S. L. Voytik-Harbin, “Tensile mechanical properties of three-dimensional type I collagen extracellular matrices with varied microstructure,” J. Biomech. Eng. 124(2), 214–222 (2002). [CrossRef] [PubMed]
11. D. Kolin, “ICS Tutorial” (Cell Migration Consortium, 2006), http://www.cellmigration.org/resource/imaging/imaging_resources.shtml#software.
12. R. Lumia, “A new three-dimensional connected components algorithm,” Comput. Vis. Graph. Image Process. 23(2), 207–217 (1983). [CrossRef]
13. N. D. Kirkpatrick, S. Andreou, J. B. Hoying, and U. Utzinger, “Live imaging of collagen remodeling during angiogenesis,” Am. J. Physiol. Heart Circ. Physiol. 292(6), H3198–H3206 (2007). [CrossRef] [PubMed]
14. J. S. Lim, Two-Dimensional Signal and Image Processing, Prentice Hall Signal Processing Series (Prentice Hall, Englewood Cliffs, N.J., 1990), pp. xvi, 694.
15. K. Zuiderveld, “Contrast limited adaptive histogram equalization,” in Graphics Gems IV (Academic Press Professional, 1994), pp. 474–485.
16. N. Otsu, “Threshold selection method from gray-level histograms,” IEEE Trans. Syst. Man Cybern. 9(1), 62–66 (1979). [CrossRef]
