Wavelet-based multifractal analysis of dynamic infrared thermograms to assist in early breast cancer diagnosis

Breast cancer is the most common type of cancer among women and despite recent advances in the medical field, there are still some inherent limitations in the currently used screening techniques. The radiological interpretation of screening X-ray mammograms often leads to over-diagnosis and, as a consequence, to unnecessary traumatic and painful biopsies. Here we propose a computer-aided multifractal analysis of dynamic infrared (IR) imaging as an efficient method for identifying women with risk of breast cancer. Using a wavelet-based multi-scale method to analyze the temporal fluctuations of breast skin temperature collected from a panel of patients with diagnosed breast cancer and some female volunteers with healthy breasts, we show that the multifractal complexity of temperature fluctuations observed in healthy breasts is lost in mammary glands with malignant tumor. Besides potential clinical impact, these results open new perspectives in the investigation of physiological changes that may precede anatomical alterations in breast cancer development.

1. Introduction

It is widely recognized that early diagnosis is the key to breast cancer survival. X-ray mammography (Nass et al., 2000; Bronzino, 2006), the golden standard for breast cancer screening detection, has a rather high false-positive rating and is not always effective in detecting cancer in young women who generally have dense breast tissue (Jorgensen and Gotzsche, 2009; Tamini et al., 2009) and this despite the increasing use of computer-aided detection/diagnosis methods (Fenton et al., 2011). Biopsy is indeed the only conclusive diagnostic test for breast cancer, but the number of unnecessary biopsies is still too high (Vinitha Sree et al., 2009). Since the original observation by Lawson (1956) that skin temperature over a malignant tumor is higher than its neighbourhood, possibly resulting from abnormal increase of metabolic activity and vascular circulation in the tissues beneath (Yahara et al., 2003; Bronzino, 2006), IR thermography has been considered as a promising non-invasive screening method of breast cancer (Ng, 2009). However the suitability of static IR imaging for routine screening has been severely questioned (Head and Elliott, 2002; Bronzino, 2006), because of insufficient sensitivity for detection of deep lesions and limited knowledge of the relationship between surface temperature distributions and thermal diseases. Renewed interest in dynamic IR imaging (Etehadtavakol and Ng, 2013) comes from the rapid development of new digital IR thermography cameras with higher temperature resolution (0.08°C or better) and faster frame rate (70 Hz) (Joro et al., 2008a), combined with increasing knowledge of tumor angiogenesis including nitric oxide production of the cancer tissue causing local disturbances in vasomotor (automatic nervous control of smooth muscles forcing blood through capillaries) and cardiogenic phenomena as compared to normal tissues (Thomsen and Miles, 1998; Anbar et al., 2001).

The basis for diagnostic application of dynamic IR imaging is the detection of intensity variations in temperature rhythms generated by the cardiogenic (1–1.5 Hz) and vasomotor (0.1–0.2 Hz) frequencies (Button et al., 2004; Joro et al., 2008b). In this study we show that beyond intensity differences in these rhythms between normal and tumor breast tissues, the complexity of temperature fluctuations about these physiological perfusion oscillations is qualitatively different. Using a wavelet-based multi-scale analysis (Muzy et al., 1991, 1994; Arneodo et al., 1995) of temperature fluctuations, we propose to characterize the multifractal properties of these temperature time-series as an effective discriminating method for early screening procedures to identify women with high risk of breast cancer.

2. Methods of Analysis

2.1. The Wavelet Transform

The wavelet transform (WT) is a mathematical microscope (Muzy et al., 1991, 1994; Arneodo et al., 1995) that is well suited for the analysis of complex non-stationary time-series such as those found in physiological systems (Ivanov et al., 1999; Goldberger et al., 2002), thanks to its ability to filter out low-frequency trends in the analyzed signal Σ (Materials and Methods). The WT is a space (or time in our study)-scale analysis which consists in expanding signals in terms of wavelets which are constructed from a single function, the “analyzing wavelet” ψ, by means of translations and dilations. The WT of a real-valued function Σ is defined as (Mallat, 1998):

$\begin{array}{cc}{W}_{\psi }[\Sigma ](t_0,a)=\frac{1}{a}{\displaystyle {\int }_{-\infty }^{+\infty }\Sigma }(t)\psi \left(\frac{t-t_0}{a}\right)dt,& (1)\end{array}$

where t₀ is the time parameter and a (> 0) the scale parameter. By choosing a wavelet ψ whose n + 1 first moments are zero [∫ tⁿ ψ (t)dt = 0, 0 ≤ m ≤ n] (Supplementary Figure 1), 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). Indeed this mathematical microscope can be seen as a singularity scanner. By increasing magnification (decreasing the scale parameter a → 0⁺) around a given point t, finer and finer details of Σ can be revealed and quantified via the estimate of the so-called Holder exponent h(t) (Muzy et al., 1991, 1994; Arneodo et al., 1995).

2.2. The Wavelet Transform Modulus Maxima Method

The WT modulus maxima (WTMM) method was originally developed to generalize box-counting techniques (Arneodo et al., 1987) and to remedy for the limitations of the structure functions method to perform multifractal analysis of one-dimensional (1D) velocity signal in fully developed turbulence (Muzy et al., 1991, 1994; Arneodo et al., 1995). It has proved very efficient to estimate scaling exponents and multifractal spectra (Muzy et al., 1994; Delour et al., 2001; Audit et al., 2002). This method has been generalized in 2D for the multifractal analysis of rough surfaces (Arneodo et al., 2003) and then for the analysis of 3D scalar and vector fields (Kestener and Arneodo, 2004; Arneodo et al., 2008). It has been successfully applied in various areas of fundamental research (Arneodo et al., 1998a,b, 2002, 2003, 2008; Khalil et al., 2006; Roland et al., 2009; Roux et al., 2009; Arneodo et al., 2011). In the context of the present study, the 1D WTMM method has proved very efficient at discriminating between healthy and sick heart beat dynamics (Ivanov et al., 1999, 2001; Goldberger et al., 2002), whereas the 2D WTMM method can be used to detect microcalcifications and has great potential to assist in cancer diagnosis from digitized mammograms (Kestener et al., 2001; Arneodo et al., 2003).

