# Combining Multifractal Analyses of Digital Mammograms and Infrared Thermograms to Assist in Early Breast Cancer Diagnosis

Evgeniya G Gerasimova-Chechkina, Brian Toner, Zach Marin, Benjamin

Audit, Stéphane G Roux, Françoise Argoul, André Khalil, Olga Gileva, Oleg

Naimark, Alain Arneodo

► To cite this version:

Evgeniya G Gerasimova-Chechkina, Brian Toner, Zach Marin, Benjamin Audit, Stéphane G Roux, et al.. Combining Multifractal Analyses of Digital Mammograms and Infrared Thermograms to Assist in Early Breast Cancer Diagnosis. Physics of cancer: Interdisciplinary problems and clinical applications (PC'16), Mar 2016, Tomsk, Russia. pp.020018, 10.1063/1.4960237. hal-01556060

HAL Id: hal-01556060

https://hal.science/hal-01556060

Submitted on 4 Jul 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Combining Multifractal Analyses of Digital Mammograms and Infrared Thermograms to Assist in Early Breast Cancer Diagnosis

E. Gerasimova-Chechkina<fnref n="1" />,ª), B. Toner<fnref n="2" />, Z. Marin<fnref n="2" />, B. Audit<fnref n="3" />, S. G. Roux<fnref n="3" />, F. Argoul<fnref n="3" />,<fnref n="4" />, A. Khalil<fnref n="2" />, O. Gileva<fnref n="5" />, O. Naimark<fnref n="1" />, and A. Arneodo<fnref n="3" />,<fnref n="4" />

**Abstract.** We used a 1D wavelet transform modulus maxima (WTMM) method to analyze the temporal fluctuations of breast skin temperature recorded with an infrared (IR) camera from a panel of patients with breast cancer. This study shows that the multifractal complexity of temperature fluctuations observed in healthy breasts, is lost in the region of the malignant tumor in cancerous breasts. Then, we applied the 2D WTMM method to analyze the spatial fluctuations of breast density in the X-ray mammograms of the same patients. Compared to the correlated roughness fluctuations observed in the healthy areas, some clear loss of correlations is detected in malignant tumor foci. These physiological and architectural changes in the environment of malignant tumors detected in both thermograms and mammograms open new perspectives in computer-aided multifractal methods to assist in early breast cancer diagnosis.

## INTRODUCTION

According to World Health Organization, breast cancer is one of the most common causes of cancer deaths among women worldwide. Clinical studies have demonstrated that patients' survivability increases if breast anomalies are detected at early stages [1]. Although X-ray mammography is currently the most effective way of early breast cancer detection, it remains difficult to interpret due to high tissue variability and to two-dimensional projection effects [2, 3]. The miss rate in mammography is increased in dense breasts of young women [4, 5]. Criticism of the use of screening mammography due to over-diagnosis made clinicians use additional modalities such as ultrasound, magnetic resonance imaging, IR imaging, etc. Since the original observation by R. Lawson [6] that skin temperature over a malignant tumor is higher than its neighborhood, possibly resulting from abnormal increase in metabolic activity and vascular circulation in the tissues beneath [7–10], IR thermography has been considered as a promising non-invasive screening method of breast cancer detection [11–13]. However, the suitability of IR imaging for routine screening has been highly questioned [14], because of insufficient sensitivity for detection of deep lesions and limited knowledge of the relationship between skin surface temperature distribution and thermal diseases. This growing wave of criticism of breast cancer screening techniques has been strengthened by a similar movement criticizing the increasing use of computer-aided detection (CAD) methods that were shown

<footnote n="1">¹Laboratory of Physical Foundation of Strength, Institute of Continuous Media Mechanics UrB RAS, Academician Korolev Str. 1, Perm, 614013 Russia</footnote>

<footnote n="2">²CompuMAINE Laboratory, Department of Mathematics and Statistics, University of Maine, Neville Hall 324, Maine, Orono, 04469-5752 USA</footnote>

<footnote n="3">³Université de Lyon, Ens de Lyon, Université Claude Bernard, CNRS (UMR 5672), Laboratoire de Physique, 46 Allée d'Italie, Lyon, F-69007 France</footnote>

<footnote n="4">⁴Université de Bordeaux, Laboratoire Ondes et Matière d'Aquitaine (UMR 5798), Cours de la Libération, TALENCE Cedex, 33405 351 France</footnote>

<footnote n="5">⁵Department of Therapeutic and Propedeutic Dentistry, Perm State Medical University named after ac. E.A. Wagner, Petropavlovskaya Str. 26, Perm, 614000 Russia</footnote>

<footnote>ª) Corresponding author: egerasimova@icmm.ru</footnote>
to decrease specificity [15], and to increase recall rates of healthy women [16–19] resulting in an increase in unnecessary stress in women. Note that most of these CAD methods [20–24] were designed for texture analysis and feature extraction with the prerequisite that background roughness fluctuations of normal breast mammograms are homogeneous and uncorrelated. In this paper, we review the results of recent and ongoing analyses of dynamic infrared thermograms and X-ray mammograms based on wavelet-based multifractal methodologies [25–27].

Application of dynamic IR imaging as a diagnostic tool is based on the detection of variations in temperature rhythms generated by the cardiogenic and vasomotor frequencies [28, 29]. Human cardiogenic frequency levels lie between 1 and 1.5 Hz and vasomotor frequencies between 0.1 and 0.2 Hz, with higher harmonic frequency rhythms. But, beyond intensity differences in these rhythms between normal and tumor tissues, there is much more to learn from dynamical IR imaging of breast cancer. Using a 1D wavelet-based multiscale analysis of the temperature intensity fluctuations about these physiologic perfusion oscillations, we propose to characterize the multifractal properties of these temperature fluctuations as a novel discriminating method that can be implemented in early screening procedures to identify women with a high risk of breast cancer [25, 26].

