2.1. SPM statistical framework

SPMs are used to evaluate the probability of change in every pixel by using decision tests based on the magnitude of the SPM values (i.e., the peak intensity of a cluster in SPM) as well as the spatial extent of these clusters formed at certain intensity thresholds⁽Reference Poline, Worsley, Evans and Friston ^{1 )}. The value of the pixel intensity reflects a parametric value of interest. Hence, these two-dimensional (2D) maps are constructed to reflect the spatial variation of a measured parameter which is important for discriminating between two given classes with regards to its intensity and area. In this way, clusters of high intensity of an SPM can correspond to a high localized parametric variation, while a large region reflects a spatially extended area of variation.

Our proposed spatial statistical learning framework relies on graph-derived parametric maps to quantify and simultaneously localize statistically significant differences across normal and disease-mimicking FN organizations. Using the pixel intensity of the maps along with the extent of the region area, these differences can be assessed both quantitatively and qualitatively, and detected as anomalies with respect to a Gaussian random field (GRF), corresponding to regions within the maps that cannot be explained by the GRF model learnt from the reference population. Hereafter, we describe the theoretical framework of GRF⁽Reference Adler ^{13 )} that enables the statistical analysis of tissue parametric maps.

GRFs, whose marginal distributions are Gaussian vectors $ X=\left({X}_{(1)},\cdots, {X}_{(n)}\right) $ , are characterized by the probability density function:

(1) $$ {f}_X(x)={\left(2\pi \right)}^{-n/2}{\left|V\right|}^{-\frac{1}{2}}\mathit{\exp}\left[-\frac{1}{2}\left(x-\mu \right){V}^{-1}{\left(x-\mu \right)}^T\right] $$

where $ \mu ={\left(E\left({X}_{(i)}\right)\right)}_{i\in \left\{1,\dots, n\right\}} $ is the expectation and $ V={\left(E\left[\left({X}_{(i)}-{\mu}_{(i)}\right)\left({X}_{(j)}-{\mu}_{(j)}\right)\right]\right)}_{i\in \left\{1,\dots, n\right\},\hskip0.35em j\in \left\{1,\dots, n\right\}} $ is the covariance matrix.

We consider clusters of pixels as connected components that are formed based on 8-pixel connectivity. Hence, upon image binarization according to a chosen threshold, the pixels (with intensity equal to 1) are grouped together in disjoint components (including single-pixel components) based on similar values of the neighboring eight pixels. It was shown in⁽Reference Adler ^{13 )}, and later adopted in fMRI-specific studies⁽Reference Poline, Worsley, Evans and Friston ^{1 ,} Reference Friston, Worsley, Frackowiak, Mazziotta and Evans ^{2 ,} Reference Worsley, Evans, Marrett and Neelin ^{14 )}, that for large thresholds $ t $ , the clusters are independent and the expectation of the number of clusters at a threshold $ t $ , of an image modeled by a zero-mean, homogeneous Gaussian field of dimension 2, is estimated by the expected Euler characteristic of the excursion set of the GRF:

(2) $$ E\left[{\mathrm{m}}_{\mathrm{t}}\right]=S{\left(2\pi \right)}^{-3/2}{\left|\Lambda \right|}^{1/2}t{\sigma}^{-3}\exp -\frac{t^2}{2{\sigma}^2} $$

where $ {\mathrm{m}}_{\mathrm{t}} $ represents the number of clusters at a certain threshold  $ t $ , S is the number of pixels of the image, $ \Lambda $ is the covariance matrix of partial derivatives of the GRF, and $ \sigma $ is the standard deviation of the GRF.

Similarly, the mean value of the number of clusters at a threshold $ t+{H}_0 $ can be written as follows:

(3) $$ E\left[{\mathrm{m}}_{\mathrm{t}+{H}_0}\right]=S{\left(2\pi \right)}^{-3/2}{\left|\Lambda \right|}^{1/2}\left(t+{H}_0\right){\sigma}^{-3}\exp -\frac{{\left(t+{H}_0\right)}^2}{2{\sigma}^2} $$