The WT modulus maxima (WTMM) method (Muzy et al., 1991, 1994; Arneodo et al., 1995) consists in computing the WT skeleton defined, at each fixed scale a, by the local maxima (a) of the WT modulus |W(t, a)|. These WTMM are disposed on curves connected across scales called maxima lines l_t (Supplementary Figure 2), along which the WTMM behave as a^h(t), where h(t) is the Hölder exponent (Muzy et al., 1991, 1994; Arneodo et al., 1995) characterizing the singularity of Σ at time t. 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 where h(t) = h. This spectrum can be obtained by investigating the scaling behavior of partition functions defined in terms of wavelet coefficients:

$\begin{array}{cc}Z(q,a)={\displaystyle \sum _{l\in ℒ(a)}|W(t,a){|}^q~a^{\tau (q)}},& (2)\end{array}$

where q ϵ ℝ. Then from the scaling function τ(q), D(h) is obtained by a Legendre transform (Muzy et al., 1991, 1994; Arneodo et al., 1995):

$\begin{array}{cc}D(h)=\underset{q}{\min }[qh-\tau (q)].& (3)\end{array}$

As originally pointed out in Muzy et al. (1994); Arneodo et al. (1995), one can avoid some practical difficulties that occur when directly performing the Legendre transform of τ(q), by computing the following expectation values:

where Ŵ_ψ[Σ](q, l, a) = $\frac{|W_{\psi }[\Sigma ](t,a){|}^q}{Z(q,a)}$ is the equivalent of Bolzmann weight in the analogy that links the multifractal formalism to thermodynamics (Arneodo et al., 1995). Then from the slopes of h(q, a) and D(q, a) vs ln a, one gets h(q) and D(q) and therefore the D(h) singularity spectrum as a curve parametrized by q.

2.3. Monofractal vs. Multifractal Functions

Homogeneous monofractal functions, i.e., functions with singularities of unique Hölder exponent H, are characterized by a linear τ(q) curve of slope H. Monofractal scaling implies that the shape of the probability distribution function (pdf) of rescaled wavelet coefficients (W(·, a)/a^H) does not depend on a, formally expressed by the self-similarity relationship (Arneodo et al., 2002):

$\begin{array}{cc}{\rho }_{W_a/a^H}(w)=\rho (w),& (5)\end{array}$

where ρ(w) is a “universal” pdf. A non-linear τ(q) is the signature of non-homogeneous multifractal functions, meaning that the Hölder exponent h(t) is a fluctuating quantity (Muzy et al., 1991, 1994; Arneodo et al., 1995) that depends on t. In this study, we fit the τ(q) data by the so-called log-normal quadratic approximation τ(q) = − c₀ + c₁q − c₂q²/2, where the coefficients c_n > 0 (Delour et al., 2001). The corresponding singularity spectrum has a characteristic single-humped shape D(h) = c₀ − (h − c₁)²/2c₂, where c₀ = − τ(0) is the fractal dimension of the support of singularities of Σ, c₁ is the value of h that maximizes D(h) and c₂, the so-called intermittency coefficient (Muzy et al., 1991, 1994; Arneodo et al., 1995; Delour et al., 2001), characterizes the width of the D(h) spectrum, i.e., the signature of a change in WT coefficient statistics across scales (Muzy et al., 1991, 1994; Arneodo et al., 1995, 2002).

3. Description of Data

3.1. Study Design and Population

Subjects were recruited for the present study from the Perm Regional Oncological Dispensary using procedures approved by the Local Ethics Committee (Gileva et al., 2012). They all gave Informed Consent to participate in this study via the recording of the IR thermograms of the both mammary glands, the cancerous one and the opposite undiagnosed one with no visible signs of pathology. Our database includes 33 females with ages between 37 and 83 (average 57 years) who all went through surgery to remove the histologically confirmed malignant tumor (invasive ductal and/or lobular breast cancer) a few weeks after thermograms were recorded. The tumors were found at different depths from 1 cm down to 12 cm with a size varying from 1.2 cm up to 6.5 cm (Table 1). As a control, we also investigated 14 women with intact mammary glands and of ages between 23 and 79 (average 49.6 years). This extensive study was preceded and encouraged by a preliminary study with only 6 patients and 3 volunteers (Gerasimova et al., 2013).

TABLE 1

Table 1. Set of (33) analyzed patients with age, cancerous breast (Right or Left), stage, size and depth of the malignant tumor and histology status.

3.2. IR Thermography Imaging

Both breasts of healthy volunteers and patients with breast cancer were imaged with a InSb photovoltaic (PV) detector camera (Joro et al., 2008a). Imaging was performed with the patient in sitting position with arms down to avoid too much discomfort during imaging. Frontal images were taken at a distance ~1 m of the patient in an environmental room temperature of 20–22°C. The image frame rate was set to 50 Hz. The image data were collected in 14 bits into the computer connected to the PV camera. Each image set comprised 30 000 256 × 320 pixel² image frames during the 10 min immobile imaging phase. To eliminate low frequency patient movements, skin surface markers were successfully used as reference points for motion correction in the analysis.

3.3. Data Sampling

