The WT is a time-scale decomposition method which consists in expanding signals in terms of wavelets constructed from a single function, the “analyzing wavelet” ψ, by means of translations and dilations. The WT of a real-valued function Σ(t) is defined as (Mallat, 1998) (1)Wψ[Σ](t0,a)=1a∫−∞+∞Σ(t)ψ(t−t0a)dt , where t₀ is a time parameter and a (> 0) a scale parameter (inverse of frequency). By choosing a wavelet ψ whose n + 1 first moments are zero [∫t^mψ(t)dt = 0, 0 ≤ m ≤ n], one makes the WT microscope blind to order-n polynomial behavior, a prerequisite for multifractal fluctuations analysis (Muzy et al., 1991, 1994; Arneodo et al., 1995, 2008).

With an adapted analyzing wavelet, one can reformulate Canny's multiscale edge detection in terms of a 2D wavelet transform (Mallat, 1998). The underlying strategy consists in smoothing the discrete image data by convolving it with a filter prior to computing the gradient of the smoothed image. Let us define two wavelets as the partial derivatives with respect to x and y of a 2D-smoothing function ϕ(x, y) (Mallat, 1998; Arneodo et al., 2003):

The WT of I(x, y) with respect to ψ₁ and ψ₂ has two components and therefore can be expressed in a vectorial form: (3)Wψ[I](x0,a)={Wψ1[I]=1a2∫f(x)ψ1(x−x0a)d2xWψ2[I]=1a2∫f(x)ψ2(x−x0a)d2x},                             =∇{ϕx0,a∗I}, where x₀ is a space parameter and a(> 0) a scale parameter. From Equation (3), we can compute the modulus and argument of the WT: (4)Mψ[I](x0,a)=[(Wψ1[I](x0,a))2+(Wψ2[I](x0,a))2]1/2, and (5)Aψ[I](x0,a)=Arg[Wψ1[I](x0,a)+iWψ2[I](x0,a)].

The wavelet transform modulus maxima (WTMM) method was originally developed to generalize box-counting techniques (Arneodo et al., 1987) and to circumvent the limitations of the structure function method to perform multifractal analysis of 1D velocity signal in fully-developed turbulence (Muzy et al., 1991, 1993, 1994; Bacry et al., 1993; Arneodo et al., 1995). It proved to be a reliable method to estimate scaling exponents and multifractal spectra (Muzy et al., 1994; Delour et al., 2001; Audit et al., 2002). The WTMM method was generalized in 2D for the multifractal analysis of rough surfaces (Arneodo et al., 2000, 2003; Decoster et al., 2000; Roux et al., 2000), and for the analysis of 3D scalar and vector fields (Kestener and Arneodo, 2003, 2004; Arneodo et al., 2008). Successful applications have already been reported in various area of fundamental research (Muzy et al., 1994; Arneodo et al., 1995, 2002, 2008, 2011). In the context of present study, the 1D WTMM method demonstrated the multifractality in physiologic dynamics and its breakdown with disease (Ivanov et al., 1999, 2001; Goldberger et al., 2002), whereas the 2D WTMM method was used to detect microcalcifications and shown to have great potential to assist in cancer diagnosis from digitized mammograms (Kestener et al., 2001; Arneodo et al., 2003; Batchelder et al., 2014).

In 1D, the 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(t, a) = |W_ψ(t, a)|. These WTMM are positioned across scales on curves called maxima lines (Figure S1). In 2D, these WTMM lie, for a given scale, on connected chains called maxima chains (Figures S2A–C). Considering the points along these maxima chains where M(x, a) is locally maximum, we define the so-called WTMMM. The WTMMM are then linked through scales to form the WT skeleton (Figure S2D). 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 Σ (resp. I) at time t (resp. position x) (Muzy et al., 1994; Arneodo et al., 2003, 2008). The multifractal formalism amounts to quantify statistically the contributions of each Hölder exponent value via the computation of the singularity spectrum defined as the fractal dimension D(h) of the set of points t (resp. x) where h(t)(resp. h(x)) = h. This spectrum can be derived from the scaling behavior of partition functions defined in terms of WT coefficients (Muzy et al., 1991, 1994; Arneodo et al., 1995, 2008, 2003): (6)Z(q,a)=∑l∈L(a)M(·,a)q~aτ(q), where q ∈ ℝ. Then, from the scaling function τ(q), D(h) is obtained by a Legendre transform: (7)D(h)=minq[qh−τ(q)].

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 (Muzy et al., 1994; Arneodo et al., 1995, 2003) (8)h(q,a)=∂∂qln(Z(q,a))=∑l∈L(a)ln(M(·,a)W^(q,l,a)), and (9)D(q,a)=q∂∂qln(Z(q,a))−Z(q,a)                 =∑l∈L(a)W^(q,l,a)ln(W^(q,l,a))​, where W^(q,l,a)=M(·,a)qZ(q,a) corresponds to the Bolzmann weight in the analogy that connects the multifractal formalism to thermodynamics (Arneodo et al., 1995). Then, from the slopes of h(q, a) and D(q, a) vs. ln a, we get h(q) and D(q), and in turn the D(h) singularity spectrum as a curve parametrized by q.