The main premise in our investigation of mammograms [27] was related to the suggestion that healthy dense and fatty tissues, indeed display persistent long-range correlations and anti-correlations in mammogram spatial fluctuations, respectively [30, 31]. Using a 2D wavelet-based multifractal methodology, we recently observed some clear loss of correlations in mammogram area surrounding the malignant tumor and which strongly correlate with regions of anomalous monofractal (as compared to healthy multifractal) temperature fluctuations previously detected in dynamic IR thermograms [27]. This loss of complexity in temperature time-series and mammographic subimages is likely to be the signature of physiological and architectural alterations of the environment of malignant breast tumors. Besides their potential clinical impact, these results open new perspectives in the investigation and understanding of the cellular mechanisms that underlie the loss of tissue homeostasis, cancer initiation and progression.

## THE WAVELET TRANSFORM MODULUS MAXIMA METHOD

The wavelet transform (WT) is a mathematical microscope [32–37] that is well suited for the analysis of complex non-stationary time-series, as those found in genomics [38–40] and physiological systems [41–44], thanks to its ability to filter out low-frequency trends in the analyzed signal. Its two-dimensional (2D) generalization has been exploited for image processing [31, 37, 45] including applications in cellular biology [46–51] and medicine [30, 31, 52, 53].

### The 1D Wavelet Transform (Time-Frequency Analysis)

The WT is a space-scale analysis which consists in expanding a signal in terms of wavelets that are constructed from a single function, the “analyzing wavelet” $\psi$, by means of translations and dilations. The WT of a function $\Sigma(t)$ is defined as [37]:

$$W_{\psi}[\Sigma](t_0, a) = \frac{1}{a} \int_{-\infty}^{+\infty} \Sigma(t) \psi\left(\frac{t-t_0}{a}\right) dt, \quad (1)$$

where $t_0$ and $a$ ($> 0$) are the time and scale (inverse of frequency) parameters, respectively. By choosing a wavelet $\psi$ whose $(n + 1)$ first moments are zero ($\int t^m \psi(t) dt = 0, 0 \le m \le n$), one makes the WT microscope blind to order-$n$ polynomial behavior, a prerequisite for multifractal fluctuations analysis [54–57].

### The 2D Wavelet Transform (Space-Scale Analysis)

With an appropriate choice of the analyzing wavelet, one can reformulate Canny’s multiscale edge detection [58] in terms of a 2D wavelet transform [37]. The general idea is to start smoothing the discrete image data by convolving it with a filter and then to compute the gradient of the smoothed image. Let us consider two wavelets that are, respectively, the partial derivatives with respect to $x$ and $y$ of a 2D-smoothing function $\phi(x, y)$:

$$\psi_1(x, y) = \partial\phi(x, y)/\partial x \quad \text{and} \quad \psi_2(x, y) = \partial\phi(x, y)/\partial y. \qquad (2)$$
The WT of I(x, y) with respect to $\psi_1$ and $\psi_2$ has two components and therefore can be expressed in a vectorial form [37]:

$$ \mathbf{W}_{\psi}[I](\mathbf{x}_0, a) = \begin{cases} W_{\psi_1}[I] = \frac{1}{a^2} \int f(\mathbf{x}) \psi_1\left(\frac{\mathbf{x}-\mathbf{x}_0}{a}\right) d^2x \\ W_{\psi_2}[I] = \frac{1}{a^2} \int f(\mathbf{x}) \psi_2\left(\frac{\mathbf{x}-\mathbf{x}_0}{a}\right) d^2x \end{cases} = \nabla\{\phi_{\mathbf{x}_0, a} * I\}, \quad (3) $$

from which we can extract the modulus and argument of WT:

$$ M_{\psi}[I](\mathbf{x}_0, a) = \left[ (W_{\psi_1}[I](\mathbf{x}_0, a))^2 + (W_{\psi_2}[I](\mathbf{x}_0, a))^2 \right]^{1/2} \quad (4) $$

and

$$ A_{\psi}[I](\mathbf{x}_0, a) = \operatorname{Arg}[W_{\psi_1}[I](\mathbf{x}_0, a) + iW_{\psi_2}[I](\mathbf{x}_0, a)]. \quad (5) $$

## The Wavelet Transform Modulus Maxima Method

The wavelet transform modulus maxima (WTMM) method consists in computing the WT skeleton defined, at each fixed scale a, by the local maxima L(a) of the WT modulus M(·, a). In 1D, these WTMM are disposed on curves connected across scales called maxima lines which defines the WT skeleton [54–56, 59, 60]. In 2D, these WTMM lie, at a given scale, on connected chains called maxima chains. When considering the points along these maxima chains where M(·, a) is locally maximum, we define the so-called WTMMM L(a). The WTMMM are then linked through scales to form the WT skeleton [31, 61–63]. Along these maxima lines the WTMM (resp. WTMMM) behave as $a^{h(t)}$ (resp. $a^{h(x)}$), where $h(t)$ (resp. $h(x)$) is the Hölder exponent characterizing the singularity of $\Sigma$ (resp. $I$) at time $t$ (resp. position $\mathbf{x}$). The multifractal formalism amounts to characterize the relative contributions of each Hölder exponent value via the estimate of the $D(h)$ singularity spectrum defined as the fractal dimension of the set of points $t$ (resp. $\mathbf{x}$) where $h(t)$ (resp. $h(\mathbf{x})$) = $h$. This spectrum can be obtained by investigating the scaling behavior of partition functions defined in terms of WT coefficients [31, 54–57]:

$$ Z(q, a) = \sum_{l \in L(a)} M(\cdot, a)^q \sim a^{\tau(q)}, \quad (6) $$

