Martin Barrere
email: barrere@inria.fr

Rémi Badonnel
email: badonnel@inria.fr

Olivier Festor
email: festor@inria.fr

Martín Barrère

Towards Vulnerability Prevention in Autonomic Networks and Systems

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction and challenges

The continuous growth and dynamics of networks, as well as the diversification of their services in the context of Future Internet has considerably increased the complexity of their management. In order to face this problem, autonomic computing [START_REF] Kephart | The Vision of Autonomic Computing[END_REF], [START_REF]An Architectural Blueprint For Autonomic Computing[END_REF], has been introduced providing new perspectives. Highly inspired on the central nervous system, this approach aims to define a strong basis for automated systems capable of managing themselves in an autonomous manner, identifying four major properties, namely, self-configuration, self-optimization, self-healing and self-protection. Despite numerous benefits have been already obtained from this new paradigm, several challenges must be addressed in order to introduce this approach into current systems and networks.

When autonomic related tasks are performed, the environment is modified in order to achieve specific objectives. Such changes may lead to potential vulnerable states, thus change management techniques for assessing change associated risks are required [START_REF] Wickboldt | Improving IT Change Management Processes with Automated Risk Assessment[END_REF]. The vulnerability management activity usually consists in checking the configurations of the system components, identifying the presence of vulnerable states and performing the required maintenance operations (typically, modification of configuration parameters and/or application of security patches). Vulnerability detection and prevention techniques not also increase systems security but also complement the change management process by providing useful information for risk assessment mechanisms.

Even though vulnerability detection techniques have been proposed [START_REF] Wang | TaintScope: A Checksum-Aware Directed Fuzzing Tool for Automatic Software Vulnerability Detection[END_REF], [START_REF] Achi | Network Security Approach For Digital Forensics Analysis[END_REF], [START_REF] Khan | Enabling Self-Configuration in Autonomic Systems Using Case-Based Reasoning with Improved Efficiency[END_REF], and mechanisms for uniformly describing vulnerabilities and exchanging related information have been provided [START_REF]OVAL Language[END_REF], there is no integration of such mechanisms within the framework of autonomic networks and systems. Such integration constitutes the target of our work as we consider that autonomic environments should exploit the knowledge provided by vulnerability repositories in order to increase their security, stability and sustainability.

In this paper we present our approach for integrating vulnerability descriptions in the autonomic management plane, considering the OVAL [START_REF]OVAL Language[END_REF] process and the autonomic maintenance system Cfengine [1]. The remainder of this paper is organized as follows. Section 2 presents the proposed approach for increasing vulnerability awareness within autonomic environments, whereas the results achieved to date are outlined in Section 3. Section 4 presents conclusions and perspectives.

Self-configuration with vulnerability prevention

Within the autonomic computing field, the self-configuration property refers to the ability of networks and systems for automatically configuring themselves according to high-level policies. When autonomous networks and systems perform changes in order to be compliant with the specified policies, collateral effects can be introduced without explicit knowledge. Such unexpected effects can vary from internal malfunction to the exposure of vulnerable states.

We propose to support the self-configuration of autonomic systems with vulnerability management mechanisms. These mechanisms can ensure safe configurations and also reduce the probability of potential attacks and failures of the involved self-managed entities. Autonomic systems must be able to perform retro-inspection, identify required changes and execute the appropriate tasks. As happens in the real world, autonomic elements coexist within dynamic environments, interacting with other autonomic and non-autonomic elements. Nevertheless, such scenarios present continuous threats that may compromise autonomic elements safety. If an autonomic element is violated in some way, its functions and abilities become untrustworthy and eventually disabled; thus autonomic elements that use services of the former become compromised as well. This inevitably leads to distrust and the failure of the autonomic system. Autonomic systems must be able to manage their own state and perform the required activities to achieve secure configurations. Autonomic elements unable to support this capability will age with time, becoming more vulnerable, insecure and useless. Automation is really possible only if autonomic networks and systems are capable of ensuring safe configurations.

We therefore argue in favor of the integration of vulnerability descriptions into the management plane of autonomic systems. Our objective is to translate these vulnerability descriptions into policies that are interpretable by an autonomic system. In particular, we propose to translate standardized OVAL 1 vulnerability descriptions into Cfengine policy rules. In this manner, the vulnerability prevention process associated to OVAL can be integrated into Cfengine devices when maintenance operations are performed as depicted in Figure 1. The OVAL language is a standard XML-based language used by vendors and security organizations for publishing security related information warning about current threats and system vulnerabilities. OVAL repositories offer a wide range of security advisories that can be used for avoiding vulnerable states as well as augmenting networks and systems security considering best practices recommendations. Autonomic maintenance systems such as Cfengine provides support for automating the management of large-scale environments based on high levelpolicies. Cfengine offers a powerful distributed agent framework that combined with the OVAL vulnerability language, provides an efficient strategy for aligning security aspects on autonomic environments. In that context, we consider a translation module, identified on Figure 1, responsible for the generation of Cfengine rules corresponding to the vulnerability descriptions defined in the OVAL repository. Generated policies are deployed by the Cfengine server into its several Cfengine agents (points in the cloud) which are in charge of managing the devices present in the target network, in order to detect and prevent vulnerable configurations when self-management activities are performed. means for extending its translation capabilities. Current plugins provide a large variety of IOS vulnerability descriptions available within the official OVAL repository. The logical data model used by the translator is automatically generated based on the OVAL specification, thus enabling a seamless declarative evolution with the OVAL language and providing support for existing and future security related knowledge. We are currently finalizing the proposed prototype and performing several experiments over a simulated environment for evaluating factors such as functionality, performance, and quality of the generated Cfengine code.

Conclusions and perspectives

Vulnerability awareness is a key challenge in autonomic networks and systems. The autonomy of such dynamic environments will really be made possible when they fully integrate support mechanisms for preventing vulnerabilities and maintaining safe configurations. In that context, we propose an approach for integrating vulnerability descriptions into the autonomic management plane. We have formalized how these descriptions can be translated into policy rules that are interpretable by an autonomic configuration tool. We have developed a first prototype based on the Cfengine configuration tool, that covers a subset of OVAL definitions and permits to generate vulnerability alerts during the self-configuration activity. For future work, we are interested in extending the coverage of our solution to a larger variety of OVAL definitions, and in investigating further the execution of treatments by the autonomic network, when a vulnerable configuration is observed.

Fig. 1 :

 1 Fig. 1: Self-configuration with vulnerability prevention

Open Vulnerability Assessment Language

Preliminary results

We have formalized our approach by designing an intermediate formal language that provides mathematical tools for supporting the translation between OVAL advisories and Cfengine policies. In order to provide a computable infrastructure for the proposed approach, we have developed a first implementation prototype of the translation engine. At this early stage we have focused on dealing with IOS vulnerabilities over Cisco devices, nevertheless we aim to provide support for managing other platforms as well. The implementation prototype has been designed over a plugin-based architecture with the purpose of enabling easy


