# Remodelling of the fibre-aggregate structure of collagen gels by cancer-
associated fibroblasts: a time-resolved grey-tone image analysis based on
stochastic modelling

Cedric J. Gommes1<EMAIL_ADDRESS>Thomas Louis2 Isabelle Bourgot2 Erik
Maquoi2 Silvia Blacher2 Agnès Noël2 1Department of Chemical Engineering,
University of Liège B6A, Allée du Six Août 3, B-4000 Liège, Belgium
2Laboratory of Tumor and Development Biology, GIGA-Cancer, University of
Liège, Liège, Belgium

###### Abstract

Solid tumors consist of tumor cells associated with stromal and immune cells,
secreted factors and extracellular matrix (ECM), that together constitute the
tumor microenvironment. Among stromal cells, cancer-associated fibroblasts
(CAFs) are of particular interest. CAFs modulate the architecture of the ECM
by exerting forces and contracting collagen fibres, creating paths that
facilitate cancer cell migration. The characterization of the collagen fibre
network and its space and time-dependent microstructural modifications is key
to investigating the interactions between cells and the ECM. Developing image
analysis tools for that purpose is still a challenge. The structural
complexity of the collagen network calls for specific statistical descriptors.
Moreover, the low signal-to-noise ratio of imaging techniques available for
time-resolved studies rules out standard image analysis methods. In this work,
we propose an innovative data analysis approach to investigate the remodelling
of fibrillar collagen in a three-dimensional spheroid model of cellular
invasion, which is based on the stochastic modelling of the gel structure and
on grey-tone image analysis. The method was used to study the reshaping of a
collagen matrix by migrating breast cancer-derived CAFs. The structure of the
gels at the scale of a few microns consists in regions with high fibre density
separated by depleted regions, which can be thought of as fibre aggregates and
pores. To characterize this structure, we developed a two-scale stochastic
model with a clipped Gaussian field model to describe the aggregates-and-pores
large-scale structure, and a homogeneous Boolean model to describe the small-
scale fibre network within the aggregates. The model parameters are identified
by fitting the grey-tone histograms and correlation functions of confocal
microscopy images. The specificity of the method is that it applies to the
unprocessed grey-tone images, and it can therefore be used with low
magnification, noisy time-lapse reflectance images. When applied to the
spheroid time-resolved images, the method revealed two different matrix
densification mechanisms for the gel in direct contact or far from the cells.

##  I Introduction

In living tissues, cells are commonly embedded within complex networks of
extracellular matrix (ECM) constituted of diverse highly cross-linked
components, including fibrous proteins, proteoglycans, glycoproteins and
polysaccharides Hynes (2009); Hynes and Naba (2012); Oskarsson (2013). Each of
the individual ECM components exhibits specific biomechanical and biochemical
properties that are related to its polymer structure, size and binding
affinities to signaling molecules Frantz  _et al._ (2010); Theocharis  _et
al._ (2016). These properties in turn modulate several cellular processes
including proliferation, differentiation, survival and motility by ligating
specialized cell surface receptors Leitinger and Hohenester (2007); Xian  _et
al._ (2010). Cell-matrix interactions also lead to ECM remodeling by the
cells, with the ECM fibres being synthesized, re-oriented, deformed and
degraded by the cells, particularly fibroblasts Shieh  _et al._ (2011). The
topology of the ECM (organization of fibres, matrix porosity and density) also
influences the mechanical properties of the matrix Doyle  _et al._ (2009);
Charras and Sahai (2014). For instance, aligned bundles of collagen increase
matrix stiffness Frantz  _et al._ (2010). Along with fibre organization, the
degree of porosity of the ECM influences its rheological properties which in
turn modulates the cell behavior Seo  _et al._ (2020).

Collagens are the most abundant components of the ECM and their structure and
composition differ across various tissue types Yue (2014); Mouw  _et al._
(2014). The collagen family is composed of 28 types consisting in fibril-
forming, fibril-associated, network-forming and other structures subfamilies
Bourgot  _et al._ (2020); Revell  _et al._ (2021). Type I collagen represents
the most common fibrillar collagen in vertebrates. Its synthesis, stiffening
and remodelling are involved in both physiological processes such as wound
healing, inflammation, tissue repair and pathological ones such as
angiogenesis or tumor fibrosis which promotes growth and invasion Provenzano
_et al._ (2019, 2008); Buchmann  _et al._ (2021); Bourgot  _et al._ (2020).
Many solid tumors such as those of the breast, lung and pancreas are
characterized by the presence of a desmoplastic stroma typified by an
accumulation of fibrillar collagen. This collagen accretion increases the
stiffness and creates discrete structural patterns called tumor-associated
collagen signatures (TACS) in the stroma. These anisotropic patterns promote
directed cell migration into and through the stroma by contact guidance
Dickinson  _et al._ (1994); Provenzano  _et al._ (2006); Ray and Provenzano
(2021); Ray  _et al._ (2022). The stroma comprises a complex ecosystem of
endothelial cells, immune cells, and cancer-associated fibroblasts (CAFs)
Hanahan and Coussens (2012). The latter make up the majority of the
desmoplastic stroma and have been shown to promote the growth of primary
tumors Kalluri (2016); Biffi and Tuveson (2021); Sahai  _et al._ (2020).

Understanding how CAFs influence the architecture of collagen is key to
improving our understanding of cancer cell invasion through a dense
collagenous stroma as observed in breast, lung and pancreatic cancers. At
macroscopic scales, this can be studied with gel-contraction assays Bell  _et
al._ (1979); Dallon and Ehrlich (2008); Mikami  _et al._ (2016) and bead
tracking Steinwachs  _et al._ (2015); Hall  _et al._ (2013); Pakshir  _et al._
(2019); Mark  _et al._ (2020), whereby the mechanical deformation of the
collagen matrix is monitored without having to explicitly consider the
underlying remodelling of the collagen network. More detailed understanding of
the process, based on local influence of individual cells on the matrix and
the biochemical interactions involved can be obtained using atomic-force and
confocal microscopy at the scale of the cells Frye  _et al._ (2018). However,
few studies have attempted to reconcile the macroscopic and cellular
approaches to characterize the influence of cells on collagen architecture.

Developing image analysis tools to characterize the collagen-fiber network and
its space- and time-dependent microstructural modifications in the context of
cellular migration is challenging for a variety of reasons. In particular, the
low signal-to-noise ratio of imaging techniques suitable for time-resolved
studies at the required resolution precludes standard image analysis methods
based on the segmentation of the fibrillar structures Bredfeldt  _et al._
(2014); Liu  _et al._ (2018). Grey-tone image analysis methods have been
developed in this context, but they focus on the analysis of fibre orientation
using a variety of methods such as pixel-wised gradient estimations Boudaoud
_et al._ (2014), Fourier transforms Pijanka  _et al._ (2019), or mathematical
morphological operations Altendorf  _et al._ (2012). In addition to fibre
orientation, however, key aspects of ECM remodelling concern the spatial
distribution of collagen fibres. Grey-tone image analysis methods that can
address that question, and statistically capture possibly non-homogeneous
fibre densification patterns, are yet to be developed.

The paper proposes such a method based on the stochastic modelling of grey-
tone structures in microscopy images, and on the identification of the model
parameters from their correlation functions and grey-tone histograms.
Acellular collagen gels with different concentrations are considered first, in
order to develop and validate the methodology. The question of ECM remodelling
is then addressed in the context of a collagen-embedded CAF spheroid model,
and the local evolution of the fibre network that accompanies cellular
migration is studied in space- and time-resolved way.

##  II Materials and Methods

Dulbecco’s modified Eagle medium (DMEM), L-glutamine, sodium pyruvate,
penicillin, streptomycin and 0.25% trypsin/EDTA solution were purchased from
ThermoFisher Scientific. Recombinant platelet-derived growth factor BB isotype
(PDGF-BB) was obtained from R&D Systems. Fetal bovine serum (FBS), sodium
bicarbonate, carboxymethylcellulose and 10$\times$ concentrated DMEM, high
viscosity carboxymethylcellulose sodium salt were obtained from Sigma-Aldrich.
High concentration acid soluble native type I rat tail collagen was purchased
from Corning and DQ-Collagen™type I from bovine skin fluorescein conjugate was
obtained from ThermoFisher. SPY650-DNA, a non-toxic, cell permeable and highly
specific live cell DNA probe was purchased from Spirochrome.

###  II.1 Preparation of collagen gels

Type I collagen gels were prepared for imaging of collagen architecture and
spheroid invasion assays. Collagen gels of 2.0 mg/mL and 3.0 mg/mL were
prepared by diluting the stock collagen solution (8-11 mg/mL) with 10x
concentrated DMEM, NaHCO3, 1N NaOH, milliQ H2O and neutralized to pH 7.2. To
visualize collagen fibers by confocal fluorescence microscopy (CFM),
fluorogenic DQ-collagen I was mixed with unlabeled diluted collagen as
previously described Del Amo  _et al._ (2018), obtaining a final concentration
of 20 $\mu$g/mL. Collagen dilutions were performed and maintained on ice until
use. Droplets (25 $\mu$L) of diluted collagens were spotted in 35mm-glass
bottom $\mu$-dishes (Ibidi) and polymerized during 45 min at 19∘C, followed by
a 30 min incubation at 37∘C. The polymerization temperature of 19∘C was chosen
to provide matrix architecture comparable to in vivo with a network consisting
of a mixture of a few thin bundles and many thick bundles Geraldo  _et al._
(2012). After completion of collagen polymerization, 1 mL of preheated DMEM
supplemented with 5% FBS was added to the dishes. Image acquisition was
performed after 24 h of incubation at 37∘C with 5% CO2.

###  II.2 Cell culture

Mouse CAFs (CAFs) have been isolated from mammary gland tumors of mammary
specific polyomavirus middle T antigen overexpression mouse model (MMTV-PyMT)
at 12 weeks as previously described Primac  _et al._ (2019) and immortalized
with the pLenti HPV16 E6-E7 RFP, expressing a cytoplasmic red fluorescent
protein (RFP). CAFs were grown in high-glucose DMEM supplemented with 10% FBS,
2mM L-glutamine, 1mM sodium pyruvate, 100IU/ml penicillin, 100$\mu$g/mL
streptomycin. Cultures were maintained at 37∘C with 5% CO2 until their
confluence reaches about 80%.

###  II.3 Spheroid invasion assay

Spheroids were prepared by seeding 1000 CAFs in 100 $\mu$L of spheroid
formation medium composed of 0.22 $\mu$m-filtered DMEM medium supplemented
with 10% FBS and 20% carboxymethylcellulose 4000 centipoise. Cells were seeded
in round-bottom non-adherent 96-well plates (CELLSTAR, Greiner Bio-One) and
centrifugated at 1000 rpm for 5 min. Plates were incubated at 37∘C with 5% CO2
for 48 hours to promote spheroid formation. The content of each well was
transferred in a petri dish with a 200 $\mu$L pipette (with cutted tip) and
individual spheroids were collected under a binocular microscope with a 10
$\mu$L pipette. Each spheroid was resuspended in 23 $\mu$L of diluted collagen
(2 mg/mL) and spotted as a 25 $\mu$l drop in a prechilled 8 well-glass bottom
chamber slide (Ibidi). The slides were then transferred immediately to 19∘C
and flipped to maintain the spheroids in the middle of the collagen drop
(preventing their sedimentation to the glass surface or to the collagen/air
interface). The extent of collagen polymerization and spheroid positioning
were carefully controlled by microscopic examination throughout the
polymerization step. After 30 min at 19∘C, the slides were transferred at 37∘C
to complete the polymerization. Preheated culture medium (300 $\mu$L/well)
supplemented with 5% FBS, 10 ng/mL PDGF-BB and SPY650-DNA (1000-fold dilution)
was added to the slides and time-lapse imaging was initiated within an hour
after collagen polymerization.

###  II.4 Confocal microscopy