Considering $ {x}_0=t+{H}_0 $ as the intensity peak of a cluster (at threshold $ t $ ), one can estimate the probability that a cluster (at a threshold $ t $ , having an intensity peak higher or equal to $ {x}_0 $ , denoted $ {C}_t^{H_0} $ ) belongs to a realization of this GRF, $ {G}_r $ . This probability, as shown in⁽Reference Poline, Worsley, Evans and Friston ^{1 )}, termed $ {P}_H $ , can be seen as the likelihood of a cluster (formed at threshold $ t $ ) of having an intensity peak higher or equal to $ t+{H}_0 $ ⁽Reference Poline, Worsley, Evans and Friston ^{1 )}:

(4) $$ P\left({C}_t^{H_0}\in {G}_r\right)=\frac{E\left[{\mathrm{m}}_{\mathrm{t}+{H}_0}\right]}{E\left[{\mathrm{m}}_{\mathrm{t}}\right]}=\frac{x_0}{t}\exp \frac{t^2-{x}_0^2}{2{\sigma}^2} $$

Next, we were interested in the estimation of the probability that a cluster (at a threshold $ t $ ) belongs to a realization of GRF, depending on its surface (spatial extent - number of pixels). To estimate the number of pixels ( $ {n}_t $ ) of a cluster at a threshold $ t $ , we use the following equation from⁽Reference Friston, Worsley, Frackowiak, Mazziotta and Evans ^{2 )}:

(5) $$ E\left[{n}_t\right]=\frac{E\left[{N}_t\right]}{E\left[{m}_t\right]} $$

where $ {N}_t $ is the number of pixels at of higher intensity than $ t $ , and $ {m}_t $ is the number of clusters at the threshold $ t $ . Since the intensity values follow a normal (zero mean value) distribution, the expectation of $ {N}_t $ is the following⁽Reference Friston, Worsley, Frackowiak, Mazziotta and Evans ^{2 )}:

(6) $$ E\left[{\mathrm{N}}_{\mathrm{t}}\right]=S{\int}_t^{\infty }{\left(2{\pi \sigma}^2\right)}^{-1/2}\exp -\frac{x^2}{2{\sigma}^2}\mathrm{dx}=S{\Phi}_{\sigma}\left(-t\right) $$

where $ {\Phi}_{\sigma}\left(-t\right) $ is the complementary cumulative distribution function. It follows, then, based on Equations (2, 5, 6) that one can approximate the mean value of $ {n}_t $ , accordingly:

(7) $$ E\left[{n}_t\right]=\frac{E\left[{\mathrm{N}}_{\mathrm{t}}\right]}{E\left[{m}_t\right]}=\frac{\Phi_{\sigma}\left(-t\right)}{{\left(2\pi \right)}^{-3/2}{\left|\Lambda \right|}^{1/2}t{\sigma}^{-3}\exp -\frac{t^2}{2{\sigma}^2}} $$

Furthermore, $ {n}_t $ follows an exponential distribution law⁽Reference Lafarge, Descombes and Zerubia ^{15 )}, which is commonly defined by a parameter $ {\lambda}_t $ , (the inverse of the mean expected value of the random variable). Consequently:

(8) $$ P\left({n}_t=x\right)={\lambda}_t\mathit{\exp}\left(-{\lambda}_tx\right) $$

where $ {\lambda}_t=\frac{{\left(2\pi \right)}^{-3/2}{\left|\Lambda \right|}^{1/2}t\hskip0.35em {\sigma}^{-3}\exp -\frac{t^2}{2{\sigma}^2}}{\Phi_{\sigma}\left(-t\right)} $

It follows then that the approximation for the probability $ {P}_S $ of a given cluster $ {C}_t^{S_0} $ having a spatial extent $ S $ greater than $ {S}_0 $ , is given by the following formulation, as shown in⁽Reference Friston, Worsley, Frackowiak, Mazziotta and Evans ^{2 )}:

(9) $$ P\left({C}_t^{S_0}\in {G}_r\right)=P\left({n}_t\ge {S}_0\right)=\mathit{\exp}\left(-{\lambda}_t{S}_0\right)=\mathit{\exp}\frac{{\left(2\pi \right)}^{-3/2}{\left|\Lambda \right|}^{1/2}{S}_0t{\sigma}^{-3}\exp -\frac{t^2}{2{\sigma}^2}}{\Phi_{\sigma}\left(-t\right)} $$