where $q \in \mathbb{R}$. Then from the scaling function $\tau(q)$, $D(h)$ is obtained by a Legendre transform:

$$ D(h) = \min_q [qh - \tau(q)]. \quad (7) $$

In practice, to avoid instabilities in the estimation of the singularity spectrum $D(h)$ through the Legendre transform, we instead compute the following expectation values [31, 54–57]:

$$ h(q, a) = \frac{\partial}{\partial q} \ln(Z(q, a)) = \sum_{l \in L(a)} \ln(M(\cdot, a)\hat{W}(q, l, a)), \quad (8) $$

$$ D(q, a) = q \frac{\partial}{\partial q} \ln(Z(q, a)) - Z(q, a) = \sum_{l \in L(a)} \hat{W}(q, l, a) \ln(\hat{W}(q, l, a)). \quad (9) $$

where $\hat{W}(q, l, a) = M(\cdot, a)/Z(q, a)$ are the equivalent of the Bolzmann weight in the analogy that links the multifractal formalism to thermodynamics [56]. Then from the slopes of $h(q, a)$ and $D(q, a)$ vs $\ln a$, we get $h(q)$ and $D(q)$ and therefore the $D(h)$ and, consequently, singularity spectrum as a curve parametrized by $q$.

## Monofractal versus Multifractal Functions

The $\tau(q)$ data are fitted by the quadratic approximation $\tau(q) = -c_0 + c_1 q - c_2 q^2/2$, where the coefficients $c_n > 0$. The corresponding singularity spectrum has a characteristic single humped shape $D(h) = c_0 - (h - c_1)^2/2c_2$, where
$c_0 = -\tau(0)$ is the fractal dimension of the support of singularities of $\Sigma$ (resp. $I$), $c_1$ is the value of $h$ that maximizes $D(h)$ and $c_2$, the so-called intermittency coefficient [64], characterizes the width of the $D(h)$ spectrum. Homogeneous (mono) fractal functions, i.e. functions with singularities of unique Hölder exponent $H$ are characterized by a linear $\tau(q)$ curve ($h = \partial\tau/\partial q = c_1 = H$). The $D(h)$ singularity then reduces to a single point ($c_2 = 0$), $D(h=H)=1$ (resp. 2) as the signature of no change in WTMM statistics across scales [31, 39, 54–57, 61–63, 65, 66]. A nonlinear $\tau(q)$ is the signature of nonhomogeneous functions exhibiting multifractal properties, meaning that the Hölder exponent $h(t)$ [54–56] (resp. $h(x)$ [31, 61–63]) is a fluctuating quantity that depends on $t$ (resp. $x$). Then $D(h)$ has a single humped shape ($c_2 \neq 0$) that traduces some change in WTMM statistics across scales [31, 55, 63, 65].

# RESULTS

## Experimental Protocols

In this work [25–27], we analyzed the thermograms of 33 patients (from 37 to 83 years old) of Perm Regional Oncological Dispensary with breast cancer (invasive ductal and/or lobular cancer) and 14 healthy volunteers (from 23 to 79 years old) [67]. Tumors (1.2 to 6.5 cm in size) were found at 1 to 12 cm depths. Both breasts of each patient and volunteer were imaged with an IR camera (FLIR SC5000) at 50 Hz frame rate during the 10 min immobile patient phase. Each image set contained 30000 image frames. The mammographic procedure involved two X-ray images of both compressed breasts of each patient, referred to as the medio-lateral oblique (MLO) and the cranio-caudal (CC) views. Spatial resolution of images was 200 pixels per 1 cm. For mammogram analysis, our database includes mammograms collected from 30 of the mentioned 33 patients with breast cancer.

## WTMM Multifractal Analysis of Breast Infrared Thermograms

We grouped single-pixel temperature time-series into $8 \times 8$ pixel² squares covering the entire breast [25, 26]. The reported results are averaged over 64 temperature time-series in these $8 \times 8$ subareas. We analyzed 1-pixel temperature time-series taken from $8 \times 8$ pixels²-squares covering patient entire breasts. As expected, these time-series generally fluctuate at a higher temperature when recorded in the tumor region of a cancerous breast (Fig. 1a) than in a symmetrically localized square on the contralateral unaffected breast (Fig. 1b), as well as on a healthy breast (Fig. 1c).

<figure><img src="image_2.png" /><figcaption>**FIGURE 1.** Multifractal analysis of cumulative IR temperature time-series in a $8 \times 8$ pixels² square in the cancerous right breast (●) and opposite left breast (°) of patient 20 and in the healthy right breast (Δ) of volunteer 14. (a–c) 1 min portions of pixel temperature time-series in the cancerous breast (top), contralateral unaffected breast (median) and healthy volunteer breast (bottom), (d) $\tau(q)$ vs $q$ estimated by linear regression fit of $\log_2 Z(q, a)$ vs $\log_2 a$ over a range of time-scales [0.3, 3] s, (e) $D(h)$ vs $h$. Adapted from Gerasimova et al. [26]</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>**FIGURE 2.** Differential multifractal signature of temperature temporal fluctuations in breasts with and without tumor. Comparative analysis of both breasts of 33 patients with breast cancer and 14 healthy volunteers. (a) Normalized histograms of $c_1$ values obtained in the cancerous breasts (red: $N = 4032$, 8×8 pixels² squares), the contralateral unaffected breasts (black: $N = 3606$) and the volunteer's healthy breasts (green: $N = 3185$); (b) normalized histograms of $c_2$ values. Adapted from Gerasimova et al. [26]</figcaption></figure>