Images of the three-dimensional (3D) collagen gels were acquired with an
inverted confocal laser scanning microscope (LSM 880 Airyscan Elyra S1, Zeiss)
with a Plan-Neofluar 10$\times$/0.30 N.A. or a Plan-Neofluar 20$\times$/0.50
N.A. objective (Zeiss). The DQ-collagen containing gels were excited with 488
nm laser. Non-fluorescent collagen gels were imaged by confocal reflectance
microscopy (CRM) in Airyscan high resolution mode with a simultaneous
excitation of the matrix by 488 nm and 633 nm lasers. To avoid edge effects,
images were acquired at least 100 $\mu$m away from the gel border, avoiding
regions close to the gel/glass and gel/medium interfaces. To visualize
collagen fibres of acellular gels (gels containing no cells), samples were
imaged both by CFM and CRM, using the 20$\times$ objective. The resulting
images have dimensions of 1000$\times$1000$\times$90 voxels with anisotropic
voxel size of 0.42$\times$0.42$\times$1.10 $\mu$m3, corresponding to a
physical volume of approximately 500$\times$500$\times$100 $\mu$m3.

Time-lapse imaging of spheroid-containing collagen gels was performed using
the 10$\times$ objective with the samples incubated at 5% CO2 and 37∘C in the
on-stage incubator (Okolab). The collagen matrix was imaged by CRM in Airyscan
high resolution mode with a 1.4$\times$ digital zoom (scaling per pixel: 0.59
$\mu$m $\times$ 0.59 $\mu$m $\times$ 3.29 $\mu$m) and CAF were imaged by CFM
in Fast Airyscan mode (scaling per pixel: 0.91$\mu$m $\times$ 0.91 $\mu$m
$\times$ 3.29 $\mu$m). CAF images were rescaled to fit the images of the
collagen matrix. CAF-derived RFP was excited by a 561 nm laser and detected at
591 nm. The SPY650 DNA nuclear stain was excited by a 633nm laser and the
signal was captured at 654 nm. The imaging procedures for the spheroids and
the matrix required a sequential imaging of the same tridimensional zone,
leading to two 3D image files. Images were recorded every 30 min up to 16 h.
3D stacks were obtained at a step size of 2-$\mu$m intervals. Type I collagen
fibrils have a diameter ranging from 20 nm to several hundred nm Christiansen
_et al._ (2000) while fibres are larger in diameter. Given that the size of
each pixel is 0.42 and 0.59 $\mu$m (for 20$\times$ and 10$\times$ objectives,
respectively), it is not possible to distinguish fibrils from fibres Lai  _et
al._ (2012), therefore the term “fibres” was used to include both fibrils and
fibres.

The raw images were converted using the Airyscan algorithm from the Zeiss Zen
Black software. The images were then subjected to a histogram stretching and
converted to tiff format through the Zeiss Zen Blue software. This operation
turned the original 16-bit images into 8-bit images, which were corrected by
histogram adjustment to maximize the conservation of valuable information.

##  III Results

###  III.1 Covariance and grey-tone correlation function

Figure 1: Images of the 2 mg/mL (a) and 3 mg/mL (b) acellular gels obtained by
confocal microscopy in fluorescence (a1 and b1) and reflectance (a2 and b2)
modes. The bottom row displays fused images with fluorescence and reflectance
data in green and magenta, respectively.

Examples of images of the dilute (2 mg/mL) and concentrated (3 mg/mL) gels
obtained through confocal fluorescence (CFM) and confocal reflectance (CRM)
microscopy are given in Fig 1. These are 2D single-$z$ images taken out of the
3D images. The structural analysis is based on 15 images such as presented in
the figure, for each gel and each imaging mode. Type I collagen fibrils have a
diameter ranging from 20 nm to several hundred nm Christiansen  _et al._
(2000) while fibres are larger in diameter. Given that the size of each pixel
is 0.42 and 0.59 $\mu$m (for 20x and 10x objectives, respectively), it is not
possible to distinguish fibrils and fibres Lai  _et al._ (2012), therefore the
term “fibres” was used to include both fibrils and fibres. The most salient
structures are the fibre aggregates that are a few tens of micrometers across.
The fused images in Fig 1a3 and b3 show that there is limited overlap between
the CRM and CFM microscopy data: the fibres and aggregates are well captured
in CFM but it is mostly the aggregates that are visible in the CRM data.

In order to follow a standard image analysis procedure, and justify further
less-standard developments, we first explored a method based on image
segmentation. In that spirit, the grey-tone images of the gels are converted
to binary images following a method described in the Supplementary Material
(see Fig. S1). Examples of segmented images of the gels, with collagen in
white and the rest in black, are provided in the insets of Figs. 2a and 2b.

The segmented gel images display complex and disordered structures, which are
also corrupted by noise in the case of the reflectance images (Figs. 2b1 and
2b2). In this context, the covariance - which describes the spatial
correlation between all pixels intensities in the images Serra (1982);
Lantuéjoul (2002); Torquato (2002); Jeulin (2021)\- was used to quantitatively
analyze the gel structures. The covariances shown in Fig. 2a and 2b were
obtained from fifteen images taken in the same gel with dimension 200 $\times$
200 $\mu m^{2}$ as illustrated in the insets, via a Fourier-transform
algorithm, and the error bars are the standard errors of the means.

Figure 2: Experimental covariance (top) and grey-tone correlation function
(bottom) of the 2 mg/mL (blue circle, a1, b1, c1, d1) and 3 mg/mL (red
diamond, a2, b2, c2, d2) gels, imaged in fluorescence (a, c) and reflectance
(b, d) modes. The insets display one of the fifteen images used for each
condition, and the error bars are the standard errors of the mean.

The covariance $C(r)$ of say, the white component of an image, has a clear
geometrical interpretation as it is defined as the probability that any couple
of points at distance $r$ from one another both belong to that phase. For very
small values of $r$, this definition coincides with the probability for a
single point to belong to the white phase, which is numerically equal to its
density $\phi_{1}$ normalized between 0 and 1. In the opposite limit, that is
for very large distances $r$ the covariance converges to the value
$\phi_{1}^{2}$, which corresponds to a horizontal asymptote in Fig 2. The
shape of the covariance curve between those two limits characterizes the
structures present in the images. In particular the progressive decrease of
$C(r)$ over distances of a few tens of microns testifies to the presence of
structures with those dimensions, which we qualitatively referred to earlier
as fibre aggregates.

The covariance of the fluorescence images of the gels (in Fig. 2a) highlights
at once the unsuitability of image segmentation in the present context. The
covariance of the 2mg/mL gel is found to be larger than that of the 3mg/mL gel
for any $r$, which means that densities of the segmented images contradict the
actual collagen concentrations of the gels. This results from the fact that
collagen-rich areas of the images display a variety of grey-tones related to
the local fibre density (see Fig. S1b), which information is lost during the
all-or-nothing segmentation procedure. This general observation calls for
grey-tone image analysis methods that preserve the structural information in
the images. In that spirit, the grey-tone correlation functions $C(r)$ of the
gels are shown in Fig. 2c and 2d. The latter are measured directly on the
unprocessed images and they characterise the statistical correlation between
grey-tones of all pixels that are at distance $r$ from one another. In the
case where the images contain only the values $0$ and $1$, this definition is
mathematically equivalent to the covariance. We postpone to a later section
the discussion of the structural significance of the grey-tone correlation
function, but we already notice at this stage that the grey-tone data in Figs.
2c and 2d scale with the actual collagen concentration of the 2mg/mL and
3mg/mL gels as they should.

###  III.2 Structural models

Covariance and grey-tone correlation functions convey indirect yet very rich
structural information Aubert and Jeulin (2000); Jiao  _et al._ (2009); Gommes
_et al._ (2012), which can notably be retrieved using structural models. We
here present two models aimed at extracting structural information from the
covariance data in Fig. 2, which we generalize later to grey-tone correlation
functions. In order to cope with the disordered structure of the gel, the
models have to be stochastic Jeulin (2000); Torquato (2002); Gommes (2018);
Jeulin (2021).

When using stochastic models, a structure is defined through probabilistic
rules and this calls for specific concepts. In particular, it is convenient to
introduce the indicator function of the structure $\mathcal{I}(\mathbf{x})$,
which takes the value 1 if the point $\mathbf{x}$ belongs to the structure and
0 otherwise Torquato (2002). With such definition, the density of the model is
calculated as

$\phi_{1}=\langle\mathcal{I}(\mathbf{x})\rangle$ (1)

where the brackets $\langle\rangle$ stand for the average value, calculated
either over $\mathbf{x}$ or over all the possible realization of the model.
Similarly, the covariance is calculated as the following two-point average

$C_{11}(r)=\langle\mathcal{I}(\mathbf{x})\mathcal{I}(\mathbf{x}+\mathbf{r})\rangle$
(2)

because the product
$\mathcal{I}(\mathbf{x})\mathcal{I}(\mathbf{x}+\mathbf{r})$ is equal to one,
only if the points $\mathbf{x}$ and $\mathbf{x}+\mathbf{r}$ belong to the
structure. In the latter equation, we have assumed statistical isotropy so
that the dependence is only through the modulus $r=|\mathbf{r}|$.

####  III.2.1 Homogeneous fibre model

The simplest model we consider to analyze the structures in Fig. 1 assumes
that the gel matrix is statistically homogeneous. The model consists in
tossing fibres (modelled as elongated rectangles) with random position and
orientation, as sketched in Fig. 3. Such model is described by three
parameters, namely: the number density of fibres $\theta$ (unit $\mu$m-2), as
well as their length $L_{F}$ and diameter $D_{F}$ (both in units of $\mu$m).
Note that the statistical homogeneity of the model does not rule out the
existence of aggregates, which form when a large number of fibres
coincidentally fall in the same region of space. The same holds for pores in
the gel matrix.

Figure 3: Homogeneous fibre model, obtained as a Boolean model of randomly
oriented rectangles. The two realisations are obtained with length and
diameter $L_{F}=20$ $\mu$m, $D_{F}=1$ $\mu$m (a and $\blacktriangle$) and
$L_{F}=100$ $\mu$m, $D_{F}=10$ $\mu$m (b and $\blacksquare$), and the
calculated covariances are compared with that of the 2mg/mL gel imaged by CFM
(same as Fig. 2a1)

The homogeneous fibre model is a particular case of a Boolean model Serra
(1982); Jeulin (2021), for which the density and covariance are known
analytically. In particular, the fibre density is

$\phi_{F}=1-\exp\left[-\eta\right]$ (3)

where $\eta=\theta D_{F}L_{F}$ is the density one would expect in absence of
fibre overlap. The covariance is given by the following expression

$C_{FF}(r)=2\phi_{F}-1+(1-\phi_{F})^{2}\exp\left[\eta K_{F}(r)\right]$ (4)

where $K_{F}(r)$ is the geometrical covariogram of the randomly-oriented
fibres, which is calculated as

$K_{F}(r)=\frac{2}{\pi}\int_{0}^{\pi/2}\left[1-\frac{r}{D_{F}}\sin(t)\right]\left[1-\frac{r}{L_{F}}\cos(t)\right]\textrm{d}t$
(5)

for $r<D_{F}$ and the upper integration bound is replaced by
$\textrm{asin}(D_{F}/r)$ for $r\geq D_{F}$.

The covariance of the homogeneous fibre model is plotted in Fig. 3 for two
fibres sizes. For the purpose of illustration, the model is compared with the
experimental covariance of the segmented image of the 2mg/mL gel (from CFM).
In the figure, the number of fibres $\theta$ is chosen to achieve a density
$\phi_{F}\simeq 0.54$ comparable to the segmented image. As a consequence, the
asymptotic values of $C_{FF}(r)$ are a close match to the gel covariance both
for $r=0$ and for $r\to\infty$. The shape of the covariance at intermediate
distances, however, cannot be captured by the homogeneous model. The model can
account for either the small- or large-$r$ data but not the two
simultaneously. In the former case, the model captures the small-scale
structure of the gel (Fig. 3a) but it is unable to account for the aggregates,
which are found to be more prevalent than what can be expected from randomness
alone. In the latter case, the large-scale structure of the gel is reasonably
reproduced, but this is done by replacing fibre aggregates by unrealistically
large rectangles (Fig. 3b).