In the following section, we illustrate our proposed approach using a simple scenario of a simulated GRF, in which the model parameters are estimated from a “normal” sample, to detect significant changes in an “abnormal” example.

2.2. Test for anomaly detection using synthetic data

To better understand the concept of detecting statistical abnormalities in a GRF realization, we considered a simple scenario which includes a synthetic example (normal), representing a simulation of a GRF of zero mean (Figure 1a). To this example we added 6 foreign objects representing ellipses with different surfaces, and intensity levels. The aim was to use our approach to detect these 6 objects within the abnormal sample, showcasing the potential to localize these anomalies at various thresholds, based on the maximum intensity of the different regions detected at each threshold, or their spatial extent (Figure 1b). The null hypothesis is that clusters of pixels computed at different thresholds in the abnormal example belong to a realization of the same GRF as the normal sample. The hypothesis is rejected if either $ {P}_H $ or $ {P}_S $ is less or equal than a p-value (pval) of 0.05.

Figure 1. (a) Realization of a GRF of zero-mean (normal case) and addition of 6 different sized ellipses (abnormal case), along with intensity-based hard thresholding for each case, respectively (at a threshold $ t $ = 10 pixels) (b). Intensity and surface-based detection at pval $ \le 0.05 $ , and at thresholds equal to (10,15,20, see colorbar). GRF sample detection will typically contain false-positive detections while within the abnormal case (GRF+ellipses), all “foreign” objects, that is, ellipses are detected at pval $ \le 0.05 $ , along with a few false-positives. (c). Clusters that are jointly detected based on the surface and intensity criteria are selected for both normal and abnormal samples.

Our method learns the GRF model parameters from the reference example, and then uses these parameters to compute the two probabilities of belonging to GRF for each region at various thresholds in the abnormal example. As shown in Figure 1a, b, the current method, compared to a naïve hard thresholding (threshold equal to 10), is better suited to localize the abnormal elements, that is, ellipses, at thresholds equal to (10,15,20). One could opt to jointly consider intensity and surface-based criteria when selecting the detected objects, which in this scenario, would lead to selecting the 6 ellipses together with one false positive (Figure 1c). It is noteworthy that the same false positive is detected on the reference image, which is consistent with the definition of the pval.

2.3. Machine-learning embedded in a GRF-based statistical parametric map framework

Providing both a quantitative and a qualitative assessment of the parameter variations is imperative for the study of spatial heterogeneity of fibers in pathological conditions. We were interested in leveraging our proposed framework for the comparison of two distinct conditions, normal and pathological ECM, given parametric maps that reflect various fiber characteristics. To do so, we learnt the GRF model’s parameters from the normal samples and using these parameters, we subsequently determined the probabilities of regions within the pathological conditions belonging to the same GRF, casting the original framework into a machine-learning setting. More specifically, we applied our proposed framework to determine whether fiber length and pore directionality can discriminate between normal and disease-mimicking states. While topographical differences between ECM of healthy and tumor tissue have been described⁽Reference Conklin, Eickhoff and Riching ^{6 )}, mainly for collagen, no current computational study can, to our knowledge, simultaneously localize and quantify them.

First, we describe the principle behind the proposed approach (Figure 2), which relies on modeling the normal maps as realizations of a GRF and testing this hypothesis on tumor-like maps⁽Reference Friston, Worsley, Frackowiak, Mazziotta and Evans ^{2 )}. We hypothesized that the tumor-like maps are realizations of the GRF learnt from the reference maps and determined a set of probabilities that characterize a degree of belonging to the GRF, for certain contiguous regions (clusters) at various intensity thresholds. In other words, the current statistical analysis identifies those foreign regions with respect to the reference GRF, within both types of maps, under the null hypothesis (i.e., at a given pval)).

Figure 2. Methodology for statistical detection of foreign regions to a GRF, in an example of a sample representing normal and tumor-like parametric maps. (a) Normal and tumor-like fiber length maps. The normal sample is modeled as a realization of a GRF, and we assume that the tumor-like sample is a realization of the same process. Clusters of regions with an intensity higher than a given threshold, t = 50 (b), t = 80 (c), $ t $ = 100 (d) are found to be statistically different to the GRF, with respect to a pval, depending on the cluster maximum intensity value or their surface.