Ali Mohammad-Djafari

Li Wang

Nicolas Gac

Folkert Bleichrodt

A Student-t based sparsity enforcing hierarchical prior for linear inverse problems and its efficient Bayesian computation for 2D and 3D Computed Tomography

In many imaging systems and in particular in X ray Computed Tomography (CT) the reconstruction problem can be written as a linear inverse problem. In these problems, one property which can often be exploited is sparsity of the solution in an appropriate basis. In this work we consider the Student-t model in its hierarchical Normal-Inverse Gamma with an appropriate dictionary based coefficient. Then, thanks to the hierarchical generative model of the observation, we derive the expression of the joint posterior law of all the unknowns and an alternate optimisation algorithm for obtaining the joint MAP solution. We then detail the implementation issues of this algorithms for parallel computation and show the results on real size 2D and 3D phantoms.

Introduction

In many imaging systems, very often the reconstruction problem can be written as a linear inverse problem. This is the case in X ray Computed Tomography (CT) [START_REF] Mohammad-Djafari | Inverse Problems in Vision and 3D Tomography[END_REF][START_REF] Batenburg | 3d imaging of nanomaterials by discrete tomography[END_REF]. In these problems, the solution often has a sparse representation in a suitable basis. This sparsity property can be exploited in the reconstruction algorithm. For example, for a piecewise continuous or constant image, using Haar Transform (HT) gives rise to very sparse coefficients. To impose sparsity, three great categories of priors have been used: i) Generalized Gaussian (GG), ii) mixture of two Gaussian and iii) the heavy tailed probability density functions such as Cauchy and its generalisation a Student-t [START_REF] Baritaux | An efficient numerical method for general lp regularization in fluorescence molecular tomography[END_REF][START_REF] Luisier | Fast Haar-wavelet denoising of multidimensional fluorescence microscopy data[END_REF][START_REF] Mohammad-Djafari | Bayesian approach with prior models which enforce sparsity in signal and image processing[END_REF][START_REF] Leary | Compressed sensing electron tomography[END_REF].

In this work we consider Student-t model in its hierarchical Normal-Inverse Gamma with an appropriate dictionary based coefficient. Then, thanks to the hierarchical generative model of the observation, we derive the expression of the joint posterior law of all the unknowns and an alternate optimisation algorithm for obtaining the joint MAP solution.

We then detail the implementation issues of this algorithms for parallel computation and show the results on real size 2D and 3D phantoms.

Proposed model

We consider a linear model

g = Hf + ✏ ( 1
)
where f is the object to be reconstructed, H is the forward model, g the observed quantities and ✏ represents the errors.

We choose an appropriate basis or dictionary D for the object f in such a way that we can write

f = Dz + ⇠ (2)
and such that z is sparse. For imposing the sparsity of z we propose to use the generalized Student-t prior for its elements z j :

St(z j |↵, ) = Z 1 0 N (z j |, 0, v zj ) IG(v zj |↵, ) dv zj (3)
which gives us the possibility to write 8 < :

p(z j |v zj ) = N (z j |0, v zj ) / v zj 1/2 exp h 1 2 P j z 2 j vz j i p(v zj ) = IG(v zj |↵ z0 , z0 ) / v zj (↵z 0 +1) exp h z 0 vz j i (4)
and where p(z|v z ) =

Q j p(z j |v zj ), p(v z |↵ z0 , z0 ) = Q j p(v zj |↵ z0 , z0 ) and p(z, v z |↵ z0 , z0 ) = p(z|v z )p(v z |↵ z0 , z0 ).
In the forward model ( 1) and the prior model ( 2) we also use Student-t for the elements of ✏ i and ⇠ j which give rise to the following relations:

⇢ p(g|f , v ✏ ) = N (g|Hf , V ✏ ) with V ✏ = diag [v ✏ ] p(v ✏ ) = Q i IG(v ✏i |↵ ✏0 , ✏0 ) (5)
and

⇢ p(f |z, v ⇠ ) = N (f |Dz, V ⇠ ) with V ⇠ = diag [v ⇠ ] p(v ⇠ ) = Q j IG(v ⇠ j |↵ ⇠0 , ⇠0 ) (6)
With these prior models, we can show the graphical generative model of the data as follows: With these relations and equations, the expression of the posterior law of all the unknowns writes:

↵ z0 , z0 ? v z ✓⌘ ◆⇣ ? z ✓⌘ ◆⇣ ? D ↵ ✏0 , ✏0 ? v ✏ ✓⌘ ◆⇣ ? ↵ ⇠0 , ⇠0 ? v ⇠ ✓⌘ ◆⇣ ? ✓⌘ ◆⇣ ⇠ @ @ R f ✓⌘ ◆⇣ ✏ ✓⌘ ◆⇣ ? H ✓⌘ ◆⇣ g ( g = Hf + ✏ f = Dz + ⇠, z sparse 8 > > > < > > > : p(g|f , v ✏ ) = N (g|Hf , V ✏ ) p(f |z) = N (f |Dz, V ⇠ ) p(z|v z ) = N (z|0, V z ) 8 > > < > > : p(v ✏ ) = Q i IG(v ✏i |↵ ✏0 , ✏0 ) p(v z ) = Q j IG(v zj |↵ z0 , z0 ) p(v ⇠ ) = Q j IG(v ⇠ j |↵ ⇠0 , ⇠0 )
p(f , z, v ✏ , v z , v ⇠ |g) / p(g|f , v ✏ ) p(f |z, v ⇠ ) p(z|v z ) p(v ✏ ) p(v z ) p(v ⇠ ) (7)
iTWIST'16: international -Traveling Workshop on Interactions between Sparse models and Technology.

(48)

We can then use it for inferring these unknowns. Classically, there are two major estimators: JMAP and Posterior Means (PM). The JMAP writes:

( b f , b z, b v ✏ , b v z , b v ⇠ ) = arg max (f ,z,v✏,vz,v⇠) {p(f , z, v ✏ , v z , v ⇠ |g)}
(8) for which, the easiest optimisation is alternate optimisation. Hopefully, with the proposed hierarchical structure and conjugate priors, in any of these steps we have analytical solutions:

8 > > > > > > < > > > > > > : f arg min f n J(f ) = kg Hf k 2 V ✏ + kf Dzk 2 V ⇠ ) o z arg min z n J(z) = kf Dzk 2 V ⇠ + kzk 2 V z o v ✏i ( ✏0 + 1 2 (g i [Hf ] i ) 2 / (↵ ✏0 + 3/2) v ⇠ j ( ⇠0 + 1 2 (f j [Dz] j ) 2 / (↵ ⇠0 + 3/2) v zj ( z0 + 1 2 |z j | 2
) / (↵ z0 + 3/2) (9) As we can see, the main part of the algorithm is two optimizations of quadratic criteria which can also be written as

8 > > > > < > > > > : J(f ) = kg Hf k 2 V ✏ + kf Dzk 2 V ⇠ = kV 1/2 ✏ (g Hf )k 2 2 + kV 1/2 ⇠ (f Dz)k 2 2 J(z) = kf Dzk 2 V ⇠ + kzk 2 V z = kV 1/2 ⇠ (f Dz)k 2 2 + kV z 1/2 zk 2 2
(10) For implementation of the gradient based optimisation algorithms we need their gradients:

⇢ rJ(f ) = 2H 0 V ✏ 1 (g Hf ) + 2V ⇠ 1 (f Dz) rJ(z) = 2D 0 V ⇠ 1 (f Dz) + 2V z 1 z (11)
The algorithm we propose here which can be implemented efficiently and scale up for real applications and Big Data is as follows

8 > > > > < > > > > : f f f rJ(f ) z z z rJ(z) v ✏i ( ✏0 + 1 2 (g i [Hf ] i ) 2 / (↵ ✏0 + 3/2) v ⇠ j ( ⇠0 + 1 2 (f j [Dz] j ) 2 / (↵ ⇠0 + 3/2) v zj ( z0 + 1 2 |z j | 2 ) / (↵ z0 + 3/2) (12)
where f and z have to be adapted and updated at each iteration.

Simulation results

To show the effectiveness of the proposed method, we show here two examples of X ray image reconstruction from 32 projections. The images have 128x128 and 512x512 pixels and in both cases we simulated 32 projections uniformly distributed between 0 and 180 degrees. A Gaussian noise is added in such a way to have a SNR of 20dB. We will show examples of 3D reconstruction results in the final paper. These results are obtained with objects of volume 256x256x256 and 32 projections of size 256x256.

Conclusion

We proposed a hierarchical Normal-Inverse Gamma prior for modelling the sparsity of both the error terms of the dataforward model and the error terms of dictionary based decomposition of the unknown images. With these priors we obtain an expression for the joint posterior law of all the unknowns (image itself f , the coefficients of the decomposition z and their respective hidden variances v ✏ , v ⇠ and v z ). An approximate Bayesian Computation (ABC) based on the alternate optimization of this joint posterior law with respect its arguments gives an algorithm which can be implemented in an efficient way which can scales up for real applications.

Figure 1 :

 1 Figure 1: The graphical model of the generative forward model for linear inverse problems and hierarchical sparsity enforcing model.

Figure 2 :

 2 Figure 2: X ray CT: Originals, Projections, Initialization by Back-projection (BP), Reconstruction results.

iTWIST'16: international -Traveling Workshop on Interactions between Sparse models and Technology.(49)


Evgeniya Gerasimova

Benjamin Audit

Stephane G Roux

André Khalil

Olga Gileva

Françoise Argoul

Oleg Naimark

Alain Arneodo

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

Keywords: breast cancer, infrared thermography, wavelet transform modulus maxima method

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

It is widely recognized that early diagnosis is the key to breast cancer survival. X-ray mammography [START_REF] Nass | Mammography and Beyond: Developing Technologies for the Early Detection of Breast Cancer[END_REF][START_REF] Bronzino | Biomedical Engeneering Handbook[END_REF], 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 [START_REF] Jorgensen | Overdiagnosis in publicly organised mammography screening programmes : systematic review of incidence trends[END_REF][START_REF] Tamini | Circulating carotenoids, mammographic density, and subsequent risk of breast cancer[END_REF] and this despite the increasing use of computer-aided detection/diagnosis methods [START_REF] Fenton | Effectiveness of computer-aided detection in community mammography practice[END_REF]. Biopsy is indeed the only conclusive diagnostic test for breast cancer, but the number of unnecessary biopsies is still too high [START_REF] Vinitha Sree | The use of skin surface electropotentials for breast cancer detection. Preliminary clinical trial results obtained using the Biofield Diagnostic System[END_REF]. Since the original observation by [START_REF] Lawson | Implications of surface temperatures in the diagnosis of breast cancer[END_REF] 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 [START_REF] Yahara | Relationship between micro vessel density and thermographic hot areas in breast cancer[END_REF][START_REF] Bronzino | Biomedical Engeneering Handbook[END_REF], IR thermography has been considered as a promising non-invasive screening method of breast cancer [START_REF] Ng | A review of thermography as promising non-invasive detection modality for breast tumor[END_REF]. However the suitability of static IR imaging for routine screening has been severely questioned [START_REF] Head | Infrared imaging: making progress in fulfilling its medical promise[END_REF][START_REF] Bronzino | Biomedical Engeneering Handbook[END_REF], 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 [START_REF] Etehadtavakol | Breast thermography as a potential noncontact method in the early detection of cancer: a review[END_REF] 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 [START_REF] Theiler | Role of nitric oxide in tumour progression: lessons from human tumours[END_REF][START_REF] Anbar | Detection of cancerous breasts by dynamic area telethermometry[END_REF].

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 [START_REF] Button | Dynamic infrared imaging for the detection of malignancy[END_REF]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 [START_REF] Muzy | Wavelets and multifractal formalism for singular signals: application to turbulence data[END_REF][START_REF] Muzy | The multifractal formalism revisited with wavelets[END_REF][START_REF] Arneodo | A wavelet-based method for multifractal image analysis: from theoretical concepts to experimental applications[END_REF] 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.

METHODS OF ANALYSIS

THE WAVELET TRANSFORM

The wavelet transform (WT) is a mathematical microscope [START_REF] Muzy | Wavelets and multifractal formalism for singular signals: application to turbulence data[END_REF][START_REF] Muzy | The multifractal formalism revisited with wavelets[END_REF][START_REF] Arneodo | A wavelet-based method for multifractal image analysis: from theoretical concepts to experimental applications[END_REF] that is well suited for the analysis of complex non-stationary time-series such as those found in physiological systems [START_REF] Ivanov | Multifractality in human heartbeat dynamics[END_REF][START_REF] Goldberger | Fractal dynamics in physiology: alterations with disease and aging[END_REF], 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 [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF]:

W ψ [ ](t 0 , a) = 1 a +∞ -∞ (t)ψ t -t 0 a dt, ( 1
)
where t 0 is the time parameter and a (> 0) the scale parameter. By choosing a wavelet ψ whose n + 1 first moments are zero [ t n ψ(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 [START_REF] Muzy | Wavelets and multifractal formalism for singular signals: application to turbulence data[END_REF][START_REF] Muzy | The multifractal formalism revisited with wavelets[END_REF][START_REF] Arneodo | A wavelet-based method for multifractal image analysis: from theoretical concepts to experimental applications[END_REF]. 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) [START_REF] Muzy | Wavelets and multifractal formalism for singular signals: application to turbulence data[END_REF][START_REF] Muzy | The multifractal formalism revisited with wavelets[END_REF][START_REF] Arneodo | A wavelet-based method for multifractal image analysis: from theoretical concepts to experimental applications[END_REF].

THE WAVELET TRANSFORM MODULUS MAXIMA METHOD

The WT modulus maxima (WTMM) method was originally developed to generalize box-counting techniques [START_REF] Arneodo | Fractal dimensions and f (α) spectrum of the Henon attractor[END_REF] 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 [START_REF] Muzy | Wavelets and multifractal formalism for singular signals: application to turbulence data[END_REF][START_REF] Muzy | The multifractal formalism revisited with wavelets[END_REF][START_REF] Arneodo | A wavelet-based method for multifractal image analysis: from theoretical concepts to experimental applications[END_REF]. It has proved very efficient to estimate scaling exponents and multifractal spectra [START_REF] Muzy | The multifractal formalism revisited with wavelets[END_REF][START_REF] Delour | Intermittency of 1D velocity spatial profiles in turbulence: a magnitude cumulant analysis[END_REF][START_REF] Audit | Wavelet-based estimators of scaling behavior[END_REF]. This method has been generalized in 2D for the multifractal analysis of rough surfaces [START_REF] Arneodo | A wavelet-based method for multifractal image analysis: from theoretical concepts to experimental applications[END_REF] and then for the analysis of 3D scalar and vector fields [START_REF] Kestener | Generalizing the wavelet-based multifractal formalism to random vector fields: application to three-dimensional turbulence velocity and vorticity data[END_REF][START_REF] Arneodo | Wavelet-based multifractal analysis[END_REF]. It has been successfully applied in various areas of fundamental research (Arneodo et al., 1998a[START_REF] Arneodo | A wavelet based multifractal formalism: application to DNA sequences, satellite images of the cloud structure and stock market data[END_REF][START_REF] Arneodo | A wavelet-based method for multifractal image analysis: from theoretical concepts to experimental applications[END_REF][START_REF] Arneodo | Wavelet-based multifractal analysis[END_REF][START_REF] Khalil | Morphological analysis of H I features. II. Wavelet-based multifractal formalism[END_REF][START_REF] Roland | Revisiting the physical processes of vapodeposited thin gold films on chemically modified glass by atomic force and surface plasmon microscopies[END_REF][START_REF] Roux | Evidence for inherent nonlinearity in temporal rainfall[END_REF][START_REF] Arneodo | Multi-scale coding of genomic information: from DNA sequence to genome structure and function[END_REF]. In the context of the present study, the 1D WTMM method has proved very efficient at discriminating between healthy and sick heart beat dynamics [START_REF] Ivanov | Multifractality in human heartbeat dynamics[END_REF][START_REF] Ivanov | From 1/f noise to multifractal cascades in heartbeat dynamics[END_REF][START_REF] Goldberger | Fractal dynamics in physiology: alterations with disease and aging[END_REF], whereas the 2D WTMM method can be used to detect microcalcifications and has great potential to assist in cancer diagnosis from digitized mammograms [START_REF] Kestener | Wavelet-based multifractal formalism to assist in diagnosis in digitized mammograms[END_REF][START_REF] Arneodo | A wavelet-based method for multifractal image analysis: from theoretical concepts to experimental applications[END_REF].

The WT modulus maxima (WTMM) method [START_REF] Muzy | Wavelets and multifractal formalism for singular signals: application to turbulence data[END_REF][START_REF] Muzy | The multifractal formalism revisited with wavelets[END_REF][START_REF] Arneodo | A wavelet-based method for multifractal image analysis: from theoretical concepts to experimental applications[END_REF] 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 [START_REF] Muzy | Wavelets and multifractal formalism for singular signals: application to turbulence data[END_REF][START_REF] Muzy | The multifractal formalism revisited with wavelets[END_REF][START_REF] Arneodo | A wavelet-based method for multifractal image analysis: from theoretical concepts to experimental applications[END_REF] 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:

Z(q, a) = l∈L(a) |W(t, a)| q ∼ a τ (q) , ( 2
)
where q ∈ R. Then from the scaling function τ (q), D(h) is obtained by a Legendre transform [START_REF] Muzy | Wavelets and multifractal formalism for singular signals: application to turbulence data[END_REF][START_REF] Muzy | The multifractal formalism revisited with wavelets[END_REF][START_REF] Arneodo | A wavelet-based method for multifractal image analysis: from theoretical concepts to experimental applications[END_REF]:

D(h) = min q [qh -τ (q)]. (3)
As originally pointed out in [START_REF] Muzy | The multifractal formalism revisited with wavelets[END_REF]; [START_REF] Arneodo | A wavelet-based method for multifractal image analysis: from theoretical concepts to experimental applications[END_REF], one can avoid some practical difficulties that occur when directly performing the Legendre transform of τ (q), by computing the following expectation values:

h(q, a) = ∂ ∂q ln (Z(q, a)) = l∈L(a) ln |W ψ [ ](t, a)| Ŵψ [ ](q, l, a), D(q, a) = q ∂ ∂q ln (Z(q, a)) -Z(q, a) = l∈L(a)
Ŵψ [ ](q, l, a) ln Ŵψ [ ](q, l, a) , (4)

where

Ŵψ [ ](q, l, a) = |W ψ [ ](t, a)| q Z(q, a)
is the equivalent of Bolzmann weight in the analogy that links the multifractal formalism to thermodynamics [START_REF] Arneodo | A wavelet-based method for multifractal image analysis: from theoretical concepts to experimental applications[END_REF]. 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 [START_REF] Arneodo | A wavelet based multifractal formalism: application to DNA sequences, satellite images of the cloud structure and stock market data[END_REF]:

ρ W a /a H (w) = ρ(w), ( 5
)
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 [START_REF] Muzy | Wavelets and multifractal formalism for singular signals: application to turbulence data[END_REF][START_REF] Muzy | The multifractal formalism revisited with wavelets[END_REF][START_REF] Arneodo | A wavelet-based method for multifractal image analysis: from theoretical concepts to experimental applications[END_REF] that depends on t.

In this study, we fit the τ (q) data by the so-called log-normal quadratic approximation τ (q) = -c 0 + c 1 qc 2 q 2 /2, where the coefficients c n > 0 [START_REF] Delour | Intermittency of 1D velocity spatial profiles in turbulence: a magnitude cumulant analysis[END_REF]. The corresponding singularity spectrum has a characteristic single-humped shape D(h) = c 0 -(hc 1 ) 2 /2c 2 , where c 0 = -τ (0) is the fractal dimension of the support of singularities of , c 1 is the value of h that maximizes D(h) and c 2 , the so-called intermittency coefficient [START_REF] Muzy | Wavelets and multifractal formalism for singular signals: application to turbulence data[END_REF][START_REF] Muzy | The multifractal formalism revisited with wavelets[END_REF][START_REF] Arneodo | A wavelet-based method for multifractal image analysis: from theoretical concepts to experimental applications[END_REF][START_REF] Delour | Intermittency of 1D velocity spatial profiles in turbulence: a magnitude cumulant analysis[END_REF], characterizes the width of the D(h) spectrum, i.e., the signature of a change in WT coefficient statistics across scales [START_REF] Muzy | Wavelets and multifractal formalism for singular signals: application to turbulence data[END_REF][START_REF] Muzy | The multifractal formalism revisited with wavelets[END_REF][START_REF] Arneodo | A wavelet-based method for multifractal image analysis: from theoretical concepts to experimental applications[END_REF][START_REF] Arneodo | A wavelet based multifractal formalism: application to DNA sequences, satellite images of the cloud structure and stock market data[END_REF].

DESCRIPTION OF DATA

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 [START_REF] Gileva | Interdisciplinary approaches to early diagnosis and screening of tumors and precancerous diseases (for example, breast cancer)[END_REF]. 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 [START_REF] Gerasimova | Multifractal analysis of dynamic infrared imaging of breast cancer[END_REF].

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.

REMARK

As commonly done for noise signals [START_REF] Muzy | The multifractal formalism revisited with wavelets[END_REF][START_REF] Audit | Wavelet-based estimators of scaling behavior[END_REF] and previously experienced when applying the 1D WTMM method to rainfall time-series [START_REF] Venugopal | Revisiting multifractality of high-resolution temporal rainfall using a waveletbased formalism[END_REF], the wavelet analysis was performed on the cumulative (or integral) of the temperature time-series (Supplementary Figure 2) In (E,F), the τ (q) and D(h) spectra were averaged over a 8 × 8 pixel 2 square.

using the second-order compactly supported version ψ (2) (3) of the Mexican hat [START_REF] Roux | Detecting vorticity filaments using wavelet analysis: about the statistical contribution of vorticity filaments to intermittency in swirling turbulent flows[END_REF] (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.

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.

RESULTS

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) vs. log 2 a, in Figures 2B,C, where the slopes h(q) = c 1 = 0.81 ± 0.01 and D(q) = 0.99 ± 0.03, do not significantly (3) (Supplementary Figure 1).

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 [START_REF] Theiler | Role of nitric oxide in tumour progression: lessons from human tumours[END_REF][START_REF] Schreiber | Improved surrogate data for nonlinear tests[END_REF] 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 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.

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 1 values extend over a rather wide range 0.6 c 1 1.8 but turn out to be quite similar with mean values c1 = 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 2 has a definite discriminatory power. The histogram for cancerous breasts (c 2 = 0.045 ± 0.001) is definitely shifted toward smaller values relative to the ones for opposite (c 2 = 0.056 ± 0.001) and healthy breasts (c 2 = 0.058 ± 0.001) (Figure 5B). The small-value leftside of the c 2 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 2 = 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%) , 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). of (white) squares were removed from our analysis because of lack of scaling (Figure 5C, Table 2 andSupplementary 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).

A common way to suspect cancer by IR thermography is to look for some dissymmetry between the two breasts of a patient [START_REF] Etehadtavakol | Breast thermography as a potential noncontact method in the early detection of cancer: a review[END_REF]. 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 4 and Supplementary Figures 4567).

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

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 [START_REF] Jorgensen | Overdiagnosis in publicly organised mammography screening programmes : systematic review of incidence trends[END_REF]. 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 [START_REF] Ivanov | Multifractality in human heartbeat dynamics[END_REF][START_REF] Ivanov | From 1/f noise to multifractal cascades in heartbeat dynamics[END_REF][START_REF] Goldberger | Fractal dynamics in physiology: alterations with disease and aging[END_REF], 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 [START_REF] Muzy | Wavelets and multifractal formalism for singular signals: application to turbulence data[END_REF][START_REF] Muzy | The multifractal formalism revisited with wavelets[END_REF][START_REF] Arneodo | A wavelet-based method for multifractal image analysis: from theoretical concepts to experimental applications[END_REF][START_REF] Arneodo | Wavelet-based multifractal analysis[END_REF]. 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 [START_REF] Ng | Computer simulation in conjunction with medical thermography as an adjunct tool for early detection of breast cancer[END_REF][START_REF] Xu | Biothermomechanical behaviour of skin tissue[END_REF][START_REF] Lin | Detecting early breast tumour by finite element thermal analysis[END_REF]. 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 [START_REF] Theiler | Role of nitric oxide in tumour progression: lessons from human tumours[END_REF][START_REF] Anbar | Detection of cancerous breasts by dynamic area telethermometry[END_REF][START_REF] Button | Dynamic infrared imaging for the detection of malignancy[END_REF]Joro et al., 2008b). More likely it can be the signature of some architectural change in the cellular microenvironment of the breast tumor [START_REF] Bissell | Why don't we get more cancer? A proposed role of the microenvironment in restraining cancer progression[END_REF] that may deeply affect heat transfer and related thermomechanics in breast tissue [START_REF] Xu | Biothermomechanical behaviour of skin tissue[END_REF][START_REF] Quail | Microenvironmental regulation of tumor progression and metastasis[END_REF]. 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.

FIGURE 1 |

 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 2 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 2 Z (q, a) vs. log 2 a over a range of time-scales [0.3, 3] s (Figure2A). (F) D(h) vs. h (Figures2B,C). The solid lines in (E,F) correspond to quadratic spectra (see text) with parameters [c 0 , c 1 , c 2 ] = [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 2 square.

FIGURE 2 |

 2 FIGURE 2 | Multifractal analysis of cumulative IR temperature time-series in a 8 × 8 pixel 2 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 2 Z (q, a) vs. log 2 (a) Equation (2). (B) h(q, a) vs. log 2 a Equation (4). (C) D(q, a) vs. log 2 a

FIGURE 3 |

 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 0 ) c1 ) of cumulative temperature of (A,A ) the cancerous right breast of patient 20: c 1 = 0.81, (B,B ) the opposite left breast of patient 20:

FIGURE 4 |

 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 2 squares covering the entire breast are colored coded (A-C) and represented as a dot in the (c 1 , c 2 ) plane (D-F). The colors have the following meaning: red: c 2 < 0.03, blue:

FIGURE 5 |

 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 1 values obtained in the breasts with a malignant tumor (red: N = 4032 8 × 8 pixel 2 squares), the opposite breasts (black: N = 3606) and the healthy breasts (green: N = 3185). (B) Normalized histograms of c 2 values. (C) Percentage of "monofractal" (red: c 2 < 0.03), "multifractal" (blue: c 2 0.03) and "no scaling" (black) 8 × 8 pixel 2 squares in cancerous, opposite and healthy breasts. (D) Percentage of "monofractal" 8 × 8 pixel 2 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 Table2).

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.

 1

		Age Breast Stage	Size (cm)	Depth (cm) Location	Histological DS
	1	L	IIIa	1.57	4	Border outer-inner upper	Invasive ductal cancer
						quadrants
	2	R	IIa	2.8	3	Border outer-inner upper	Invasive ductal cancer
						quadrants
	3	R	IIa	2.5 × 2.85	5	Upper outer quadrant	Invasive ductal cancer
	4	L	IIa	2.6 × 3.3	2	Border lower-upper outer	Cystadenocarcinoma
						quadrant
	5	L	IIa	4.15 × 3.45	1	Border inner-outer upper	Cystadenocarcinoma
						quadrants
	6	L	III	4.5	4	Upper outer quadrant	Invasive lobular cancer
	7	R	IIa	2.7 × 2.2	3	Border inner-outer upper	Invasive lobular cancer
						quadrants
	8	L	III	1.6 × 1.8	3	Upper inner quadrant	Invasive ductal cancer
	9	L	IIa	1.8 × 3	6	Upper outer quadrant	Invasive ductal cancer
	10	R	IIa	-	-		Invasive ductal cancer
	11	R	IIb	5.3	1	Upper quadrant	Ductal cancer
	12	R	Ib	1.8 and 1.26	12	Border between upper quadrants	Invasive ductal and lobular cancer
	13	L	IIa	1.5	1	Underarm area	Invasive lobular cancer
	14	R	I	1.2	5	Border between upper quadrants	Invasive ductal cancer
	15	R	I	-	-		Invasive ductal cancer
	16	R	IIa	3.4	7	Upper inner quadrant	Invasive lobular cancer
	17	R	IIa	3 × 1.5 × 1.4	4	Upper outer quadrant	Invasive lobular and ductal cancer
	18	L	IIa	3.49 × 2.39	6	Upper outer quadrant	Invasive ductal cancer
	19	L	-	-	-		Invasive ductal cancer
	20	R	IIa	3.1 × 2.5 × 3.6	3	Upper outer quadrant	Invasive lobular and ductal cancer
	21	R	IIIb	6.5 × 4	2	Upper outer qadrant	Invasive cancer
	22	R	IIa	-	-		Invasive ductal cancer
	23	L	-	2.9 × 3.4	1	Border between upper quadrants	Invasive lobular cancer
	24	L	IIa	1.7 × 2.2	2	Border between outer quadrants	Paget disease of the nipple
	25	L	IIa	2.5 and 1.8 × 2, 2.7 × 1.7	3	Border between upper quadrants,	Invasive ductal cancer
						upper outer quadrant
	26	R	Susp. ms	1.9	1	Underarm area	No metastasis
	27	L	IIa	2.2 × 2.1	2	Border between inner quadrants	Invasive ductal and lobular cancer
	28	L	IIb	3.49	8	Border between lower quadrants	Invasive ductal cancer
	29	R	II	2.1 × 1.6	6	Border between lower quadrants	Invasive ductal cancer
	30	R	IIa	3.5	3	Upper outer quadrant	Invasive ductal cancer
	31	R	IIb	-	-		Invasive ductal cancer
	32	R	IIb	4.4 × 3 × 2.5	9	Border between upper quadrants	Invasive ductal cancer
	33	R	IIa	1.57	1.5	Border Upper-lower outer	Invasive lobular cancer
						quadrants, close to the nipple

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.

 2

			Cancerous breast			Opposite breast
		N mono	N multi	N noscaling	Total	N mono	N multi	N noscaling	Total
	1	2 (2.4%)	20 (24.4%)	60 (73.2%)	82 (100%)	3 (4.7%)	46 (71.9%)	15 (23.4%)	64 (100%)
	2	39 (42.9%)	23 (25.3%)	29 (31.9%)	91 (100%)	5 (4.1%)	61 (50.0%)	56 (45.9%)	122 (100%)
	3	16 (16.3%)	60 (61.2%)	22 (22.4%)	98 (100%)	31 (25.8%)	77 (64.2%)	12 (10.0%)	120 (100%)
	4	56 (33.1%)	90 (53.3%)	23 (13.6%)	169 (100%)	7 (8.3%)	62 (73.8%)	15 (17.9%)	84 (100%)
	5	1 (0.8%)	108 (82.4%)	22 (16.8%)	131 (100%)	0 (0.0%)	122 (91.7%)	11 (8.3%)	133 (100%)
	6	54 (30.5%)	101 (57.1%)	22 (12.4%)	177 (100%)
	7	44 (32.4%)	59 (43.4%)	33 (24.3%)	136 (100%)	26 (23.0%)	75 (66.4%)	12 (10.6%)	113 (100%)
	8	20 (27.0%)	46 (62.2%)	8 (10.8%)	74 (100%)	1 (1.4%)	29 (39.2%)	44 (59.5%)	74 (100%)
	9	65 (40.1%)	85 (52.5%)	12 (7.4%)	162 (100%)	39 (25.7%)	37 (24.3%)	76 (50.0%)	152 (100%)
	10	77 (51.3%)	42 (28.0%)	31 (20.7%)	150 (100%)	40 (24.2%)	96 (58.2%)	29 (17.6%)	165 (100%)
	11	62 (46.3%)	54 (40.3%)	18 (13.4%)	134 (100%)	23 (17.6%)	87 (66.4%)	21 (16.0%)	131 (100%)
	12	5 (4.9%)	97 (94.2%)	1 (1.0%)	103 (100%)	1 (0.9%)	114 (99.1%)	0 (0.0%)	115 (100%)
	13	48 (37.8%)	18 (14.2%)	61 (48.0%)	127 (100%)	50 (40.0%)	18 (14.4%)	57 (45.6%)	125 (100%)
	14	17 (18.9%)	65 (72.2%)	8 (8.9%)	90 (100%)	8 (8.9%)	72 (80.0%)	10 (11.1%)	90 (100%)
	15	1 (1.3%)	70 (90.9%)	6 (7.8%)	77 (100%)	0 (0.0%)	26 (89.7%)	3 (10.3%)	29 (100%)
	16	0 (0.0%)	79 (54.1%)	67 (45.9%)	146 (100%)	52 (29.4%)	111 (62.7%)	14 (7.9%)	177 (100%)
	17	106 (57.6%)	64 (34.8%)	14 (7.6%)	184 (100%)	3 (2.1%)	104 (73.8%)	34 (24.1%)	141 (100%)
	18	11 (7.7%)	111 (77.6%)	21 (14.7%)	143 (100%)	0 (0.0%)	111 (93.3%)	8 (6.7%)	119 (100%)
	19	60 (34.1%)	106 (60.2%)	10 (5.7%)	176 (100%)	29 (19.1%)	113 (74.3%)	10 (6.6%)	152 (100%)
	20	76 (49.7%)	66 (43.1%)	11 (7.2%)	153 (100%)	11 (7.7%)	127 (89.4%)	4 (2.8%)	142 (100%)
	21	25 (14.0%)	148 (82.7%)	6 (3.4%)	179 (100%)	17 (9.3%)	156 (85.7%)	9 (4.9%)	182 (100%)
	22	50 (29.8%)	41 (24.4%)	77 (45.8%)	168 (100%)	86 (44.1%)	82 (42.1%)	27 (13.8%)	195 (100%)
	23	7 (4.6%)	143 (93.5%)	3 (2.0%)	153 (100%)	27 (16.2%)	111 (66.5%)	29 (17.4%)	167 (100%)
	24	8 (4.9%)	144 (88.3%)	11 (6.7%)	163 (100%)	0 (0.0%)	163 (94.2%)	10 (5.8%)	173 (100%)
	25	77 (57.0%)	55 (40.7%)	3 (2.2%)	135 (100%)	24 (18.0%)	101 (75.9%)	8 (6.0%)	133 (100%)
	26	87 (47.3%)	31 (16.8%)	66 (35.9%)	184 (100%)	18 (14.4%)	37 (29.6%)	70 (56.0%)	125 (100%)
	27	7 (4.3%)	148 (91.4%)	7 (4.3%)	162 (100%)	0 (0.0%)	139 (91.4%)	13 (8.6%)	152 (100%)
	28	11 (6.6%)	140 (84.3%)	15 (9.0%)	166 (100%)	27 (16.5%)	84 (51.2%)	53 (32.3%)	164 (100%)
	29	65 (38.9%)	67 (40.1%)	35 (21.0%)	167 (100%)	7 (4.0%)	144 (82.8%)	23 (13.2%)	174 (100%)
	30	116 (70.3%)	40 (24.2%)	9 (5.5%)	165 (100%)	65 (42.8%)	72 (47.4%)	15 (9.9%)	152 (100%)
	31	13 (6.8%)	175 (91.1%)	4 (2.1%)	192 (100%)	1 (0.6%)	140 (82.8%)	28 (16.6%)	169 (100%)
	32	88 (48.1%)	90 (49.2%)	5 (2.7%)	183 (100%)	11 (5.5%)	160 (80.0%)	29 (14.5%)	200 (100%)
	33	23 (17.0%)	109 (80.7%)	3 (2.2%)	135 (100%)	5 (3.6%)	112 (80.0%)	23 (16.4%)	140 (100%)
	Number and percentage of monofractal (c 2 < 0.03), multifractal (c 2	0.03), no-scaling (and total) 8 × 8 pixel 2 squares in the cancerous breast and in the opposite
	breast (see Figure

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

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.


Grégoire Burgé

Florian Chemarin

Claire Saulou-Berion

Henry-Eric Spinnler

Violaine Athes-Dutour

Marwen Moussa

Biobased 3-hydroxypropionic acid through a new integrated process of glycerol bioconversion and membrane-assisted reactive extraction

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Material and methods

Integrated process Principle of hollow fiber membrane contactor

-

  Chemarin a,b , Claire Saulou-Bérion a , Henry-Eric Spinnler a , Violaine Athès a , Marwen Moussa *a a) UMR GMPA, AgroParisTech, INRA, Université Paris-Saclay, 78850, Thiverval-Grignon, France b) Chaire Agro-Biotechnologies Industrielles (ABI) -AgroParisTech, 247 rue Paul Vaillant Couturier, F-51100 Reims, France * email: marwen.moussa@agroparistech.fr Tremendous growth of biodiesel manufacturing industries  glycerol as a main bypropduct Development of biotechnological processes to convert glycerol into high-added value chemicals  3-Hydroxypropionic acid (3-HP): significant platform chemical from which various specialty chemicals can be synthesized (Werpy and Petersen, 2004) Currently produced by chemical methods, but biotechnological production not well established  Until now, among lactic acid bacteria, only bacteria of the Lactobacillus genus and specially L. reuteri have been shown to produce 3-HP from glycerol, although at low productivity  3-HP and its metabolic intermediate 3-hydroxypropionaldehyde (3-HPA) are suspected to exhibit inhibitory or toxic effects on the producing microorganisms (Burgé et al., 2015) ISPR (In Situ Product Recovery) = potential strategy to relieve the stress, increase the performance of microbial cells and recover the molecule of interest  Study of the impact of the integrated process on bioconversion productivity and cell physiological state -Higher glycerol consumption and 3-HPA production with increasing the pH and 3-HP/1,3-PDO molar ratio > 1 -Lower impact on cell physiological state with increasing the pH Biobased 3-hydroxypropionic acid through a new integrated process of glycerol bioconversion and membrane-assisted reactive extraction  Distribution coefficient: integrated process coupling glycerol bioconversion and 3-HP reactive extraction assisted by hollow fiber membrane contactor Bioconversion at pH 5, 250 rpm, 5 x 10 9 cells/mL Reactive extraction with TOA 20 % (v/v) + decanol 80 % (v/v) cell circulation -Toxic effect of 3-HP and 3-HPA Harmful impact of Aliquat 336 compared to extractant phase with only TOA / decanol -Better 3-HP extraction at low pH but possible and favorable in all the range of pH tested (K D > 1) -Lower 3-HP extraction from the real bioconversion medium than from the solutions of pure 3-HP in water Effect of the soluble molecules (proteins, phospholipids, salts) -Good selectivity of the reactive extraction (Moussa et al., 2015; Burgé et al., 2016) growth 5 L of MRS medium + 20 g/L of glucose, 37°C, anaerobic conditions, pH 6 regulated with KOH (10 N) -Harvesting and washing 3 x Centrifugation 10 min, 5000 g, 4°C, in Potassium Phosphate Buffer, pH 6.5 Glycerol bioconversion into 3-HP 2.5 L of Glycerol (18 g/L in distilled water), 37°C, microaerobic conditions, 5 x 10 9 Cells/mL, pH regulated with KOH (10 N) and HCl (5 N) -In situ product recovery of 3-HP Liquid-liquid reactive extraction assisted by membrane contactor in an organic phase containing decanol and 20 % v/v amines (TOA with or without Aliquat 336), 25 °C -3-HP recovery in aqueous phase Back-extraction of the 3-HP in aqueous phase Petersen G (2004) Top value added chemicals from biomass, vol 1: results of screening for potential candidates from sugars and synthesis gas. US Department of Energy. Burgé G, Saulou-Bérion C, Moussa M, Pollet B, Flourat A, Allais F, Athès V, Spinnler H E (2015), Diversity of Lactobacillus reuteri in converting glycerol into 3-hydroxypropionic acid. Applied Biochemistry and Biotechnology. Moussa M, Burgé G, Chemarin F, Bounader R, Saulou-Bérion C, Allais F, Spinnler H E, Athès V (2015), Reactive extraction of 3-hydroxypropionic acid from model aqueous solutions and real bioconversion media. Comparison with its isomer 2-hydroxypropionic (lactic) acid. Journal of Chemical Technology and Biotechnology Burgé G, Chemarin F, Moussa M, Saulou-Bérion C, Allais F, Spinnler H E, Athès V (2016), Reactive extraction of bio-based 3-hydroxypropionic acid assisted by hollow-fiber membrane contactor using TOA and Aliquat 336 in n-decanol. Journal of Chemical Technology and Biotechnology