####  III.2.2 Fibre aggregates model

The inability of the homogeneous fibre model to reproduce the experimental
covariance of gels proves that structures larger than individual fibres are
more frequent than what can be accounted for by statistical fluctuations
alone. The existence of such large pores and aggregates is notably apparent
when comparing the realization of the homogeneous model in Fig. 3a with the
insets of Fig. 2.

To address this issue, we introduce a second model that builds on two
distinctly different structures. At the smallest scale the structure is
assumed to be that of homogeneously-distributed fibres, corresponding to
indicator function $\mathcal{I}_{F}(\mathbf{x})$.The larger-scale structure,
however, is accounted for by creating the indicator function of the entire
structure $\mathcal{I}(\mathbf{x})$ through the following multiplication

$\mathcal{I}(\mathbf{x})=\mathcal{I}_{F}(\mathbf{x})\times\mathcal{I}_{A}(\mathbf{x})$
(6)

where $\mathcal{I}_{A}(\mathbf{x})$ is the indicator function of the
aggregates, equal to one if $\mathbf{x}$ is inside an aggregate.
Mathematically, Eq. (6) is equivalent to starting with a homogeneous fibre
structure and subsequently carving pores out of it, using
$\mathcal{I}_{A}(\mathbf{x})$ as a mathematical cookie-cutter.

Independently of the specific models chosen for $\mathcal{I}_{F}(\mathbf{x})$
and $\mathcal{I}_{A}(\mathbf{x})$, evaluating the average value of Eq. (6)
yields the following density for the two-scale structure

$\phi_{1}=\phi_{F}\phi_{A}$ (7)

where $\phi_{A}$ is the density of the aggregates, and $\phi_{F}$ is the
density of the fibres within the aggregates. Equation (7) results from the
general definition of the density in Eq. (1), with the assumption that the
fibre and aggregate models are statistically independent from one another. The
same assumption provides the following relation for the covariance of the
solid phase

$C_{11}(r)=C_{FF}(r)C_{AA}(r)$ (8)

as a consequence of Eq. (2).

For the small-scale structure, we assume the same fibre model as considered
earlier, with density $\phi_{F}$ and covariance $C_{FF}(r)$ given in Eqs. (3)
and (4). As the aggregates are very disordered with no well-defined shape, we
model them with a clipped Gaussian-field approach Quiblier (1984); Berk
(1987); Gommes (2018) as described in the Supplementary Material. The two
parameters of the large-scale model are the density of the aggregates
$\phi_{A}$ and a single characteristic length $L_{A}$ that controls their
size.

Figure 4: Covariance of the fibre-aggregate model, with homogeneous fibre
model at small scale (a and $\diamond$) and clipped Gaussian-field model for
the large-scale aggregates (b and $\square$), combined to yield a two-scale
structure (c and $\blacktriangle$). The calculated covariance is compared with
that of the 2mg/mL gel imaged by fluorescence microscopy (same as Fig. 2a1)

As illustrated in Fig. 4 for the segmented image of the 2mg/mL gel in
fluorescence mode, the fibre-aggregate model captures well the experimental
covariance of the gels. For the fitting of the data, the parameter $L_{F}$ in
the fibre model is irrelevant as the actual length of the fibres is controlled
by the size of the aggregates. Equation 5 was therefore simplified to its
limit $L_{F}\to\infty$. The values of the remaining parameters are gathered in
Tab. 1 for the two gels and the two imaging modes. From the fitted parameters
$\phi_{F}$ and $\phi_{A}$ of the fibre-aggregate model, the total density of
the fibres $\phi_{1}$ was calculated through Eq. (7) and is also reported in
Tab. 1.

Table 1: Structural parameters of the gels, obtained from fitting the
covariance (binary images) or the correlation function and histogram (grey-
tone images) with the fibre-aggregate model.

Imaging | Processing | Gel |  $\phi_{F}$ (-) |  $\phi_{A}$ (-) |  $D_{F}$ ($\mu$m) |  $L_{A}$ ($\mu$m) |  $\phi_{1}$ (-)
---|---|---|---|---|---|---|---
CFM | Binary | 2 mg/mL | $0.76\pm 0.03$ | $0.72\pm 0.01$ | $1.9\pm 0.2$ | $8.2\pm 1$ | $0.55\pm 0.01$
3 mg/mL | $0.76\pm 0.01$ | $0.69\pm 0.01$ | $1.5\pm 0.1$ | $6.3\pm 0.4$ | $0.52\pm 0.01$
Grey-tone | 2 mg/mL | $0.73\pm 0.04$ | $0.32\pm 0.05$ | $1.9\pm 0.1$ | $6.6\pm 1$ | $0.23\pm 0.03$
3 mg/mL | $0.78\pm 0.06$ | $0.49\pm 0.04$ | $2.3\pm 0.1$ | $7.7\pm 1$ | $0.38\pm 0.06$
CRM | Binary | 2 mg/mL | $0.60\pm 0.02$ | $0.88\pm 0.04$ | $0.7\pm 0.09$ | $15\pm 1$ | $0.53\pm 0.04$
3 mg/mL | $0.67\pm 0.01$ | $0.90\pm 0.02$ | $0.90\pm 0.03$ | $13\pm 1$ | $0.61\pm 0.01$
Grey-tone | 2 mg/mL | $0.78\pm 0.05$ | $0.41\pm 0.03$ | $0.98\pm 0.2$ | $7.6\pm 0.5$ | $0.32\pm 0.01$
3 mg/mL | $0.61\pm 0.06$ | $0.56\pm 0.01$ | $0.62\pm 0.1$ | $14\pm 7$ | $0.34\pm 0.02$

$\phi_{F}$, $D_{F}$: density and diameter of the fibres; $\phi_{A}$, $L_{A}$:
density and size of the aggregates; $\phi_{1}$: total fibre density. The error
bars are the standard deviations observed from the fits of the fifteen images
in each condition.

Realizations of the fibre aggregate model for the two gels and the two imaging
modes are given in Fig. 5. These realizations illustrate the structures
captured by the covariance in the segmented images of the gels, and they
largely illustrate the drawbacks of a data-analysis based on image
segmentation. As mentioned earlier, the segmentation of the fluorescence data
leads to inconsistencies between the density of the segmented images and the
gel concentrations. In addition, the significant noise in the reflectance data
lead to unrealistically large values for the aggregate density (see Figs.
2b1-b2 and 5b1-b2).

Figure 5: Realizations of the fibre aggregate model, for the 2mg/mL (a1, a2)
and 3mg/mL (b1, b2) gels, based on the average values of the parameters fitted
from the covariance of the CFM (a1 and b1) and CRM (a2 and b2) segmented
images.

####  III.2.3 Greytone model

In order to overcome the limitations of image segmentation, we expand here the
binary fibre-aggregate model so as to allow the identification of its
parameters directly from grey-tone images. The binary models ignore the fact
that large image intensities are associated with large local fibre
concentrations. In the grey-tone model, we therefore assume that every fibre
contributes additively (by a quantity $\Delta$) to the local grey-tone of the
image. Namely, the image intensity where 2, 3 or more fibres overlap is
$2\Delta$, $3\Delta$, etc. The model also accounts explicitly for the
noisiness of the data, through uncorrelated Gaussian noise $n(\mathbf{x})$
with variance $\sigma^{2}$, and for a background intensity $b$.

These assumptions can be formally written by expressing the local image
intensity at point $\mathbf{x}$ as

$I(\mathbf{x})=b+\mathcal{I}_{A}(\mathbf{x})\times\Delta\sum_{i}\mathcal{I}_{F}^{(i)}(\mathbf{x})+n(\mathbf{x})$
(9)

where $\mathcal{I}_{A}(\mathbf{x})$ is the indicator function of the
aggregates as before, and the sum accounts for the overlapping of the fibres.
In the sum, each term $\mathcal{I}_{F}^{(i)}(\mathbf{x})$ is the indicator
function of a homogeneous fibre model with vanishingly small number density
$\eta^{(i)}$ so as to the overlapping of fibres. The fibre density $\phi_{F}$
is then obtained again through Eq. (3)., with the total number density
$\eta=\sum_{i}\eta^{(i)}$.

The statistical independence of all contributions in Eq. (9) -namely,
$\mathcal{I}_{A}(\mathbf{x})$, $n(\mathbf{x})$ and
$\mathcal{I}_{F}^{(i)}(\mathbf{x})$ for all $i$’s- enables one to calculate
the corresponding grey-tone correlation function. This is done through to
application of Eq. (2) to the image intensity in Eq. (9), and the result is

$C(r)=\langle
I\rangle^{2}+\eta\Delta^{2}K_{F}(r)C_{AA}(r)+(\eta\Delta)^{2}\left(C_{AA}(r)-\phi_{A}^{2}\right)+\sigma^{2}\delta(r)$
(10)

where $K_{F}(r)$ is the geometrical covariogram of the fibres, $C_{AA}(r)$ is
the covariance of the aggregates, the last term is the correlation function of
uncorrelated Gaussian noise where $\delta(r)$ is equal to 1 if $r=0$ and to 0
otherwise, and

$\langle I\rangle=b+\eta\phi_{A}\Delta$ (11)

is the average intensity of the image.

Figure 6: Greytone fibre-aggregates model of the 2mg/mL type-I collagen gel,
imaged by fluorescence microscopy, with (a1) histogram of grey levels (bars:
experimental values; blue shade: model), (b) the correlation function (dots:
experimental value; solid line: model), as well as realizations of the fibres
(c), aggregates (d) and complete model (e). In (a2) the values inferred from
the segmented images in the pores (back bars) and fibre areas (white bars) are
compared with the model contributions to the background (dark shade) and
increasing number of overlapping fibres (bright shades).

In addition to correlation functions, the grey-tone modelling enables one to
extract structural information from the grey-tone distribution itself. Because
the fibres in the model are distributed according to a Boolean process, the
probability for exactly $k$ fibres to overlap at any point of space is a
Poisson variable, namely

$\textrm{Prob}\Big{\\{}k\textrm{ overlapping
fibres}\Big{\\}}=\frac{\eta^{k}\exp[-\eta]}{k!}$ (12)

which contains Eq. (3) as a particular case, when all probabilities for $k$
larger or equal to 1 are added. As each fibre contributes a quantity $\Delta$
to the local intensity, Eq. (12) can also be interpreted as the probability
for observing the intensity $k\times\Delta$ for any point inside an aggregate.
Accounting also for the background intensity $b$ and for the Gaussian noise,
the intensity distribution is therefore

$f(I)=(1-\phi_{A}\phi_{F})g_{\sigma}[I-b]+\phi_{A}\exp[-\eta]\sum_{k=1}^{\infty}\frac{\eta^{k}}{k!}g_{\sigma}[I-(b+k\Delta)]$
(13)

where

$g_{\sigma}[x]=\frac{1}{\sqrt{2\pi}\sigma}\exp[-\frac{x^{2}}{2\sigma^{2}}]$
(14)

is the centred Gaussian probability density. The quantity $f(I)\textrm{d}I$ is
the probability for a pixel to have intensity in the interval
$[I,I+\textrm{d}I]$. In Eq. (13) the first term accounts for grey-tone
distribution of the background, outside the aggregates or between fibres
within the aggregates. In the second term, the sum over $k$ is on the
increasing number of overlapping fibres within the aggregates.

Figure 7: Realizations of the grey-level fibre aggregate model, for the 2mg/mL
(a1, a2) and 3mg/mL (b1, b2) gels, based on the average values of the
parameters fitted from the fluorescence (a1 and b1) and reflectance (a2 and
b2) images.