In this context, the parametric maps are described by the union of two classes of pixels: those representing a realization of a GRF modeling the normal case, and those that are foreign to the GRF. We expect these foreign elements to occur in regions with very high pixel intensity and/or in larger clusters taken at a specific threshold.

Modeling the parametric maps with GRF is only possible upon gaussianization (i.e., conversion of the empirical distributions into normal distributions) of the GRF marginal distributions. In practice, we only performed the gaussianization of the first-order marginal distribution of the GRF, that is, the image intensity histogram, considering that the parameter maps under study were smooth enough. Therefore, the image intensity histogram was the only distribution to be gaussianized using an approach based on optimal transport⁽Reference Peyré and Cuturi ^{16 ,} Reference Peyré ^{17 )} (Supplementary Figure S1). Thereby, the resulting intensity histogram follows a normal distribution of zero mean with identical variance as the empirical native histogram.

To estimate the likelihood of a certain cluster formed at an intensity threshold $ t $ to belong to a GRF, depending on the maximal intensity of this cluster, we relied on the formulations taken from the theory of random fields⁽Reference Adler ^{13 )}, as previously described. Considering $ {x}_0=t+{H}_0 $ as the intensity peak of a cluster (at intensity threshold $ t $ ), one can estimate the probability that a cluster having an intensity peak higher or equal to $ {x}_0 $ , belongs to a realization of GRF. This probability can be seen as the likelihood of a cluster (taken at threshold $ t $ ) of having an intensity peak higher or equal to $ t+{H}_0 $ (Equation 4). Furthermore, as previously shown, the approximation for the probability of a given cluster having a spatial extent S greater than $ {S}_0 $ is given by Equation 9. At pval $ \le 0.05 $ , the clusters of pixels identified at $ t $ are considered significantly different from the normal GRF model.

In our experiments, we focused on two different FN parametric maps that could potentially discriminate between normal and pathological conditions, fiber length and pore directionality maps, in both reference and disease-mimicking states. We embedded the SPM framework, initially developed to analyze single datasets independently, into a machine learning paradigm. To evaluate the differences between two given groups of maps (e.g., FN in normal state vs disease-like state), we considered one of the groups as the normal realization of the GRF which we divided into a learning set and a smaller test set. The second group was tested for anomalous regions, therefore all the images belonging to this group were considered part of the test set. The proposed method learns the normal GRF model specific parameters, that is, the average value of $ {\lambda}_j,{\sigma}_j $ , from the training set. These learnt parameters during the learning phase were subsequently used to compute the two relevant probabilities, $ {P}_H $ and $ {P}_S $ , as described hereafter:

For all the images $ {\mathrm{I}}_j $ (previously Gaussianized) in the learning set:

• Computation of $ {\Lambda}_j $ (Equation 2, empirical estimator of the covariance of partial derivatives of $ {I}_j $ ). If for an image function $ f\in {\mathrm{\mathbb{R}}}^2 $ , we consider its gradient vector $ \mathrm{\nabla}f=(\hskip2.5pt {f}_x,{f}_y)=(\frac{\mathrm{\partial}f}{\mathrm{\partial}x},\frac{\mathrm{\partial}f}{\mathrm{\partial}y}) $ , then $ {\Lambda}_j=\mathrm{cov}[\hskip2.5pt {f}_x,{f}_y]=E[(\hskip2.5pt {f}_x-E[\hskip2.5pt {f}_x])(\hskip2.5pt {f}_y-E[\hskip2.5pt {f}_y])] $ .

• Computation of $ {\sigma}_j $ , as the $ {I}_j $ ’s sample standard deviation.

The last step involves storing the average $ {\Lambda}_m,{\sigma}_m $ of the learning dataset.

During the test phase, for all the clusters identified at a threshold $ t $ within the test set, $ {P}_H $ and $ {P}_S $ are evaluated using the model’s previously learnt parameters. At pval $ \le 0.05 $ , the clusters are significantly different from the normal GRF model, and can be considered for subsequent analysis (e.g., quantification).

For all the images $ {\mathrm{I}}_j $ in the test set:

• Gaussianization of each sample image $ {I}_j $ . The result is a new image $ {I}_g $ , whose histogram is Gaussian with identical variance to that of $ {I}_j $ .