As commonly done for noise signals [38, 55] and previously experienced when analyzing rainfall time-series [68], we applied the WTMM method to the cumulative of these temperature time-series using the second-order compactly supported version $\psi_{(3)}^{(2)}$ of the Mexican hat [26, 69] (hence the singularities with a possible negative Hölder exponent $-1 < h < 0$, become singularities with $0 < h_{\text{cum}} = h+1 < 1$ in the cumulative). We found that the partition functions $Z(q, a)$ (Eq. (6)) display nice scaling properties for $q = -1$ to 5, over a range of time-scales limited to [0.43 s, 2.30 s] for linear regression fit estimates in a logarithmic representation. The $\tau(q)$ so-obtained are well approximated by quadratic spectra (Fig. 1d). For the cancerous breast of patient 20, $\tau(q)$ is nearly linear as quantified by a very small value of the intermittency coefficient $c_2 = (4.4 \pm 0.6) \cdot 10^{-3}$. This signature of monofractality is confirmed, when respectively plotting $h(q, a)$ (Eq. (8)) and $D(q, a)$ (Eq. (9)) vs $\log_2 a$ (see Fig. 2 in Ref. [26]) where the slopes $h(q) = c_1 = 0.81 \pm 0.01$ and $D(q) = 0.99 \pm 0.03$, do not significantly depend on $q$, meaning that the $D(h)$ singularity spectrum nearly reduces to a single point $D(h = c_1 = 0.81) = 1$ (Fig. 1e). In contrast, the $\tau(q)$ spectra obtained for the contralateral unaffected breast of patient 20 and the right healthy breast of volunteer 14 (Fig. 1d) are definitely nonlinear with a no longer negligible quadratic term $c_2 = 0.080 \pm 0.001$ and $c_2 = 0.069 \pm 0.002$, respectively. This hallmark of multifractal scaling is confirmed by the computation of the corresponding $D(h)$ spectra in Fig. 1e that both display a single-humped shape [26].

The results of our comparative wavelet-based multifractal analysis of (cumulative) temperature fluctuations over the 33 cancerous breasts, the 32 contralateral unaffected breasts (no right breast for patient 6) and the 28 volunteer healthy breasts are shown in Fig. 2 [25, 26]. In Fig. 2a, the corresponding histograms of $c_1$ values extend over a rather wide range $0.6 \le c_1 \le 1.8$, but turn out to be quite similar with the mean values $\bar{c}_1 = 1.066 \pm 0.002$ (cancerous), $1.104 \pm 0.002$ (contralateral unaffected) and $1.103 \pm 0.002$ (healthy). In contrast, the intermittency parameter $c_2$ has a definite discriminatory power. The histogram for cancerous breasts ($\bar{c}_2 = 0.045 \pm 0.001$) is definitely shifted towards smaller values relative to the ones for contralateral unaffected ($\bar{c}_2 = 0.056 \pm 0.001$) and normal breasts of healthy volunteers ($\bar{c}_2 = 0.058 \pm 0.001$) (Fig. 2b). The small values of $c_2$ dominate in cancerous breasts compared to normal breasts, confirming that the cancerous breasts are enriched in squares where temperature fluctuations display significantly reduced multifractal properties. This justifies that we considered $c_2 = 0.03$ as the threshold below (resp. above) which a square was qualified as monofractal (resp. multifractal) (Figs. 3d–3f).

According to the proposed classification criteria, 8×8 pixel² squares covering breast thermograms were color coded in the following way: squares with $c_2 < 0.03$ were shown in black, squares with $c_2 \ge 0.03$ in white and squares without scaling in gray. A high percentage of squares in the cancerous breast displays monofractal temperature fluctuations with small intermittency coefficient values (49.7% in Fig. 3a), while only few of those squares are found in the contralateral unaffected breast (7.7% in Fig. 3b) and in the healthy volunteer breast (11% in Fig. 3c). Consistently, in the thermograms of the contralateral unaffected and volunteer's normal breasts, a majority
of squares (89.4% and 65%, respectively) manifests multifractal scaling with an intermittency coefficient $c_2 \ge 0.03$ (Figs. 3b and 3c). Note that 43.1% of the squares in the cancerous breast also display multifractal temperature fluctuations, as observed for normal breasts (Fig. 3a). These squares indeed cover regions of the breast that are far from the tumor area (left upper quadrant) mostly covered by monofractal square. Let us point out that for each breast, a small percentage (<20%) of (gray) squares was removed from our analysis because of lack of scaling.

A comparative statistical analysis of monofractal and multifractal square distributions in breast thermograms confirmed that the percentage of monofractal (black) squares covering cancerous breasts (26.8 ± 3.5%) is about twice higher than the ones covering contralateral unaffected breasts (13.1 ± 2.3%) and breasts of healthy volunteers (11.3 ± 2.2%). This excess is compensated by a smaller percentage of multifractal (white) squares in cancerous breasts (56.9 ± 4.4%) than in contralateral unaffected breasts (68.5 ± 3.8%) and normal breasts (69.4 ± 4.3%) in healthy volunteers (Table 1).

We have compared the percentages of monofractal squares in cancerous and contralateral breasts [26] and found that 25 patients have higher quantity of monofractal squares in the cancerous breast. 18 cancerous breasts have a percentage of monofractal squares greater than the mean value (26%). Among the other 15 cancerous breasts, 7 have a smaller percentage of monofractal squares but these squares are well localized in the tumor region. The remaining 8 cancerous breasts were detected as false negatives due to small percentage of monofractal squares and their distant location with regard to the tumor region. Four of them correspond to rather deep tumors (6 cm, 7 cm, 8 cm and 12 cm depth) in the breasts with massive adipose tissue that can explain the absence of significant changes in skin temperature dynamics. Importantly, in five of the contralateral unaffected breasts a high (>26%) percentage of monofractal squares was detected which could be an early sign of pathology and require thus further check-ups. As a control, we reproduced this comparative analysis on both normal breasts of the 14 healthy volunteers. Only 4 volunteers manifested high (>26%) percentage of monofractal squares, while in most of them this value did not exceed 10% [26].