Ispra (VA); Italy;dejanbrkic0611@gmail.com,dejan.brkic@jrc.ec.europa.eu This discussion communicates comparisons, practical applications and accuracy check related to the closed-form solutions of the Colebrook-White equation [START_REF] Colebrook | Turbulent flow in pipes with particular reference to the transition region between the smooth and rough pipe laws[END_REF][START_REF] Colebrook | Experiments with fluid friction in roughened pipes[END_REF] showed in the discussed paper and provides comments and references to already available similar works.

The first exact closed-form solution of the Colebrook-White equation in terms of Lambert Wfunction appears in [START_REF] Keady | Colebrook-White formula for pipe flow[END_REF] [Eq. ( 9) of the discussed paper]. [START_REF] Sonnad | Constraints for using Lambert W function-based explicit Colebrook-White equation[END_REF], confirmed also by Brkić (2012a), show constraints for the use of Keady's equation because it contains exponential form as the argument of the Lambert W-function (as defined in [START_REF] Hayes | Why W?[END_REF], "W" represents the Lambert function). Very high numerical values caused by exponential form in Keady's case very often make impossible calculation using computers (these values are too high to be processed by registers of memory and processor unit). These cases are associated in engineering practice with high speed of flow through very rough pipes (higher numerical values of the Reynolds (R) and higher values the relative roughness of inner pipe surface (ε/D) occurred simultaneously in pairs). Rollman and Spindler (2015), referred to general properties of the Lambert W-function [START_REF] Corless | On the Lambert Wfunction[END_REF], give their solution in order to avoid such exponential form in the transformed Colebrook-White equation (similar as Eq. ( 14) of the discussed paper). To overwhelm such inconvenience Brkić (2011a) also transforms the Colebrook-White equation using Lambert Wfunction avoiding exponential term (1):

                                                    D D R R W f R W 71 .

 