Pixel based and windowed regional power spectra and wavelet-based multifractal analysis of normal and cancer breasts were tested to define the best procedure to minimize the effect of the camera noise and to ensure statistical convergence in the multifractal spectra estimation. We grouped single-pixel temperature time-series (Figures 1A–C) into 8 × 8 pixel² squares spanning 10 × 10 mm² and covering the entire breast (see Figures 4A–C). The results reported correspond to averaged power spectra, partition functions, singularity spectra and WT pdfs over 64 temperature time-series in these 8 × 8 subareas.

FIGURE 1

Figure 1. Multifractal analysis of IR temperature time-series. Comparative analysis of the cancerous right breast (red) and healthy left breast (black) of patient 20 (age 56) and of the healthy right breast (green) of volunteer 14 (age 60). (A–C) 1 min portions of pixel temperature time-series. (D) Averaged temperature power spectra in a 8 × 8 pixel² square. The straight lines correspond to power-law scaling 1/f^β with exponent β = τ(2) = 0.62 (red), 1.32 (black), and 1.22 (green) as estimated with the WTMM method in (E). (E) τ(q) vs. q estimated by linear regression fit of log₂Z(q, a) vs. log₂ a over a range of time-scales [0.3, 3] s (Figure 2A). (F) D(h) vs. h (Figures 2B,C). The solid lines in (E,F) correspond to quadratic spectra (see text) with parameters [c₀, c₁, c₂] = [0.99, 0.81, 0.0044] (red), [0.99, 1.23, 0.080] (black), and [0.99, 1.171, 0.069] (green). In (E,F), the τ(q) and D(h) spectra were averaged over a 8 × 8 pixel² square.

3.4. Remark

As commonly done for noise signals (Muzy et al., 1994; Audit et al., 2002) and previously experienced when applying the 1D WTMM method to rainfall time-series (Venugopal et al., 2006), the wavelet analysis was performed on the cumulative (or integral) Σ of the temperature time-series (Supplementary Figure 2) using the second-order compactly supported version ψ⁽²⁾₍₃₎ of the Mexican hat (Roux et al., 1999) (Supplementary Figure 1). Hence the singularities with possible negative Hölder exponent −1 < h < 0, become singularities with 0 < h_c = h + 1 < 1 in the cumulative.

3.5. Software and Documentation

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

LastWave is an open source software written in C. We recommend interested users to 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.

4. Results

4.1. Spectral Analysis Reveals Scale-Invariance Properties in Skin Temperature Dynamics in Both Cancerous and Healthy Breasts

We analyzed individual 1-pixel temperature time-series taken from 8 × 8 pixels²-squares covering the patients entire breasts. As expected, these time-series generally fluctuate at a higher temperature when recorded in the tumor region of a malignant breast (Figure 1A) than in a symmetrically positioned square on the opposite breast (Figure 1B) as well as on a healthy breast (Figure 1C). When averaging the corresponding power-spectra over the 64 pixels of the considered squares, we observed for the two breasts of patient 20 (age 56) and the healthy breast of volunteer 14 (age 60), a rather convincing 1/f^β power-law scaling over a range of frequencies that extends from the characteristic human respiratory frequency (≳0.3 Hz) up to the cross-over frequency (≲4 Hz) toward (instrumental) white noise (Figure 1D). As confirmed by the WTMM method (Figure 1E), the exponent β = τ (q = 2) = 0.62 ± 0.19 found in the malignant breast is smaller than in the opposite breast of patient 20, β = 1.32 ± 0.11, and in the volunteer 14 healthy breast, β = 1.22 ± 0.11. This difference looks quite significant and very promising in a discriminatory perspective. Unfortunately, the histograms of β values obtained for all 8 × 8 pixel² squares covering 33 cancerous breasts, the 32 opposite breasts (patient 6 had mastectomy of right breast) and the 28 volunteer healthy breasts are quite similar (Supplementary Figure 3) with mean values β = 1.09 ± 0.01 (cancer), 1.14 ± 0.01 (opposite) and 1.14 ± 0.01 (healthy). Indeed these histograms extend over a rather wide range of β values: 0.5 ≲ β ≲ 1.9.

4.2. WTMM Analysis Discriminates between Monofractal (Tumor Area) and Multifractal (Healthy Area) Temperature Temporal Fluctuations

When applying the WTMM method to the cumulative of these temperature time-series, we confirmed that the partition functions Z(q, a) Equation (2) display nice scaling properties for q = − 1 to 5, over a range of time-scales that we strictly limited to (0.43, 2.30 s) for linear regression fit estimates in a logarithmic representation (Figure 2A). The τ(q) so-obtained are well approximated by quadratic spectra (Figure 1E). For the malignant breast of patient 20, τ(q) is nearly linear as quantified by a very small value of the intermittency coefficient c₂ = (4.4 ± 0.6) · 10⁻³. This signature of monofractality is confirmed, when respectively plotting h(q, a) and D(q, a) Equation (4) vs. log₂ a, in Figures 2B,C, where the slopes h(q) = c₁ = 0.81 ± 0.01 and D(q) = 0.99 ± 0.03, do not significantly depend on q, meaning that the D(h) singularity spectrum nearly reduces to a single point D(h = c₁ = 0.81) = 1 (Figure 1F). This monofractal diagnosis is confirmed when comparing the WT pdfs obtained at different time-scales (Figure 3A); according to Equation (5), they all collapse on a single curve when using the exponent H = c₁ (Figure 3A′).

FIGURE 2

