Anca-Ioana Grapa

Georgios Efthymi

Ellen Van Obberghen-Schilling

Laure Blanc- Féraud

Xavier Descombes
email: xavier.descombes@inria.fr

A spatial statistical framework for the parametric study of fiber networks: application to fibronectin deposition by normal and activated fibroblasts

Keywords: Fibronectin, Extracellular matrix, Oncofetal isoforms, Statistical parametric maps, Graph networks

Due to the complex architectural diversity of biological networks, there is an increasing need to complement statistical analyses with a qualitative and local description of their spatial properties. One such network is the extracellular matrix (ECM), a biological scaffold for which changes in its spatial organization significantly impact tissue functions in health and disease. Quantifying variations in the fibrillar architecture of major ECM proteins should considerably advance our understanding of the link between tissue structure and function. Inspired from the analysis of functional magnetic resonance imaging (fMRI) images, we propose a novel statistical analysis approach embedded into a machine learning paradigm, to measure and detect local variations of meaningful ECM parameters. We show that parametric maps representing fiber length and pore directionality can be analyzed within the proposed framework to differentiate among various tissue states. The parametric maps are derived from graph-based representations that reflect the network architecture of fibronectin (FN) fibers in a normal, or disease-mimicking in vitro setting. Such tools can potentially lead to a better characterization of dynamic matrix networks within fibrotic tumor microenvironments and contribute to the development of better imaging modalities for monitoring their remodeling and normalization following therapeutic intervention.

Impact Statement

Quantification of phenotypic variation during tissue development and/or disease progression is essential for the understanding of different pathologies. All organs and tissues contain a non-cellular core component known as the ECM, composed of a network of macromolecules whose architecture depends on the pathophysiological state of the tissue. To derive a meaningful comparison of ECM between healthy and diseased tissues, computational frameworks that account for the localization of areas of phenotypic variation are needed. Here we introduce a novel framework for the statistical analysis of parametric maps calculated from graph-based representations of fibers composed of FN, a provisional ECM component that guides ECM organization. Our framework is inspired from the statistical analysis of fMRI parametric maps and is embedded in a machinelearning model to compare distinct states of ECM fiber networks, both quantitatively and qualitatively. These methods may be further developed and implemented in ECM profiling of tumor/fibrotic tissue to provide both valuable insight into specific roles of ECM landscapes and their remodeling in disease, and more specific diagnostic, prognostic, and predictive companion biomarkers in the clinic.

Introduction

During normal development and disease progression, tissues undergo various remodeling processes, which can, in turn, affect their physical characteristics, yielding heterogeneous morphologies. Automated detection and quantification of these phenotypic changes in the tissue landscape are essential for an accurate characterization of a given pathology. Statistical tests that are commonly used for the comparison of two different conditions based on the distributions of morphological properties, are applied at a global scale, and do not account for any explicit spatial information. Here we were interested in exploiting a known spatial statistical approach (historically proposed for functional imaging (fMRI) analysis) and recasting it into a machine learning framework to facilitate the comparison of two tissue conditions. Within the proposed framework relying on statistical parametric mapping (SPM) [START_REF] Poline | Combining Spatial Extent and Peak Intensity to Test for Activations in Functional Imaging[END_REF][START_REF] Friston | Assessing the Significance of Focal Activations Using Their Spatial Extent[END_REF] , the comparison of various physical tissue characteristics is thus achieved both at a quantitative and qualitative level. It does so by enabling the localization and quantification of local variations of certain morphological properties in the sample that are significantly different and relevant to a given pathology.

SPM is a long-established methodology, specifically developed in fMRI for the detection of significantly activated regions of the brain in a given image sample. Mapping of activated regions is achieved by assessing the probability of random occurrences of activated regions with pixel intensities higher than a given threshold or having a larger spatial extent at lower intensities.

To address the need for taking spatial localization into account when designing frameworks that can discriminate between two different conditions of a given tissue, we recast the SPM paradigm into a data-driven machine learning framework for detection of significant parametric differences between the two classes. Thus, we train the model on a given population (for example the control case) and we detect local areas in the second population of samples that deviate from this model. In our work, we applied this approach to the characterization of two distinct states of the extracellular matrix (ECM), a non-cellular component of organs and tissues.

The ECM is a biological scaffold with multiple forms and functions. It acts as a biomechanical and structural support ensuring tissue integrity, it relays chemical and physical signals to the residing cells through cell surface receptors and it sequesters growth factors and regulates their bioavailability.

The composition and architecture of the ECM is tissue-and organ-specific, and depends on the pathophysiological state of the tissue (i.e. normal vs diseased) [START_REF] Frantz | The extracellular matrix at a glance[END_REF] . For example, while a healthy connective tissue displays a loose meshwork-like ECM, a fibrotic or cancerous stroma is characterized by the presence of dense, aligned ECM fibrils. Thus, the physical and structural traits of the tumor matrix have recently drawn much attention as cancer hallmarks and potential therapeutic targets [START_REF] Afratis | Novel Approaches for Extracellular Matrix Targeting in Disease Treatment[END_REF][START_REF] Nia | Physical traits of cancer Science[END_REF] . Collagen, the most abundant matrix component, has been extensively investigated in this context and several studies addressing its structural features and their association to cancer progression, metastasis and treatment have been published [START_REF] Conklin | Aligned collagen is a prognostic signature for survival in human breast carcinoma[END_REF][START_REF] Provenzano | Collagen reorganization at the tumorstromal interface facilitates local invasion[END_REF][START_REF] Burke | Using second harmonic generation to predict patient outcome in solid tumors[END_REF] . Collagen deposition, however, depends on the presence of fibronectin (FN), a dimeric glycoprotein that forms a provisional matrix framework to which other ECM components integrate to generate a mature ECM [START_REF] Hynes | The extracellular matrix: Not just pretty fibrils[END_REF][START_REF] Barker | The provisional matrix: setting the stage for tissue repair outcomes[END_REF] .

During inflammation, wound healing, or tumor development, the expression of FN is induced and assembled into an insoluble matrix. This FN produced primarily by fibroblasts, corresponds to cellular FN, as opposed to plasma FN. At the molecular level, cellular FN differs from plasma FN by the presence of one or two 90-amino acid-long alternatively spliced sequences, termed Extra Domains (EDB and EDA). Extra Domain-containing FN displays enhanced assembly, making it the most prevalent FN isoform in diseased tissue. This enhanced FN deposition results in a highly modified ECM architecture with increased fiber density, directionality and stiffness, that together tune cellular responses and impact tissue homeostasis [START_REF] Efthymiou | Fibronectin Extra Domains tune cellular responses and confer topographically distinct features to fibril networks[END_REF] .