(1)

To develop (1), Brkić (2011a) use limiting case of the Colebrook-White equation under smooth pipe law of flow as showed in Appendix I of the discussed paper, where ε/D→0 [START_REF] Goudar | Explicit friction factor correlation for turbulent flow in smooth pipes[END_REF][START_REF] Sonnad | Turbulent flow friction factor calculation using a mathematically exact alternative to the Colebrook-White equation[END_REF][START_REF] Brkić | New explicit correlations for turbulent flow friction factor[END_REF][START_REF] Brkić | Can pipes be actually really that smooth?[END_REF]. To make Eq. ( 1) more applicable for engineering practice, [START_REF] Brkić | An explicit approximation of the Colebrook equation for fluid flow friction factor[END_REF] replaces Lambert W-function with approximate calculus. Further for the purpose of this discussion, numerical values of the parameters from this explicit approximation are optimized using genetic algorithms developed by [START_REF] Ćojbašić | Very accurate explicit approximations for calculation of the Colebrook friction factor[END_REF] where the relative error, δ(%)=(|f-f0|/f0)•100%, decreases to 1.28% (before optimization it was 2.2% and 3.16%, respectively referred to Eq. ( 2)): (2)

                                                 R . R . . R A D ε . . D ε . R A . . f A .
Regarding soft computation techniques such as optimization through genetic algorithms, it should be noted that the Colebrook-White equation can be simulated very accurately using Artificial Neural Networks as showed in Brkić and Ćojbašić (2016).