Figure 6 illustrates the fitting of the microscopy images with the grey-tone
model. The model captures nicely both the distribution of grey levels (Fig.
6a1) and the correlation function (Fig. 6b). As an independent check, Fig. 6a2
compares the modelled distribution of grey tones with the values measured in
the regions of the image that were classified as pores or fibres in the
segmented images (see Fig. S1). The comparison should be considered with
caution, as it is the unsuitability of the segmentation that justified the
present grey-tone approach. The agreement of the two methods, however, seem
reasonable. The various contributions to the fibre intensity in Fig. 6a2 also
testify to the importance of fibre overlap for the image analysis. The values
of the fitted parameters on the two gels and two imaging modes are reported in
Tab. 1. Realizations obtained from the average values of the fitted parameters
are given in Fig. 7.

##  IV Time- and space-resolved analysis of the gel structure in the spheroid
model

Time-resolved reflectance images of the platelet-derived growth factor-BB
(PDGF-BB)-treated CAF spheroid and surrounding gel are shown in Fig. 8, for
times ranging from 30 min to 750 min. In order to characterize the
spatiotemporal dynamics of cell interaction with the surrounding collagen
matrix during cell migration, the system was imaged at lower resolution than
the acellular gels considered so far. At the scale of the spheroid, individual
fibres are not visible and the only structures detected in the gels are the
aggregates. The images testify to a progressive local densification of the gel
structure close to the spheroid and migrating cells but it is unclear whether
other structural characteristics of the gel evolve as well. Moreover, as the
gel heterogeneity seems to develop at the same scale as the spheroid, it is
difficult to ascertain whether the effect of the cells is through direct
contact with the gel, or if long-range effects are also operative.

Figure 8: Images of the CAF spheroid and surrounding gel at various times:
$t_{1}=30$ min, $t_{3}=210$ min, $t_{9}$=750 min. The images correspond to one
height $z$ cutting through the spheroid, with the cells shown in red and the
density of the gel coded from blue to yellow.

In order to quantitatively investigate the modification of the gel structure
in relation with cell migration, the gel area surrounding the spheroid was
decomposed into layers corresponding to increasing distances to the closest
cells (shown as different colors in Fig. 9 left). At each time step, the
layers were recalculated according to the updated position of the cells. The
shapes of the equal-distance layers become increasingly complex with time, as
a consequence of the irregularity of the cell invasion pattern. The histogram
of grey tones and the correlation function of the collagen were measured
within each layer at each time step (Fig. 9 middle and right).

Figure 9: Collagen structure analysis in the neighbourhood of the CAF
spheroid, at different times ($t_{1}=30$ min, $t_{3}=210$ min $t_{9}=750$ min)
and distances from the cells ($d_{1}<100\ \mu$m, $100\
\mu\textrm{m}<d_{2}<200\ \mu\textrm{m}$ and $200\ \mu\textrm{m}<d_{3}<300\
\mu\textrm{m}$). In the original images (left) the gel is shown in grey and
the cells in bright red, and the same colour code is used for the distances in
the grey-tone distributions (middle) and correlation functions (right). The
middle panel displays the distributions of grey-tone intensity $I$, with the
bars being the experimental distributions and the 90 % percentile shown as a
vertical line. The shaded areas are the grey-tone model, with the background
contribution (darker) and increasing number of fibre overlaps. In the right
panel the grey-tone correlation functions $C(r)$ are shown, with the dots
being the experimental values and the solid lines being the grey-tone model.

The densification of the fibres in the vicinity of the spheroid as well as of
the migrating cells is manifest in Fig. 9 through the progressive shifting of
the grey-tone distributions towards brighter values (middle panel), and
notably the 90 % percentile (vertical lines). A similar evolution is observed
in the correlation functions that progressively shift towards larger values
with time. The evolution is quite rapid for the collagen in direct contact
with the cells (d1) and much slower for larger distance (d3).

In order to analyse the structural modifications corresponding to the changes
in grey-tone distributions and correlation functions, all the data measured at
9 successive time steps and 5 different distances from the cells were fitted
simultaneously with the grey-tone model (see Fig. 9 for 3 times and 3
distances). Among the parameters of the grey-tone model, the densities and
sizes bear a structural meaning (see Tab. 1), but other parameters merely
characterize the imaging mode. This is notably the case for the fibre contrast
$\Delta$, the background intensity $b$, and the noise intensity $\sigma$ (see
Eqs. 10 and 13). For the fitting of the data, different values of the
structural parameters were allowed for each time and distance, but a unique
value of the imaging parameters was allowed for all images. This is justified
as the time series images of the gel and spheroid were measured through a
time-laps protocol with unchanged imaging parameters. The values of the
imaging parameters from the fit are $b\simeq 9$, $\Delta\simeq 15$ and
$\sigma\simeq 5$. The structural parameters are $D_{F}\simeq 2$ $\mu$m and
$L_{A}\simeq 19$ $\mu$m for all distances and times, and the volume fractions
$\phi_{A}$ and $\phi_{F}$ are as shown in Fig. 10a and 10b.

Figure 10: Fitted collagen structure for various times and distances to the
CAF cells, with aggregate density $\phi_{A}$ (a) and fibre density within the
aggregates $\phi_{F}$ (b). The grey surface is a guide to the eye. Specific
realizations are shown on the right panel to illustrate the significance of
the parameters. The colour code for distances is the same as in Fig. 9; the
non-colored symbols point to intermediate distances not shown in Fig. 9.

Fig. 10 testifies to qualitatively different structural changes in the gel in
direct contact with the CAF cells and distant from them. Close to the cells
(blue series), a large number of aggregates have already formed before the
first measurement time ($t<30$ min). It remains relatively constant afterward
with density $\phi_{A}\simeq 80$ %, which can also be described in terms of a
matrix porosity that remains constant around 20 %. In the meantime, the fibre
density within the aggregates increases steadily all throughout the
experiment. This contrasts with the situation far from the cells (yellow
series), where the aggregate density $\phi_{A}$ increases far more
progressively and the fibre density within the aggregates $\phi_{F}$ remains
altogether constant. These qualitatively different structural changes close
and far from the cells are illustrated with right panel of Fig. 10 for various
distances and times.

##  V Discussion

The structure of collagen gels is a highly hierarchical one Sun (2021).
Starting form molecular dimensions, the smallest structures are the
tropocollagen triple helices (diameter 1.5 nm), assembled into protofibrils
and fibrils (6 to 100 nm), woven into fibres (1 $\mu$m), which are eventually
connected to form the gel matrix. The matrix, however, is not homogeneous at a
scale larger than 10 $\mu$m, as visible in all images of Fig. 1. A central
structure of collagen at large scale is he presence of fibre aggregates, which
can also be thought of in terms of pores in the matrix. The latter are the
most salient structures in the gels at the scale of the cells (Fig. 8), and it
is therefore natural to enquire how they are modified by neighboring cells. We
addressed this question by developing first a novel image analysis methodology
to characterize the large-scale structure of collagen, by validating it on two
acellular gels, and finally applying it to characterise the collagen structure
surrounding a spheroid of CAF cells.

The characterisation of the disordered pore-and-aggregate structure of the
gels is challenging, and calls for a statistical description. The approach we
explored here is based on a mathematical modelling, which aims at
statistically capturing the main structural characteristics of the gels
through a small number of meaningful parameters. It is important to stress
that the goal of the models is not to reproduce a local gel structure, as one
would typically observe in a single image, but also to capture the structure
variability. This feature is essential to assign a single value to the
structural parameters based on many different images of one and a single gel.

In that spirit we first considered a homogeneous fibre model (Fig. 3). This
simple model was found to be incompatible the experimental covariances of the
collagen, which statistically confirmed the presence of fibre aggregates and
pores. We then proposed a two-scale model, where the small-scale structure is
described as a homogeneous model of random fibres, and the large-scale
aggregate structure is obtained by using a clipped Gaussian-field to carve
pores out of it Fig. 4). In this two-scale model, two parameters $\phi_{A}$
and $L_{A}$ describe the density and size of the aggregates, and two
parameters $\phi_{F}$ and $D_{F}$ describe the density and the diameter of the
fibres inside the aggregates. Analytical expressions were derived to fit
experimental characteristics of the gels measured from the images, and infer
the values of the mentioned structural parameters.

A classical image analysis method consists in segmenting first the images to
identify objects of interest and subsequently measure them. This approach was
proved unsuitable in the present context, as segmentation would result in the
loss of valuable structural information. The collagen-rich areas of the gels
are indeed characterized by broad grey-tone distributions (see Fig. S1b), with
large intensities being associated with the close proximity or overlapping of
many fibres. Image segmentation would amount to treating all these regions on
the same footing, and would therefore bias any measurement.

We developed a method to identify the model parameters directly from grey-tone
microscopy images, through the fitting of the covariance function and
histogram of intensities. The covariance -or correlation function, in the case
of grey-tone images- is a very informative structural descriptor Jiao  _et
al._ (2011); Gommes  _et al._ (2012) that can easily be measured on any image
through fast Fourier transforms. Covariance analysis offers a variety of
advantages. First, it provides a robust and fully objective statistical
description of the spatial distribution of the objects that make up an image.
This is particularly useful for structures as complex and disordered as
fibrillar collagen gels. Second, it reduces the impact and possible bias of
image preprocessing, which is absent altogether in the case of grey-tone
correlation functions. Finally, it also provides an efficient way to cope with
noisy images, as noise is typically defined as the non-correlated contribution
to an image.

Globally, the procedure based on grey-tone measurements (correlation functions
and intensity distributions) and on their modelling proved to be quite
consistent with respect to the two types of imaging modes considered in the
paper. The absolute values of the parameters obtained through reflectance and
fluorescence data differ slightly, but identical trends are detected when
comparing the 2mg/mL and 3mg/mL acellular gels (Tab. 1). The overall fibre
density $\phi_{1}$ is larger in the 3 mg/mL gel, as it should. Our analysis
shows that this is due largely to an increased density of aggregates
$\phi_{A}$ while the fibre density within the aggregates $\phi_{F}$ is almost
identical in the two gels. The characteristic size of the aggregates $L_{A}$
is larger in the 3mg/mL gel because the aggregates are more numerous and form
a larger connected structure. To further highlight the importance of using
grey-tone images, we also performed the same analysis on segmented images. In
that case, segmentation leads to the undesirable situation where the density
of the segmented images is opposite to the known collagen concentration of the
gels in the case of fluorescence data (Tab. 1 and Fig. 2a).

From a methodological point of view, the application of the grey-tone model to
analyze the structure of the 2mg/mL and 3 mg/mL gels in both reflectance and
fluorescence modes, validated the overall image analysis approach and enabled
us to apply it to the dynamic model of CAF spheroid invasion (Figs. 9 and 10).
Matrix remodelling has been widely documented in earlier works Balcioglu  _et
al._ (2016); Kopanska  _et al._ (2016); Piotrowski-Daspit  _et al._ (2017);
Winkler  _et al._ (2020); Mark  _et al._ (2020). Contractility of collagen-
embedded cells generates strains on the collagen fibres, leading to their
local densification Doyle  _et al._ (2015, 2021) and reorganization Kopanska
_et al._ (2016); Piotrowski-Daspit  _et al._ (2017); Chen  _et al._ (2019).
This process is especially evident for cells of mesenchymal origin such as
CAFs, which are characterized by high expression levels of the contractile
myosin protein and collagen-binding integrins Doyle  _et al._ (2021). In our
experimental model, CAF spheroids were treated with PDGF-BB, a growth factor
known to stimulate collagen gel contraction Gullberg  _et al._ (1990); Reyhani
_et al._ (2017). By pulling on collagen fibers, cells stiffen their
microenvironment and induce irreversible collagen deformations, thereby
altering locally the mechanical properties of the ECM. This remodeled
microenvironment in turn provides guidance cues for the nearby cells through
durotactic and/or topotactic migrations SenGupta  _et al._ (2021); Lo  _et
al._ (2000); DuChez  _et al._ (2019). Besides these guidance cues, high
concentrations of fibrillar collagen also promote the local invasion of cells
by inducing the formation of specialized cellular extensions termed
invadopodia implicated in the local proteolytic remodeling of the ECM Artym
_et al._ (2015).