<figure><img src="image_4.png" /><figcaption>FIGURE 3. Breast-wide multifractal analysis of skin temperature temporal fluctuations. As estimated from the τ(q) spectra computed with the WTMM method (Fig. 1d), the 8×8 pixels² squares covering the entire breast are coded (a-c) and represented as dots in the (c₁, c₂) plane (d-f). The coding has the following meaning: black: $c_2 < 0.03$, white: $c_2 \ge 0.03$ and gray: no scaling (see text). (a, d) Cancerous right breast of patient 20: the tumor is located in the upper outer quadrant, (b, e) contralateral unaffected left breast of patient 20, (c, f) healthy right breast of volunteer 14. Adapted from Gerasimova et al. [26]</figcaption></figure>
<table><caption>TABLE 1. Percentage of monofractal ($c_2 < 0.03$), multifractal ($c_2 \ge 0.03$) and no-scaling 8x8 pixels² squares in cancerous, contralateral unaffected and healthy breasts [26]</caption><thead><tr><th></th><th>Cancerous breasts</th><th>Opposite breasts</th><th>Healthy breasts</th></tr></thead><tbody><tr><td>Monofractal</td><td>26.8 ± 3.5</td><td>13.1 ± 2.3</td><td>11.3 ± 2.2</td></tr><tr><td>Multifractal</td><td>56.9 ± 4.4</td><td>68.5 ± 3.8</td><td>69.4 ± 4.3</td></tr><tr><td>No-scaling</td><td>16.3 ± 2.9</td><td>18.4 ± 2.8</td><td>20.0 ± 5.0</td></tr></tbody></table>

To evaluate the efficiency of the test based on the multifractal analysis of breast infrared thermograms, we computed sensitivity (Se) and specificity (Sp). The results reported in Table 2 for our rather small set of patients, namely 76% sensitivity and 86% specificity, are quite encouraging for future applications to larger cohorts of females with breast cancer.

## WTMM Multifractal Analysis of X-ray Mammograms

We divided each mammogram view into 360x360 pixel² squares spanning 18x18 mm² [27]. From the central 256x256 pixels² part of the WT skeleton, we computed the partition functions $Z(q, a)$ (Eq. (6)) that were found to display scaling properties over a range of scales $1 < \log_2 a < 3$, corresponding to [0.7 mm, 2.8 mm], for linear regression fit estimates in a logarithmic representation. Within the range $-1 < q < 3$, statistical convergence was achieved and the so-obtained $\tau(q)$-spectrum was found remarkably linear (Fig. 4e) for a great majority of squares in both the cancerous and contralateral unaffected breasts, as the signature of monofractal roughness fluctuations [31, 61]. We systematically run our 2D WTMM algorithm for all mammographic views in our database. We confirmed the result of our preliminary study [30] that normal tissue regions display monofractal scaling properties ($c_1 = H, c_2 = 0$) as characterized by the Hurst exponent values $H < 0.5$ in fatty areas which look like anti-persistent self-similar random surfaces, while $H > 0.5$ in dense areas which exhibit persistent long-range correlations in mammographic roughness fluctuations [31, 61]. Importantly, our multiscale method provides an objective way to identify areas with some clear loss of correlations in roughness fluctuations ($H \approx 0.5$). We extended our wavelet-based multifractal analysis of CC and MLO mammographic images to the entire sets of 256x256 pixels² squares that cover each breast of the 30 patients. As illustrated in Figs. 4a–4d for patient 20, we color-coded each of these squares according to the monofractal diagnosis obtained with the 2D WTMM method. Besides confirming the existence of areas of correlated spatial roughness fluctuations, our methodology reveals that peculiar areas of uncorrelated roughness fluctuations are mainly found in the neighborhood of malignant tumors and strongly correlate with the regions of anomalous monofractal (as compared to healthy multifractal) temperature fluctuations previously detected in dynamic IR thermograms (Fig. 3a).

We computed the percentages of monofractal uncorrelated $H = 0.5$ squares in the mammograms of the two breasts of 30 patients with breast cancer. Good agreement was observed between the percentages of uncorrelated $H = 0.5$ squares in the CC and the corresponding MLO mammograms of both breasts of each patient (Fig. 5a). Obviously, the percentage of uncorrelated $H = 0.5$ squares in the cancerous mammograms is much larger than in the opposite breast mammograms. When comparing the two breasts of patients with breast cancer (Fig. 5b), we definitely evidenced some important asymmetry with a significant excess of $H = 0.5$ squares in the mammograms of the cancerous breast compared to the contralateral unaffected breast. This is a very promising indication that the loss of correlations in the mammographic spatial density fluctuations can be used as an indicator of malignancy. We correlated our findings in IR thermograms and X-ray mammograms of the 30 patients (Fig. 5c).