The nth formula is developed (Eq. ( 26) of the discussed paper) using similar approach as in Brkić (2011aBrkić ( , 2012a)). In this case, the Boyd's "shifted" Lambert W-function [START_REF] Boyd | Global approximations to the principal real-valued branch of the Lambert Wfunction[END_REF]) is used which is noted as the Y function the discussed paper. In addition to the solution based on Boyd's function, Brkić (2011aBrkić ( , 2012a,c) ,c) presents some further solutions based on works of [START_REF] Barry | Analytical approximations for real values of the Lambert W-function[END_REF] and [START_REF] Winitzki | Uniformapproximations for transcendental functions[END_REF]. Also, Brkić (2011aBrkić ( , 2012a,c) ,c) uses the series expansion of the Lambert W-function in a similar way as it is done for of Y function (Eq. ( 44) of the discussed paper).

Approximations of the Colebrook-White equation based on the nth formula, Eq (21) of the discussed paper, are given by (3).

                                                                                                                                                   D R R x x x x x R f x x x x R f x x x R f x x R f n n n n n 71 . 3
1 1 1 1 1 1 1 1 1 1 1 3 1 1 1 2 1 1 1        (3)
Using methodology from [START_REF] Brkić | Review of explicit approximations to the Colebrook relation for flow friction[END_REF], the maximal percentage relative error is evaluated. It is 1.1877% for n=1, 0.1826% for n=2, 0.0278% for n=3, 0.004249% for n=4, 6.48•10 -4 % for n=4, 9.89•10 - 5 % for n=5, etc. The error aggregates in the zone of small values of the Reynolds number (R) and relative roughness (ε/D) and it decreases with geometric progression. On the other hand, these approximations are computationally demanded since they contain many logarithmic terms [START_REF] Clamond | Efficient resolution of the Colebrook equation[END_REF][START_REF] Giustolisi | Some explicit formulations of Colebrook-White friction factor considering accuracy vs. computational speed[END_REF].