Our use of a two-scale stochastic model to describe the collagen structure and
the identification of its parameters from grey-tone images, enables us to
describe the ECM remodelling in more detail than simply the collagen density.
The pristine state of the ECM in the CAF spheroid model is the one observed at
early times $t$ and large distances $d$ from the cells. The relevant values of
the aggregate and fibre densities are $\phi_{A}\simeq 0.6$ and $\phi_{F}\simeq
0.5$ based on Figs 10a and 10b. These values are reasonably consistent with
those of the 2mg/mL acellular gels in Tab. 1, and the slight differences can
be attributed to the lower resolution used in the time-resolved analysis of
the spheroid. Starting from that initial state, our space- and time-dependent
analysis of the collagen structure shows that the ECM densification happens
via two distinct mechanisms, namely: the increased density of fibre aggregates
$\phi_{A}$ and the increased fibre density $\phi_{F}$ within the aggregates.
The two mechanisms are both at work in the gel in close contact with the
cells. The aggregate density increases very fast, from about $\phi_{A}\simeq
60$% to 80% in less than 30 mins (Fig. 10a), similar to the cancer cell
spheroids investigated by Chen  _et al._ (2019) which exert a strong
contraction of the surrounding collagen immediately after embedding in the
matrix. The fibre densification also takes place inside the aggregates is a
much more progressive way, as it occurs over hours (Fig. 10b).

The global densification of the gel resulting from the two processes is quite
significant, as the overall fibre density in contact with the cells -
estimated as $\phi_{1}=\phi_{A}\times\phi_{F}$ \- approximately passes from
$\phi_{1}\simeq$ 30% to 50%. It is important to stress that this only concerns
the fraction of fibres that are visible in our experiments, which is
determined by the imaging technique (CRM vs CFM) and by the spatial resolution
of the imaging system ($\sim$0.59$\mu$m/pixel in our case). As smaller fibrils
and protofibrils are undetected at the considered scale, the apparent
densification does not contradict the fact that the total collagen
concentration has to remain locally constant. We also cannot exclude that a
fraction of the aggregate densification observed by CRM in our spheroid model
results from a CAF-mediated partial reorientation of the more vertical
collagen fibres relative to the imaging plane Jawerth  _et al._ (2010). In any
event, the analysis testifies to significant yet different remodelling of the
gel structure at two scales simultaneously.

Interestingly, the remodelling of the collagen is not limited to the regions
in direct contact with the cells. Significant structural changes are observed
also in regions as far as 300 $\mu$m away from the closest CAF cell, but the
phenomenology is distinctly different there. In those regions, the density of
fibre aggregates increases slowly, over hours, while the fibre density within
the aggregates remains constant. Similar long-range (up to 1300 $\mu$m)
mechanical signals generated by fibroblasts have been shown to propagate
through fibrillar collagen networks. Such mechanical cues can be sensed far
beyond the signal source by cells sharing the same substrate, which is
responsible for the ability of contractile fibroblasts to induce the migration
of macrophages towards the force source Pakshir  _et al._ (2019).

From a methodological point of view, it is important to stress that the
discrimination between the two types of densification processes is robust.
Considering Fig. 6a2, an increase of aggregate density $\phi_{A}$ with
constant fibre density $\phi_{F}$ would be manifest through a reduction of the
background contribution (dark blue) in favor of the fibre contribution (bright
blue) but the shapes of the two contributions would remain unchanged. By
contrast, any increase in the fibre density $\phi_{F}$ is accompanied by an
increasing number of pixels where the fibres overlap, which would necessarily
extend the grey-tone distribution towards higher intensities. With that in
mind, the observed shifting of the 90%-percentile of the image intensity in
Fig. 9 can only be interpreted as an increase of $\phi_{F}$, independently of
any evolution of the aggregate density $\phi_{A}$. As for the grey-tone
correlation function, its two asymptotic values for $r=0$ and $r\to\infty$
depend only on the total fibre density $\phi_{1}=\phi_{A}\times\phi_{F}$. The
combination of the grey-tone distribution and correlation function is
therefore key to discriminate the two different densification processes.

##  VI Conclusion

In this study, we developed an innovative image analysis approach to
investigate the remodelling of fibrillar collagen in a 3D spheroid model of
cellular invasion. Unlike existing work, most of which focus on fibre the
densification of the collagen network and on small-scale fibre reorientation,
we focused here on the structural modification of the collagen matrix at the
scale of a few microns, comparable with that of the cells. This was achieved
by first developing a novel image analysis method based on the stochastic
modelling of the acellular gel structure, and applying it afterwards to study
the space- and time-dependent reshaping of the collagen matrix by migrating
CAFs.

The analysis of acellular collagen gels (without embedded cells) investigated
by confocal microscopy in both reflectance and fluorescence modes, shows that
the structure of the gels is not homogeneous at the scale of about 10 microns.
The structure consists in regions with high fibre density separated by
depleted regions, which can be thought of as fibre aggregates and pores. In
order to mathematically describe this structure, we developed a two-scale
stochastic model with a clipped Gaussian-field model for the aggregates and
pores, and a homogeneous Boolean model to describe the fibre network within
the aggregates. We also developed a method to identify the model parameters
from the grey-tone distributions and correlation functions of the gel images.
The specificity of the method is that it applies to the unprocessed grey-tone
images, and it can therefore be used with noisy time-lapse reflectance images
of non-fluorescent collagen.

When applied to the collagen-embedded CAF spheroid images, the developed
method testifies that the invasion of PDGF-BB-treated CAFs is accompanied by
an overall densification of the collagen gel while the sizes of the pores and
aggregates, as well as that of the fibres, remains largely unchanged.
Interestingly, the densification occurs differently for the gel in direct
contact with the cells or far away from them. The gel in close contact with
the invading cells densifies through the rapid increase of the number of
aggregates, over less than 30 min, followed by the slow increase of fibre
density within the aggregates, over hours. By contrast, the densification
occurring in the gel located farther away from the cells occurs via the slow
increase of the aggregate density, while the density of fibres within the
aggregates remains constant.

At the present stage, one can only speculate on the biomechanical mechanisms
responsible for the two-scale densification. The very observation of two
distinct phenomenologies hint at diverse mechanisms, which presumably involve
both biochemical and mechanical effects.

## Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any
commercial or financial relationships that could be construed as a potential
conflict of interest.

## Author Contributions

AN and EM designed the biological experiments; CJG and SB designed the
mathematical analysis; IB prepared the biological samples; TL performed the
confocal microscopy imaging; CJG developed the mathematical models, analyzed
the data and prepared all figures; CJG, EM, SB and AN wrote the manuscript.

## Funding

CJG and EM are grateful to the Funds for Scientific Research (F.R.S.-FNRS,
Belgium) for Research Associate positions. This work was supported by FNRS-
Télévie grants 7.4589.16 and 7.6527.18.

## Acknowledgments

Technical support from the GIGA Imaging platform of the Université de Liège is
gratefully acknowledged.

## Supplemental Data

A Supplementary Material file is available with details about the method used
to segment the gel images into fibre and background pixels, and about the
Gaussian Random Field approach used to model the fibre aggregates.