Despite the pivotal role of FN in health and disease [START_REF] Efthymiou | Shaping Up the Tumor Microenvironment With Cellular Fibronectin[END_REF] , comprehensive studies of FN fiber features are lacking. In our previous work [START_REF] Efthymiou | Fibronectin Extra Domains tune cellular responses and confer topographically distinct features to fibril networks[END_REF] , we set out to develop a robust pipeline of numerical analyses for https://doi.org/10.1017/S2633903X23000247 Published online by Cambridge University Press the extraction of biologically relevant metrics to discriminate among different isoforms of cellular FN from confocal microscopy images of FN matrices. The goal of the present work was to capture and analyze physiologically relevant ECM fibrillar features that discriminate between normal and diseased states. To that end, we generated FN-rich ECM using an in vitro model of cell derived matrices (CDMs) produced by normal fibroblasts, or fibroblasts activated with transforming growth factor beta 1 (TGF-𝛽1), a fibrosis-promoting cytokine known to induce a tumor-like state.

Herein we show that the proposed SPM-based machine learning methodology can be used to distinguish between normal and disease-mimicking FN fiber networks. To capture and localize significant differences of fiber properties, we built parametric maps, such as fiber length and pore directionality relying on a graph-based fiber representation that recapitulates the FN fiber networks from confocal microscopy images. In the following sections, we will provide an outline of the proposed methodology applied to SPMs and describe how our proposed machine-learning embedding can yield significant localized parametric variations between different tissue states.

Statistical parametric maps (SPM) and Gaussian Random Fields

SPM statistical framework