Notation

The following symbols are used in this discussion:

  Discussion of "Exact Analytical Solutions of the Colebrook-White Equation" by Yozo Mikata and Walter S. Walczak, J. Hydraul. Eng. 04015050; doi. 10.1061/(ASCE)HY.1943-7900.0001074 Dejan Brkić, PhD, Research Fellow; European Commission; DG Joint Research Center (JRC); Institute for Energy and Transport (IET); Energy Security, Systems and Market Unit; Via Enrico Fermi 2749;

R

  -Reynolds number (dimensionless) ε/D -Relative roughness of inner pipe surface (dimensionless) f -Darcy (Moody) flow friction factor (dimensionless) "=" -exactly equivalent to the Colebrook-White equation "≈"-approximately equivalent to the Colebrook-White equation f0 -Darcy (Moody) flow friction factor (dimensionless); obtained from the Colebrook-White equation using iterative procedure and hence treated as accurate δ(%)=(|f-f0|/f0)•100% -relative error (%) W -Lambert function X1-parameter defined by Eq. (14) of the discussed paper (here in Eq. 3)A-auxiliary term used in Eq. 2 'ln' denotes the natural log function Disclaimer: Neither the European Commission nor any person acting on behalf of it is responsible for the use which might be made of this publication.

Answers to the reviewer's comment:

The text is thoroughly checked and grammatical errors are hopefully corrected.