## References

  * Hynes (2009) R. O. Hynes, Science 326, 1216 (2009).
  * Hynes and Naba (2012) R. O. Hynes and A. Naba, Cold Spring Harbor Perspect. Biol. 4 (2012).
  * Oskarsson (2013) T. Oskarsson, Breast 22, S66 (2013).
  * Frantz  _et al._ (2010) C. Frantz, K. M. Stewart, and V. M. Weaver, J. Cell Sci. 123, 4195 (2010).
  * Theocharis  _et al._ (2016) A. D. Theocharis, S. S. Skandalis, C. Gialeli, and N. K. Karamanos, Adv. Drug Delivery Rev. 97, 4 (2016).
  * Leitinger and Hohenester (2007) B. Leitinger and E. Hohenester, Matrix Biol. 26, 146 (2007).
  * Xian  _et al._ (2010) X. Xian, S. Gopal, and J. R. Couchman, Cell Tissue Res. 339, 31 (2010).
  * Shieh  _et al._ (2011) A. C. Shieh, H. A. Rozansky, B. Hinz, and M. A. Swartz, Cancer Res. 71, 790 (2011).
  * Doyle  _et al._ (2009) A. D. Doyle, F. W. Wang, K. Matsumoto, and K. M. Yamada, J. Cell. Biol. 184, 481 (2009).
  * Charras and Sahai (2014) G. Charras and E. Sahai, Nat. Rev. Mol. Cell Biol. 15, 813 (2014).
  * Seo  _et al._ (2020) B. R. Seo, X. Chen, L. Ling, Y. H. Song, A. A. Shimpi, S. Choi, J. Gonzalez, J. Sapudom, K. Wang, R. C. A. Eguiluz, D. Gourdon, V. B. Shenoy, and C. Fischbach, Proc. Natl. Acad. Sci. U.S.A. 117 (2020).
  * Yue (2014) B. Yue, J. Glaucoma 23, S20 (2014).
  * Mouw  _et al._ (2014) J. K. Mouw, G. Ou, and V. M. Weaver, Nat. Rev. Mol. Cell Biol. 15, 771 (2014).
  * Bourgot  _et al._ (2020) I. Bourgot, I. Primac, T. Louis, A. Noël, and E. Maquoi, Front. Oncol. 10 (2020).
  * Revell  _et al._ (2021) C. K. Revell, O. E. Jensen, T. Shearer, Y. Lu, D. F. Holmes, and K. E. Kadler, Matrix Biol. Plus 12 (2021).
  * Provenzano  _et al._ (2019) P. Provenzano, D. Inman, K. Eliceiri, and P. J. Keely, Oncogene 28, 4326 (2019).
  * Provenzano  _et al._ (2008) P. P. Provenzano, D. R. Inman, K. W. Eliceiri, J. G. Knittel, L. Yan, C. T. Rueden, J. G. White, and P. J. Keely, BMC Med. 6 (2008).
  * Buchmann  _et al._ (2021) B. Buchmann, L. Meixner, P. Fernandez, M. Raich, F. Hutterer, C. Scheel, and A. Bausch, “Dataset for ”mechanical plasticity of collagen directs branch elongation in human mammary gland organoids”,”  (2021).
  * Dickinson  _et al._ (1994) R. B. Dickinson, S. Guido, and R. T. Tranquillo, Ann. Biomed. Eng. 22, 342 (1994).
  * Provenzano  _et al._ (2006) P. P. Provenzano, K. W. Eliceiri, J. M. Campbell, D. R. Inman, J. G. White, and P. J. Keely, BMC Med. 4 (2006).
  * Ray and Provenzano (2021) A. Ray and P. P. Provenzano, Curr. Opin. Cell Biol. 72, 63 (2021).
  * Ray  _et al._ (2022) A. Ray, M. K. Callaway, N. J. Rodríguez-Merced, A. L. Crampton, M. Carlson, K. B. Emme, E. A. Ensminger, A. A. Kinne, J. H. Schrope, H. R. Rasmussen, H. Jiang, D. G. DeNardo, D. K. Wood, and P. P. Provenzano, JCI Insight 7 (2022).
  * Hanahan and Coussens (2012) D. Hanahan and L. M. Coussens, Cancer Cell 21, 309 (2012).
  * Kalluri (2016) R. Kalluri, Nat. Rev. Cancer 16, 582 (2016).
  * Biffi and Tuveson (2021) G. Biffi and D. A. Tuveson, Physiol. Rev. 101, 147 (2021).
  * Sahai  _et al._ (2020) E. Sahai, I. Astsaturov, E. Cukierman, D. G. DeNardo, M. Egeblad, R. M. Evans, D. Fearon, F. R. Greten, S. R. Hingorani, T. Hunter, R. O. Hynes, R. K. Jain, T. Janowitz, C. Jorgensen, A. C. Kimmelman, M. G. Kolonin, R. G. Maki, R. S. Powers, E. Puré, D. C. Ramirez, R. Scherz-Shouval, M. H. Sherman, S. Stewart, T. D. Tlsty, D. A. Tuveson, F. M. Watt, V. Weaver, A. T. Weeraratna, and Z. Werb, Nat. Rev. Cancer 20, 174 (2020).
  * Bell  _et al._ (1979) E. Bell, B. Ivarsson, and C. Merrill, Proc. Natl. Acad. Sci. U.S.A. 76, 1274 (1979).
  * Dallon and Ehrlich (2008) J. C. Dallon and H. P. Ehrlich, Wound Repair Regen. 16, 472 (2008).
  * Mikami  _et al._ (2016) Y. Mikami, H. Matsuzaki, H. Takeshima, K. Makita, Y. Yamauchi, and T. Nagase, J. Visualized Exp. (2016).
  * Steinwachs  _et al._ (2015) J. Steinwachs, C. Metzner, K. Skodzek, N. Lang, I. Thievessen, C. Mark, S. Münster, K. E. Aifantis, and B. Fabry, Nat. Methods 13, 171 (2015).
  * Hall  _et al._ (2013) M. S. Hall, R. Long, X. Feng, Y. Huang, C.-Y. Hui, and M. Wu, Exp. Cell Res. 319, 2396 (2013).
  * Pakshir  _et al._ (2019) P. Pakshir, M. Alizadehgiashi, B. Wong, N. M. Coelho, X. Chen, Z. Gong, V. B. Shenoy, C. A. McCulloch, and B. Hinz, Nat. Commun. 10 (2019).
  * Mark  _et al._ (2020) C. Mark, T. J. Grundy, P. L. Strissel, D. Böhringer, N. Grummel, R. Gerum, J. Steinwachs, C. C. Hack, M. W. Beckmann, M. Eckstein, R. Strick, G. M. O’Neill, and B. Fabry, Elife 9 (2020).
  * Frye  _et al._ (2018) M. Frye, A. Taddei, C. Dierkes, I. Martinez-Corral, M. Fielden, H. Ortsäter, J. Kazenwadel, D. P. Calado, P. Ostergaard, M. Salminen, L. He, N. L. Harvey, F. Kiefer, and T. Mäkinen, Nat. Commun. 9 (2018).
  * Bredfeldt  _et al._ (2014) J. S. Bredfeldt, Y. Liu, C. A. Pehlke, M. W. Conklin, J. M. Szulczewski, D. R. Inman, P. J. Keely, R. D. Nowak, T. R. MacKie, and K. W. Eliceiri, J. Biomed. Opt. 19 (2014).
  * Liu  _et al._ (2018) Z. Liu, L. Speroni, K. P. Quinn, C. Alonzo, D. Pouli, Y. Zhang, E. Stuntz, C. Sonnenschein, A. M. Soto, and I. Georgakoudi, Biomaterials 179, 96 (2018).
  * Boudaoud  _et al._ (2014) A. Boudaoud, A. Burian, D. Borowska-Wykret, M. Uyttewaal, R. Wrzalik, D. Kwiatkowska, and O. Hamant, Nat. Protoc. 9, 457 (2014).
  * Pijanka  _et al._ (2019) J. K. Pijanka, P. P. Markov, D. Midgett, N. G. Paterson, N. White, E. J. Blain, T. D. Nguyen, H. A. Quigley, and C. Boote, J. Biophotonics 12 (2019).
  * Altendorf  _et al._ (2012) H. Altendorf, E. Decenciere, D. Jeulin, P. D. S. Peixoto, A. Deniset-Besseau, E. Angelini, G. Mosser, and M. . Schanne-Klein, J. Microscopy 247, 161 (2012).
  * Del Amo  _et al._ (2018) C. Del Amo, V. Olivares, M. Cóndor, A. Blanco, J. Santolaria, J. Asín, C. Borau, and J. M. García-Aznar, J. Mech. Behav. Biomed. Mater. 83, 52 (2018).
  * Geraldo  _et al._ (2012) S. Geraldo, A. Simon, N. Elkhatib, D. Louvard, L. Fetler, and D. M. Vignjevic, Eur. J. Cell. Biol. 91, 930 (2012).
  * Primac  _et al._ (2019) I. Primac, E. Maquoi, S. Blacher, R. Heljasvaara, J. V. Deun, H. Y. Smeland, A. Canale, T. Louis, L. Stuhr, N. E. Sounni, D. Cataldo, T. Pihlajaniemi, C. Pequeux, O. D. Wever, D. Gullberg, and A. Noel, J. Clin. Invest. 129, 4609 (2019).
  * Christiansen  _et al._ (2000) D. L. Christiansen, E. K. Huang, and F. H. Silver, Matrix Biol. 19, 409 (2000).
  * Lai  _et al._ (2012) V. K. Lai, C. R. Frey, A. M. Kerandi, S. P. Lake, R. T. Tranquillo, and V. H. Barocas, Acta Biomater. 8, 4031 (2012).
  * Serra (1982) J. Serra,  _Image Analysis and Mathematical Morphology_ , Vol. 1 (Academic Press, London, 1982).
  * Lantuéjoul (2002) C. Lantuéjoul,  _Geostatistical Simulations_ (Springer, Berlin, 2002).
  * Torquato (2002) S. Torquato,  _Random Heterogeneous Materials_ (Springer, New York, 2002).
  * Jeulin (2021) D. Jeulin,  _Morphological Models of Random Structures_ (Springer, Cham, Switzerland, 2021).
  * Aubert and Jeulin (2000) A. Aubert and D. Jeulin, Patt. Rec. 33, 1083 (2000).
  * Jiao  _et al._ (2009) Y. Jiao, F. H. Stillinger, and S. Torquato, Proc. Natl. Acad. Sci. U.S.A. 106, 17634 (2009).
  * Gommes  _et al._ (2012) C. J. Gommes, Y. Jiao, and S. Torquato, Phys. Rev. E 85 (2012).
  * Jeulin (2000) D. Jeulin, Stat. Comput. 10, 121 (2000).
  * Gommes (2018) C. J. Gommes, Microp. Mesop. Mater. 257, 62 (2018).
  * Quiblier (1984) J. A. Quiblier, J. Coll. Interf. Sci. 98, 84 (1984).
  * Berk (1987) N. Berk, Phys. Rev. Lett. 58, 2718 (1987).
  * Sun (2021) B. Sun, Cell Rep. Phys. Sci. 2, 100515 (2021).
  * Jiao  _et al._ (2011) Y. Jiao, H. Berman, T. Kiehl, and S. Torquato, PLoS ONE 6 (2011).
  * Balcioglu  _et al._ (2016) H. E. Balcioglu, B. van de Water, and E. H. J. Danen, Sci. Rep. 6, 22580 (2016).
  * Kopanska  _et al._ (2016) K. S. Kopanska, Y. Alcheikh, R. Staneva, D. Vignjevic, and T. Betz, PLoS One 11, e0156442 (2016).
  * Piotrowski-Daspit  _et al._ (2017) A. S. Piotrowski-Daspit, B. A. Nerger, A. E. Wolf, S. Sundaresan, and C. M. Nelson, Biophys. J. 113, 702 (2017).
  * Winkler  _et al._ (2020) J. Winkler, A. Abisoye-Ogunniyan, K. J. Metcalf, and Z. Werb, Nat. Commun. 11, 5120 (2020).
  * Doyle  _et al._ (2015) A. D. Doyle, N. Carvajal, A. Jin, K. Matsumoto, and K. M. Yamada, Nat. Commun. 6 (2015).
  * Doyle  _et al._ (2021) A. D. Doyle, D. J. Sykora, G. G. Pacheco, M. L. Kutys, and K. M. Yamada, Dev. Cell 56, 826 (2021).
  * Chen  _et al._ (2019) Y.-Q. Chen, J.-C. Kuo, M.-T. Wei, Y.-C. Chen, M.-H. Yang, and A. Chiou, Acta Biomater. 84, 280 (2019).
  * Gullberg  _et al._ (1990) D. Gullberg, A. Tingström, A.-C. Thuresson, L. Olsson, L. Terracio, T. K. Borg, and K. Rubin, Exp. Cell Res. 186, 264 (1990).
  * Reyhani  _et al._ (2017) V. Reyhani, M. Tsioumpekou, T. van Wieringen, L. Rask, J. Lennartsson, and K. Rubin, Sci. Rep. 7, 8924 (2017).
  * SenGupta  _et al._ (2021) S. SenGupta, C. A. Parent, and J. E. Bear, Nat. Rev. Mol. Cell Biol. 22, 529 (2021).
  * Lo  _et al._ (2000) C.-M. Lo, H.-B. Wang, M. Dembo, and Y. li Wang, Biophys. J. 79, 144 (2000).
  * DuChez  _et al._ (2019) B. J. DuChez, A. D. Doyle, E. K. Dimitriadis, and K. M. Yamada, Biophys. J. 116, 670 (2019).
  * Artym  _et al._ (2015) V. V. Artym, S. Swatkoski, K. Matsumoto, C. B. Campbell, R. J. Petrie, E. K. Dimitriadis, X. Li, S. C. Mueller, T. H. Bugge, M. Gucek, and K. M. Yamada, J. Cell. Biol. 208, 331 (2015).
  * Jawerth  _et al._ (2010) L. M. Jawerth, S. Münster, D. A. Vader, B. Fabry, and D. A. Weitz, Biophys. J. 98, L1 (2010).


# Electromagnetic Induction Imaging: Signal Detection Based on Tuned-Dressed
Optical Magnetometry

Giuseppe Bevilacqua  Valerio Biancalana<EMAIL_ADDRESS>Dept. of
Information Engineering and Mathematics - DIISM, University of Siena – Via
Roma 56, 53100 Siena, Italy  Yordanka Dancheva  Dept. of Physical Sciences,
Earth and Environment - DSFTA, University of Siena – Via Roma 56, 53100 Siena,
Italy  currently with Aerospazoio Tecnologie, Rapolano (SI), Italy  Alessandro
Fregosi  Dept. of Physical Sciences, Earth and Environment - DSFTA, University
of Siena – Via Roma 56, 53100 Siena, Italy  CNR Istituto Nazionale di Ottica,
via Moruzzi 1, 56124 Pisa, Italy  Gaetano Napoli  Dipartimento di Matematica e
Fisica "Ennio De Giorgi" Università del Salento, Lecce, Italy  Antonio
Vigilante  Dept. of Physical Sciences, Earth and Environment - DSFTA,
University of Siena – Via Roma 56, 53100 Siena, Italy

(August 27, 2024)

###### Abstract

A recently introduced tuning-dressed scheme makes a Bell and Bloom
magnetometer suited to detect weak variations of a radio-frequency (RF)
magnetic field. We envisage the application of such innovative detection
scheme as an alternative (or rather as a complement) to RF atomic
magnetometers in electromagnetic-induction-imaging apparatuses.

##  I Introduction

In 2001, Griffiths griffiths_mst_01  proposed an imaging technique based on
inferring one or more of the three passive electromagnetic properties
(conductivity $\sigma$, permittivity $\epsilon$ and permeability $\mu$) to
produce images on the basis of the response to a position dependent
oscillating magnetic field. Several denominations are used to identify this
kind of methodologies, among which electromagnetic induction imaging (EII)
deans_apl_20 , electro-magnetic tomograhy (EMT) yu_el_93 , magnetic induction
tomography (MIT) Korjenevsky_pm_00 ; bevington_as_20 , and also mutual
inductance tomography (same acronym) soleimani_ieee_07 , the latter three
stressing the potential of the technique to provide 3D mapping. A review on
the subject was recently authored by Ma and Soleimani ma_mst_17 .

As a general feature, the technique uses an AC magnetic field to excite eddy
currents in the specimen. The secondary magnetic field generated by those
currents is then detected and analyzed.