<table><caption>TABLE 2. Sensitivity Se = TP/(TP+FN) and specificity Sp = TN/(TN+FP) of the test based on multifractal analysis of breast infrared thermograms</caption><thead><tr><th></th><th>Has the disease</th><th>Does not have the disease</th></tr></thead><tbody><tr><td>Positive test score</td><td>True Positives (TP)<br>25</td><td>False Positives (FP)<br>4</td></tr><tr><td>Negative test score</td><td>False Negatives (FN)<br>8</td><td>True Negatives (TN)<br>24</td></tr></tbody></table>
<figure><img src="image_5.png" /><figcaption>**FIGURE 4.** Wavelet-based multifractal segmentation of X-ray mammograms. Patient 20: cancerous right breast (a, b) and contralateral unaffected left breast (c, d). (a, c) As estimated from the $\tau(q)$ spectra of CC mammographic view computed with the 2D WTMM method, 360×360 pixel$^2$ squares spanning 18×18 mm$^2$ were colored coded according to monofractal $H < 0.45$ (blue), $0.45 \le H \le 0.55$ (yellow), $H > 0.55$ (red) and no scaling (pink). (b, d) Same as (a, c) for MLO mammographic view. (e) $\tau(q)$ vs. $q$ estimated by linear regression fit of $\log_2 Z(q,a)$ vs. $\log_2 a$ over a range of space-scales [0.7, 2.8] mm for a yellow (●), a red (Δ) and a blue (□) squares. (f) $D(h)$ vs. $h$ for a yellow (●), a red (Δ) and a blue (□) squares. Adapted from Gerasimova et al. [27]</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**FIGURE 5.** Analysis of the percentage of “risky” monofractal uncorrelated $H = 0.5$ squares in mammograms and “risky” monofractal squares in thermograms of the breasts of 30 patients with breast cancer. (a) Percentage of $H = 0.5$ squares in cancerous (●) and contralateral unaffected (○) breasts: MLO view vs. CC view. (b) Mean percentage (MLO + CC)/2 of $H = 0.5$ squares in MLO and CC views: contralateral unaffected breast vs. cancerous breast. (c) Percentage of monofractal squares in thermograms vs. percentage of $H = 0.5$ squares in mammograms of the cancerous (▲) and contralateral unaffected (Δ) breasts of the 30 patients with breast cancer. Adapted from Gerasimova et al. [27]</figcaption></figure>
<table><caption>TABLE 3. Mean values of *H* = 0.5 square clusters characteristics</caption><thead><tr><th></th><th>Cancerous<br>breasts</th><th>Contralateral<br>unaffected breasts</th></tr></thead><tbody><tr><td>Size of a cluster</td><td>3.7</td><td>2.3</td></tr><tr><td>Size of the largest cluster</td><td>8.9</td><td>4.0</td></tr></tbody></table>

It is clear that the majority of patients have consistent large numbers of “risky” uncorrelated *H* = 0.5 squares in mammograms and of “risky” monofractal squares in thermograms of their cancerous breast as compared to their contralateral unaffected breast. We further investigated the spatial distribution of uncorrelated *H* = 0.5 squares over the entire breast defining a *H* = 0.5 cluster as an aggregation of *H* = 0.5 squares sharing a common edge [27]. Isolated *H* = 0.5 squares were considered as cluster-singlet. As reported in Table 3, the mean size of a cluster in a cancerous breast is larger than in the contralateral unaffected breast. This clustering effect is even more evident for the size of the largest *H* = 0.5 square that is twice bigger for the cancerous breasts than for the opposite ones. These findings suggest that uncorrelated *H* = 0.5 squares conglomerate in rather big clusters in cancerous breasts likely defining “risky” area that can be used as a new indicator of malignancy.

CONCLUSIONS

By combining a 1D WTMM multifractal analysis of IR temperature time-series and a 2D WTMM multifractal analysis of X-ray mammograms, we have revealed some change in both temporal and spatial scaling properties in the region of the malignant tumor in cancerous breasts. This is a clear indication of some physiological alterations and architectural disorganization of the “niche” surrounding breast tumor that may be the signature of a transition from an antitumorigenic to a tumorigenic microenvironment underlying the process of oncogenic transformation, tissue invasion and metastasis invasion during cancer progression [70–73]. Despite the exploratory nature of this preliminary study involving a rather small cohort of females who all went through surgery to remove the histologically confirmed malignant tumor, the obtained results are quite encouraging in the perspective of generalizing our approach to the analysis of longitudinal studies. This will allow us to quantify the statistical performance (sensitivity, specificity) of our computational approach and its potentiality to assist in early breast cancer diagnosis, progression and treatment.

ACKNOWLEDGMENTS

This work was supported by INSERM, ITMO Cancer for its financial support under contract PC201201-084862 “Physiques, mathématiques ou sciences de l’ingénieur appliqués au Cancer”, the Russian Foundation for Basic Research (grant No. 16-41-590235) and the Maine Cancer Foundation.

The study reported in this article was conducted according to accepted ethical guidelines involving research in humans and/or animals and was approved by an appropriate institution or national research organization.

The study is compliant with the ethical standards as currently outlined in the Declaration of Helsinki.

All individual participants discussed in this study, or for whom any identifying information or image has been presented, have freely given their informed written consent for such information and/or image to be included in the published article.

REFERENCES

1. C. H. Lee, Radiol. Clin. N. Am. **40**, 395–407 (2002).

2. R. E. Bird, T. Wallace, and B. Yankaskas, Radiology **184**, 613–617 (1992).

3. M. L. Giger, Acad. Radiol. **9**, 1–3 (2002).

4. K. J. Jorgensen and P. C. Gotzsche, BMJ **339**, b2587 (2008).

5. R. M. Tamini, G. A. Golditz, and S. E. Hankinson, Cancer Res. **69**, 9323 (2009).

6. R. N. Lawson, Can. Med. Assoc. J. **75**, 309–310 (1956).

7. M. Gautherie, Ann. N. Y. Acad. Sci. **335**, 383–415 (1980).
8. J. R. Keyserlingk, P. D. Ahlgren, E. Yu, N. Belliveau, and M. Yassa, Functional infrared imaging of the breast: Historical perspectives, current application and future considerations, in *Biomedical Engineering Handbook*, edited by J. D. Bronzino (CRC Press, 2006), pp. 97–158.

9. E. E. Sterns, B. Zee, S. SenGupta, and F. W. Saunders, Cancer **77**, 1324–1328 (1996).

10. T. Yahara, T. Koga, S. Yoshida, S. Nakagawa, H. Deguchi, and K. Shirouzou, Surg. Today. **33**, 243–248 (2008).

11. J. F. Head, F. Wang, C. Lipari, and R. L. Elliott, IEEE Eng. Med. Biol. **19**, 52–57 (2000).