Figure 2. Multifractal analysis of cumulative IR temperature time-series in a 8 × 8 pixel² square in malignant right breast (red) and opposite left breast (black) of patient 20 and of the healthy right breast (green) of volunteer 14. (A) log₂Z(q, a) vs. log₂(a) Equation (2). (B) h(q, a) vs. log₂a Equation (4). (C) D(q, a) vs. log₂ a Equation (4). The τ(q) spectra shown in Figure 1E were estimated by linear regression fit of the data in (A) over the range 2.8 ⩽ log₂ a ⩽ 5.2. The D(h) spectra in Figure 1F were obtained by linear regression fit in (B,C) over the same range of time-scales. The analyzing wavelet is ψ⁽²⁾₍₃₎ (Supplementary Figure 1).

FIGURE 3

Figure 3. Evolution of the temperature WT coefficient pdfs across time-scales. Pdfs of WT coefficients W and rescaled WT coefficients (W/(a/a₀)^c₁) of cumulative temperature of (A,A′) the cancerous right breast of patient 20: c₁ = 0.81, (B,B′) the opposite left breast of patient 20: c₁ = 1.23, and (C,C′) the healthy right breast of volunteer 14: c₁ = 1.17. The different curves correspond to seven different scales from a = 0.43 s to 2.30 s, larger scales are in darker shades; a₀ = 0.43 s. The pdfs were averaged over a 8 × 8 pixel² square.

In contrast, the τ(q) spectrum obtained for the opposite breast of patient 20, is definitely non-linear with a no longer negligible (one order of magnitude larger) quadratic term c₂ = 0.080 ± 0.001 (Figure 1E), the hallmark of multifractal scaling. As shown in Figures 2B,C, the slopes h(q) and D(q) of h(q, a) and D(q, a) vs log₂ a now depend on q. From the estimate of h(q) and D(q), we get the single-humped D(h) spectrum shown in Figure 1F, which is well approximated by a quadratic spectrum with parameters c₀ = 0.99 ± 0.05, c₁ = 1.23 ± 0.01 and c₂ = 0.080 ± 0.001. Because there is no longer a unique scaling exponent c₁, the self-similarity Equation (5) is not verified meaning that the shape of the WT coefficient pdf now evolves across scales (Figure 3B), with fatter tails appearing at small scales (Figure 3B′). Interestingly, the τ(q) (Figure 1E) and D(h) (Figure 1F) spectra obtained for the healthy breast of volunteer 14 are quite similar quadratic spectra with parameter values c₀ = 0.99 ± 0.03, c₁ = 1.17 ± 0.01 and c₂ = 0.069 ± 0.002. Again this multifractal diagnosis is strengthened by the observation that the WT coefficient pdf has a shape that evolves across time-scales (Figures 3C,C′).

To check the statistical relevance of our multifractal spectra estimates, we have generated so-called surrogate series (Theiler et al., 1992; Schreiber and Schmitz, 1996) that can be created with an identical pdf and optimally similar power spectrum to the original series (Supplementary Method). Examples of surrogate series for the IR temperature time-series shown in Figures 1A–C can be seen in Supplementary Figures 4A–C respectively. Visually, there is an obvious resemblance with the original time-series, but when computing the corresponding τ(q) (Supplementary Figure 4E) and D(h) (Supplementary Figure 4F), we find now for the three breasts a τ(q) spectrum that is quite linear and a D(h) singularity spectrum that almost reduces to a single point D(h) = c₁ with a very small width c₂ ≤ 0.01. This monofractal diagnostic is confirmed when reproducing this analysis for 100 surrogate series; the histograms of intermittency coefficients c₂ obtained for the three breasts are very similar and mainly confined to very small values with means c₂ = 0.012 (cancer), 0.005 (opposite) and 0.006 (healthy) that are all much smaller than the threshold c₂ = 0.03 we will further use to discriminate between monofractal (c₂ ≤ 0.03) and multifractal (c₂ > 0.03) cumulative temperature time-series (Supplementary Figure 5). These results indicate that the cumulative temperature time-series of healthy breasts are not generated by an underlying linear Gaussian process, but have an inherently non-linear structure that is apparently lost in the presence of a malignant tumor.

4.3. Multifractal-Based Segmentation of Breast Thermograms into Physiologically Altered (Risky) and Normal Regions

When extending our wavelet-based multifractal analysis of cumulative temperature time-series to the entire set of 8 × 8 pixels²-squares that cover the right breast with malignant tumor of patient 20 (Figures 4A,D), her opposite left breast (Figures 4B,E) and the healthy right breast of volunteer 14 (Figures 4C,F), we confirmed, except in a minority of squares, the existence of scaling. In the cancerous breast, a large proportion of squares (49.7%) display monofractal temperature fluctuations with small intermittency coefficient values (c₂ < 0.03 in Figure 4D), whereas only few of those squares are found in the opposite breast (7.7% in Figure 4E) and in the volunteer 14 healthy breast (11% in Figure 4F). Both these healthy breasts have a large majority of squares where multifractal scaling is observed (c₂ ⩾ 0.03), namely 89.4% for the former (Figure 4B) and 65% for the latter (Figure 4C). Note that 43.1% of the squares in the cancerous breast also display multifractal temperature fluctuations as observed for healthy breasts (Figure 4A). These squares indeed cover regions of the breast that are far from the tumor area (left upper quadrant) mostly covered by monofractal squares. These results are indeed quite representative of the outcome of the statistical analysis of our entire data set.

FIGURE 4

Figure 4. Breast-wide multifractal analysis of skin temperature temporal fluctuations. As estimated from the τ(q) spectra computed with the WTMM method (Figure 1E), 8 × 8 pixel² squares covering the entire breast are colored coded (A–C) and represented as a dot in the (c₁, c₂) plane (D–F). The colors have the following meaning: red: c₂ < 0.03, blue: c₂ ⩾ 0.03, and white: no scaling (see text). (A,D) Cancerous right breast of patient 20: the tumor is located in the upper outer quadrant (Table 1); (B,E) healthy left breast of patient 20; (C,F) healthy right breast of volunteer 14. The blue diamond shapes in (A–C) are skin surface markers used as references for motion correction.