Statistical parametric maps 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 [START_REF] Poline | Combining Spatial Extent and Peak Intensity to Test for Activations in Functional Imaging[END_REF] . The value of the pixel intensity reflects a parametric value of interest. Hence these (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 a 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 diseasemimicking FN organization. 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 [START_REF] Adler | The Geometry of Random Fields[END_REF] that enables the statistical analysis of tissue parametric maps. GRF, whose marginal distributions are Gaussian vectors 𝑋 = (𝑋 (1) , • • • , 𝑋 (𝑛) ), are characterized by the probability density function:

𝑓 𝑋 (𝑥) = (2𝜋) -𝑛/2 |𝑉 | -1 2 𝑒𝑥 𝑝 [- 1 2 (𝑥 -𝜇)𝑉 -1 (𝑥 -𝜇) 𝑇 ] (1)
where 𝜇 = (𝐸 (𝑋 (𝑖) )) 𝑖 ∈ {1,...,𝑛} is the expectation and 𝑉 = (𝐸 [(𝑋 (𝑖) -𝜇 (𝑖) ) (𝑋 ( 𝑗 ) -𝜇 ( 𝑗 ) )]) 𝑖 ∈ {1,...,𝑛}, 𝑗 ∈ {1,...,𝑛} 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 neighbouring 8 pixels. It was shown in [START_REF] Adler | The Geometry of Random Fields[END_REF] , and later adopted in fMRI-specific studies [START_REF] Poline | Combining Spatial Extent and Peak Intensity to Test for Activations in Functional Imaging[END_REF][START_REF] Friston | Assessing the Significance of Focal Activations Using Their Spatial Extent[END_REF][START_REF] Worsley | A Three-Dimensional Statistical Analysis for CBF Activation Studies in Human Brain[END_REF] , that for large thresholds 𝑡, the clusters are independent and the expectation of the number of clusters at a threshold 𝑡, 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:

𝐸 [m t ] = 𝑆(2𝜋) -3/2 |Λ| 1/2 𝑡𝜎 -3 exp - 𝑡 2 2𝜎 2
(2)

https://doi.org/10.1017/S2633903X23000247 Published online by Cambridge University Press where:

• m t represents the number of clusters at a certain threshold 𝑡

• S is the number of pixels of the image

• Λ is the covariance matrix of partial derivatives of the GRF

• 𝜎 is the standard deviation of the GRF Similarly, the mean value of the number of clusters at a threshold 𝑡 + 𝐻 0 can be written as such:

𝐸 [m t+𝐻 0 ] = 𝑆(2𝜋) -3/2 |Λ| 1/2 (𝑡 + 𝐻 0 )𝜎 -3 exp - (𝑡 + 𝐻 0 ) 2 2𝜎 2 (3)
Considering 𝑥 0 = 𝑡+𝐻 0 as the intensity peak of a cluster (at threshold 𝑡), one can estimate the probability that a cluster (at a threshold 𝑡, having an intensity peak higher or equal to 𝑥 0 , denoted 𝐶 𝐻 0 𝑡 ) belongs to a realization of this GRF, 𝐺 𝑟 . This probability, as shown in [START_REF] Poline | Combining Spatial Extent and Peak Intensity to Test for Activations in Functional Imaging[END_REF] , termed 𝑃 𝐻 , can be seen as the likelihood of a cluster (formed at threshold 𝑡) of having an intensity peak higher or equal to 𝑡 + 𝐻 0 [START_REF] Poline | Combining Spatial Extent and Peak Intensity to Test for Activations in Functional Imaging[END_REF] :

𝑃(𝐶 𝐻 0 𝑡 ∈ 𝐺 𝑟 ) = 𝐸 [m t+𝐻 0 ] 𝐸 [m t ] = 𝑥 0 𝑡 exp 𝑡 2 -𝑥 2 0 2𝜎 2 (4)
Next, we were interested in the estimation of the probability that a cluster (at a threshold 𝑡) belongs to a realization of GRF, depending on its surface (spatial extent -number of pixels). To estimate the number of pixels (𝑛 𝑡 ) of a cluster at a threshold 𝑡, we use the following equation from [START_REF] Friston | Assessing the Significance of Focal Activations Using Their Spatial Extent[END_REF] :

𝐸 [𝑛 𝑡 ] = 𝐸 [𝑁 𝑡 ] 𝐸 [𝑚 𝑡 ] (5)
where 𝑁 𝑡 is the number of pixels at of higher intensity than 𝑡, and 𝑚 𝑡 is the number of clusters at the threshold 𝑡. Since the intensity values follow a normal (zero mean value) distribution, the expectation of 𝑁 𝑡 is the following [START_REF] Friston | Assessing the Significance of Focal Activations Using Their Spatial Extent[END_REF] :

𝐸 [N t ] = 𝑆 ∫ ∞ 𝑡 (2𝜋𝜎 2 ) -1/2 exp - 𝑥 2 2𝜎 2 dx = 𝑆Φ 𝜎 (-𝑡) (6)
where Φ 𝜎 (-𝑡) is the complementary cumulative distribution function. It follows, then, based on equations (2,5,6) that one can approximate the mean value of 𝑛 𝑡 , accordingly:

𝐸 [𝑛 𝑡 ] = 𝐸 [N t ] 𝐸 [𝑚 𝑡 ] = Φ 𝜎 (-𝑡) (2𝜋) -3/2 |Λ| 1/2 𝑡𝜎 -3 exp -𝑡 2 2𝜎 2 (7)
Furthermore, 𝑛 𝑡 follows an exponential distribution law [START_REF] Lafarge | Détection de feux de forêt par analyse statistique de la radiométrie d'images satellitaires HAL[END_REF] , which is commonly defined by a parameter 𝜆 𝑡 , (the inverse of the mean expected value of the random variable). Consequently:

𝑃(𝑛 𝑡 = 𝑥) = 𝜆 𝑡 𝑒𝑥 𝑝(-𝜆 𝑡 𝑥) (8)
where

𝜆 𝑡 = (2 𝜋 ) -3/2 |Λ| 1/2 𝑡 𝜎 -3 exp -𝑡 2 2𝜎 2
Φ 𝜎 (-𝑡 )

https://doi.org/10.1017/S2633903X23000247 Published online by Cambridge University Press

It follows then that the approximation for the probability 𝑃 𝑆 of a given cluster 𝐶 𝑆 0 𝑡 having a spatial extent 𝑆 greater than 𝑆 0 , is given by the following formulation, as shown in [START_REF] Friston | Assessing the Significance of Focal Activations Using Their Spatial Extent[END_REF] :

𝑃(𝐶 𝑆 0 𝑡 ∈ 𝐺 𝑟 ) = 𝑃(𝑛 𝑡 ≥ 𝑆 0 ) = 𝑒𝑥 𝑝(-𝜆 𝑡 𝑆 0 ) = 𝑒𝑥 𝑝 (2𝜋) -3/2 |Λ| 1/2 𝑆 0 𝑡𝜎 -3 exp -𝑡 2 2𝜎 2 Φ 𝜎 (-𝑡) (9)
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.

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 𝑃 𝐻 or 𝑃 𝑆 is less or equal than a p-value (pval) of 0.05.

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 and Figure 1B, the current method, compared to a naïve hard thresholding (threshold equal to 10), is better suited to localize the abnormal elements, i.e.

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.

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 [START_REF] Conklin | Aligned collagen is a prognostic signature for survival in human breast carcinoma[END_REF] , 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 [START_REF] Friston | Assessing the Significance of Focal Activations Using Their Spatial Extent[END_REF] . 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)). 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, i.e. 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 [START_REF] Peyré | Computational Optimal Transport: With Applications to Data Science[END_REF][START_REF] Peyré | The Numerical Tours of Signal Processing -Advanced Computational Signal and Image Processing[END_REF] (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 𝑡 to belong to a GRF, depending on the maximal intensity of this cluster, we relied on the formulations taken from the theory 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, 𝑡 = 50 (B), 𝑡 = 80 (C), 𝑡 = 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 .

of random fields [START_REF] Adler | The Geometry of Random Fields[END_REF] , as previously described. Considering 𝑥 0 = 𝑡 + 𝐻 0 as the intensity peak of a cluster (at intensity threshold 𝑡), one can estimate the probability that a cluster having an intensity peak higher or equal to 𝑥 0 , belongs to a realization of GRF. This probability can be seen as the likelihood of a cluster (taken at threshold 𝑡) of having an intensity peak higher or equal to 𝑡 + 𝐻 0 (Equation 4). Furthermore, as previously shown, the approximation for the probability of a given cluster having a spatial extent S greater than 𝑆 0 is given by Equation 9. At pval ≤ 0.05, the clusters of pixels identified at 𝑡 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 https://doi.org/10.1017/S2633903X23000247 Published online by Cambridge University Press 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, i.e. the average value of 𝜆 𝑗 , 𝜎 𝑗 , from the training set. These learnt parameters during the learning phase were subsequently used to compute the two relevant probabilities, 𝑃 𝐻 and 𝑃 𝑆 , as described hereafter:

For all the images I 𝑗 (previously Gaussianized) in the learning set:

• Computation of Λ 𝑗 (Equation 2, empirical estimator of the covariance of partial derivatives of 𝐼 𝑗 ). If for an image function 𝑓 ∈ R 2 , we consider its gradient vector

∇ 𝑓 = ( 𝑓 𝑥 , 𝑓 𝑦 ) = ( 𝜕 𝑓 𝜕𝑥 , 𝜕 𝑓 𝜕𝑦 ), then Λ 𝑗 = cov[ 𝑓 𝑥 , 𝑓 𝑦 ] = 𝐸 [( 𝑓 𝑥 -𝐸 [ 𝑓 𝑥 ]) ( 𝑓 𝑦 -𝐸 [ 𝑓 𝑦 ])].
• Computation of 𝜎 𝑗 , as the 𝐼 𝑗 's sample standard deviation.

The last step involves storing the average Λ 𝑚 , 𝜎 𝑚 of the learning dataset.

During the test phase, for all the clusters identified at a threshold 𝑡 within the test set, 𝑃 𝐻 and 𝑃 𝑆 are evaluated using the model's previously learnt parameters. At pval ≤ 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 I 𝑗 in the test set:

• Gaussianization of each sample image 𝐼 𝑗 . The result is a new image 𝐼 𝑔 , whose histogram is Gaussian with identical variance to that of 𝐼 𝑗 .

• For a given list of thresholds

𝑇 = (𝑡 1 , 𝑡 2 , • • • , 𝑡 𝑛 ) :
-Binarization of the image I 𝑔 according to the threshold 𝑡 𝑖 -Once the list of connected components in the binary image resulted from thresholding is achieved, then for every (labeled) connected-component

(𝑙 1 , 𝑙 2 , • • • , 𝑙 𝑝 ) :
* Evaluation of 𝑃 𝐻 using the learnt model parameter 𝜎 𝑚 .

* Evaluation of 𝑃 𝑆 using the learnt model parameters Λ 𝑚 , 𝜎 𝑚 .

Generation of FN variant-specific fibroblast-derived matrices and induction of a tumor-like state

To test our method, we utilized a previously established in vitro system to generate normal and diseasemimicking ECMs by normal mouse fibroblasts. Fibroblasts are the major ECM-producing cells of tissues. In pathological conditions (e.g. tumors), quiescent fibroblasts become activated by environmental cues that induce their phenotypic conversion to 'myofibroblasts' with a pro-tumoral phenotype (Figure 3A, top). This process is characterized by the upregulation of cellular FN expression, actin reorganization and increased cell contractility that result in their elongation and the deposition of a highly anisotropic FN-rich ECM [START_REF] Efthymiou | Fibronectin Extra Domains tune cellular responses and confer topographically distinct features to fibril networks[END_REF] . In vitro, these changes can be mimicked by treating normal resting fibroblasts with TGF-𝛽1, a potent cytokine involved in fibroblast activation in the tumor microenvironment [START_REF] Efthymiou | Fibronectin Extra Domains tune cellular responses and confer topographically distinct features to fibril networks[END_REF] . For our analyses, FN-rich normal or tumor-like matrices were generated by presenting FN-null mouse embryo fibroblasts with recombinant cFN (prepared as previously described [START_REF] Efthymiou | Fibronectin Extra Domains tune cellular responses and confer topographically distinct features to fibril networks[END_REF] ), as schematized in Figure 3B. For the induction of a tumor-like phenotype, fibroblasts were incubated with TGF-𝛽1 (Figure 3C). Cultures were decellularized after 7 days, and the resulting cell-derived matrices https://doi.org/10.1017/S2633903X23000247 Published online by Cambridge University Press

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 [START_REF] Efthymiou | Fibronectin Extra Domains tune cellular responses and confer topographically distinct features to fibril networks[END_REF] . Existing ECM analysis methods focusing on collagen exploit alternative fiber detection techniques, such as Fast Fourier Transform bandpass filters [START_REF] Morrill | A validated software application to measure fiber organization in soft tissue[END_REF] , ridge detection [START_REF] Wershof | A FIJI macro for quantifying pattern in extracellular matrix[END_REF] , or fast discrete curvelet transform [START_REF] Liu | Fibrillar Collagen Quantification With Curvelet Transform Based Computational Methods[END_REF] , 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 [START_REF] Efthymiou | Fibronectin Extra Domains tune cellular responses and confer topographically distinct features to fibril networks[END_REF] , 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 [START_REF] Grapa | Characterization of fibronectin networks using graph-based representations of the fibers from 2D confocal images[END_REF] . 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 [START_REF] Kollmannsberger | The small world of osteocytes: Connectomics of the lacuno-canalicular network in bone[END_REF][START_REF] Aguilar | Graph-Based Methods for Retinal Mosaicing and Vascular Characterization[END_REF] .

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 [START_REF] Efthymiou | Fibronectin Extra Domains tune cellular responses and confer topographically distinct features to fibril networks[END_REF] .

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.

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

https://doi.org/10.1017/S2633903X23000247 Published online by Cambridge University Press 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 [START_REF] Bresenham | Algorithm for computer control of a digital plotter[END_REF] 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 [START_REF] Errico | Interpolation and extrapolation of elements in a 2D array[END_REF] and smoothing of this result with a Gaussian kernel. Concerning the pore directionality parametric maps (Figure 5B), the starting point was the skeleton graph.

https://doi.org/10.1017/S2633903X23000247 Published online by Cambridge University Press 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 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.

Test for anomaly detection using fiber network simulations

To simulate fiber networks, we considered a set of scattered point patterns 𝑃 𝑑 , where 𝑃 𝑑 is the set of points (𝑥 + 𝑛 𝑥 * 𝑑, 𝑦 + 𝑛 𝑦 * 𝑑), for 𝑛 𝑦 odd, and (𝑥 + 𝑛 𝑥 * 𝑑, 𝑦 + (𝑛 𝑦 + 0.5) * 𝑑) for 𝑛 𝑦 even, 𝑛 𝑥 , 𝑛 𝑦 , 𝑑 ∈ N. Random noise was added at the points' location :

∀𝑝 = (𝑋, 𝑌 ) ∈ 𝑃 𝑑 , 𝑝 ′ = (𝑋 +𝜂 𝑥 , 𝑌 +𝜂 𝑦 ),
where 𝜂 𝑥 and 𝜂 𝑦 follow a uniform law between 0 and 𝑁, 𝑁 ∈ N. In this way, we generated two patterns 𝑃 1 and 𝑃 2 using different values for d, as well as a mask M (for the second scenario), https://doi.org/10.1017/S2633903X23000247 Published online by Cambridge University Press

which is an image equal to zero except for specific areas (e.g. ellipses Figure 6A). We consider P as:

{(𝑥, 𝑦) ∈ 𝑃 1 : 𝑀 (𝑥, 𝑦) = 0} {(𝑥, 𝑦) ∈ 𝑃 2 : 𝑀 (𝑥, 𝑦) ≠ 0}.
The fibers were subsequently defined by the edges of the Delaunay graph of P.

The first isotropic fiber network corresponds to a 'normal' example (𝑑 = 30, 𝑁 = 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 (𝑑 ∈ {40, 50}, 𝑁 = 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 ≤ 0.05. By only keeping the regions which were detected at a certain threshold (surface-based detection) having a nonnull 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.

Results -statistical analysis of fiber parametric maps

The graph-based representation of FN networks enabled the subsequent design of a novel framework to perform a spatial statistical analysis of ECM patterns, using graph-derived statistical parametric maps. This methodology was applied for a quantitative and qualitative analysis of fiber length and pore directionality differences, across all FN variant networks in normal and tumor-like states. We were thus interested in determining whether the proposed SPM analysis of the selected spatial fiber features could reveal significant variant-specific differences between the FN variant networks in normal (N) and tumor-like (T) states.

To apply our framework to the available data, we first divided the available sets of confocal images (1024 x 1024 pixels, 0.27 𝜇m/pixel; 70 images/variant for normal FN (N) and 65 images/variant for tumor-mimicking FN networks) as follows. For comparison of (N) vs (T) FN networks, we considered 50 (N) samples as the learning dataset, 20 (N) as a test set for normal, and 65 (T) as a test set for disease-like networks. In all scenarios, a cluster is considered significantly different from the normal GRF model at pval ≤ 0.05. Anomalies in fiber length (Figure 7A) detected using either intensity or surface-based criteria (Figure 7B,C) and pore directionality maps (Figure 8A,B,C) were detected at a few intensity thresholds (e.g. 70,80,90) and (10,12,14), respectively (Figure 7, Figure 8). Thus, using our approach, differences in fiber length and pore directionality could be localized in regions formed at different intensity thresholds. This property is very useful for obtaining a qualitative analysis of parametric maps, where clusters of pixels that are statistically different from a normal model can be localized.

For the quantitative analysis of tissue parametric variations, we set out to determine significant differences between FN variant networks through the average number of identified foreign clusters, as well as the average cluster area per image. It is noteworthy that the group for which anomalous clusters were found at superior thresholds, had higher parametric values than in the normal model. with respect to fiber length, then fibers are statistically more elongated within former networks than normal counterparts. As shown in Figure 7,8, we found both fiber length and pore directionality to be significantly different for pairwise comparisons of normal and tumor-like FN variant networks. Essentially, the latter type of FN architecture relative to normal ECM, is represented by statistically longer fibers (Table 1) with a more pronounced pore directionality (Table 2). The increased length of FN fibers is consistent with the elongated phenotype of the TGF-𝛽1-treated fibroblasts that assemble them. The statistically significant increase in pore directionality of tumor-like matrices compared to normal FN matrices is in agreement with published reports [START_REF] Park | Extracellular matrix anisotropy is determined by TFAP2C-dependent regulation of cell collisions[END_REF] of higher FN alignment in cell-derived matrices from cancer-associated fibroblasts and in tumor tissue. Detailed results, including the average area and number of identified anomalous clusters, at multiple thresholds, for four FN variants, are presented in Supplementary TableS1-TableS8. Table 1. Quantitative analysis for detection of differences in fiber length -Anomalous cluster quantification (with respect to the normal statistical model, at pval ≤ 0.05), for the comparison of normal (N) and tumor-like (T) FN (1024x1024 pixels, or 276.48 𝜇𝑚 x 276.48 𝜇𝑚). The average number of significant clusters per test database for each variant is shown here, for either surface or intensity criteria, if higher than for the normal model, at any selected threshold). Average number of clusters detected in fiber length maps, for normal vs tumor-mimicking FN, for one of the test thresholds (70,80,90). Table 2. Quantitative analysis for detection of differences in pore directionality -Anomalous cluster quantification (with respect to the normal statistical model, at pval ≤ 0.05), for the comparison of normal (N) and tumor-like (T) FN ( 1024x 1024 pixels, or 276.48 𝜇𝑚 x 276.48 𝜇𝑚). The average number of significant clusters per test database for each variant is shown here, for either surface or intensity criteria, if higher than for the normal model, at any selected threshold). Average number of clusters detected in pore directionality maps, for normal vs tumor-mimicking FN, for one of the test thresholds (10,12,14). '-' is recorded if no significant detection was present.

FN variant Average number of clusters for (N) vs (T) FN B-A-0.12 (intensity) B+A-0.22 (surface) B-A+ 0.14 (intensity); 0.34 (surface) B+A+ -

Discussion and conclusion

352

The proposed methodology was designed to quantify the differences in terms of spatial organization 353 between normal and disease-like architectures of variant-specific FN matrices generated in vitro. We have previously been able to discriminate the matrix patterns of four alternatively spliced FN variants deposited by cultured fibroblasts using different learning approaches and relying on graph-based feature analysis. The pipeline which includes steps for fiber detection and representation [START_REF] Efthymiou | Fibronectin Extra Domains tune cellular responses and confer topographically distinct features to fibril networks[END_REF] , and generation of parametric fiber maps is available with a MATLAB GUI: both graph and Gabor filter-based fiber features can be extracted from different images representing fiber networks (e.g. ECM-specific proteins) for downstream analysis. In principle, the steps required for the characterization of ECM (i.e. fiber detection and representation using graphs) in cell-derived matrices generated by cells of different origin remain the same. However, in the case of tissue samples, the ECM is more complex and heterogeneous, and further pre-processing steps may be needed to filter structures that are normally found in tissues (e.g. blood vessels). The number and type of additional steps, however, is dependent on the type of organ/tissue, the pathology under evaluation, and the staining procedure.

Here, we developed a statistical parametric map framework for the quantitative and qualitative analysis of fibers, capable of simultaneously detecting and measuring variations of specific ECM features within two different tissue conditions. Importantly, our approach can be evaluated on parametric maps at different thresholds, producing results that are more reliable and statistically relevant (by providing a pval) than a simple hard thresholding of the parametric maps. Our framework was tested using two relevant fiber features, fiber length and pore directionality, whose parametric maps revealed significant differences between normal and disease-mimicking states. However, parametric maps can be extended to include other fiber or pore-specific parameters (e.g. fiber density, width, length, orientation, waviness, and straightness), which could be useful for differentiating among various biological networks in normal and pathological states.

Computational analyses of ECM structures can provide essential information about their role in shaping the cellular microenvironment topology in health and during disease progression. Prognostic ECM-specific signatures have already been inferred in cancer-related studies, and in diseases with prominent fibrosis [START_REF] Parker | Extracellular matrix profiles determine risk and prognosis of the squamous cell carcinoma subtype of non-small cell lung carcinoma[END_REF][START_REF] Feng | Characterization of Matricellular Protein Expression Signatures in Mechanistically Diverse Mouse Models of Kidney Injury[END_REF][START_REF] Basak | Quantitative proteomic profiling of extracellular matrix and site-specific collagen post-translational modifications in an in vitro model of lung fibrosis Matrix[END_REF] . There is also a growing interest in the integration of cell and ECM analyses in a spatially resolved manner to further understand the interactions between cells and their matrix microenvironment [START_REF] Vasiukov | Integrated Cells and Collagen Fibers Spatial Image Analysis Front Bioinform[END_REF] . Indeed, the present work proposing a versatile pipeline for the analysis of ECM produced by cultured fibroblasts, is being extended to studies of ECM organization in human tumor tissue and aims to integrate the phenotypes of cellular components. Hence, a combined local analysis of parametric maps and metrics describing the organization/morphology of adjacent cells (e.g. tumor, immune, vascular cells) will potentially help elucidate the complex interplay between cellular and noncellular components of the tumor microenvironment.

Materials and methods

Materials and FN preparations

Recombinant human TGF-𝛽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 [START_REF] Efthymiou | Fibronectin Extra Domains tune cellular responses and confer topographically distinct features to fibril networks[END_REF] .

Cells and culture conditions

Fn1 -/mouse kidney fibroblasts were generated and cultured as previously described [START_REF] Efthymiou | Fibronectin Extra Domains tune cellular responses and confer topographically distinct features to fibril networks[END_REF] . 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-𝛽1 (5 ng/ml). Absence of Mycoplasma sp. contamination was routinely verified by PCR as described elsewhere [START_REF] Kong | Species-Specific PCR for Identification of Common Contaminant Mollicutes in Cell Culture[END_REF] .

https://doi.org/10.1017/S2633903X23000247 Published online by Cambridge University Press

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). Fluorescentlylabeled (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 [START_REF] Schindelin | Fiji: an open-source platform for biological-image analysis[END_REF] .

Figure 1 .

 1 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 𝑡 = 10 pixels) (B). Intensity and surface-based detection at pval ≤ 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, i.e ellipses are detected at pval ≤ 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 .