Homogeneous monofractal functions are functions with singularities of unique Hölder exponent H. Their τ(q) spectrum is a linear curve of slope H. Monofractal scaling indeed means that the shape of the probability distribution function (pdf) of rescaled wavelet coefficients (M(·, a)/a^H) is independent of the scale a: ρMa/aH(w)=ρ(w), where ρ(w) is a “universal” pdf (Arneodo et al., 2002, 2011). A nonlinear τ(q) is the signature of nonhomogeneous multifractal functions with Hölder exponent h(t) (resp. h(x)) fluctuating over time t (Muzy et al., 1991, 1994; Arneodo et al., 1995, 2008) [resp. over space x (Arneodo et al., 2000, 2003; Decoster et al., 2000; Roux et al., 2000)]. τ(q) data are generally well approximated by the log-normal quadratic spectrum: (10)τ(q)=−c0+c1q−c2q2/2, where the coefficients c_n > 0. The corresponding singularity spectrum has a quadratic single-humped shape (11)D(h)=c0−(h−c1)2/2c2, where c₀ = −τ(0) is the fractal dimension of the support of singularities of Σ (resp. I), c₁ is the value of h that maximizes D(h) and the intermittency coefficient c₂ (Delour et al., 2001) characterizes the width of the D(h) spectrum as an indication of a change in WT coefficient statistics across scales (Muzy et al., 1991, 1994; Arneodo et al., 1995, 2002, 2003, 2008).

Subjects were recruited for the present study from Perm Region Oncological Dispensary using procedures approved by the Local Ethics Committe (Gileva et al., 2012). Patients gave informed consent to participate in this study via recording IR thermograms and X-ray mammograms of their two mammary glands: the cancerous breast and the contralateral unaffected breast. The thermography (resp. mammography) database included 33 (resp. 30 among these 33) females, aged 37–83 (average 57 years), who all went through surgery to remove a histologically-confirmed malignant tumor (invasive ductal and/or lobular cancer) a few weeks after thermo-(mammo-)grams were recorded. The tumors were found at depths varying from 1 to 12 cm into the tissue, with a size (diameter) varying from 1.2 cm up to 6.5 cm (Table S1).

Recording protocol was described in our preliminary study (Gerasimova et al., 2014). The patient's two breasts were imaged with an InSb photovoltaic (PV) detector (Joro et al., 2008). We performed imaging with the patient sitted with arms down to avoid too much discomfort. Frontal images were taken at a distance ~1 m from the patient, with an environmental room temperature ~20−22 °C. The image frame rate was 50 Hz. During the 10 min immobile imaging period, we collected 30,000 256 × 320 pixel² image frames at 14 bits on a computer connected to the PV camera. We used skin surface markers as reference points for low-frequency motion correction in the analysis.

The mammographic procedure involved taking two X-ray images of the (two) compressed breasts of each patient, the standard medio-lateral oblique (MLO) and cranio-caudal (CC) views. Spatial resolution of the images was 50 μm per pixel. The ionizing radiation dose was 0.4 mSv per patient.