Magnetometers can detect and record in-phase ($\epsilon$ dependent) and out-
of-phase ($\sigma$ dependent) response of specimens subjected to an excitation
radio-frequency field $B_{\mathrm{rf}}$. The response obtained while scanning
the position of a sample with respect to the field generator and detector
enables the registration of 2D maps. The possibility of varying the
$B_{\mathrm{rf}}$ frequency (and correspondingly the skin depth in the
specimen) enables 3D (tomographic) capabilities.

Optical magnetometers, in the so-called radio-frequency (RF) implementation
deans_apl_16 ; wickenbrock_apl_16  are excellent detectors of weak (variation
of) magnetic fields oscillating at a resonant frequency and hence constitute
favourite detectors for EII, particularly in the case of low-conductivity
and/or small size specimens.

This led recently to important steps toward high resolution EII of weakly
conductive materials. Appositely developed RF magnetometers with opportune
field-specimen arrangement, jensen_prr_19 ; marmugi_apl_19 , demonstrated a
sensitivity sufficient to detect and characterize sub S m-1 conductive
material in small size ($5\text{\,}\mathrm{m}\mathrm{l}$) specimens, despite
the operation in unshielded environment) deans_apl_20 .

Such specification level makes the technique suited to develop tools for
medical diagnostics scharfetter_ieee_03 , where the negligible invasivity of
EII represents a valuable attractive, but the low-conductivity and the need of
high-spatial resolution constitute a severe requirement.

This work considers the potential of a peculiar configuration of a Bell&Bloom
(BB) magnetometer as a detector of eddy currents induced in small and/or
weakly conductive specimens, which could be eventually used to produce maps in
2D or 3D scans.

While the basic principle of operation (detection of a secondary magnetic
field) is shared with the above mentioned research, the proposed sensor works
with a particular arrangement that does not require resonant conditions for
the time-dependent excitation field. The proposed arrangement can constitute
an alternative as well as a complement to the commonly used RF magnetometers.

The standard BB implementation of atomic magnetometers is based on optically
pumping atoms into a given Zeeman sublevel (typically with a maximum magnetic
number along the quantization axis defined by the pump light wavevector) that
is not stationary due to the presence of a transverse magnetic field.

In the usual picture of the BB operation, the macroscopic magnetization
precesses around the magnetic field direction and it is periodically
reinforced when its orientation is along the optical axis of the pump beam.
The periodic reinforcement is obtained by modulating the pump radiation –its
intensity, polarization or wavelength– synchronously with the precession. A
weak and unmodulated probe radiation interrogates the evolution of the atomic
state and produces an output signal whose dynamics is driven by the field
under measurement.

The presence of time dependent field can modify the spin dynamics also in non-
resonant conditions. As an example, in conventional BB magnetometry, slow-
varying magnetic field variations, oriented along the bias field direction,
are detected with a high sensitivity. Field variations along perpendicular
directions produce some signal, as well, but with a second-order response
biancalana_arx_21 .

Interesting dynamic responses of precessing spins are observed in a different
regime, namely in the case of an intense fast-varying magnetic field. In
particular, a strong magnetic field oscillating along the probe-beam axis
(i.e. transversely to the static one) at a frequency much above the resonance
may effectively freeze the atomic precession, according with a phenomenon
originally studied in the late Sixties haroche_prl_1970  and generally known
as magnetic dressing. .

In a recent work bevilacqua_prl_20 , we have demonstrated that the dressing
phenomenon can be deeply tuned by the application of a secondary (and much
weaker) field that oscillates at the same frequency –or at a low-order
harmonic– of the dressing one, perpendicularly to the latter. Among the
peculiarities of such tuned-dressed configuration, we pointed out that the
system shows a remarkable response to small variations of the amplitude and/or
of the phase of the tuning field. This work investigates with a proof-of-
principle experiment the potential of such a tuning-dressing apparatus to
detect tiny field variations due to eddy currents induced by the tuning field.

##  II Experimental setup

The experimental setup is built around a BB magnetometer making use of cesium
thermal vapor in a centimetric, gas buffered cell. The atomic vapour is
synchronously pumped by $D_{1}$ radiation (at
$\text{\,}\mathrm{mW}\text{/}{\mathrm{cm}}^{2}$ irradiance) and
polarimetrically probed by a weak (at $\text{\,}\mathrm{\SIUnitSymbolMicro
W}\text{/}{\mathrm{cm}}^{2}$ irradiance) $D_{2}$ radiation, co-propagating
with the pump one along the $x$ direction, further details can be found in
Ref. biancalana_apb_16 .

The atomic sample is merged in a constant field at $\mu$T level (let its
direction be $z$) oriented perpendicularly to the laser beams. This static
field $B_{0}$ is obtained by partially compensating the environmental magnetic
field. The task is accomplished by means of three large size (180 cm) mutually
orthogonal Helmholtz pairs. Additional adjustable quadrupoles let improve the
field homogeneity.

Figure 1: Magnetic field and laser beams geometry. Circularly polarized
$D_{1}$ pump radiation and linearly polarized $D_{2}$ probe radiation
propagate along $x$, which is the same direction of a strong oscillating
(dressing) field $B_{D}$. A weaker (tuning) field $B_{t}$ oscillates along $y$
synchronously with $B_{D}$; a weak static field $B_{0}$ is oriented along $z$.

The tuning-dressing effect bevilacqua_prl_20  occurs when two phase-related
oscillating fields ($B_{D}$ and $B_{t}$) are applied along $x$ to dress the
atoms, and along $y$ to tune the dressing effect, respectively. A solenoidal
coil ($11\text{\,}\mathrm{c}\mathrm{m}$ in length and
$4\text{\,}\mathrm{c}\mathrm{m}$ in diameter) surrounding the cell generates
$B_{D}$ and a small (millimetric) solenoid wound on a ferrite nucleus
generates $B_{t}$. The interaction geometry is sketched in Fig.1.

The $B_{D}$ and $B_{t}$ coils are supplied by two waveform generators (Agilent
33250A) phase-locked to each other. Series capacitors help to adapt the
impedances, and a linear amplifier can be used to enhance the $B_{t}$
amplitude. Each of these coils has $10\text{\,}\mathrm{O}\mathrm{h}\mathrm{m}$
series resistor to precisely monitor the current phases and amplitudes via a
$16\text{\,}\mathrm{b}\mathrm{i}\mathrm{t}$,
$500\text{\,}\mathrm{k}\mathrm{S}\mathrm{/}\mathrm{s}$ DAQ card (NI 6346).

The tuning-dressed phenomenon occurs and is well modeled and characterized for
a tuning-field frequency $\omega_{t}$ that is an integer multiple of the
dressing frequency $\omega_{D}$ ($\omega_{t}=p\omega_{D}$). It is worth
stressing that, in the considered application, switching among different $p$
values may enable fast and relevant variations of the skin depth, which in EII
applications may constitute an interesting feature.

The coil-specimen-sensor geometry is sketched in Fig.2. The used specimens are
one or more Aluminum disks of assigned thickness $d$, diverse radii $R$,
centered on the $y$ axis, and lying on a $xz$ plane between the $B_{t}$
generator and the atomic sensor. The experimental results reported in this
work are obtained with $a=$$40\text{\,}\mathrm{m}\mathrm{m}$,
$b=$$130\text{\,}\mathrm{m}\mathrm{m}$,
$d=$$30\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$, and $D$ ranging from
$24\text{\,}\mathrm{m}\mathrm{m}$ to $350\text{\,}\mathrm{m}\mathrm{m}$. In
this proof-of-concept experiment, these values are chosen on the basis of
fortuitous constraints of the available setup. In real imaging applications,
all of them could be selected in view of optimizing the detection performance,
in terms of sensitivity and/or spatial resolution.

Figure 2: Relative positions of specimen, sensor, and tuning/dressing field
sources. The specimen is a circular disk of radius R (thick grey line), which
is parallel to the $xz$ plane and is centered on the $y$ axis, at a quote
$y=b$ over the sensor $S$. A solenoidal coil $C_{D}$ surrounds $S$ and
produces a homogeneous field $B_{D}$ along the $x$ axis. The latter coincides
with the optical axis of the pump and probe lasers $D_{1}$ and $D_{2}$. A
small dipolar source $C_{t}$ located on the y axis at distance $a$ from the
specimen (and $a+b$ from the sensor) produces a tuning field $B_{t}$ which is
oriented along $y$ on the sensor, and excites tangential eddy currents in the
disk. Those currents modify both the amplitude and the phase of the tuning
field on $S$.

The eddy currents induced in the conductive disks modify the amplitude and the
phase of $B_{t}$ in the sensor location. As summarized in Sec.III, both these
parameters play a role in shifting the effective Larmor frequency, and this is
the key feature at the basis of the proposed detection technique.

##  III Effective Larmor Frequency

As shown and discussed in the Ref.bevilacqua_prl_20 , a magnetic field that
oscillates in the direction perpendicular to both the transverse static field
$B_{0}$ and the longitudinal dressing field $B_{D}$ at a frequency $p$ times
larger than that of the dressing field, modifies the effective Larmor
frequency according to

$\Omega_{L}=\gamma\left[B_{0}J_{0}(\xi)+B_{t}J_{p}(\xi)\sin(\phi)\right],$ (1)

if $p$ is an odd integer, and to

$\Omega_{L}=\gamma\sqrt{\left[B_{0}J_{0}(\xi)\right]^{2}+\left[B_{t}J_{p}(\xi)\cos(\phi)\right]^{2}},$
(2)

if it is an even one. Here, $\gamma$ is the atomic gyromagnetic factor
($\gamma\approx$$3.5\text{\,}\mathrm{Hz}\text{/}\mathrm{nT}$ for Cesium), and
the argument $\xi$ is set by the amplitude of the dressing field $B_{D}$ and
by its frequency $\omega_{D}$ according to

$\xi=\frac{\gamma B_{D}}{\omega_{D}}.$ (3)

The validity of the eqs.(1) and (2) requires that

$\xi\gg\frac{\gamma B_{t}}{\omega_{D}},\frac{\omega_{0}}{\omega_{D}}$

In the measurements shown in this paper, the tuning field oscillates at the
same frequency of $B_{D}$ (i.e. $p=1$), while both its amplitude and phase
vary either due to different waveform-generator settings or to the presence of
conductive samples placed in the proximity of the field source and of the
sensor. In particular assuming that $\phi=\varphi-\theta$, where $\varphi$ is
the relative phase of the two RF generators, and $\theta$ is the dephasing
caused by the eddy currents in the specimen, the effective Larmor frequency
can be expressed as

$f=\frac{\gamma}{2\pi}\left[J_{0}(\xi)B_{0}+J_{1}(\xi)B_{t}\sin{\phi}\right]=f_{0}+A\sin(\varphi-\theta),$
(4)

where $f_{0}$ is the dressed frequency in absence of the tuning field. Beside
causing the dephasing $\theta$, the eddy currents attenuate $B_{t}$, and this
shielding effect determines a reduction of the parameter $A$.

The eq.(4) lets define the conditions in which the system response to either
the attenuation or the dephasing is maximized:

$\left|\frac{\partial f}{\partial A}\right|=\left|\sin(\varphi-\theta)\right|$
(5)

is maximal for $\phi=\varphi-\theta=\pi/2$, while

$\left|\frac{\partial
f}{\partial\theta}\right|=\left|A\cos(\varphi-\theta)\right|$ (6)

is maximal when $J_{1}$ is maximum, that is at $\xi\approx 1.84$ and
$\phi=0,\pi$. These relations provide an indication of good experimental
working conditions.

##  IV Results

The $B_{t}$ variation due to the eddy currents induced in the specimen can be
evaluated analytically for large (infinite radius) disks dodd_jap_68 , or
numerically for finite radius ones. The variation is monotonically dependent
on the disk radius, and concerns both the phase and the amplitude of $B_{t}$.
This is verified experimentally as shown in Fig.3