12. W. C. Amalu, W. B. Hobbins, J. F. Head, and R. L. Elliott, Infrared Imaging of the Breast–An Over View, in *Biomedical Engineering Handbook*, edited by J. D. Bronzino (CRC Press, 2006), pp. 1–21.

13. E. Y. K. Ng, Int. J. Therm. Sci. **49**, 849–859 (2009).

14. J. F. Head and R. L. Elliot, IEEE Eng. Med. Biol. **21**, 80–85 (2002).

15. J. J. Fenton, L. Abraham, S. H. Taplin, B. M. Geller, P. A. Carney, C. D'Orsi, J. G. Elmore, and W. E. Barlow, J. Natl. Cancer Inst. **103**, 1152–1161 (2011).

16. J. J. Fenton, G. Xing, J. G. Elmore, H. Bang, S. L. Chen, K. K. Lindfors, and L. M. Baldwin, Ann. Int. Med. **158**, 580–587 (2013).

17. S. J. Kim, W. K. Moon, M. H. Seong, N. Cho, and J. M. Chang, Acta Radiol. **50**, 999–1004 (2009).

18. M. Noble, W. Bruening, S. Uhl, and K. Schoelles, Arch. Gynecol. Obstet. **279**, 881–890 (2009).

19. L. H. Eadie, P. Taylor, and A. P. Gibson, Eur. J. Radiol. **81**, e70–e76 (2012).

20. A. N. Karahaliou, I. S. Boniatis, S. G. Skiadopoulos, F. N. Sakellaropoulos, N. S. Arikidis, E. A. Likaki, G.S. Panayiotakis, and L. I. Costaridou, IEEE Trans. Inf. Technol. Biomed. **12**, 731–738 (2008).

21. T. Ayer, M. U. Ayvaci, Z. X. Liu, O. Alagoz, and E. S. Burnside, Imaging Med. **2**, 313–323 (2010).

22. N. C. Tsai, H. W. Chen, and S. L. Hsu, Comp. Med. Imaging and Graph. **35**, 1–8 (2011).

23. M. Meselhy Eltoukhy, I. Faye, and B. Belhaouari Samir, Comp. Biol. Med. **42**, 123–128 (2012).

24. L. Häberle, F. Wagner, P. A. Fasching, S. M. Jud, K. Heusinger, C. R. Loehberg, A. Hein, C. M. Bayer, C. C. Hack, M. P. Lux, K. Binder, M. Elter, C. Münzenmayer, R. Schulz-Wendtland, M. Meier-Meitinger, B. R. Adamietz, M. Uder, M. W. Beckmann, and T. Wittenberg, Breast Cancer Res. **14**, R59 (2012).

25. E. Gerasimova, B. Audit, S.-G. Roux, A. Khalil, F. Argoul, O. Naimark, and A. Arneodo, Europhys. Lett. **104**, 68001 (2013).

26. E. Gerasimova, B. Audit, S.-G. Roux, A. Khalil, O. Gileva, F. Argoul, O. Naimark, and A. Arneodo, Front. Physiol. **5**, 176 (2014).

27. E. Gerasimova Chechkina, B. Toner, Z. Marin, B. Audit, S.-G. Roux, F. Argoul, A. Khalil, O. Gileva, O. Naimark, and A. Arneodo, Cancer Res. (2016), submitted.

28. T. M. Button, H. Li, P. Fisher, R. Rosenblatt, K. Dulaimy, S. Li, B. O'Hea, M. Salvitti, V. Geronimo, S. Jambawalikar, and R. Weiss, Phys. Med. Biol. **49**, 3105–3116 (2004).

29. R. Joro, A.-L. Lääperi, S. Soimakallio, R. Järvenpää, T. Kuukasjä, T. Toivonen, R. Saaristo, and P. Dastida, J. Med. Eng. Technol. **32**, 325–335 (2008).

30. P. Kestener, J.M. Lina, P. Saint-Jean, and A. Arneodo, Image Anal. Stereol. **20**, 169–174 (2001).

31. A. Arneodo, N. Decoster, P. Kestener, and S.-G. Roux, Adv. Imaging Electr. Phys. **126**, 2–93 (2003).

32. A. Arneodo, G. Grasseau, and M. Holschneider, Phys. Rev. Lett. **61**, 2281–2284 (1988).

33. F. Argoul, A. Arneodo, J. Elezgaray, G. Grasseau, and R. Murenzi, Phys. Rev. A **41**, 5537–5560 (1990).

34. E. Freysz, B. Pouligny, F. Argoul, and A. Arneodo, Phys. Rev. Lett. **64**, 745–748 (1990).

35. A. Arneodo, F. Argoul, J. F. Muzy, M. Tabard, and E. Bacry, Fractals **1**, 629–649 (1993).

36. A. Khun, F. Argoul, J. F. Muzy, and A. Arneodo, Phys. Rev. Lett. **73**, 2998–3001 (1994).

37. S. Mallat, *A Wavelet Tour of Signal Processing* (Academic press, New York, 1999).

38. S. Nicolay, E. B. Brodie of Brodie, M. Touchon, B. Audit, Y. d'Aubenton-Carafa, C. Thermes, and A. Arneodo, Phys. Rev. E **75**, 032902 (2007).

39. A. Arneodo, C. Vaillant, B. Audit, F. Argoul, Y. d'Aubenton-Carafa, and C. Thermes, Phys. Rep. **498**, 45–188 (2011).

40. B. Audit, A. Baker, C. L. Chen, A. Rappailles, G. Guilbaud, H. Julienne, A. Goldar, Y. d'Aubenton-Carafa, O. Hyrien, C. Thermes, and A. Arneodo, Nat. Protoc. **8**, 98–110 (2013).