Statistical analyses were performed using the R statistical package (https://cran.r-project.org/). A nonparametric Wilcoxon rank-sum test was used to calculate p-values when comparing CC and MLO mammograms, cancer breasts and contralateral unaffected breasts, and thermograms and mammograms. Linear correlation of data sets were analyzed via a Pearson test.

The numerical procedure to perform the WTMM analysis of 1D signals is avaliable at http://perso.ens-lyon.fr/benjamin.audit/LastWave

LastWave is an open source software written in C. We recommend interested users read the LastWave C-Application Programming Interface documentation and to contact the corresponding author to be directed to the part of the code of most relevance to them. For 2D WTMM analysis, we adapted home-made Xsmurf software written in C, for specific application to mammograms.

We analyzed individual 360 × 360 pixel² mammographic subimages of patient breasts with malignant tumor using the 2D WTMM method (Materials and Methods). From the central 256 × 256 pixels² part of the WT skeleton we computed the partition functions Z(q, a) (Equation 6) that were found to display nice scaling properties over a range of scales 1 ≲ log₂ a ≲ 3 corresponding to [0.7 mm, 2.8 mm] for linear regression fit estimates in a logarithmic representation (Figure 2A). The scaling deteriorates when considering larger scales due to finite size effects. In the range −1 ≲ q ≲ 3, statistical convergence is achieved; the so-obtained τ(q)-spectrum is remarkably linear, as the signature of monofractal roughness fluctuations and this wherever the location of the mammographic subimages over the cancerous breast (Figure 2B). Indeed the data are well-fitted by the theoretical spectrum τ(q) = qH − 2 of monofractal rough surfaces, which are almost everywhere singular with a unique hölder exponent h = H (Arneodo et al., 2000, 2003, 2008; Decoster et al., 2000). This is confirmed when computing the D(h) singularity spectrum from the slopes h(q) and D(q) of the partition functions h(q, a) (Equation 8) and D(q, a) (Equation 9) vs. ln a (Figure S3). h(q) and D(q) do not significantly depend on q, meaning that the D(h) singularity spectrum reduces to a single point D(h = H) = 2 (Figure 2C). Actually, what possibly changes when investigating different regions of cancerous breasts is the Hurst exponent H determined by the WTMM method (Figures 2B,C). As previously noted in a preliminary study (Kestener et al., 2001; Arneodo et al., 2003), we recovered H ≃ 1/3 in regions of fatty breast tissues as an indication of antipersistent (anti-correlated) roughness fluctuations, while dense tissues more likely displayed monofractal scaling with H ≃ 2/3 as an indication of persistent long-range correlations. But what was undetermined in previous analysis, was the existence of breast regions where the estimated Hurst exponent H ≃ 1/2, the hallmark of uncorrelated monofractal rough surfaces (Figures 2B,C). Importantly, the loss of correlations in roughness fluctuations is mainly observed in the region of the breast where the malignant tumor is located.

We extended our wavelet-based multifractal analysis of CC and MLO mammographic images to the entire sets of 256 × 256 pixels² squares that cover every breast of the 30 (among 33) patients that proceeded through X-ray mammography prior to surgery. As shown in Figures 1B,B′,C,C′ for patient 20, we color-coded each of these squares according to the parameter H estimated with the 2D WTMM method: (i) H < 0.45, blue (“fatty”) anti-correlated rough surface, (ii) H > 0.55, red (“dense”) long-range correlated rough surface, (iii) 0.45 ≤ H ≤ 0.55, yellow uncorrelated rough surface, and pink when no convincing scaling was observed in the considered square. For the two mammographic views of each breast, we calculated the numbers of blue (N_b), red (N_r), yellow (N_y) and pink (N_n) squares and corresponding percentages of breast coverage (Table S2). We first confirmed that, except in a minority of squares, monofractal scaling was clear. For the cancerous breasts, the mean percentage of no-scaling pink squares was 6.5% (resp. 5.7%) in CC (resp. MLO) mammograms, which was consistent with the mean percentage 5.9% (resp. 6.1%) in CC (resp. MLO) mammograms obtained for the (a priori healthy) contralateral unaffected breasts. Most of the patients had breasts consisting of mostly fatty tissues, with a majority of anti-correlated blue squares covering on average 70.9% (resp. 72.4%) of CC (resp. MLO) mammograms of the cancerous breasts and 79.1% (resp. 80.1%) in CC (resp. MLO) mammograms of the contralateral unaffected breasts. This excess was compensated by a rather small percentage of long-range correlated red squares 6.4% (resp. 6.5%) in CC (resp. MLO) mammograms of the cancerous breasts, as well as in the CC (resp. MLO) mammograms of the contralateral unaffected breasts [6.5% (resp. 6.5%)]. A similar agreement was observed between the percentages of uncorrelated yellow squares (further called H = 0.5 squares) in the CC and corresponding MLO mammograms of each breast of each patient (Figure 3A and Table S2). However, the number (Figures S4, S5A,A′) and percentage (Figures 3A, 4A,A′) of uncorrelated H = 0.5 squares in the cancerous breast mammograms were both larger than in the contralateral unaffected breast mammograms. In the CC (resp. MLO) view, N_y/N_Tot ranges in the interval [4.5%, 41.4%] (resp. [4.9%, 32%]) with a mean percentage 16.1% (resp. 15.4%) for the breasts with the malignant tumor, as compared to the interval [0.6%, 27%] (resp. [0%, 16.7%]) and mean percentage 8.5% (resp. 7.3%) for the opposite breasts. A statistical comparison between cancerous breasts and contralateral unaffected breasts was performed using a Wilcoxon rank-sum test which yielded p < 10⁻⁴ and 10⁻⁵ for the moduli of the values in Figure 3A and Figure S4A, respectively. Moreover, the moduli of values for both cancerous breasts and contralateral unaffected breasts in CC and MLO mammograms exhibited linear correlation. In Figures 4A,A′, cancerous breasts and contralateral unaffected breasts were different with p < 10⁻⁸. This statistically significant excess of H = 0.5 squares in the mammograms of cancerous breasts is evidence that the loss of correlations in the mammographic spatial density fluctuations can be used as an indicator of malignancy.

Asymmetry of patient's breast can be a sign of breast cancer. When comparing the percentages of uncorrelated H = 0.5 squares on both the cancerous and contralateral unaffected breast mammograms of each patient, we found that a majority (25/30) have more H = 0.5 squares (average over CC and MLO views) on the cancerous breast (Figure 3B and Figure S4B and Table S2), including patients 20 (Figure 1), 17 (Figure S6), 24 (Figure S7), and 29 (Figure S8). For the 5 patients left including patients 3 (Figure S9), 7 (Figure 5) and 30 (Figure S10), an important and slightly larger (or equal) percentage of H = 0.5 squares was found in the second breast as a probable indication of some physiological and architectural changes in the undiagnosed opposite breast. When looking at the excess of H = 0.5 squares in the cancerous breasts we realized that these uncorrelated squares were not sparsely distributed all over the breast, but were concentrated in specific regions likely surrounding the underlying malignant tumor (Figures 1B,C and Figures S6B,C, S7B,C, S8B,C). To quantify this inhomogeneous distribution, we investigated the size distribution of H = 0.5 clusters defined as H = 0.5 squares sharing a common edge or corner (Table 1). If cancerous and contralateral unaffected breasts display about the same number of isolated H = 0.5 squares (singlet), the former are comprised of clusters with sizes larger than 2 and up to 26 (patient 13). The corresponding histogram of H = 0.5 cluster size displays a long tail in cancerous breasts that has no counterpart in the contralateral unaffected breasts (Figure 4B). As seen on the cumulative histogram (Figure 4B′), more than 20% of H = 0.5 clusters have a size larger than 5 in the cancerous breast as compared to only 8.7% in the contralateral unaffected breasts. This clustering effect is even more pronounced on the histogram of the size of the largest H = 0.5 square cluster (Figure 4C), with 33.3% of the cancerous breasts with a largest H = 0.5 square cluster bigger than 10 as compared to only 3.4% for the contralateral unaffected breasts (Figure 4C′). Using Wilcoxon rank-sum test to compare cancerous breasts and contralateral unaffected breasts in Figures 4A–C yielded p < 10⁻⁸, 10⁻³, and 10⁻⁸, respectively. Note that the only contralateral unaffected breast with a H = 0.5 square cluster of size 22 is the second breast of patient 3 that was previously shown to have more H = 0.5 squares that the breast with the malignant tumor (Figure S9). Quite similar statistical results were obtained when defining H = 0.5 clusters from squares sharing a common edge only (Figures S5B,B′,C,C′ and Table S3).

For the same cohort of patients, the analysis of dynamic IR thermograms put into light some significant change in the nature of temporal fluctuations of breast skin temperature in cancerous breasts (Gerasimova et al., 2013, 2014). When using the 1D WTMM method, the computation of the partition functions Z(q, a) (Equation 6), h(q, a) (Equation 8) and D(q, a) (Equation 9) revealed that, besides the cardiogenic and vasomotor perfusion oscillations, the temperature dynamics display some scale invariant background component that extends from the characteristic human respiratory frequency (≳ 0.3 Hz) up to the cross-over frequency (≲ 4 Hz) toward instrumental white noise (Figure 2D). For normal breasts as well as in healthy regions of cancerous breasts, the τ(q) spectrum is definitely non-linear (Figure 2E) and accordingly the D(h) spectrum has a single humped shape (Figure 2F), the hallmark of multifractal scaling. By contrast, in the region of the malignant tumor in cancerous breasts the τ(q) spectrum is nearly linear and the D(h) spectrum reduces to a single point as the signature of monofractal scaling. As illustrated in Figures 1A,A′ for patient 20, we have segmented the two breasts of each patient into 8 × 8 pixels squares that were further color-coded according to the monofractal (red), multifractal (blue) or no-scaling (white) diagnosis obtained with the 1D WTMM method. In good agreement with our previous segmentation of CC (Figures 1B,B′) and MLO (Figures 1C,C′) mammograms, more “risky” monofractal squares (49.7% coverage) were found in the cancerous breast (Figure 1A) than in the contralateral unaffected breast (7.7%) (Figure 1A′). In addition to their number, the location and clustering of these “risky” squares, where the multifractal complexity of temperature fluctuations is lost, correlate with the clustering of the “risky” H = 0.5 squares in CC and MLO mammograms indicating some additional evidence of loss of correlations in spatial density fluctuations in the upper outer quadrant of the right breast of patient 20 where the underlying tumor is located.

The results obtained for patient 20 are quite representative of the outcome of our overall comparative analysis of IR thermograms and X-ray mammograms for our set of 30 patients (Figure 6, Figure S11). The use of Wilcoxon rank-sum test to compare cancerous breasts and contralateral unaffected breasts in Figure 6 and Figure S11 yielded p < 10⁻⁴ and 10⁻³, respectively. But both types of breasts did not exhibit linear correlation between mammogram and thermogram data. Most patients have consistent large numbers (breast coverage ≳ 10%) of “risky” uncorrelated H = 0.5 (yellow) squares in mammograms and monofractal (red) squares in thermograms of their cancerous breast as compared to their contralateral unaffected breast (breast coverage ≲ 10%), including patients 17 (Figure S6) and 29 (Figure S8). Note that for the subset of (5) patients, including patients 3 (Figure S9), 7 (Figure 5), and 30 (Figure S10), who were shown to have also a lot of “risky” H = 0.5 squares on their contralateral unaffected breast mammograms, they also have a lot of risky monofractal squares on the thermograms of their two breasts which might be an additional sign of the possible extension of cancer to their second breast. Among the patients with few “risky” monofractal squares (≲ 10%) on their IR thermograms of their cancerous breast, 4 correspond to rather deep tumors, namely patients 12 (size 1.8 cm, depth 12 cm), 16 (3.4 cm, 7 cm), 18 (3.49 cm, 6 cm) (Figure 7) and 28 (3.49 cm, 8 cm). For these 4 patients the percentage of “risky” H = 0.5 squares on their corresponding cancerous breast mammograms is significantly high (≳ 10%), meaning that, although imperceptible in temperature dynamics at the skin surface, the change in the microenvironment of these deep tumors turns out to be detectable with X-ray mammography. Let us also mention that if the Paget disease of the nipple of the left breast of patient 24 (Figure S7) was hardly identified in the IR thermograms with only few (breast coverage 4.9%), but well localized, “risky” monofractal squares (Figure S7A), it could not be missed on the CC (Figure S7B) and MLO (Figure S7C) mammograms, with respective 41.4% and 24% coverages by “risky” uncorrelated H = 0.5 squares.

As a control, we reproduced our 2D WTMM multifractal analysis on the mammograms of the two breasts of patient 30 that were recorded 11 months after surgical treatment of invasive ductal carcinoma of the right breast. For both mammographic views, the number of “risky” H = 0.5 squares in the cancerous breast was reduced from N_y = 8 (Figure 8A) to 3 (Figure 8C) in CC mammograms, and from N_y = 13 (Figure 8B) to 1 (Figure 8D) in MLO mammograms. However, less expected was the quite similar quantitative reduction observed in the other breast of patient 30, who was one of the 5 patients having slightly more H = 0.5 squares on the contralateral unaffected breast than on the cancerous breast (Figure S10): N_y = 12 (Figure 8A′) to 0 (Figure 8C′) in CC mammograms, and N_y = 9 (Figure 8B′) to 3 (Figure 8D′) in MLO mammograms. These observations suggest that the surgery and associated therapeutic treatment were efficient with regards to removing the malignant tumor, as well as of its tumorigenic microenvironment including cancer stem cells. They also attest to some regeneration of the breast tissue, probably via the activation of quiescent epithelial stem cells housed in the ducts, and which were shown to exhibit a high proliferative, self-renewal and morphogenic capacity in culture (Villadsen et al., 2007). The underlying recovery mechanisms appear to be engaged and operating not only in the cancerous breast but also in the other breast which contained an above average number of “risky” H = 0.5 squares on both mammographic views.

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.