Figure 3: The resonance frequency as a function of the different relative
phase $\varphi$, for various disk diameters whose size in mm is reported in
the legend. The vertical error bar ($\approx$3\text{\,}\mathrm{Hz}\text{/}$$)
is smaller than the symbols used in the graph. The solid lines represent the
best fit profile obtained with the target function eq(4), with $A$ and
$\theta$ as fitting parameters.

The atomic resonance is experimentally analyzed by scanning the frequency of
the pump-laser modulation around the effective Larmor frequency expressed by
the eq.(4). The resonance line-width is about 25 Hz, however the high S/N
ratio let a best fit targeted to a Lorentzian profile determine the peak
frequency with a sub-Hz accuracy. On the other hand, the ambient field
fluctuations (the BB magnetometer is operated in an unshielded environment,
with deactivated field-stabilization system biancalana_prapplml_19 ) are at nT
level, and this leads to a few Hz uncertainty in the resonant frequency
estimation.

The measurements shown here are performed at
$B_{0}=$$1.7\text{\,}\mathrm{\SIUnitSymbolMicro T}\text{/}$ corresponding to
an undressed Larmor frequency of $6\text{\,}\mathrm{kHz}\text{/}$. The
dressing field frequency is set at $15\text{\,}\mathrm{kHz}\text{/}$ or
$30\text{\,}\mathrm{kHz}\text{/}$, and its amplitude $B_{D}$ is set to get a
dressing parameter $\xi=1.84$, according to the indication provided by the
eq.(6). In these conditions, the dressed frequency with no tuning field, is
$f_{0}=J_{0}(1.84)\gamma B_{0}/2\pi\approx$
$1900\text{\,}\mathrm{Hz}\text{/}$.

When the tuning field is activated, the measured peak frequency has a
remarkable dependence on both the amplitude and the phase of $B_{t}$ at the
sensor location, according to eq.(4). Both $A$ and $\theta$ can be determined
by scanning the relative phase of the RF generators in the $[0,2\pi]$
interval.

The Fig.3 shows clearly that an increase in disk radius causes both a
reduction of the amplitude $A$ and a phase shift of the sinusoidal profiles.
It is also clear that larger variations of the effective Larmor frequency
occur when $\phi=\pm\pi/2$, while operating near the zero-crossings
($\phi=0,\pi$) improves the sensitivity to $\phi$ (i.e to $\theta$)
variations, consistently with the eqs.(5), (6).

Figure 4: Dependence of the $B_{t}$ amplitude on the Al disk radius inferred
from experimental measurements (dots) and evaluated from a simplified
numerical simulation (line, see text for more details). Both the cases of
$\omega_{D}=2\pi\cdot$$15\text{\,}\mathrm{k}\mathrm{H}\mathrm{z}$ and
$2\pi\cdot$$30\text{\,}\mathrm{k}\mathrm{H}\mathrm{z}$ are considered. The
thinner skin depth at $30\text{\,}\mathrm{k}\mathrm{H}\mathrm{z}$ produces a
faster decay to a lower level, when the disk radius increases, as expected.

A clearer visualization is facilitated by a best fit procedure targeted to the
function expressed in eq.4. Plots of the best-fit parameter $A$ versus the
disk radius are shown with dots in Fig.4. The solid lines in that figure are
results from a finite-elements numerical simulation. The latter is based on
the geometry of Fig.2, without the $C_{D}$ solenoid and considering $B_{t}$ at
the center of the sensor as produced by a point dipolar source and attenuated
by on axis Al disks.

This analysis has been performed with the set of Al disks described in Sec.II
and we report results obtained at both the tested tuning-dressing frequencies.
Selecting different values of the oscillating fields changes the skin depth
and hence modifies the shielding properties of the Al samples. As known, the
skin depth depends on the absolute permeability $\mu$ of the medium and on its
resistivity $\rho$, according to:

$\delta=\sqrt{\frac{2\rho}{p\omega_{D}\mu}}$ (7)

In the case of Al, this formula gives $\delta=$
$460\text{\,}\mathrm{\SIUnitSymbolMicro m}\text{/}$ and
$650\text{\,}\mathrm{\SIUnitSymbolMicro m}\text{/}$, for the two considered
frequencies, respectively. Both values are much larger than the thickness of
the foil from which the Aluminum disks ($30\text{\,}\mathrm{\SIUnitSymbolMicro
m}\text{/}$) are cropped.

The general trend of the experimental results reported in Fig.4 is
qualitatively reproduced by numerical simulation. In particular both the
experiment and the simulation show that the curve width and the asymptotic
values decrease at higher $\omega_{D}$ (smaller $\delta$). However, some
quantitative discrepancies emerge: double thickness Al disks are considered in
the simulation to produce the well matching profiles reported in Fig.4.
Indeed, several details are neglected in the simulation, such as the presence
of the solenoid $C_{D}$ and of an electric heater surrounding the Cs cell, and
the finite size of the $B_{t}$ source and of the sensor volume. Very likely
one of more of these factors are responsible for the mentioned discrepancies.
Similar issues are encountered in the phase estimations $\theta(r)$.

##  V Discussion

The proposed tuning-dressing technique, extends to all-optical atomic
magnetometric sensors applicability as detectors in EII setups. Interestingly,
the proposed arrangement could be easily implemented in pre-existing EII
apparatuses, and dual modality (RF and tuning-dressing) setups could be built,
with incremental complexity of the available instrumentation.

In arrangements where the pump and probe beams are not parallel, a dual-mode
operation would make possible to apply the RF field along different
orientation, in such a way to induce variously distributed eddy currents. In
fact, while in RF magnetometry the RF field is applied perpendicularly to the
pump beam, in the tuning-dressing BB arrangement $B_{t}$ is perpendicular to
the probe beam.

Compared to conventional RF apparatuses, the tuning-dressing implementation
besides not requiring that the detection field oscillates at a resonant
frequency. Moreover, fast selection of different $B_{t}$ frequencies is
feasible, via an appropriate selection of the harmonic parameter $p$. The
possibility of selecting arbitrary values for both $\omega$ and $p$ is
relevant for the skin depth dependence on $p\omega_{D}$ (eq. 7), which is of
interest to improve the tomographic (three-dimensional) potential of the EII
apparatuses.

Concerning the detection efficiency, the proposed methodology takes advantage
from the fact that in the tuning-dressed configuration the weak tuning field
borrows strength from the primary dressing field. The linear dependence on
$B_{t}$ occurring for odd $p$ (which can be maximized with appropriate choices
of $\xi$ and selecting $\phi=\pm\pi/2$) makes possible and efficient the use
of approaches as that described in ref.deans_apl_20 . Specifically it would be
possible to apply a vanishing $B_{t}$ obtained as a superposition of two
opposite contributions, in such way that the specimen causes a $B_{t}$
variation via an unbalance of the two applied terms.

Improved detection techniques could be developed on the basis of the harmonic
dependence on $\phi$. As an example, if $B_{t}$ oscillates at frequency
$\omega_{t}$ slightly different from $p\omega$, the signal will appear as a
low-frequency beating term oscillating at $\omega_{t}-p\omega$: a feature that
would enable the application of phase-sensitive detection techniques.

##  VI Conclusion

We have presented a proof of principle experiment proposing a new kind of
magnetometric detector of small or weakly conductive specimens. Specifically,
we have used a recently studied dressing-field configuration to make an all-
optical Bell-and-Bloom magnetometer suited to detect faint oscillating fields
as occurring in magnetic induction tomography setups. The experiment
demonstrated that the variation of a non-resonantly oscillating probe field
can be effectively measured to detect the presence of conductive objects. The
described approach can be regarded as non-resonant alternative to the well
known methodologies based on RF magnetometers. It constitutes indeed a
complementary approach to the problem and dual-mode (RF and tuning-dressing)
apparatuses can be envisaged, with interesting novel features and enhanced
flexibility.

## References

  * (1) H. Griffiths, “Magnetic induction tomography,” Measurement Science and Technology, vol. 12, pp. 1126–1131, jul 2001.
  * (2) C. Deans, L. Marmugi, and F. Renzoni, “Sub-Sm–1 electromagnetic induction imaging with an unshielded atomic magnetometer,” Applied Physics Letters, vol. 116, no. 13, p. 133501, 2020.
  * (3) Z. Z. Yu, A. T. Peyton, M. S. Beck, W. F. Conway, and L. A. Xu, “Imaging system based on electromagnetic tomography (EMT),” Electronics Letters, vol. 29, no. 7, pp. 625–626, 1993.
  * (4) A. Korjenevsky, V. Cherepenin, and S. Sapetsky, “Magnetic induction tomography: experimental realization,” Physiological Measurement, vol. 21, pp. 89–94, feb 2000.
  * (5) P. Bevington, R. Gartman, and W. Chalupczak, “Inductive imaging of the concealed defects with radio-frequency atomic magnetometers,” Applied Sciences, vol. 10, no. 19, 2020.
  * (6) M. Soleimani, W. R. B. Lionheart, and A. J. Peyton, “Image reconstruction for high-contrast conductivity imaging in mutual induction tomography for industrial applications,” IEEE Transactions on Instrumentation and Measurement, vol. 56, no. 5, pp. 2024–2032, 2007.
  * (7) L. Ma and M. Soleimani, “Magnetic induction tomography methods and applications: a review,” Measurement Science and Technology, vol. 28, p. 072001, jun 2017.
  * (8) C. Deans, L. Marmugi, S. Hussain, and F. Renzoni, “Electromagnetic induction imaging with a radio-frequency atomic magnetometer,” Applied Physics Letters, vol. 108, no. 10, p. 103503, 2016.
  * (9) A. Wickenbrock, N. Leefer, J. W. Blanchard, and D. Budker, “Eddy current imaging with an atomic radio-frequency magnetometer,” Applied Physics Letters, vol. 108, no. 18, p. 183507, 2016.
  * (10) K. Jensen, M. Zugenmaier, J. Arnbak, H. Stærkind, M. V. Balabas, and E. S. Polzik, “Detection of low-conductivity objects using eddy current measurements with an optical magnetometer,” Phys. Rev. Research, vol. 1, p. 033087, Nov 2019.
  * (11) L. Marmugi, C. Deans, and F. Renzoni, “Electromagnetic induction imaging with atomic magnetometers: Unlocking the low-conductivity regime,” Applied Physics Letters, vol. 115, no. 8, p. 083503, 2019.
  * (12) H. Scharfetter, R. Casanas, and J. Rosell, “Biological tissue characterization by magnetic induction spectroscopy (mis): requirements and limitations,” IEEE Transactions on Biomedical Engineering, vol. 50, no. 7, pp. 870–880, 2003.
  * (13) G. Bevilacqua, V. Biancalana, Y. Dancheva, A. Fregosi, and A. Vigilante, “Spin dynamic response to a time dependent field,” 2021.
  * (14) S. Haroche, C. Cohen-Tannoudji, C. Audoin, and J. P. Schermann, “Modified Zeeman hyperfine spectra observed in 1H and 87Rb ground states interacting with a nonresonant RF field,” Phys. Rev. Lett., vol. 24, pp. 861–864, 1970.
  * (15) G. Bevilacqua, V. Biancalana, A. Vigilante, T. Zanon-Willette, and E. Arimondo, “Harmonic fine tuning and triaxial spatial anisotropy of dressed atomic spins,” Phys. Rev. Lett., vol. 125, p. 093203, Aug 2020.
  * (16) G. Bevilacqua, V. Biancalana, P. Chessa, and Y. Dancheva, “Multichannel optical atomic magnetometer operating in unshielded environment,” Applied Physics B, vol. 122, no. 4, p. 103, 2016.
  * (17) C. V. Dodd and W. E. Deeds, “Analytical solutions to eddy-current probe-coil problems,” Journal of Applied Physics, vol. 39, no. 6, pp. 2829–2838, 1968.
  * (18) G. Bevilacqua, V. Biancalana, Y. Dancheva, and A. Vigilante, “Self-adaptive loop for external-disturbance reduction in a differential measurement setup,” Phys. Rev. Applied, vol. 11, p. 014029, Jan 2019.