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.


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(Muzy et al., , 1994Arneodo 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.

THE WAVELET TRANSFORM
The wavelet transform (WT) is a mathematical microscope (Muzy et al., 1991(Muzy et al., , 1994Arneodo 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): where t 0 is the time parameter and a (> 0) the scale parameter. By choosing a wavelet ψ whose n + 1 first moments are Figure 1), one makes the WT microscope blind to order-n polynomial behavior, a prerequisite for multifractal fluctuations analysis (Muzy et al., 1991(Muzy et al., , 1994Arneodo 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(Muzy et al., , 1994Arneodo et al., 1995).

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(Muzy et al., , 1994Arneodo 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(Arneodo et al., , 2003(Arneodo et al., , 2008Khalil 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(Ivanov et al., , 2001Goldberger 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(Muzy et al., , 1994Arneodo et al., 1995) consists in computing the WT skeleton defined, at each fixed scale a, by the local maxima L(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(Muzy et al., , 1994Arneodo 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: where q ∈ R. Then from the scaling function τ (q), D(h) is obtained by a Legendre transform (Muzy et al., 1991(Muzy et al., , 1994Arneodo et al., 1995): 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: 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.

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 : 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(Muzy et al., , 1994Arneodo et al., 1995) that depends on t.

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

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

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 2 squares spanning 10 × 10 mm 2 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.

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.

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

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 2 = (4.4 ± 0.6) · 10 −3 . This signature of monofractality is confirmed, when respectively plotting h(q, a) and D(q, a) Equation (4)  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 2 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 ψ (2) (Supplementary Figure 1).

www.frontiersin.org
May 2014 | Volume 5 | Article 176 | 5 depend on q, meaning that the D(h) singularity spectrum nearly reduces to a single point D(h = c 1 = 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 1 (Figure 3A ). 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 2 = 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 2 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 = 0.99 ± 0.05, c 1 = 1.23 ± 0.01 and c 2 = 0.080 ± 0.001. Because there is no longer a unique scaling exponent c 1 , 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 = 0.99 ± 0.03, c 1 = 1.17 ± 0.01 and c 2 = 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) 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 1 with a very small width c 2 ≤ 0.01. This monofractal diagnostic is confirmed when reproducing this analysis for 100 surrogate series; the histograms of intermittency coefficients c 2 obtained for the three breasts are very similar and mainly confined to very small values with means c 2 = 0.012 (cancer), 0.005 (opposite) and 0.006 (healthy) that are all much smaller than the threshold c 2 = 0.03 we will further use to discriminate between monofractal (c 2 ≤ 0.03) and multifractal (c 2 > 0.03) cumulative temperature time-series (Supplementary Figure 5). These results indicate that the cumulative temperature timeseries 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.

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 2squares 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 2 < 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 2 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.
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   Figure 4 and Supplementary Figures 4-7).

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 1 ), and healthy breasts with multifractal temperature fluctuations requiring a wide range of singularity exponents as quantified by the intermittency coefficient c 2 0. This is strikingly analogous to the results of a similar waveletbased analysis of human heart beat dynamics (Ivanov et al., 1999(Ivanov et al., , 2001Goldberger 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 powerspectrum exponent β = τ (q = 2) = −c 0 + 2c 1 − 2c 2 , and not to the full τ (q) spectrum required for multifractal diagnosis (Muzy et al., 1991(Muzy et al., , 1994Arneodo et al., 1995Arneodo et al., , 2008. Furthermore the fact that c 1 (∼ 1) is about one order of magnitude larger than the intermittency coefficient c 2 ( 0.1), explains why very much like c 1 (Figure 5A), the spectral exponent β ∼ −c 0 + 2c 1 ∼ 1 c 2 (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.