• For a given list of thresholds $ T=\left({t}_1,{t}_2,\cdots, {t}_n\right): $

– Binarization of the image $ {\mathrm{I}}_g $ according to the threshold $ {t}_i $

– Once the list of connected components in the binary image resulted from thresholding is achieved, then for every (labeled) connected-component $ \left({l}_1,{l}_2,\cdots, {l}_p\right): $

* Evaluation of $ {P}_H $ using the learnt model parameter $ {\sigma}_m $ .

* Evaluation of $ {P}_S $ using the learnt model parameters $ {\Lambda}_m,{\sigma}_m $ .

4.1. Fiber detection using Gabor filters and graph extraction

To detect and quantify fiber-specific properties of FN networks (Figure 4a), we developed a pipeline that was primarily utilized for the extraction of local topological fiber properties from 2D confocal microscopy images, as described previously⁽Reference Efthymiou, Radwanska and Grapa ^{11 )}. Existing ECM analysis methods focusing on collagen exploit alternative fiber detection techniques, such as Fast Fourier Transform bandpass filters⁽Reference Morrill, Tulepbergenov, Stender, Lamichhane, Brown and Lujan ^{18 )}, ridge detection⁽Reference Wershof, Park and Barry ^{19 )}, or fast discrete curvelet transform⁽Reference Liu, Keikhosravi and Pehlke ^{20 )}, which is arguably the best suited method to detect curvilinear anisotropic objects, among the previous options. The latter option was used in conjunction with a fiber extraction algorithm. Our method relies on a more flexible detection scheme using Gabor filters, thereby avoiding translation/rotation errors, and unlike other methods, associates graph networks to fiber morphological skeletons, enabling diverse fiber analyses. The current study builds on our previously described fiber enhancement approach⁽Reference Efthymiou, Radwanska and Grapa ^{11 )}, for which the key steps are summarized as follows. Fibers in confocal images (Figure 4a) were detected and enhanced using Gabor filters (Supplementary S1) tailored to capture a range of different fiber elements that occur at multiple frequencies and orientations (Figure 4b). Subsequently, we opted for a graph-based framework to construct morphological fiber skeletons (Figure 4c) that would ultimately provide a geometrical characterization of fiber patterns. Further steps for improving fiber representation (e.g., fiber pruning, post-processing fiber reconnection) were implemented as described in previous work⁽Reference Grapa ^{21 )}. Graphs (i.e., collections of nodes connected by edges) are powerful tools for the structural and pattern analysis of objects, which can be utilized for the mathematical study of relations between entities, including fiber-like object detection⁽Reference Kollmannsberger, Kerschnitzki, Repp, Wagermaier, Weinkamer and Fratzl ^{22 ,} Reference Aguilar, Martinez-Perez, Frauel, Escolano, Lozano and Espinosa-Romero ^{23 )}.

Figure 4. Fiber enhancement and graph-based representation starting from confocal 2D images. (a) Representative region (512x512 pixels) of a sample image (FN B-A+) at a resolution of 0.27 $ \mu $ m/pixel. (b) Fiber enhancement with Gabor filters (c) Morphological fiber skeleton extraction (d) Skeleton-based graph (left) and simplified graph representation (right) which is derived from the latter.

Within our current analysis, we employ two different graph types to measure fiber-specific properties. First, graphs are used to depict a morphological skeleton representation (Figure 4d, left). Here, the nodes represent either fiber crosslinks (actual fiber junctions or junctions due to the 2D projection of the network onto the image plane) or fiber ends, and the edges capture the fiber length between two given nodes. We previously showed how such a representation can provide a description of distinct local features among four FN variant networks, in a normal state⁽Reference Efthymiou, Radwanska and Grapa ^{11 )}.

The second type of graph-based representation introduced in this work (Figure 4d, right) is meant to simplify fiber delineation, as described hereafter. Starting from the skeleton graph, we kept all nodes corresponding to fiber extremities, and connected all pairs of nodes with a straight line, if a fiber had previously been identified. For the sake of simplicity, we refer to fiber length as the length of any straight line connecting a pair of graph nodes.