Figure 2 .

 2 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, 𝑡 = 50 (B), 𝑡 = 80 (C), 𝑡 = 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 .

Figure 3 . 10 4.

 310 Figure 3. (A) Schematic representation of a simple cuboidal epithelium displaying the different architectures of the underlining ECM in normal (left) and pathological conditions (right). (B) Workflow diagram featuring the linear structure of the purified recombinant FN (rFN) variants and the relative positions of the alternatively spliced Extra Domains, the generation of fibroblast-derived matrices, and image acquisition and analysis. (C) Phase contrast images (top row) of FN-null mouse fibroblasts presented with FN B+A+ variant (15 𝜇g/ml) in the presence or absence of TGF-𝛽1 (5 ng/ml) to mimic the changes that take place in the tumor/fibrotic microenvironment. After removal of the cells, matrices were stained with a rabbit-anti-FN polyclonal antibody and visualized with confocal microscopy. Scale bars: phase, 100 𝜇m; IF, 50 𝜇m .

Figure 4 .

 4 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 𝜇m/pixel. (B) Fiber enhancement with Gabor filters (C) Morphological fiber skeleton extraction (D) Skeletonbased graph (left) and simplified graph representation (right) which is derived from the latter .

Figure 5 .

 5 Figure 5. Computation of fiber parametric maps: (A) Starting from the skeleton graph (FN B-A+ disease-like sample, 1024x1024 pixels, 0.27 𝜇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 .

Figure 6 .

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

Figure 7 .

 7 Figure 7. Qualitative analysis -Anomalous cluster detection (with respect to the normal statistical model), applied to two samples of fiber length map (FN B-A+ tumor-like), 1024x1024 pixels, 0.27 𝜇m/pixel (A), (B) and (C) depict the anomalous clusters (pval ≤ 0.05) at various intensity thresholds (70,80,90) .

Figure 8 .

 8 Figure 8. Qualitative analysis -Anomalous cluster detection (with respect to the normal statistical model), applied to two samples of pore directionality map (FN B-A+ tumor-like), 1024x1024 pixels, 0.27 𝜇m/pixel (A), (B) and (C) depict the detected clusters (pval ≤ 0.05) at various intensity thresholds (10,12,14) .

https://doi.org/10.1017/S2633903X23000247 Published online by Cambridge University Press

Acknowledgements. We thank the members of the Adhesion Signaling and Stromal Reprogramming in the Tumor Microenvironment team for critical discussion on the manuscript, iBV PRISM imaging platform for the use of their machines and support, in particular, Sébastien Schaub. Funding Statement. This work was supported by the French Government (National Research Agency, ANR) through the 'Investments for the Future' LABEX SIGNALIFE: program reference ANR-11-LABX-0028-01, ANR-16-CE93-0005 (ANGIO-FIB), ANR-19-P3IA-0002 (3IA Côte d'Azur) and La Fondation ARC (PJA20151203207).

Data Availability Statement. Sample image data and code can be found on the GitHub platform: github.com/aigrapa/ ECM-fiber-graph.

Ethical Standards. The research meets all ethical guidelines, including adherence to the legal requirements of the study country.

Supplementary Material. A supplementary file has been provided for submission.


Arnaud Allais

Jean-Maxime Saugrain

Beate West

Nicolas Lallouet

Hervé Caron

David Ferandelle

Louise Terrien

Grégory Bouvier

Ghazi Hajiri

Kévin Berger

Loïc Quéval

SuperRail -World-first HTS cable to be installed on a railway network in France

Keywords: Cooling system, railway network, superconducting cables, Turbo Brayton technology

teaching and research institutions in France or abroad, or from public or private research centers.

I. INTRODUCTION

HE French rail network supplied with direct current at a voltage of 1500 V is 5904 km long. This represents around a quarter of the national rail electrification [START_REF] Steimel | Under Europe's incompatible catenary voltages a review of multi-system traction technology[END_REF]. A major disadvantage of this low-voltage configuration is the voltage drops along the lines, which necessitate the installation of substations close together. The French railway network is also subject to a number of constraints, notably linked to the 60% increase in traffic between 1998 and 2008, which is set to

The SuperRail project is supported in part by the French government via the Banque Publique d'Investissement (BPI) as part of the France 2030 plan. Corresponding author: Kévin Berger.

Arnaud Allais, Jean-Maxime Saugrain and Nicolas Lallouet are with Nexans,

France (e-mail: arnaud.allais@nexans.com, jean_maxime.saugrain@nexans.com, nicolas.lallouet@nexans.com).

Beate West is with Nexans, Germany (e-mail: beate.west@nexans.com). Hervé Caron and David Ferandelle are with Département de la Traction Électrique, SNCF Réseau, 93418 La Plaine Saint-Denis, France (e-mail: herve.caron@reseau.sncf.fr, david.ferandelle@reseau.sncf.fr).

multiply by a factor of 4 by 2030. In addition, it is becoming increasingly difficult to install substations in densely populated urban areas. Superconducting power cables make it possible to avoid major construction work, particularly on civil engineering structures, while improving the quality of the power supply and reducing energy losses. This is the choice made by SNCF as part of the SuperRail project, launched in February 2022, which involves the development, manufacture, installation and long-term operation of two 60-meter long 1.5 kV-3.5 kA high-temperature superconducting (HTS) DC cables to reinforce the power supply to the Montparnasse station in Paris, France.

This project follows on from several superconducting cable projects and demonstrators around the world. Since 2000, we have pinpointed no less than 35 superconducting cable projects: 26 for AC [START_REF] Stovall | Installation and operation of the Southwire 30-meter high-temperature superconducting power cable[END_REF]- [START_REF]Nexans installs and commissions superconducting cable for Chicago's Resilient Electric Grid project[END_REF] and 9 for DC [START_REF] Hamabe | Critical Current and Its Magnetic Field Effect Measurement of HTS Tapes Forming DC Superconducting Cable[END_REF]- [START_REF] Ballarino | The BEST PATHS Project on MgB2 Superconducting Cables for Very High Power Transmission[END_REF].

II. MONTPARNASSE-VOUILLÉ DEMONSTRATOR

A. Motivation

The Montparnasse railway station in Paris, shown in Fig. 1, offers 28 tracks and a daily throughput of 200,000 passengers on 750 trains. But recent electrical failures [START_REF] Julien | Hauts-de-Seine: un impressionnant incendie s'est déclenché sur un poste électrique à Issy-les-Moulineaux[END_REF], [START_REF] Beaugé | Pourquoi la gare Montparnasse, touchée par une nouvelle panne dimanche, est particulièrement vulnérable aux perturbations[END_REF] underscore the urgent need for reinforced electrical infrastructure. In this context, SNCF needs to bring additional power from the Ouest-Ceinture Vouillé traction substation to a railway tracks platform 60 meters away at Montparnasse station (Fig. 2). The requirement is for 4.5 MW of nominal power (3000 A), and 10.5 MW of inrush power (7000 A) during the train's acceleration phases.

The existing connection between the substation and the tracks relies on 400 mm 2 copper cables, installed in a cable duct bank integrated into an aging bridge. However, with only two unoccupied conduits remaining (Fig. 3), there is insufficient capacity to accommodate the six copper cables required to carry the nominal 3000 A at 1500 V. Reinforcing with copper cable would have required the installation of a new cable duct bank, which would have necessitated substantial civil engineering work, estimated at over 10 million euros.

Louise Terrien and Grégory Bouvier are with Absolut System, 38170 Seyssinet-Pariset, France (e-mail: louise.terrien@absolut-system.com, gregory.bouvier@absolut-system.com).

Ghazi Hajiri and Kévin Berger are with Université de Lorraine, GREEN, 54000 Nancy, France (e-mail: ghazi.hajiri@univ-lorraine.fr, kevin.berger@univ-lorraine.fr).

Loïc Quéval is with Université Paris-Saclay, CentraleSupélec, CNRS, GeePs, 91190 Gif-sur-Yvette, France. And with Sorbonne Université, CNRS, GeePs, 75252 Paris, France. (e-mail: loic.queval@centralesupelec.fr).

T

The SuperRail project aims to perform the electrical reinforcement using two HTS cables laid in the two existing conduits, thus avoiding major civil engineering work. This could have been done with only one HTS cable but, for operational reasons, it has been decided to install two 1500 A DC HTS cables. This is an excellent example of how the high energy density of HTS cables can be leveraged to reinforce the electric grid in dense urban areas.

B. Cable and termination layout

A 3D digital twin of the installation site was created using LiDAR (Light Detection And Ranging) analysis. It enables the virtual integration of the HTS cable system into its future environment. This is useful to estimate the length of the HTS cable to be manufactured, confirm the design of the various components, and anticipate problems.

Fig. 4 shows the planned cable path. Each of the two 60meter HTS cables will be pulled inside one empty 100 mm diameter conduit. Cable pulling requires adapting to the geometry of the existing conduits, installed several decades ago, and to the presence of the copper cables. One of the first tasks carried out was to check that the HTS cables could be correctly pulled into the conduits, taking advantage of two newly built manholes. The investigations identified strong gradients and two critical areas with bending radii of 1.7 m and 2.1 m. The pulling of a dummy cable was then carried out onsite to ensure that the pulling of the real cable would be successful (Fig. 5). It gives an insight into the installation of an HTS cable on a railway network.

Two terminations will be installed in the substation and two on a platform as shown in Fig. 6. Each termination has two current leads: one for the conductor layer and one for the screen layer. The conductor is connected to the positive pole of the AC/DC converter on the substation side and to the feeder on the platform side. The screen is connected to the neutral pole of the AC/DC converter on the substation side only. The screens of both cables are connected to each other on the platform side.

III. SUPERRAIL PROJECT INNOVATIONS

The HTS cable system is subjected to several operating constraints which have influenced its design. Table I summarizes all the figures provided in this section.

A. HTS cables

Each cable must carry a nominal current of 1500 A (maximum 3% harmonics below 5 kHz) and sustain repetitive inrush currents of up to 3500 A during train acceleration. Even more demanding, each cable must withstand a short-circuit current of 67 kA for 200 ms, corresponding to a fault power of 100 MW. Fig. 7 illustrates the cable that has been developed by Nexans. The cable core is made of several layers. The copper former has a diameter of approximately 18 mm. This figure is relatively high to withstand the fault current. The HTS conductor consists of 25 4.4 mm wide second-generation (2G) HTS tapes arranged in two layers. It has a minimum critical current of 4000 A to withstand inrush currents of 3500 A for a few seconds. Two different manufacturers of HTS tapes have been selected: SuNAM [START_REF] Sunam | 2G HTS Wire[END_REF] and AMSC [START_REF]Amperium® HTS Wire[END_REF]. The electric insulation is a 1-mm thick Polypropylene Laminated Paper (PPLP) that is impregnated by liquid nitrogen which contributes to the dielectric strength. The cable insulation is designed for 3000 V, but it will operate at 1500 V for the SuperRail installation. The 1.8-mm thick copper screen is made up of several copper tapes. The cable core is inserted into a linear cryostat. This is a vacuum-insulated flexible cryostat made of corrugated stainless-steel tube that provides the thermal insulation to maintain the cable core around 70 K while minimizing thermal losses. The outer diameter of the cryostat is 74 mm, and the cable has a minimum bending radius of 1.5 m, in line with the site constraints.

As part of the SuperRail project, an HTS cable production line has been assembled at the Nexans' Bourg-en-Bresse factory in France.

B. Terminations

In a project such as SuperRail, where the cable length is short, the design of the terminations plays an important role, as it is the element that contributes most to the losses that must be extracted by the cooling system. In this project, terminations account for around 70% of the system total losses at nominal current. The terminations developed by Nexans for this project were adapted from earlier designs made for the REG project in Chicago [START_REF] Ross | Resilient electric grid networks: achieving enhanced resiliency in today's grid[END_REF]. These are specialized components that still need to be standardized, so that they can be supplied and installed easily and reliably.

C. Cooling system

Absolut System designed the cooling system, working with all the partners to define the system requirements:

▪ Cooling capacity of 1700 W at 67 K for the cables and terminations (1200 W) and other cryogenic losses (500 W); ▪ Closed loop cooling system; ▪ Low maintenance. These requirements rule out Pulse Tube and Gifford-McMahon technologies. A Stirling machine was initially envisaged, but it has three major drawbacks: no cold power regulation to adapt to load variations, high maintenance needs due to the piston, and need to import key elements. Internal R&D was unable to find solutions to these problems. The decision was therefore taken to develop a Reverse Turbo-Brayton (RTB) machine, the main sub-assemblies of which are manufactured in France.

The 5 stages of the RTB closed cycle are shown in Fig. 8:  Compression of helium gas with centrifugal compressors operating at ambient temperature: increasing the pressure and temperature of the gas;  Pre-cooling in the first branch of the recuperative heat exchanger;  Expansion in the centripetal turbine: pressure and temperature are reduced to the target cryogenic temperature;  Extraction of cold power from the helium gas in the cold generation cycle and recovery by the liquid nitrogen in the cables loop in the application exchanger;  Heating of the fluid in the recuperative heat exchanger.

The RTB uses high-speed turbomachines and high-efficiency recuperative heat exchangers. Hydrodynamic gas bearings are used for the turbomachine rotor, which operate at speeds ranging from 30,000 rpm for the largest to over 300,000 rpm for the smallest. These non-contact bearings guarantee reduced maintenance and high reliability. This technology is scalable for future applications (not limited to railways) and is notably proposed for future European observation satellites. The consumption of an RTB machine is typically 10 kW of electrical power for 1 kW of cold thermal power at 80 K. Moreover, the cooling capacity can be controlled from 30% to 100% of the maximum power. The overall architecture of the cooling system connected to the HTS cables is shown in Fig. 9. Note that the cooling system is one-sided: the LN2 will flow in one cable and return in the other one. It consists of:

▪ A LN2 tank used for the initial cooling of the cable and refills after maintenance; ▪ A skid that integrates the circulation and pressurization of LN2 and overall cooling management. A LN2 pump in a vacuum cryostat to circulate liquid nitrogen in the cables loop; ▪ A Reverse Turbo-Brayton (RTB) machine to remove the heat from the cryogenic circuit; ▪ A set of valves, sensors and a programmable logic controller to control the system. The skid has been assembled and is being tested. It is pictured in Fig. 10. A 3D CAD drawing of the RTB machine is presented in Fig. 11. It is now being manufactured by Absolut System. This cooling system will be installed right next to the Vouillé substation. Thanks to the closed loop operation, there will be almost no LN2 consumption, except for filling or maintenance operations. This limits the number of times the LN2 tank has to be refilled each year, greatly simplifying operation.

IV. TYPE TEST LOOP

Cable and termination designs are complete. A 35-meter cable type test loop shown in Fig. 12 is being installed at the SNCF Railway Test Agency (SNCF-AEF), starting in summer 2023. The test loop includes two terminations, one joint, and the cooling system. Although it is not needed for the installation in the Vouillé substation, a joint has been designed for future projects and will be tested. These tests will qualify cables, terminations and joints for SuperRail and similar projects, whatever the length. results of these tests will be the subject of a specific report.

V. ANTICIPATION OF FUTURE PROJECTS

In addition to being a technology demonstrator, SuperRail aims to assess the competitiveness of using a superconducting cable versus a resistive cable solution for future railway installations. This is essential to gauge the market potential and economic prospects of superconducting cable systems. In this regard, the GREEN laboratory at the University of Lorraine is actively developing computational tools to simulate superconducting cables in railway networks. The models take into account electromagnetism, heat transfers, hydraulics and circuit equations, including the structure of the electrical network [START_REF] Hajiri | Design and modelling tools for DC HTS cables for the future railway network in France[END_REF]- [START_REF] Hajiri | Optimization of the terminations of an HTS cable operating on a DC railway network[END_REF]. These models will be validated using data from the SNCF-AEF type test loop and the Montparnasse-Fig. 10. Cryogenic skid, being tested. Vouillé demonstrator. In parallel, the GeePs laboratory at University Paris-Saclay is carrying out feasibility studies to integrate a superconducting transformer and a cryogenic power electronic converter into the cable system, as shown in Fig. 13 [45], [START_REF] Ferreira | Feasibility study of a cryogenic power supply for superconducting DC devices[END_REF]. This architecture would eliminate the need for two terminations while taking advantage of cryogenic temperatures to improve the converter's efficiency and power density. These efforts collectively contribute to estimating potential and forecasting future projects, thus playing an essential role in the advancement of the technology.

VI. CONCLUSION

The SuperRail project responds to the need to improve railway power supply in urban environments, such as the Montparnasse station in Paris. Through the use of a hightemperature superconducting DC cable system, it introduces an innovative approach to reinforcing electrical infrastructure, avoiding the need for extensive civil engineering works. SuperRail is the world-first HTS cable to be installed on a railway network, and the first HTS cable commercial installation in France.

This article described the technology adopted for the demonstrator and the first steps towards the installation. In parallel with the demonstrator, SuperRail aims to evaluate the technology for other installations and anticipate the future.

The cable will be installed and commissioned in 2024. It will then be operated by SNCF for its entire lifetime. Feedback will be provided on the commissioning and operation of this HTS cable system.

Fig. 1 .

 1 Fig. 1. Montparnasse railway station in Paris, 28 tracks and a daily throughput of 200,000 passengers on 750 trains.

Fig. 2 .

 2 Fig. 2. The Vouillé substation (in orange) in the heart of Paris and the platform (in green) that will house the HTS cable terminations right next to the Montparnasse rail tracks.

Fig. 3 .

 3 Fig.3. Existing cable duct bank between the substation and the platform. Each conduit is filled with one 400 mm 2 copper cable (500 A each).

Fig. 2 .

 2 Fig. 2. The Vouillé traction substation (in orange) in the heart of Paris and the railway tracks platform (in green) that will house the HTS cable terminations at Montparnasse station.

Fig. 3 .Fig. 4 .

 34 Fig.3. Existing cable duct bank between the substation and the platform. Each conduit is filled with one 400 mm 2 copper cable (500 A each).

Fig. 5 .Fig. 6 .Fig. 7 .

 567 Fig. 5. Pulling of a dummy cable in Montparnasse-Vouillé.

Fig. 8 .

 8 Fig. 8. RTB cooling machine operating principle.

Fig. 11 .

 11 Fig. 11. 3D CAD drawing of the RTB machine, being assembled.

Fig. 9 .

 9 Fig. 9. Architecture of the cooling system connected to the HTS cables, in nominal operation.

Fig. 12 .

 12 Fig. 12. Cable type test loop installed at the SNCF Railway Test Agency (SNCF-AEF), including two terminations and a joint.

Fig. 13 .Fig. 11 .

 1311 Fig. 13. a) Conventional solution, b) Montparnasse-Vouillé demonstrator (TRL 6-7), c) Lab prototype (TRL 4) with a superconducting transformer and a power electronic converter operating at cryogenic temperature.

