Measuring the anisotropy of features in biological images is of increasing interest as the degree of alignment can inform on both the underlying cellular behaviors and material properties of the sample. For example, different cell types have unique emergent capacities to align in culture, controlled by cellular packing and geometry (Duclos et al., 2014, Duclos et al., 2017). When cells are exposed to external forces such as cyclical stretch, they tend to align their cytoskeleton perpendicular to the axis of stretch (Standley et al., 2002; Yoshigi et al., 2005; Livne et al., 2014). Subcellular cytoskeletal networks can also spontaneously organize in response to the stresses of their environment (Gupta et al., 2015, 2019) or changes in biochemical signalling (Ridley and Hall, 1992). Similar changes in cytoskeletal architecture can be seen in reconstituted protein systems (Falzone et al., 2013; Linsmeier et al., 2016). Additionally, the extracellular matrix (ECM) also has an inherent capacity to align in different tissues or pathologies, such as cancer (Ouellette et al., 2021) or tissue fibrosis (Park et al., 2020; Mascharak et al., 2021), which is thought to alter the ECM network mechanical properties. These examples illustrate the variety of environments where alignment of features reveals important biological properties. As such, it is necessary to develop approaches to efficiently quantify anisotropy of features across a range of length scales, from subcellular organization to tissue level alignment.
A number of different approaches have been developed to quantify the alignment of image features. These methods can be categorized based on the type of algorithm used to highlight features and subsequently quantify anisotropy. Fiber segmentation tools, such as the Fiji Ridge Detection plugin (Lindeberg, 1998), the curvelet-based CurveAlign/CT-FIRE suite (Bredfeldt JS. et al., 2014; Bredfeldt J. S. et al., 2014; Liu et al., 2017, 2020), and Filament Sensor (Eltzner et al., 2015), provide individual fiber information. Tools based on representation/transformation of the image, such as the Fiji plugins OrientationJ (Rezakhaniha et al., 2012) and FibrilTool (Boudaoud et al., 2014) or the CytoSpectre suite (Kartasalo et al., 2015), supply overall fiber alignment information. Hybrid tools, such as TWOMBLI, which exploit a combination of approaches (i.e., fiber segmentation followed by matrix representation of the image) have also been recently made available (Wershof et al., 2021).
Many image transformation algorithms rely on Fourier transformation of the image and exploit its frequency space representation to obtain alignment information (Chaudhuri et al., 1987; Pourdeyhimi et al., 1997; Marquez, 2006; Sander and Barocas, 2009; Goldyn et al., 2010; Kartasalo et al., 2015; Clemons et al., 2018). These algorithms tend to be computationally efficient (Sander and Barocas, 2009) and insensitive to the image background noise (Liu et al., 2020), which is desirable for rapidly processing large volumes of experimental data (Püspöki et al., 2016). Another desirable trait for an alignment detection tool is its flexibility in terms of size and location of alignment patterns (Püspöki et al., 2016), which is handled well by the only open-source software available in this category (CytoSpectre (Kartasalo et al., 2015)) due to a power spectrum analysis of the image’s discrete Fourier transform.
One feature missing from the available tools in the literature is an analysis of the alignment length scale. Depending on the scientific question of interest, one may want to measure alignment over a small region of the image (e.g., local cytoskeletal features inside cells) or over a broad length scale (e.g., ECM organization in whole tissues). Additionally, when comparing experimental datasets, it is important to understand the precise length scale over which the observed alignment is significant and determine the spatial decay of the anisotropy, which will allow for better interpretation and analysis of results (e.g., is the alignment of a region of interest a subcellular or supracellular phenomenon). Here we explain the implementation of an open-source alignment quantification algorithm, AFT − Alignment by Fourier Transform, that can be run on a variety of biological images and has a number of advantages over pre-existing approaches. This FFT-based quantification is rapid and computationally efficient, flexible to different types of microscopy images and samples, easy to use with a low number of parameters, and importantly allows for a user-defined length scale over which to measure feature anisotropy.
2 Materials and Methods
2.1 Biological Images
Sample datasets of microscopy images were used to illustrate the analysis and are provided as Supplementary Material. In particular, cultured fibroblasts transfected with GFP-actin, fibroblasts fixed and stained for filamentous actin (F-actin, phalloidin), fixed cell-derived matrices immunostained for actin (before decellularization) and fibronectin (after decellularization), and second harmonic generation imaging of tissue samples were utilized. Fixed actin images were acquired at different magnifications using either single cells or cells at confluence, to allow for evaluation of the length scale of alignment. Y-27632 treated cells were incubated in the indicated concentration of the drug for 30 min and then fixed and stained with phalloidin.
Fixed samples were imaged using a Zeiss LSM 880 equipped with a 40x NA 1.3 Plan-Apochromat oil objective or 63x NA 1.4 Plan-Apochromat oil objective. Decellularized matrices were imaged using a Zeiss LSM 880 equipped with 63x NA 1.4 Plan-Apochromat oil objective. For second harmonic generation imaging, tissue sections were imaged using a Zeiss LSM 7MP equipped with a 20x NA 1.0 water immersion objective. Y-27632 treated cells were imaged on a 3i Marianas Imaging System consisting of a Zeiss Axio Observer 7 inverted microscope attached to a Yokogawa W1 Confocal Spinning Disk using a 63x NA 1.4 Plan-Apochromat oil objective.
2.2 Synthetic Data
Images of fibrillar features were simulated in MATLAB (Mathworks, v2018b) to create sets of synthetic data with specific length scales. To this aim, squared units with a pre-defined size were randomly distributed in the field of view. Each unit was filled with parallel lines with 5 px spacing, and then rotated with a random angle between 0° and 180°. Three 10-image sets were generated, with unit sizes equal to 100, 200, or 400 px.
Image pre-processing is not required for this analysis, but can be used to additionally highlight the fibrillar features. In the context of this work, pre-processing was performed only for the cell-derived matrix images, as their surface is often not flat, resulting in uneven signal across the field of view. To account for this, contrast was enhanced in Fiji (0.35% saturated pixels) and subsequently a local contrast adjustment was performed (CLAHE plugin, default parameters). All images shown are maximum projections of Z-stack acquisitions, except for the live cell images which represent a single confocal slice.
2.4 Comparison With Other Available Alignment Tools
The performance of the algorithm presented in this work was compared with OrientationJ (Rezakhaniha et al., 2012) and TWOMBLI (Wershof et al., 2021), by using two 10-image samples of fibronectin-stained cell-derived matrices that were either isotropic or anisotropic. To access the alignment vector field for OrientationJ, a Fiji macro was designed to batch run the OrientationJ Vector Field plugin. This operation was either performed on its own or preceded by fiber segmentation by the Ridge Detection plugin (Lindeberg, 1998), this second approach mimicking the TWOMBLI workflow (hybrid method). To obtain a comparable spatial resolution of the vector field, a window of 100 px with 50% overlap was used to run the FFT analysis, and a matching local window σ of 100 px with a grid size of 50 was employed in OrientationJ. The Ridge Detection parameters were set to: line width = 20 px, high contrast = 120, low contrast = 0, sigma = 6.27, lower threshold = 0, upper threshold = 0.17, minimum line length = 10 px. These parameters were selected to match the ones assigned by the TWOMBLI macro when calling the Ridge Detection plugin (see below for the user input parameters that were chosen in TWOMBLI). From the acquired vector fields, an order parameter over neighborhoods of 5x vectors was calculated and compared for the three approaches.
A second comparison was based on coherency, the metric traditionally used for global alignment output by both OrientationJ and TWOMBLI. Coherency of the FFT was measured on windows of 100 px with 50% overlap and averaged to obtain a median value representing the global alignment score for each image. To obtain the coherency score from OrientationJ, a macro was designed to batch process the images using the OrientationJ Dominant Direction plugin. TWOMBLI was run on the sample data as per the developer instructions. A parameter optimization step was performed on test images and the obtained parameter file was used for batch processing (line width = 20 px, curvature window = 150 px, branch length = 10).
Computational time was measured on a 3.6 GHz Quad-Core Intel Core i7 machine with 32 GB memory, with suitable functions in MATLAB and ImageJ macro language. Only the effective analysis time for a 10-image sample (1 MB per image) was taken into account, excluding the input and parameter upload. Timings for OrientationJ refer to the designed macro, no manual handling of data was involved.
2.5 Statistical Analysis
Statistical analysis was performed in Prism (Graphpad, v9). Mann-Whitney two-tailed tests were used to compare metrics on isotropic vs. anisotropic samples. Significance is reported in figure panels as follows: ‘****’ for p-values lower than 0.0001, ‘***’ lower than 0.001, ‘**’ lower than 0.01, ‘*’ lower than 0.05, ‘ns’ otherwise.
This approach builds upon previous work employed in (Aratyn-Schaus et al., 2011; Cetera et al., 2014; Fernandes et al., 2020) and adds features such as periodic decomposition to improve the accuracy of the FFT angle determination along with the ability to search for the optimal length scale for comparison. The current implementation is available at https://github.com/OakesLab/AFT-Alignment_by_Fourier_Transform, both in MATLAB and Python languages. The MATLAB version is provided with a simple user interface for inputting parameters, while the Python version is presented as a set of Jupyter notebooks. The code is divided into two separate suites. The first suite can be used to perform the alignment quantification on a sample of images containing fibrillar structures. The second suite can be used to run a search to optimise the analysis parameters, with the aim of maximizing differences between two samples, i.e., locating the length scale for which the difference in alignment between two samples is greatest. Further code documentation is provided in the repository.
3.1 Measuring Alignment
The algorithm is inspired by the general approach of Particle Image Velocimetry (Willert and Gharib, 1991), where the image is broken down into a series of windows that are analyzed independently to create a vector field that represents the entire image. Windows are analyzed in frequency space to reveal information about both the fibrillar structure and local alignment as illustrated in Figure 1 using filamentous actin images. If the image contains aligned features in the real space, the corresponding FFT in the frequency domain will be asymmetrically skewed, with the direction of skew orthogonal to the original feature orientation (Figure 1A,B). This process can then be repeated on each successive window across the image, resulting in a vector field that represents the local alignment in the image. A detailed protocol for determining the local alignment of each window and the order parameter of global alignment follows.
FIGURE 1. Using FFTs to measure local alignment of biological images. (A) Small windows of a microscopic image are considered, containing either aligned filamentous actin fibers, or isotropic signal. The 2D Fast Fourier Transform (FFT) shows an elongated (i.e., skewed) shape if the image in the real space contains aligned fibrillar features, and displays a more round shape otherwise. Scale bar, 5 μm. (B) A small window of a microscopic image containing aligned filamentous actin fibers is considered. Its FFT shows a predominant skewness (black line) of an angle that is orthogonal to the direction of the fibers in the real space image. The obtained alignment vector (yellow line) can be overlaid to the original image to highlight the measured fiber orientation. Scale bar, 5 μm. (C) The alignment vector is calculated for an image containing aligned actin fibers by exploiting the characteristics of the FFT. The same image, rotated by 45°, is analyzed, to demonstrate robustness of the method to rotation. Scale bar, 5 μm. (D) A square representing the window size for analysis (250 px, light blue) is overlaid on a filamentous actin image. The black squares represent subsequent instances of such a window with the defined overlap, as the image is scanned during the analysis to calculate local alignment vectors (yellow). The heatmap shows a different representation of the alignment vector field, with a wrapped colour scale ranging from 0° to 180°. Scale bar, 20 μm.
3.1.1 Local Alignment
Square windows of n x n pixels are chosen from the original real space image and a 2D FFT is performed. Typically, n is chosen as odd to ensure that the zero-order frequency component, which is always the greatest in magnitude, is situated at the centre of the FFT. FFTs of non-periodic signals are often plagued by strong horizontal and vertical components, appearing as a cross, due to the mismatch of intensity in the images at the edges. To avoid this effect, we break down the image into its smooth and periodic components following the approach of Moisan (2011). We use the periodic component of this decomposition to take the FFT and then take its norm to only deal with real numbers. The resulting image is then multiplied by a mask of diameter n/2 to capture all the relevant high frequency components and to ensure the sample is symmetric. The central image moments of the masked FFT are then calculated as
where x and y represent the position and I is the magnitude of the FFT norm. We next construct the covariance matrix of the image as
which has eigenvalues of
The orientation of the FFT can thus be calculated by
and the eccentricity, a measure of how oblong the FFT is, as
The orientation of the features in the real image is orthogonal to the orientation of the FFT, and thus we apply a 90°rotation (Figure 1B). Orientation vectors do not have a polarity, and thus angles of 0° or 180° are considered equivalent. Therefore, all angles are mapped to a range of 0°–180°(Figure 1D). This methodology is robust to rotational transformations and does not suffer from the inherited bias of FFT for vertical and horizontal components (Figure 1C).
3.1.2 Order Parameter Calculation
Both the size of the window and the degree of overlap between the windows is customizable (Figure 1D). To measure how aligned fibrillar features are within adjacent windows, we define an order parameter by correlating the directionality of vectors within the neighborhood of customizable size. Specifically, the order parameter can be calculated as
with θij representing the angle between the orientation of a central reference vector and its neighbors (Figure 2A). The order parameter can take values between −1 and 1 (Figure 2A), with 1 representing perfect alignment (i.e., all vectors in the neighborhood have the same orientation as the reference vector), 0 representing random orientation, and −1 representing opposite alignment (i.e., all vectors in the neighborhood are pointing in the orthogonal direction compared to the reference vector). While a number of different formulations of the order parameter have been used previously, in the context of the current analysis, the chosen order parameter normally ranges between 0 (random) and 1 (perfectly aligned), as fiber polarity is not taken into account and the neighborhood area is kept relatively local (Figure 2B). Similar results were obtained for real and simulated data (and Figure 3).
FIGURE 2. Calculating an image order parameter. (A) The equation for the order parameter used to assign an alignment score for each neighborhood is shown, together with examples. The order parameter takes a value of 1 for complete alignment, 0 for random alignment, and −1 for orthogonal alignment. (B) A square representing a neighborhood of 5x vectors (magenta, neighborhood radius of 2x vectors around the central reference) is overlaid on an actin image and its corresponding alignment vector field. The black squares represent subsequent instances of such neighborhood, as the vector field is scanned during the analysis to obtain local alignment scores. The resulting order parameters for each neighborhood are reported in the heatmap, together with their median value, representing the output of the analysis.
FIGURE 3. Examples of image neighborhoods and corresponding order parameter. (A) Examples of neighborhoods in simulated data displaying different degrees of alignment are shown together with the calculated order parameter (window size 35 px). Scale bar, 25 px. (B) Examples of neighborhoods in actin images displaying different degrees of alignment are shown together with the calculated order parameter (window size 100 μm). Scale bar, 50 μm.
To calculate an overall alignment score for each image, we define a neighborhood size. The previously calculated vector field is then split into all possible overlapping neighborhoods of such size. The order parameter is subsequently calculated for each neighborhood, and the overall order parameter for the image is defined as the median value (Figure 2B).
3.2 Investigating the Optimal Alignment Length Scale Between Two Samples
When comparing two samples displaying different degrees of alignment, it may be of interest to evaluate the length scale for which this difference is greatest. The approach described above allows for this, as it is possible to evaluate the order parameter using a range of window and neighborhood sizes which correspond to different length scales (i.e., each permutation of window and neighborhood size define a specific length scale). By comparing the obtained alignment scores between the two experimental conditions for the different parameter permutations, it is possible to investigate the precise length scale at which the difference between the two samples is most pronounced. The comparison can be carried out by looking at the order parameter difference between the two samples or by running statistical tests and analyzing the resulting p-values.
We tested this parameter search approach on actin images of cultured fibroblasts with different degree of isotropy (Figure 4A, 10 images for each sample), for window sizes ranging from 25 to 325 px and neighborhood radii ranging from 1x to 38x vectors (Figures 4B,C). We performed the comparison in alignment scores for each window/neighborhood size permutation by quantifying the difference between the sample order parameters (Figure 4B) or by calculating the p-value of a Mann-Whitney statistical test between the two populations (Figure 4C).
FIGURE 4. Determining the optimal length scale difference. (A) Example actin images taken from two separate samples and displaying different degrees of alignment (higher on the left, more anisotropic; lower on the right, more isotropic). Scale bar, 20 μm. (B) The parameter search runs the analysis for a range of window sizes (y-axis, from 25 px to 325 px) and a range of neighborhood sizes (x-axis, from 3x to 77x vectors, where the neighborhood size is defined as 2*neighborhood radius + 1). Each permutation of window and neighborhood represents a length scale. Each colored square shows the difference in the median order parameter between the anisotropic sample and the isotropic sample (mean values over 10 images for each sample) for a specific pair of window and neighborhood sizes. Small window sizes paired with small neighborhoods (upper left corner) display noisy output, due to the limited amount of image information available to process. Three example vector fields are shown for a window size of 25 px and neighborhood of 55x vectors (corresponding to 92 μm), 75 px and 17x vectors (89 μm), and 200 px and 5x vectors (79 μm), with the light blue and magenta squares showing the size of the window and neighborhood respectively. Despite relatively similar neighborhood sizes, the output varies depending on the size of the examined window. Reading the graph horizontally from left to right for each window size up to ∼ 150 px, it is possible to note how an optimal length scale of ∼ 90 μm is detected (i.e., a peak difference value is displayed), likely correlated to relevant biological dimensions (e.g., cell size). (C) A similar graph to the one in (B) is shown for the p-value of Mann-Whitney tests between the two samples for all the combinations of window and neighborhood sizes. In this case, greater differences are represented by lower values (yellow).
It should be noted that small windows paired with small neighborhoods display noisy output (Figures 4B,C), as the corresponding FFT for such small windows will tend towards low eccentricity values (i.e., increasing the likelihood of spurious vectors, Figure 1A). Interestingly, when comparing two samples with known differences in alignment (Figure 4A), it can be observed how the difference in order parameter between the two samples displays a peak when keeping the window size constant and increasing the neighborhood (Figure 4B). This suggests a well-defined optimum length scale where a clear difference between the two samples occurs, and it is likely to be correlated with relevant dimensions in the sample (e.g., cell size). It is important to note that the optimal length scale will likely exist for a range of different window/neighborhood pairs, which together will define a similar local area. In the example shown in Figure 4B, the order parameter difference peaks for a length scale of ∼ 90 μm, as calculated by multiplying the window size by the neighborhood size for pairs displaying high difference values. This suggests that this length scale should be investigated further to evaluate the biological relevance of features of such size.
3.3 Evaluating the Length Scale Decay in Alignment
A similar approach to the parameter search can be employed to look at how the alignment decays over increasingly larger length scales, by progressively enlarging the size of the region over which the alignment is measured. To this aim, the window size over which the angle vector field is calculated is kept constant, while progressively increasing the neighborhood over which the order parameter is evaluated.
To confirm the performance of the algorithm in locating the alignment length scale, we created three 10-image sets of synthetic data. These images show a number of squared units of pre-defined size (i.e., 100, 200, or 400 px), randomly arranged in the field of view. Each unit is filled with parallel lines representing fibers (Figure 5A). We first measured the median order parameter for each image in each set with a window size of 35 px and a neighborhood of 5x vectors, corresponding to 100 px (Figure 5B). The order parameter increases with the selected unit size, as expected as neighborhoods of a given size will represent a progressively more local point of view for increasing unit sizes. Moreover, we analyzed the decay in the order parameter by keeping a fixed 10 px window size and progressively increasing the neighborhood size (Figure 5C). The half-life of the fitted decay curve shows values very close to the pre-defined unit size for each set. This result highlights the ability of our approach to correctly identify the alignment length scale within an image.
FIGURE 5. Evaluating length scale on synthetic data. (A) Synthetic images were created by randomly distributing units of parallel lines within the field of view. The size of each unit was pre-defined to be either 100 px (left, black), 200 px (centre, orange) or 400 px (right, green). Each sample contains 10 images. The blue and magenta boxes represent the window and neighborhood sizes used for the analysis in (B), set to 35 px and 100 px, respectively. (B) Measuring the order parameter using the same parameters for the three sets of images shows an increase in its median value with the unit size. This can be expected, as the information within this length scale will be more aligned the bigger the unit size. (C) Keeping the window size constant (10 px) and progressively increasing the neighborhood shows a decay in the order parameter for all three samples (dots). By fitting a one-phase decay to each data set (shaded lines) and measuring its half-life, it is possible to appreciate how the algorithm is able to locate the pre-defined length scale (i.e., unit size).
We then used the alignment decay analysis to quantify the length scale on real biological data, as exemplified by using ∼ 1 mm2 tilescan images of cell-derived matrices (i.e., matrices produced in situ by the cells) stained for F-actin and fibronectin (before and after decellularization, respectively, Figure 6A). We first calculated the angle vector field for both actin and fibronectin using the FFT approach for a window size of 250 px and overlap of 50% (Figure 6B). Subsequently, we calculated the median order parameter over neighborhoods of increasing radii, ranging from 1x to 21x vectors, corresponding to 50 x 50 μm to 1050 x 1050 μm in the original image (Figures 6A,C,D). By plotting the median order parameter over the neighborhood size, it is possible to observe the alignment decaying from the local to the global length scale (Figure 6D). As expected, actin and fibronectin displayed a similar decay in alignment as cells are responsible for depositing and remodelling the ECM in this experiment.
FIGURE 6. Measuring the decay in order parameter as a function of length scale. (A) Example tilescan images about 1 mm2 in size for actin and fibronectin on the same sample (fibroblast-derived matrix). Scale bar, 100 μm. (B) Angle heatmap showing the alignment vector field for the images in (A), using a window size of 250 px and 50% overlap. Actin and fibronectin display similar alignment. (C) Schematic of the alignment decay analysis, where the window size is maintained constant while increasing the size of the neighborhood (i.e., including more and more vectors in the order parameter calculation to evaluate alignment from the local to the global scale). The solid magenta square represents the starting neighborhood, which then increases in size at each iteration (dashed magenta squares). (D) The calculated order parameter for actin and fibronectin for the images in (A) is graphed for a window size of 250 px, for each neighborhood ranging from 50 μm (3x vectors, neighborhood radius of 1x vector) to 1050 μm (43x vectors, neighborhood radius of 21x vectors). Alignment decays in a similar manner when going from local (small neighborhood size) to global (large neighborhood size) length scales.
3.4 Filtering Input for Alignment Analysis
Some biological images might display regions that are either blank (i.e., little to no signal, lack of image features), or isotropic (i.e., no obvious fibrillar features). In certain applications, such regions should be excluded from the alignment analysis, as images containing noisy or blank areas will affect the output order parameter. By exploiting our window-based algorithm, it is possible to implement optional filtering of features to automatically exclude regions with such characteristics from the analysis.
To filter out blank regions, a threshold can be set on the mean intensity for each window, as demonstrated on second harmonic generation imaging of tissue samples (Figure 7A). In this example, windows that are considered background have a mean intensity signal lower than 20 (with intensity values ranging from 0-black to 255-white). This value can be set as a threshold: windows displaying mean intensity values lower than the threshold will be ignored during the calculation of the angle vector field (and subsequent alignment score, Figure 7B).
FIGURE 7. Filtering vector fields based on intensity and eccentricity. (A) A tissue section second harmonic generation (SHG) image is shown. Two regions with the same size of the window used for the alignment analysis (100 px) are enlarged and their mean intensity value measured. The intensity of an 8-bit image can vary from 0 to 255. The area on the top left contains little to no signal and it is considered background (mean intensity of 11); the area on the right contains many dark pixels, for a mean intensity of 84. A threshold value of 20 is set, meaning that all windows with mean intensity lower than such values will be considered as blank (i.e., background) and excluded from the analysis. Scale bar, 50 μm. (B) The alignment vector field (yellow) is overlaid on the raw SHG image, either without filtering for blank regions (left) or filtering with a threshold on the mean intensity of 20 (right). Magenta asterisks highlight regions excluded from the analysis. (C) An immunostaining for fibronectin in cell-derived matrices is shown. Two regions with the same size of the window used for the alignment analysis (100 px) are enlarged and their 2D Fast Fourier Transform (FFT) displayed. The FFT eccentricity can range from 0 (circular) to 1 (highly elongated). The region on the top lacks clear aligned fibers, resulting in a more uniform FFT (eccentricity of 0.40); the area on the bottom has aligned fibrillar features, leading to a highly skewed FFT (eccentricity of 0.84). A threshold value of 0.73 is set, meaning that all windows with eccentricity lower than such values will be considered as isotropic and excluded from the analysis. Scale bar, 20 μm. (D) The alignment vector field (yellow) is overlaid on the raw fibronectin image, either without filtering for isotropic regions (left) or filtering with a threshold on the eccentricity of 0.73 (right). Magenta asterisks highlight the regions excluded from the analysis. (E) The actin in a single cell is shown. If one wants to exclude the background from the analysis, a binary mask can be provided highlighting the regions of the image to be analyzed (window size 50 px with 50% overlap). Scale bar, 50 μm.
Regions in which image features display isotropic organization can also be excluded by setting a threshold on the eccentricity of the FFT calculated for a given window, as shown for cell derived matrices immunostained for fibronectin in Figure 7C. The FFT eccentricity can be used as a measure of its skewness (Figure 1A): regions containing oriented fibers will display a more elongated FFT, hence higher eccentricity (values closer to 1); regions where no fibers can be detected, or where there are multiple fibers without a clear orientation, will display a more homogeneously shaped (i.e., round) FFT, hence lower eccentricity (values closer to 0). In the example shown in Figure 7C, threshold on the eccentricity was set to 0.73, allowing for the exclusion of isotropic regions from the angle vector field calculation (Figure 7D).
Finally, it is possible that only a specific region within an image is to be analyzed, as in the case of actin in single cells (Figure 7E). A binary image can be used to mask the angle vector field, in order to only consider a given area of the original input.
3.5 Extracting Kinetics of Organization From Live Imaging
Many biological processes involve the evolution of architectures and organizations over time. By measuring an order parameter for each frame of a time series it is possible to extract the kinetics of alignment in response to perturbations. As an example, we transfected an NIH 3T3 fibroblast with GFP-actin and imaged the cell before, during, and after treatment with the Rho kinase inhibitor Y-27632 (Figure 8). Inhibition of Rho kinase leads to a reduction in myosin activity and thus a dissolution of stress fibers (e.g., a loss of fibrillar structures; Figure 8A). As a test case, treatment with Y-27632 is convenient because it acts rapidly and is easily washed out by replacing the media in the sample chamber with fresh media. Upon removal of Y-27632, the actin cytoskeleton begins to reform stress fibers and again takes on an aligned appearance (Figure 8B). As evidenced in the plot of the order parameter, we can see a loss of alignment immediately after addition of Y-27632 and an immediate recovery following washout (Figure 8C).
FIGURE 8. Measuring kinetics of alignment. (A) An NIH 3T3 fibroblast transfected with GFP-actin is shown before, during, and after washout of 20 μM Y-27632. (B) Local alignment was calculated for each frame in the time series using a window size of 17 px and an overlap of 50%. Yellow vectors indicate the local direction of alignment in the image while insets show the angle of orientation for each window. (C) Plot of the calculated order parameter over time for a neighborhood radius of 2x vectors. Arrows indicate the time points of addition of Y-27632 and wash out of the drug. Red dashed line indicates the fit to Eq. 7. (D) Boxplots of the average order parameter for fibroblasts treated for 30 min with different concentrations of Y-27632 and then fixed and stained with phalloidin to label actin. Each dot represents the median order parameter of a single cell with a window size of 33 px, overlap of 50%, and neighborhood radius of 2x.
To measure the kinetics of this interaction we can fit the data to any relevant model. For the washout of Y-27632 we fit the post-washout data to a simple exponential recovery curve of the form
We can then define a t1/2 value, or the time it takes to reach half its maximum recovery value, as
This results in a value of t1/2 ≈ 7 min for the Y-27632 washout experiment, consistent with previous reports (Aratyn-Schaus et al., 2011). To further illustrate the robustness of this approach, we imaged numerous cells treated with varying concentrations of Y-27632 (Figure 8D). Cells were treated with concentrations ranging from 0–20 μM for 30 min, and then fixed and stained with phalloidin. Analysis of these data sets reveals a steady decrease in order parameter as a function of Y-27632 concentration, as expected. The difference in magnitude of the order parameter between the fixed and live imaging is the product of superior signal to noise and better labelling that is achieved with phalloidin compared to fluorescently expressed actin.
3.6 Comparing With Other Available Algorithms
The performance of the algorithm presented in this work was tested against two other open-source tools widely used by the biology community, OrientationJ (Rezakhaniha et al., 2012) and TWOMBLI (Wershof et al., 2021). OrientationJ is a Fiji plugin based on image representation by structure tensors using real space, and evaluates orientation by the coherency metric, that indicates the degree of feature orientation with values ranging from 0 (random) to 1 (perfectly aligned) (Rezakhaniha et al., 2012; Clemons et al., 2018). Using this approach, it is possible to output either an alignment vector field (OrientationJ Vector Field) or a global coherency score for the image (OrientationJ Dominant Direction). TWOMBLI is a Fiji macro envisioned as a hybrid approach performing fiber segmentation followed by alignment analysis. Fiber segmentation is carried out with the Ridge Detection plugin (Lindeberg, 1998), while the alignment of the binarized image is subsequently measured via OrientationJ. In addition to the coherency metric as a measure of alignment and due to the segmentation step, TWOMBLI also allows the user to obtain further metrics that characterize the fibrillar features, such as number of branch and end points, density, curvature, and fractal dimensions.
We first compared the alignment vector field calculated with the three methods (FFT, structure tensor, and hybrid approach) for cell-derived matrices immunostained for fibronectin displaying different degrees of isotropy (Figure 9). When comparing the alignment vector field for images with strongly anisotropic fibrillar features, similar output was observed for the three methods (Figure 9A). When working with less aligned input, however, the FFT was able to better resolve the feature alignment in the isotropic regions compared to the other two algorithms (Figure 9B).
FIGURE 9. Comparing approaches to measure organization in biological images. (A) A fibronectin image showing fibrillar features with strong alignment is analyzed with the FFT approach presented in this work, with the Fiji OrientationJ Vector Field plugin, and with a combined approach using the Fiji Ridge Detection plugin followed by OrientationJ (mimicking the TWOMBLI hybrid approach, the binary mask calculated by the Ridge Detection is shown). A window size of 100 px with 50% overlap was used, and the order parameter was calculated over neighborhoods of 5x vectors (neighborhood radius of 2x vectors). Both the alignment vector field and the angle heatmap are shown for each case; zooms for the areas highlighted by the white boxes are displayed. The three algorithms behave similarly. Scale bar, 20 μm. (B) The analysis of a more randomly organized fibronectin image is compared as in (A). It can be noted how the FFT returns more accurate alignment vectors for the more isotropic regions of the image (insets).
Two alignment metrics were analyzed, the order parameter S (Eq. 6) and the coherency; the first metric is the one used in this work, while the coherency is the standard output of both OrientationJ and TWOMBLI. The comparison dataset contained two samples of images, with isotropic (Figure 9B) and anisotropic (Figure 9A) features. All the tested algorithms were able to recognize the alignment difference in the data, however, the order parameter scores calculated with the TWOMBLI approach (fiber segmentation via Ridge Detection followed by alignment with OrientationJ) did not reach statistical significance (Figure 10A). Moreover, the order parameter values calculated with both OrientationJ and the hybrid approach were very close to the upper limit of the range (values equal to 1 signify perfect alignment). This suggests that these methods might not be reliable in detecting subtle differences between samples, due to their limited ability to resolve local isotropy (Figure 9). Global coherency obtained with the three methods was also compared, showing that all of them could successfully distinguish the global alignment differences in the two samples (Figure 10B).
FIGURE 10. Comparing output metrics of different approaches to measure organization in biological images. (A) Boxplots of the order parameters for the samples shown in Figure 9 calculated with the three methods, the FFT (as presented here), OrientationJ (Vector Field plugin), and the hybrid approach (Ridge Detection followed by OrientationJ Vector Field). Local windows of 100 px with 50% overlap and 5x vector neighborhoods were considered. Insets represent zooms of the y-axis to better display distributions. Mann-Whitney two-tailed test (****p < 0.0001, **p = 0.0039, ns = 0.063). (B) Boxplots of the global coherency calculated with the three tools, AFT − Alignment by Fourier Transform (as presented here), OrientationJ (Dominant Direction plugin), and TWOMBLI. Mann-Whitney two-tailed test (***p = 0.0002, ***p = 0.0002, *p = 0.0147).
Computational time is comparable between the FFT and OrientationJ (with a dedicated macro) with each image being evaluated in less than half a second, while it increases for TWOMBLI with about 7 s per image (due to the fiber segmentation step). While this difference might appear insignificant for small samples as the ones used in this context, it might become important should batch processing be required for larger sample sizes or more high-throughput applications. The FFT approach presented here can be easily used to address length scale analysis. While OrientationJ also offers the option of a length scale analysis (Vector Field plugin), this requires ad hoc macros and further post-processing to obtain a unique alignment score for each image for a given window size. TWOMBLI allows the user to obtain a wider range of metrics characterizing the fibrillar features, thanks to its fiber segmentation step, but this comes at the cost of setting more parameters to start the analysis. The FFT performs best at resolving local alignment differences (Figure 9), and the available filters on window intensity and eccentricity make it straightforward to automatically exclude regions of the images. This avoids having to manually discard entire images containing unsuitable regions, as would happen with the other two methods. The pros and cons of each methodology should be taken into account depending on the application (Table 1); overall, the AFT approach demonstrated efficient performance and broad flexibility to different types of input images.
3.7 Implementing a Streamlined Software for Measuring Alignment
Finally, we have developed a workflow named AFT − Alignment by Fourier Transform to automate image analysis using the approach described above. This includes the pipeline described in Section 3.1, the parameter search in Section 3.2, the analysis of length scale decay in Section 3.3, the filtering options in Section 3.4, and the extraction of organization kinetics in Section 3.5. Open-source implementation is made available in both MATLAB and Python at https://github.com/OakesLab/AFT-Alignment_by_Fourier_Transform.
Some considerations regarding the choice of parameters follows. The user is requested to set three mandatory parameters in order to run the alignment analysis: the window size, the overlap, and the neighborhood size (Section 3.1, Figures 1D, 2B). Each of these parameters will affect the calculated order parameter by changing the length scale of the alignment analysis. While default parameters are suggested, it is important to test different values on individual images to gain some insight on the optimal length scale over which to carry out this analysis. In addition to the mandatory parameters, a number of optional filtering features are available as described in Section 3.4 (Figure 7).
3.7.1 Mandatory Parameters
• Window size. The size of the local regions in the image to analyse. This value depends on the image resolution, the size of the fibrillar features to be examined, and the length scale of interest (Figure 1D). In practice, within a window one should be able to observe enough information to discern by eye an orientation of features. There is no hard-set lowest limit for the window size, but it should be noted that smaller windows will lead to more spurious vectors as the calculated FFT will be nearly circular (i.e., not skewed, Figure 1A) for regions where fibrillar features are not recognized. A visual check of the output images is always recommended during the parameter optimization phase.
• Overlap. The overlap between adjacent windows (Figure 1D). Increasing the overlap, in conjunction with choosing smaller window sizes, allows for increasing the resolution of the analysis, with more sampled areas (i.e., vectors) within the field of view. A starting value of 50% is recommended.
• Neighborhood radius. The length scale over which to compare neighboring measures of alignment. This should be set depending on the length scale of interest. It is defined as the number of nearest neighbor vectors to be compared to a central reference in order to obtain an alignment score (Figure 2). Neighborhood size is calculated as 2*neighborhood radius + 1.
3.7.2 Optional Parameters
• Output images. This parameter enables the ability to save the analyzed images, consisting of the original image overlaid with the alignment vector field and the angle heatmap (Figure 1D). The color scale for the angle heatmap is wrapped, with similar colors for 0° and 180°, as fibers are not considered polar (Section 3.1, Figure 1D).
• Filter blank spaces. Windows that contain little or no signal can be excluded from the analysis. When enabled a threshold on the mean window intensity can be set (Figures 7A,B). Windows below this mean intensity are not analyzed. Threshold values can be estimated by opening a sample image in Fiji, selecting regions considered background with similar size to the window size, and measuring their mean intensity. Expected values range from 0 (black) to 255 (white) for an 8-bit image.
• Filter isotropic regions. It is possible to exclude windows from the analysis should the eccentricity of their FFT be lower than a threshold (Figures 7C,D). Values for the FFT eccentricity range from 0 (circular, isotropic) to 1 (elongated, anisotropic). To estimate this threshold, the analysis can be run iteratively for increasing values until the desired output is achieved (i.e., until the regions considered isotropic do not contain alignment vectors in the output images).
• Masking. By default, the whole image is being analyzed. Alternatively, the user can input a folder containing binary masks of selected areas to be analyzed for each input image (Figure 7E). A similar output as the one shown in Figure 7E could have also been obtained by filtering out the blank spaces, as the regions of interest have brighter signal than the background. However, if regions devoid of signal would also be present inside the cell outline, they would be disregarded as well. Therefore, a binary masking approach might be preferable in this case.
Quantifying alignment of fibrillar features in biological images has gained wide interest as a possible biomarker in disease aetiology and progression (Ouellette et al., 2021). While tools are available for this aim (Püspöki et al., 2016; Liu et al., 2017), many are difficult to implement and fairly rigid on the length scale over which anisotropy can be interrogated. The method presented in this work exploits the representation of the image in the frequency domain paired with a custom window-based approach, offering a computationally efficient algorithm that can be easily applied to investigate local to global alignment. Decay in alignment scores with increasing distance can be evaluated to reveal the length scale at which fibrillar features remain anisotropic.
An open-source suite, AFT − Alignment by Fourier Transform, is presented: its flexibility to a wide range of biological images is showcased through the examples, and its performance tested against synthetic data. The AFT − Alignment by Fourier Transform tool is demonstrated to be robust for different imaging modes (fixed, live-cell, SHG, fluorescence), fluorescent probes (F-actin, fibronectin), and image resolutions (from single cells to tissue tilescans). Further applications for which this method could be used that were not shown in the present work could entail the evaluation of alignment across the depth of a sample (measuring alignment for each slice of a Z-stack), or the investigation of cytoskeletal architectures in reconstituted protein systems.
The FFT approach performance was comparable to both OrientationJ and TWOMBLI, with increased accuracy in resolving local isotropic regions. Its computational efficiency makes it a good candidate for analyzing high volumes of data. Moreover, the availability of filters and masking options allow for this methodology to be applied to a wide range of microscopy images.
Excitingly, more advanced 3-dimensional imaging modalities are becoming available, with techniques being developed to achieve isotropic resolution. While currently the FFT approach is designed to work with 2D data, the algorithm could be adapted to evaluate alignment of 3D fibrillar features, and its low computational cost could be exploited for these more computationally intensive data.
Here we presented a methodology to measure fibrillar feature alignment in biological images. We created an open-source tool that can be used on a wide range of biological images. Its performance was compared against other common workflows, and it was shown to be computationally efficient without compromising on accuracy. It is easy to implement with a small number of analysis parameters, and allows for interrogating the length scale of fiber anisotropy.
Data Availability Statement
All the experimental datasets used as examples are included in the Supplementary Material (https://doi.org/10.6084/m9.figshare.15326472.v1), further inquiries can be directed to the corresponding authors. The software implementations can be found at https://github.com/OakesLab/AFT-Alignment_by_Fourier_Transform, together with further documentation on how to run the code.
Conceptualization and Methodology, PO; Software, SM, DF, PO; Data acquisition, LT, FK; Formal analysis and Data curation, SM, PO; Writing SM, DF, LT, FK, TS, BS, PO; Funding Acquisition and Supervision, TS, BS, PO.
This project has been funded from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement no. 681808 – SM, BS), the Wellcome Trust (grant no. 107859/Z/15/Z – FK, BS), the London Interdisciplinary Doctoral Programme (BB/J014567/1 – DBDF), the National Institutes of Health (NIH) National Institute of Allergy and Infectious Disease (NIAID) (Award # P01 AI02851 – LDT, PWO), and the National Science Foundation (NSF) (CAREER award # 2000554 – PWO). For the purpose of open access, the author has applied a CC BY public copyright licence to any Author Accepted Manuscript version arising from this submission.
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
The Supplementary Material for this article can be found online at: https://doi.org/10.6084/m9.figshare.15326472.v1
Boudaoud, A., Burian, A., Borowska-Wykręt, D., Uyttewaal, M., Wrzalik, R., Kwiatkowska, D., et al. (2014). FibrilTool, an ImageJ Plug-In to Quantify Fibrillar Structures in Raw Microscopy Images. Nat. Protoc. 9, 457–463. doi:10.1038/nprot.2014.024
Bredfeldt, J. S., Liu, Y., Conklin, M. W., Keely, P. J., Mackie, T. R., and Eliceiri, K. W. (2014a). Automated Quantification of Aligned Collagen for Human Breast Carcinoma Prognosis. J. Pathol. Inform. 5, 28. doi:10.4103/2153-3539.139707
Bredfeldt, J. S., Liu, Y., Pehlke, C. A., Conklin, M. W., Szulczewski, J. M., Inman, D. R., et al. (2014b). Computational Segmentation of Collagen Fibers from Second-Harmonic Generation Images of Breast Cancer. J. Biomed. Opt. 19, 016007. doi:10.1117/1.jbo.19.1.016007
Cetera, M., Ramirez-San Juan, G. R., Oakes, P. W., Lewellyn, L., Fairchild, M. J., Tanentzapf, G., et al. (2014). Epithelial Rotation Promotes the Global Alignment of Contractile Actin Bundles during Drosophila Egg Chamber Elongation. Nat. Commun. 5, 1. doi:10.1038/ncomms6511
Chaudhuri, S., Nguyen, H., Rangayyan, R. M., Walsh, S., and Frank, C. B. (1987). A Fourier Domain Directional Filterng Method for Analysis of Collagen Alignment in Ligaments. IEEE Trans. Biomed. Eng. BME-34, 509–518. doi:10.1109/TBME.1987.325980
Clemons, T. D., Bradshaw, M., Toshniwal, P., Chaudhari, N., Stevenson, A. W., Lynch, J., et al. (2018). Coherency Image Analysis to Quantify Collagen Architecture: Implications in Scar Assessment. RSC Adv. 8, 9661–9669. doi:10.1039/c7ra12693j
Eltzner, B., Wollnik, C., Gottschlich, C., Huckemann, S., and Rehfeldt, F. (2015). The Filament Sensor for Near Real-Time Detection of Cytoskeletal Fiber Structures. PLoS ONE 10, e0126346–28. doi:10.1371/journal.pone.0126346
Falzone, T. T., Oakes, P. W., Sees, J., Kovar, D. R., and Gardel, M. L. (2013). Actin Assembly Factors Regulate the Gelation Kinetics and Architecture of F-Actin Networks. Biophysical J. 104, 1709–1719. doi:10.1016/j.bpj.2013.01.017
Fernandes, N. R. J., Reilly, N. S., Schrock, D. C., Hocking, D. C., Oakes, P. W., and Fowell, D. J. (2020). CD4+ T Cell Interstitial Migration Controlled by Fibronectin in the Inflamed Skin. Front. Immunol. 11, 1501–1514. doi:10.3389/fimmu.2020.01501
Goldyn, A. M., Kaiser, P., Spatz, J. P., Ballestrem, C., and Kemkemer, R. (2010). The Kinetics of Force-Induced Cell Reorganization Depend on Microtubules and Actin. Cytoskeleton 67, NA. doi:10.1002/cm.20439
Gupta, M., Doss, B. L., Kocgozlu, L., Pan, M., Mège, R.-M., Callan-Jones, A., et al. (2019). Cell Shape and Substrate Stiffness Drive Actin-Based Cell Polarity. Phys. Rev. E 99, 1–16. doi:10.1103/PhysRevE.99.012412
Gupta, M., Sarangi, B. R., Deschamps, J., Nematbakhsh, Y., Callan-Jones, A., Margadant, F., et al. (2015). Adaptive Rheology and Ordering of Cell Cytoskeleton Govern Matrix Rigidity Sensing. Nat. Commun. 6, 1. doi:10.1038/ncomms8525
Kartasalo, K., Pölönen, R. P., Ojala, M., Rasku, J., Lekkala, J., Aalto-Setälä, K., et al. (2015). CytoSpectre: A Tool for Spectral Analysis of Oriented Structures on Cellular and Subcellular Levels. BMC Bioinformatics 16, 1–23. doi:10.1186/s12859-015-0782-y
Linsmeier, I., Banerjee, S., Oakes, P. W., Jung, W., Kim, T., and Murrell, M. P. (2016). Disordered Actomyosin Networks Are Sufficient to Produce Cooperative and Telescopic Contractility. Nat. Commun. 7, 12615. doi:10.1038/ncomms12615
Liu, Y., Keikhosravi, A., Pehlke, C. A., Bredfeldt, J. A., Dutson, M., Liu, H., et al. (2020). Fibrillar Collagen Quantification with Curvelet Transform Based Computational Methods. Front. Bioeng. Biotechnol. 8, 1–14. doi:10.3389/fbioe.2020.00198
Mascharak, S., des Jardins-Park, H. E., Davitt, M. F., Griffin, M., Borrelli, M. R., Moore, A. L., et al. (2021). Preventing Engrailed-1 Activation in Fibroblasts Yields Wound Regeneration without Scarring. Science 372, 1. doi:10.1126/science.aba2374
Ouellette, J. N., Drifka, C. R., Pointer, K. B., Liu, Y., Lieberthal, T. J., Kao, W. J., et al. (2021). Navigating the Collagen Jungle: The Biomedical Potential of Fiber Organization in Cancer. Bioengineering 8, 1–19. doi:10.3390/bioengineering8020017
Park, D., Wershof, E., Boeing, S., Labernadie, A., Jenkins, R. P., George, S., et al. (2020). Extracellular Matrix Anisotropy Is Determined by TFAP2C-dependent Regulation of Cell Collisions. Nat. Mater. 19, 227–238. doi:10.1038/s41563-019-0504-3
Püspöki, Z., Storath, M., Sage, D., and Unser, M. (2016). “Transforms and Operators for Directional Bioimage Analysis: A Survey,” in Focus On Bio-Image Informatics. Advances in Anatomy Embryology and Cell Biology. Editors W. De Vos, S. Munck, and J. Timmermans (Cham: Springer), 69–93. doi:10.1007/978-3-319-28549-8_3
Rezakhaniha, R., Agianniotis, A., Schrauwen, J. T., Griffa, A., Sage, D., Bouten, C. V. C., et al. (2012). Experimental Investigation of Collagen Waviness and Orientation in the Arterial Adventitia Using Confocal Laser Scanning Microscopy. Biomech. Model. Mechanobiology 11, 461–473. doi:10.1007/s10237-011-0325-z
Ridley, A. J., and Hall, A. (1992). The Small Gtp-Binding Protein Rho Regulates the Assembly of Focal Adhesions and Actin Stress Fibers in Response to Growth Factors. Cell 70, 389–399. doi:10.1016/0092-8674(92)90163-7
Standley, P. R., Camaratta, A., Nolan, B. P., Purgason, C. T., and Stanley, M. A. (2002). Cyclic Stretch Induces Vascular Smooth Muscle Cell Alignment via NO Signaling. Am. J. Physiol. – Heart Circulatory Physiol. 283, 1907–1914. doi:10.1152/ajpheart.01043.2001
Wershof, E., Park, D., Barry, D. J., Jenkins, R. P., Rullan, A., Wilkins, A., et al. (2021). A FIJI Macro for Quantifying Pattern in Extracellular Matrix. Life Sci. Alliance 4, e202000880. doi:10.26508/lsa.202000880