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.

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.

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.

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.

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.

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).

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.

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.

We then used the alignment decay analysis to quantify the length scale on real biological data, as exemplified by using ∼ 1 mm ² 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.

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).

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.

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).

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 S ( t ) = A − B ⁡ exp ( − t / C ) . (7) We can then define a t _1/2 value, or the time it takes to reach half its maximum recovery value, as t 1 / 2 = − C ⁡ ln A 2 B . (8) This results in a value of t _1/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.

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).

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).

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.

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).