We note that both representations are useful to extract different local or global fiber properties. The graph-based skeleton fiber delineation faithfully represents (according to a visual assessment performed by a trained biologist) the geometrical and topological properties of the fibers from the 2D confocal images, while the Gabor-specific (e.g., fiber local orientation, thickness) and graph-derived parameters (e.g., fiber length, number of nodes, etc.) are linked to meaningful physical fiber attributes. This biologically relevant representation enabled us to develop here a statistical analysis of the variation of certain fiber parameters for both the normal and a tumor-like state of the FN networks.

4.2. Generation of fiber parametric maps from graph-derived fiber representations

We next sought to develop a statistical framework for differentiating between parametric maps of activated and non-activated FN network configurations, computed from graph-based fiber representations. To create tissue variation maps (Figure 5), we considered different fiber attributes computed from graphs, representing either morphological skeletons (Figure 5a) or simplified graph depictions (Figure 5c). This framework is exemplified on two different types of parametric maps (Figure 5b,d) reflecting discriminative features, namely the individual fiber lengths (i.e., the length of the connecting line between two graph nodes), and the fiber pore (“gap”) directionality (i.e., the inverse value of the absolute difference between the median and individual pore orientation).

Figure 5. Computation of fiber parametric maps: (a) Starting from the skeleton graph (FN B-A+ disease-like sample, 1024x1024 pixels, 0.27 $ \mu $ m/pixel), a pore directionality map is derived (b), as the inverse value of the difference between the median pore angle and each individual one. (c) Starting from the fiber skeleton associated graph, a parametric map (fiber length, (d)) associates the fiber length, in pixels, to each corresponding connecting line.

To build a fiber length map (Figure 5d), the starting point was the simplified graph-based representation, where nodes depict fiber ends or crosslinks, and fibers are represented by the straight connecting line between these nodes. During the next step, we identified the 2D pixels coordinates that approximate the straight line between the nodes⁽Reference Bresenham ^{24 )} and assigned the length value of the connecting line to each one of the corresponding pixels. The last step for generating dense fiber length maps consists of including the extrapolation of the fiber length values⁽Reference D’Errico ^{25 )} and smoothing of this result with a Gaussian kernel. Concerning the pore directionality parametric maps (Figure 5b), the starting point was the skeleton graph. Pore orientation was computed by first fitting ellipses to each pore and was subsequently obtained by measuring the angle i between the horizontal axis and the major ellipse axis. Pore median value was then subtracted from each i to remove image rotations from the analysis. Finally, the inverse absolute value of the resulting individual score per pore (indicative of the directionality) was assigned to all pixels filling its corresponding surface, and subsequently smoothed out with a Gaussian kernel. High values within the pore directionality maps correspond to those regions in which pores are oriented similarly to the median pore angle, ultimately indicating regions characterized by a predominant pore orientation.

To complement these analyses, we developed a graphical user interface (GUI) for the analysis of a single image/batch displaying fiber networks. Fibers are enhanced using Gabor filters and represented by graphs. Parametric maps, such as fiber length and pore directionality maps can subsequently be derived. The output results can be written into .png image files, while the fiber specific graph/Gabor-derived features are collected in .csv files. The MATLAB source code and sample images for testing can be found on the GitHub platform, at github.com/aigrapa/ECM-fiber-graph.

4.3. Test for anomaly detection using fiber network simulations

To simulate fiber networks, we considered a set of scattered point patterns $ {P}_d $ , where $ {P}_d $ is the set of points $ (x+{n}_x\ast d,y+{n}_y\ast d) $ , for $ {n}_y $ odd, and $ (x+{n}_x\ast d,y+({n}_y+0.5)\ast d) $ for $ {n}_y $ even, $ {n}_x,{n}_y,d\in \mathrm{\mathbb{N}} $ . Random noise was added at the points’ location : $ \forall p=\left(X,Y\right)\in {P}_d,{p}^{\prime }=\left(X+{\eta}_x,Y+{\eta}_y\right) $ , where $ {\eta}_x $ and $ {\eta}_y $ follow a uniform law between $ 0 $ and $ N $ , $ N\in \mathrm{\mathbb{N}} $ . In this way, we generated two patterns $ {P}_1 $ and $ {P}_2 $ using different values for d, as well as a mask M (for the second scenario), which is an image equal to zero except for specific areas (e.g., ellipses Figure 6a). We consider P as: $ \left\{\left(x,y\right)\in {P}_1:M\left(x,y\right)=0\right\}\cup \left\{\left(x,y\right)\in {P}_2:M\left(x,y\right)\ne 0\right\} $ . The fibers were subsequently defined by the edges of the Delaunay graph of P.