4.4. Comparative Statistical Analysis of Skin Temperature Dynamics in Women Breasts with and without Malignant Tumor

The results of our comparative wavelet-based multifractal analysis of (cumulative) temperature fluctuations over the 33 cancerous breasts, the 32 opposite breasts (no right breast for patient 6) and the 28 volunteer healthy breasts are reported in Figure 5 and Table 2. In Figure 5A, the corresponding histograms of c₁ values extend over a rather wide range 0.6 ≲ c₁ ≲ 1.8 but turn out to be quite similar with mean values c₁ = 1.066 ± 0.002 (cancer), 1.104 ± 0.002 (opposite) and 1.103 ± 0.002 (healthy). This is reminiscent to the similar histograms that were also obtained for the power-spectrum exponent β (Supplementary Figure 3). In contrast, the intermittency parameter c₂ has a definite discriminatory power. The histogram for cancerous breasts (c₂ = 0.045 ± 0.001) is definitely shifted toward smaller values relative to the ones for opposite (c₂ = 0.056 ± 0.001) and healthy breasts (c₂ = 0.058 ± 0.001) (Figure 5B). The small-value left-side of the c₂ histogram is much more populated in cancerous than in healthy breasts, confirming that cancerous breasts are enriched in squares where temperature fluctuations display significantly reduced multifractal properties. This justifies that we considered c₂ = 0.03 as the threshold below (resp. above) which a square was qualified as monofractal (resp. multifractal) (Figure 4). Note that for each breast, a small percentage (≲20%) of (white) squares were removed from our analysis because of lack of scaling (Figure 5C, Table 2 and Supplementary Table 1). Importantly, the percentage of monofractal (red) squares covering cancerous breasts (26.8 ± 3.5%) is about twice larger than the ones covering opposite (13.1 ± 2.3%) and healthy (11.3 ± 2.2%) breasts. This excess is compensated by a smaller percentage of multifractal (blue) squares in cancerous breasts (56.9 ± 4.4%) than in opposite (68.5 ± 3.8%) and healthy (69.4 ± 4.3%) breasts (Figure 5C).

FIGURE 5

Figure 5. Differential multifractal signature of temperature temporal fluctuations on 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₁ values obtained in the breasts with a malignant tumor (red: N = 4032 8 × 8 pixel² squares), the opposite breasts (black: N = 3606) and the healthy breasts (green: N = 3185). (B) Normalized histograms of c₂ values. (C) Percentage of “monofractal” (red: c₂ < 0.03), “multifractal” (blue: c₂ ⩾ 0.03) and “no scaling” (black) 8 × 8 pixel² squares in cancerous, opposite and healthy breasts. (D) Percentage of “monofractal” 8 × 8 pixel² squares in the cancerous breast vs. the percentage in the opposite breast of each patient; the symbols have the following meaning: (•) 17 patients with a percentage of monofractal squares on the cancerous breast larger than the mean value 26.8%, namely patients 2, 4, 7, 8, 9, 10, 11, 13, 17, 19, 20, 22, 25, 26, 29, 30, and 32 (patient 6 is not taken into account because of mastectomy of right breast); (*) 7 patients with the percentage of monofractal squares on the cancerous breast smaller than 26.8% but well localized on the tumor region, namely patients 3 (16.3%), 14 (18.9%), 21 (14.0%), 23 (4.6%), 24 (4.9%), 31 (6.8%), and 33 (17.0%); (°) 8 false negatives, namely patients 1 (2.4%), 5 (0.8%), 12 (4.9%), 15 (1.3%), 16 (0.0%), 18 (7.7%), 27 (4.3%), and 28 (6.6%) (see Table 2).

TABLE 2

Table 2. Results of our WTMM multifractal analysis of the skin (cumulative) temperature temporal fluctuations of the two breasts of our 33 patients with breast cancer.

A common way to suspect cancer by IR thermography is to look for some dissymmetry between the two breasts of a patient (Etehadtavakol and Ng, 2013). When comparing the percentages of monofractal squares on both the cancer and opposite breasts of the 33 patients (except patient 6), we found that 25 (/32) have more monofractal squares on the cancerous breast (Table 2). Indeed, we found 18 (/33) malignant breasts that have a percentage of monofractal squares greater than the mean value 26.8 ± 3.5 (Figure 5D). Among the other 15 cancerous breasts, 7 have a smaller percentage of monofractal squares but well localized on the tumor region (Figure 4A and Supplementary Figure 6). The remaining 8 cancerous breasts correspond to false negatives for which not only the percentage of monofractal squares is small but their location is far from the tumor region (Supplementary Figure 7). Among these false negatives, 4 correspond to rather deep tumors in fatty breasts which can explain that they do not manifest in a qualitative change in temperature dynamics at the skin surface in patients 12 (size 1.8 cm, depth 12 cm) (Supplementary Figure 8), 16 (3.4 cm, 7 cm), 18 (3.49 cm, 6 cm) and 28 (3.49 cm, 8 cm). When investigating the 32 opposite breasts, 5 of them have a large percentage of monofractal squares (Figure 5D, Table 2 and Supplementary Figure 9). These important percentages similar to those obtained in malignant breasts are probable indication of some physiological changes in the opposite breast that may announce the possible extension of cancer to the second breast. As a control, we reproduced this comparative analysis on the two breasts of the 14 healthy volunteers (Supplementary Figure 10 and Table 1). Among the 28 breasts analyzed, only 4 have a large (⩾26.8%) percentage of monofractal squares, whereas most of them have a percentage ≲10% (Supplementary Figures 11, 12). Overall, we thus obtained 25 (/33) true positives and 4 (/28) false positives, i.e., a sensitivity of 76% and a specificity of 86%, respectively.