TABLE I SUPERRAIL

 I CABLE SYSTEM IN SUMMARY

	Component	Parameters	Value
		Electrical parameters
		Number of cables	2
		Nominal voltage	1500 V
		Inrush current	3500 A
		Nominal current	1500 A
		Critical current	4000 A
		Fault current	67 kA -200 ms
		Tape parameters
	HTS Cable	Tape manufacturers	SuNAM & AMSC
		Tape width	4.4 mm
		Geometric parameters
		Copper core diameter	18 mm
		Number of tapes	12 + 13
		PPLP thickness	1 mm
		Copper screen thickness	1.8 mm
		Cryostat outer diameter	74 mm
		Length	60 m
		Minimum bending radius	1.5 m
		Cooling system parameters
	Reverse	Overall cooling capacity	1700 W @ 67 K
	Turbo-
	Brayton	Available power for the cables	1200 W @ 67 K
		Cable inlet temperature	> 64 K
	Cooling	Cable outlet temperature	< 74 K
	system	Inlet pressure	< 14 bar
		Outlet pressure	> 7.5 bar

ACKNOWLEDGMENT

We would like to thank our colleagues at SNCF, Nexans, Absolut System, University Paris-Saclay and the University of Lorraine for the valuable insights they provided and the indepth discussions we had.