Figure 6. Anomaly detection within parametric maps of simulated fiber networks (a) Simulations of fiber networks (1024x1024 pixels), isotropic (left) and with local defects (center), ground-truth mask (right). (b) Graph-based representations of fiber networks, and corresponding fiber length parametric maps for both samples. (c) Detection of anomalous clusters with respect to the normal GRF model (at pval $ \le 0.05 $ ), at various thresholds, on the parametric map containing defects, for intensity-based (left) and surface-based criteria (center). The regions detected at a threshold of 20, based on a surface-based criterion having a non-null intersection with those detected at a threshold of 35, according to an intensity-based criterion (right).

The first isotropic fiber network corresponds to a “normal” example ( $ d=30 $ , $ N=20 $ ), while the fiber network with local defects (i.e., fibers are more elongated in the regions containing “defects,” corresponding to the three regions within the mask) is considered here an “abnormal” example ( $ d\in \left\{40,50\right\} $ , $ N=20 $ ) (Figure 6a). Fibers were detected as explained in 4.1, and fiber graphs were correspondingly derived for both images (Figure 6b). Starting from the graph-based fiber representation, fiber length parametric maps were generated accordingly, as described in 4.2 (Figure 6b), and subsequently “gaussianised,” as explained in 2.3. We were interested in applying the same principle described in 2.2, in order to detect the three regions of fiber length variation corresponding to the proposed ground-truth mask (Figure 6a). According to this principle, the null hypothesis is that clusters of pixels in the parametric map, computed at different thresholds in the abnormal example, belong to a realization of the same GRF as the normal sample. The method learns the GRF model parameters from the parametric map corresponding to the reference-normal example, and then uses these parameters to compute the two probabilities of belonging to GRF, for each region at various thresholds, in the abnormal sample’s parametric map, using an intensity or surface-based criterion. Different regions were identified at various thresholds, for both intensity and surface-based detection (Figure 6c), at a pval $ \le 0.05 $ . By only keeping the regions which were detected at a certain threshold (surface-based detection) having a non-null intersection with the clusters detected according to the intensity-based criterion, we were able to accurately detect the three regions of fiber length variation, as well as two additional false positive smaller regions within the parametric map of the fiber network with local defects.

7.1. Materials and FN preparations

Recombinant human TGF- $ \beta $ 1 was from R $ \& $ D Systems Inc. (Minneapolis, MN, USA). All other chemicals and reagents were purchased from Sigma Aldrich (St Louis, MO, USA) unless otherwise stated. Purified recombinant FN variants were produced as previously described⁽Reference Efthymiou, Radwanska and Grapa ^{11 )}.

7.2. Cells and culture conditions

Fn1 -/- mouse kidney fibroblasts were generated and cultured as previously described⁽Reference Efthymiou, Radwanska and Grapa ^{11 )}. For experiments, FN was depleted from fetal calf serum using gelatin sepharose-4B columns (GE Healthcare, Uppsala, Sweden), and the culture medium was supplemented with Penicillin-Streptomycin 100 U/ml and, where indicated, TGF- $ \beta $ 1 (5 ng/ml). Absence of Mycoplasma sp. contamination was routinely verified by PCR as described elsewhere⁽Reference Kong, James, Gordon, Zelynski and Gilbert ^{31 )}.

7.3. Generation of fibroblast-derived matrices, immunofluorescence staining and microscopy

Fibroblast-derived matrices were generated as described previously. For FN immunostaining, primary antibody (rabbit polyclonal anti-FN) was from Merck-Millipore (Darmstadt, Germany). Fluorescently-labeled (Alexa Fluor 488-conjugated) secondary antibody was from Thermo Fisher Scientific (Waltham, Massachusetts). After staining, the coverslips were mounted in ProLong Gold antifade reagent (Thermo Fischer Scientific). Confocal imaging was performed on a Zeiss LSM710 confocal system equipped with a 10X/0.45 NA objective. For visual representation, image treatment was performed using Fiji⁽Reference Schindelin, Arganda-Carreras and Frise ^{32 )}.