5. Conclusions

Over the course of a lifetime, 1 in 8 women will be diagnosed with breast cancer. There are no well-established ways to avoid breast cancer (as opposed to lung cancer for example) and in the context of breast cancer screening, abnormalities should be detected at an early stage to improve prognosis. Criticism of the use of screening mammography due to over-diagnosis led some researchers to show that one in three breast cancers identified by mammography would not cause symptoms in a patient's lifetime (Jorgensen and Gotzsche, 2009). Therefore, alternative and accurate screening technologies must be developed. The functional and technical background of dynamic IR imaging has the potential for early detection of breast cancer and treatment response evaluation if optimal diagnostic algorithms are developed. We have shown that the wavelet-based multifractal analysis of dynamic IR thermograms is able to discriminate between cancerous breasts with monofractal (cumulative) temperature temporal fluctuations characterized by a unique singularity exponent (h = c₁), and healthy breasts with multifractal temperature fluctuations requiring a wide range of singularity exponents as quantified by the intermittency coefficient c₂ ≫ 0. This is strikingly analogous to the results of a similar wavelet-based analysis of human heart beat dynamics (Ivanov et al., 1999, 2001; Goldberger et al., 2002), where the multifractal character and non-linear properties of the healthy heart rate were shown to be lost in pathological condition, congestive heart failure. Indeed, this distinction was intrinsically beyond the capability of spectral (Fourier) analysis which only gives access to the power-spectrum exponent β = τ(q = 2) = −c₀ + 2c₁ − 2c₂, and not to the full τ(q) spectrum required for multifractal diagnosis (Muzy et al., 1991, 1994; Arneodo et al., 1995, 2008). Furthermore the fact that c₁ (~ 1) is about one order of magnitude larger than the intermittency coefficient c₂ (≲ 0.1), explains why very much like c₁ (Figure 5A), the spectral exponent β ~ −c₀ + 2c₁ ~ 1 ≫ c₂ (Supplementary Figure 3) has no discriminatory power.

Interdisciplinary effort revealing specific fractal characteristics for healthy and cancerous breast tissues definitely challenges current knowledge in physical, physiological and clinical fundamentals of oncogenesis. Fundamentally, our results indicate that skin temperature fluctuations of healthy breasts are more complex (multifractal) than previously suspected. They definitely raise new challenging questions to ongoing efforts to develop physiological 3D breast models that account for the skin surface temperature distribution in the presence (or absence) of an internal tumor (Ng and Sudharsan, 2004; Xu et al., 2008; Lin et al., 2009). The observed drastic simplification from multifractal to monofractal skin temperature dynamics may result from some increase in blood flow and cellular activity associated with the presence of a tumor (Thomsen and Miles, 1998; Anbar et al., 2001; Button et al., 2004; Joro et al., 2008b). More likely it can be the signature of some architectural change in the cellular microenvironment of the breast tumor (Bissell and Hines, 2011) that may deeply affect heat transfer and related thermomechanics in breast tissue (Xu et al., 2008; Quail and Joyce, 2013). Identifying the regulation mechanisms that originate in a loss of multifractal temperature dynamics will be an important step toward understanding breast cancer development, tumor growth and progression. Dynamic IR thermography is a non-invasive and objective screening method that is inexpensive, quick and painless for the patient. Future use of wavelet-based multifractal processing of dynamic IR thermography, could help identifying women with high risk of breast cancer, prior to more traumatic and painful examination such as mammography and biopsy. It can also prove to be a valuable and reliable adjunct tool for early detection of tumors in other locations than in mammary glands.

Grant Support

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 Perm Regional Government (Russia) for the contract “Multiscale approaches in mechanobiology for early cancer diagnosis” and the Maine Cancer Foundation.

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.

Acknowledgments

We are thankful to Prof. O. Orlov for his support in IR thermal imaging at Perm Region Cancer Hospital and to I. Panteleev, Y. Bayandin, O. Plekhov, and S. Uvarov for fruitful discussions.

Supplementary Material

The Supplementary Material for this article can be found online at: http://www.frontiersin.org/journal/10.3389/fphys.2014.00176/abstract

References

Anbar, M., Milescu, L., Naumov, A., Brown, C., Button, T., Carty, C., et al. (2001). Detection of cancerous breasts by dynamic area telethermometry. IEEE Eng. Med. Biol. 20, 80–91. doi: 10.1109/51.956823

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Arneodo, A., Audit, B., Decoster, N., Muzy, J. F., and Vaillant, C. (2002). “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, eds A. Bunde, J. Kropp, and H. J. Schellnhuber (Berlin: Springer Verlag), 26–102.

Arneodo, A., Audit, B., Kestener, P., and Roux, S. G. (2008). Wavelet-based multifractal analysis. Scholarpedia 3: 4103. doi: 10.4249/scholarpedia.4103

CrossRef Full Text

Arneodo, A., Bacry, E., and Muzy, J. F. (1995). The thermodynamics of fractals revisited with wavelets. Physica A 213, 232–275. doi: 10.1016/0378-4371(94)00163-N

CrossRef Full Text

Arneodo, A., Decoster, N., Kestener, P., and Roux, S. G. (2003). A wavelet-based method for multifractal image analysis: from theoretical concepts to experimental applications. Adv. Imaging Electr. Phys. 126, 1–92. doi: 10.1016/S1076-5670(03)80014-9