41. P. C. Ivanov, L. A. Amaral, A. L. Goldberger, S. Havlin, M. G. Rosenblum, Z. R. Struzik, and H. E. Stanley, Nature **399**, 461–465 (1999).

42. P. C. Ivanov, L. A. N. Amaral, A. L. Goldberger, S. Havlin, M. G. Rosenblum, H. E. Stanley, and Z. R. Struzik, Chaos **11**, 641–652 (2001).

43. A. L. Goldberger, L. A. Amaral, J. M. Hausdorff, P. C. Ivanov., C. K. Peng, and H. E. Stanley, Proc. Natl. Acad. Sci. USA **99**, 2466–2472 (2002).
44. C. D. Richard, A. Tanenbaum, B. Audit, A. Arneodo, A. Khalil, W. N. Frankel, and J. Neurosci, Methods **242**, 127–140 (2015).

45. J. P. Antoine, R. Murenzi, P. Vandergheynst, and S. T. Ali, *Two-Dimensional Wavelets and their Relatives* (Cambridge University Press, Cambridge, UK, 2008).

46. A. Khalil, J. L. Grant, L. B. Caddle, E. Atzema, K. D. Mills, and A. Arneodo, Chrom. Res. **15**, 899–916 (2007).

47. C. J. Snow, M. Goody, M. W. Kelly, E. C. Oster, R. Jones, A. Khalil, and C. A. Henry, PLoS Genet. **4**, e1000219 (2008).

48. M. F. Goody, M. W. Kelly, K. N. Lessard, A. Khalil, and C. A. Henry, Dev. Biol. **344**, 809–826 (2010).

49. J. Grant, C. Verrill, V. Coustham, A. Arneodo, F. Palladino, K. Monier, and A. Khalil, Chrom. Res. **18**, 873–885 (2010).

50. C. Martinez-Torres, L. Berguiga, L. Streppa, E. Boyer-Provera, L. Schaeffer, J. Elezgaray, A. Arneodo, and F. Argoul, J. Biomed. Opt. **19**, 036007 (2014).

51. C. Martinez-Torres, B. Laperrousaz, L. Berguiga, E. Boyer-Provera, J. Elezgaray, F. E. Nicolini, V. Maguer-Satta, A. Arneodo, and F. Argoul, J. Biomed. Opt. **20**, 096005 (2015).

52. A. Khalil, C. Aponte, R. Zhang, T. Davisson, I. Dickey, D. Engelman, and M. Mason, Med. Eng. Phys. **31**, 775–783 (2009).

53. K. A. Batchelder, A. B. Tanenbaum, S. Albert, L. Guimond, P. Kestener, A. Arneodo, and A. Khalil, PLoS One **9**, 107580 (2014).

54. J.-F. Muzy, E. Bacry, and A. Arneodo, Phys. Rev. Lett. **67**, 3515–3518 (1991).

55. J.-F. Muzy, E. Bacry, and A. Arneodo, Int. J. Bifurc. Chaos **4**, 242–302 (1994).

56. A. Arneodo, E. Bacry, and J.-F. Muzy, Physica A **213**, 232–275 (1995).

57. A. Arneodo, B. Audit, P. Kestener, and S.-G. Roux, Scholarpedia **3**, 4103 (2008).

58. J. Canny, IEEE Trans. Pattern Anal. Machine Intelligence **8**, 679–698 (1986).

59. J. F. Muzy, E. Bacry, and A. Arneodo, Phys. Rev. E **47**, 875–884 (1993).

60. E. Bacry, J. F. Muzy, and A. Arneodo, J. Stat. Phys. **70**, 635–674 (1993).

61. A. Arneodo, N. Decoster, and S.-G. Roux, Eur. Phys. J. B **15**, 567–600 (2000).

62. N. Decoster, S.-G. Roux, and A. Arneodo, Eur. Phys. J. B **15**, 739–764 (2000).

63. S.-G. Roux, A. Arneodo, and N. Decoster, Eur. Phys. J. B **15**, 765–786 (2000).

64. J. Delour, J. F. Muzy, and A. Arneodo, Eur. Phys. J. B **23**, 243–248 (2001).

65. A. Arneodo, B. Audit, N. Decoster, J.-F. Muzy, and C. Vaillant, A wavelet based multifractal formalism: Application to DNA sequences, satellite images of the cloud structure and stock market data, in *The Science of Disasters: Climate Disruptions, Heart Attacks, and Market Crashes*, edited by A. Bunde, J. Kropp, and H.J. Schellnhuber (Springer Verlag, Berlin, 2002), pp. 26–102.

66. B. Audit, E. Bacry, J. F. Muzy, and A. Arneodo, IEEE Trans. Info. Theory **48**, 2938–2954 (2002).

67. O. S. Gileva, G. G. Freynd, O. A. Orlov, T. V. Libik, E. I. Gerasimova, O. A. Plekhov, Y. V. Bayandin, and I. A. Panteleev, RFBR J. **74–75**, 93–99 (2012).

68. V. Venugopal, S. G. Roux, E. Foufoula-Georgiou, and A. Arneodo, Water Resour. Res. **42**, W06D14 (2006).

69. S. Roux, J.F. Muzy, and A. Arneodo, Eur. Phys. J. B **8**, 301–322 (1999).

70. M. J. Bissel, P. A. Kenny, and D. C. Radisky, Cold Spring Harb. Symp. Quant. Biol. **70**, 343–356 (2005).

71. M. J. Bissell and W. C. Hines, Nat. Med. **17**, 320–329 (2011).

72. V. Maguer-Satta, The stem cell niche: the black master of cancer, in *Cancer Stem Cells: Theories and Practice*, edited by S. Shostak (Tech, 2011), pp. 215–240.

73. P. Lu, V. M. Weaver, and Z. Werb, J. Cell Biol. **196**, 395–406 (2012).