CrossRef Full Text

Arneodo, A., Grasseau, G., and Kostelich, E. J. (1987). Fractal dimensions and f(α) spectrum of the Henon attractor. Phys. Lett. A 124, 426–432. doi: 10.1016/0375-9601(87)90546-9

CrossRef Full Text

Arneodo, A., Manneville, S., and Muzy, J. F. (1998). Towards log-normal statistics in high Reynolds number turbulence. Eur. Phys. J. B 1, 129–140. doi: 10.1007/s100510050162

CrossRef Full Text

Arneodo, A., Muzy, J. F., and Sornette, D. (1998). “Direct” causal cascade in the stock market. Eur. Phys. J. B 2, 277–282. doi: 10.1007/s100510050250

CrossRef Full Text

Arneodo, A., Vaillant, C., Audit, B., Argoul, F., d'Aubenton-Carafa, Y., and Thermes, C. (2011). Multi-scale coding of genomic information: from DNA sequence to genome structure and function. Phys. Rep. 498, 45–188. doi: 10.1016/j.physrep.2010.10.001

CrossRef Full Text

Audit, B., Bacry, E., Muzy, J. F., and Arneodo, A. (2002). Wavelet-based estimators of scaling behavior. IEEE Trans. Info. Theory 48, 2938–2954. doi: 10.1109/TIT.2002.802631

CrossRef Full Text

Bissell, M. J., and Hines, W. C. (2011). Why don't we get more cancer? A proposed role of the microenvironment in restraining cancer progression. Nat. Med. 17, 320–329. doi: 10.1038/nm.2328

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Bronzino, J. D. (2006). Biomedical Engeneering Handbook. Boca Raton, FL: CRC Press.

Button, T. M., Li, H., Fisher, P., Rosenblatt, R., Dulaimy, K., Li, S., et al. (2004). Dynamic infrared imaging for the detection of malignancy. Phys. Med. Biol. 49, 3105–3116. doi: 10.1088/0031-9155/49/14/005

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Delour, J., Muzy, J. F., and Arneodo, A. (2001). Intermittency of 1D velocity spatial profiles in turbulence: a magnitude cumulant analysis. Eur. Phys. J. B 23, 243–248. doi: 10.1007/s100510170074

CrossRef Full Text

Etehadtavakol, M., and Ng, E. Y. K. (2013). Breast thermography as a potential non-contact method in the early detection of cancer: a review. J. Mech. Med. Biol. 13, 1330001. doi: 10.1142/S0219519413300019

CrossRef Full Text

Fenton, J. J., Abraham, L., Taplin, S. H., Geller, B. M., Carney, P. A., D'Orsi, C., et al. (2011). Effectiveness of computer-aided detection in community mammography practice. J. Natl. Cancer Inst. 103, 1152–1161. doi: 10.1093/jnci/djr206

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Gerasimova, E., Audit, B., Roux, S., Khalil, A., Argoul, F., Naimark, O., et al. (2013). Multifractal analysis of dynamic infrared imaging of breast cancer. Eur. Phys. Lett. 104, 68001. doi: 10.1209/0295-5075/104/68001

CrossRef Full Text

Gileva, O. S., Freynd, G. G., Orlov, O. A., Libik, T. V., Gerasimova, E. I., Plekhov, O. A., et al. (2012). Interdisciplinary approaches to early diagnosis and screening of tumors and precancerous diseases (for example, breast cancer). RFBR J. 74–75, 93–99.

Goldberger, A. L., Amaral, L. A. N., Hausdorff, J. M., Ivanov, P. C., Peng, C. K., and Stanley, H. E. (2002). Fractal dynamics in physiology: alterations with disease and aging. Proc. Natl. Acad. Sci. U.S.A. 99, 2466–2472. doi: 10.1073/pnas.012579499

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Head, J. F., and Elliott, R. L. (2002). Infrared imaging: making progress in fulfilling its medical promise. IEEE Eng. Med. Biol. 21, 80–85. doi: 10.1109/MEMB.2002.1175142

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Ivanov, P. C., Amaral, L. A., Goldberger, A. L., Havlin, S., Rosenblum, M. G., Struzik, Z. R., et al. (1999). Multifractality in human heartbeat dynamics. Nature 399, 461–465. doi: 10.1038/20924

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Ivanov, P. C., Amaral, L. A. N., Goldberger, A. L., Havlin, S., Rosenblum, M. G., Stanley, H. E., et al. (2001). From 1/f noise to multifractal cascades in heartbeat dynamics. Chaos 11, 641–652. doi: 10.1063/1.1395631

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Jorgensen, K. J., and Gotzsche, P. C. (2009). Overdiagnosis in publicly organised mammography screening programmes : systematic review of incidence trends. BMJ 339:b2587. doi: 10.1136/bmj.b2587

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Joro, R., Lääperi, A. L., Dastidar, P., Soimakallio, S., Kuukasjäri, T., Toivonen, T., et al. (2008). Imaging of breast cancer with mid- and long-wave infrared camera. J. Med. Eng. Technol. 32, 189–197. doi: 10.1080/03091900701234358

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Joro, R., Lääperi, A. L., Soimakallio, S., Järvenpää, R., Kuukasjäri, T., Toivonen, T., et al. (2008). Dynamic infrared imaging in identification of breast cancer tissue with combined image processing and frequency analysis. J. Med. Eng. Technol. 32, 325–335. doi: 10.1080/03091900701541240

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Kestener, P., and Arneodo, A. (2004). Generalizing the wavelet-based multifractal formalism to random vector fields: application to three-dimensional turbulence velocity and vorticity data. Phys. Rev. Lett. 93, 044501. doi: 10.1103/PhysRevLett.93.044501

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Kestener, P., Lina, J. M., Saint-Jean, P., and Arneodo, A. (2001). Wavelet-based multifractal formalism to assist in diagnosis in digitized mammograms. Image Anal. Stereol. 20, 169–174. doi: 10.5566/ias.v20.p169-174

CrossRef Full Text

Khalil, A., Joncas, G., Nekka, F., Kestener, P., and Arneodo, A. (2006). Morphological analysis of H_I features. II. Wavelet-based multifractal formalism. Astrophys. J. Suppl. Ser. 165, 512–550. doi: 10.1086/505144

CrossRef Full Text

Lawson, R. N. (1956). Implications of surface temperatures in the diagnosis of breast cancer. Can. Med. Assoc. J. 75, 309–310.

Pubmed Abstract | Pubmed Full Text

Lin, Q. Y., Yang, H. Q., Xie, S. S., Wang, Y. H., Ye, Z., and Chen, S. Q. (2009). Detecting early breast tumour by finite element thermal analysis. J. Med. Eng. Technol. 33, 274–280. doi: 10.1080/03091900802106638

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Mallat, S. (1998). A Wavelet Tour of Signal Processing. New York, NY: Academic Press.

Muzy, J. F., Bacry, E., and Arneodo, A. (1991). Wavelets and multifractal formalism for singular signals: application to turbulence data. Phys. Rev. Lett. 67, 3515–3518. doi: 10.1103/PhysRevLett.67.3515

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Muzy, J. F., Bacry, E., and Arneodo, A. (1994). The multifractal formalism revisited with wavelets. Int. J. Bifurc. Chaos 4, 245–302. doi: 10.1142/S0218127494000204

CrossRef Full Text

Nass, S. J., Henderson, I. C., and Lashof, J. C. (2000). Mammography and Beyond: Developing Technologies for the Early Detection of Breast Cancer. Washington, DC: National Academy Press.

Ng, E. Y. K. (2009). A review of thermography as promising non-invasive detection modality for breast tumor. Int. J. Therm. Sci. 49, 849–859. doi: 10.1016/j.ijthermalsci.2008.06.015

CrossRef Full Text

Ng, E. Y. K., and Sudharsan, N. M. (2004). Computer simulation in conjunction with medical thermography as an adjunct tool for early detection of breast cancer. BMC Cancer 4:17. doi: 10.1186/1471-2407-4-17

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Quail, D. F., and Joyce, J. A. (2013). Microenvironmental regulation of tumor progression and metastasis. Nat. Med. 19, 1423–1437. doi: 10.1038/nm.3394

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Roland, T., Khalil, A., Tanenbaum, A., Berguiga, L., Delichère, P., Bonneviot, L., et al. (2009). Revisiting the physical processes of vapodeposited thin gold films on chemically modified glass by atomic force and surface plasmon microscopies. Surf. Sci. 603, 3307–3320. doi: 10.1016/j.susc.2009.09.021

CrossRef Full Text

Roux, S. G., Venugopal, V., Fienberg, K., Arneodo, A., and Foufoula-Georgiou, E. (2009). Evidence for inherent nonlinearity in temporal rainfall. Adv. Water Resourc. 32, 41–48. doi: 10.1016/j.advwatres.2008.09.007

CrossRef Full Text

Roux, S., Muzy, J. F, and Arneodo, A. (1999). Detecting vorticity filaments using wavelet analysis: about the statistical contribution of vorticity filaments to intermittency in swirling turbulent flows. Eur. Phys. J. B 8, 301–322. doi: 10.1007/s100510050694

CrossRef Full Text

Schreiber, T., and Schmitz, A. (1996). Improved surrogate data for nonlinear tests. Phys. Rev. Lett. 77, 635–638. doi: 10.1103/PhysRevLett.77.635

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Tamini, R. M., Colditz, G. A., and Hankinson, S. E. (2009). Circulating carotenoids, mammographic density, and subsequent risk of breast cancer. Cancer Res. 69, 9323. doi: 10.1158/0008-5472.CAN-09-1018

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Theiler, J., Eubank, S., Longtin, A., Galdikian, B., and Farmer, J. D. (1992). Testing for nonlinearity in time series: the method of surrogate data. Physica D 58, 77–94. doi: 10.1016/0167-2789(92)90102-S

CrossRef Full Text

Thomsen, L. L., and Miles, D. W. (1998). Role of nitric oxide in tumour progression: lessons from human tumours. Cancer Metast. Rev. 17, 107–118. doi: 10.1023/A:1005912906436

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Venugopal, V., Roux, S. G., Foufoula-Georgiou, E., and Arneodo, A. (2006). Revisiting multifractality of high-resolution temporal rainfall using a wavelet-based formalism. Water Resour. Res. 42, W06D14. doi: 10.1029/2005WR004489

CrossRef Full Text

Vinitha Sree, S., Ng, E. Y. K., Kaw, G., Acharya, U. R., and Chong, B. K. (2009). The use of skin surface electropotentials for breast cancer detection. Preliminary clinical trial results obtained using the Biofield Diagnostic System. J. Med. Syst. 35, 79–86. doi: 10.1007/s10916-009-9343-0

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Xu, F., Lu, T. J., and Steffen, K. A. (2008). Biothermomechanical behaviour of skin tissue. Acta Mech. Sin. 24, 1–23. doi: 10.1007/s10409-007-0128-8

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Yahara, T., Koga, T., Yoshida, S., Nakagawa, S., Deguchi, H., and Shirouzou, K. (2003). Relationship between micro vessel density and thermographic hot areas in breast cancer. Surg. Today 33, 243–248. doi: 10.1007/s005950300055

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text