Anne-Sophie Vignion-Dewalle
email: anne-sophie.vignion@inserm.fr

Grégory Baert

Elise Thecua

Claire Vicentini

Laurent Mortier

Serge Mordon

Gregory Baert

Photodynamic therapy for actinic keratosis: Is the European consensus protocol for daylight PDT superior to conventional protocol for Aktilite CL 128 PDT?

Keywords: Photodynamic therapy, local damage comparison, mathematical modeling, red light, daylight, Aktilite CL 128. I

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

Photodynamic therapy for actinic keratosis: Is the European consensus protocol for daylight PDT superior to conventional protocol for Aktilite CL 128 PDT? Introduction

Photodynamic therapy (PDT) is a cancer treatment combining a light of an appropriate wavelength, a photosensitizer (PS), and sufficient molecular oxygen to generate reactive oxygen species and destroy (pre-) malignant cells [START_REF] Plaetzer | Photophysics and photochemistry of photodynamic therapy: fundamental aspects[END_REF]. Over the last 15 years, topical PDT using 5-aminolevulinic acid (ALA) (ALA-PDT) or methyl aminolevulinate (MAL) (MAL-PDT) has proven to be successful in the treatment of various dermatological conditions including actinic keratoses (AK) [START_REF] Bissonette | Large surface photodynamic therapy with aminolevulinic acid: treatment of actinic keratoses and beyond[END_REF][START_REF] Braathen | Guidelines on the use of photodynamic therapy for nonmelanoma skin cancer: an international consensus[END_REF][START_REF] Morton | Guidelines for topical photodynamic therapy: update[END_REF][START_REF] Morton | European Dermatology Forum Guidelines on topical photodynamic therapy[END_REF][START_REF] Wiegell | Update on photodynamic treatment for actinic keratosis[END_REF]. Topical administration of ALA or MAL leads to the selective accumulation of the photosensitizer protoporphyrin IX (PpIX) in the AK lesions and subsequent light irradiation leads to the destruction of the lesions.

Red light irradiation with the Aktilite CL 128 lamp (Galderma SA, Lausanne, Switzerland) using a total light dose of 37 J/cm 2 after three hours of incubation with MAL, is a conventional protocol that is approved and widely used in Europe for the PDT treatment of AK. This protocol, referred to in this paper as the conventional protocol for Aktilite CL 128 PDT, has been demonstrated to be an effective treatment with similar efficacy and better cosmetic results compared with standard therapies [START_REF] Morton | Intraindividual, right-left comparison of topical methyl aminolaevulinate-photodynamic therapy and cryotherapy in subjects with actinic keratoses: a multicentre, randomized controlled study[END_REF]. However, due in particular to high pain scores reported by patient during the treatment [START_REF] Tyrrell | The effect of air cooling pain relief on protoporphyrin IX photobleaching and clinical efficacy during dermatological photodynamic therapy[END_REF] and high room occupancy for dermatologists, many other protocols for PDT involving either shorter incubation times [START_REF] Braathen | Short incubation with methyl aminolevulinate for photodynamic therapy of actinic keratoses[END_REF], lower fluence rates [START_REF] Apalla | The impact of different fluence rates on pain and clinical outcome in patients with actinic keratoses treated with photodynamic therapy[END_REF], irradiation with other light sources [START_REF] Wiegell | Continuous activation of PpIX by daylight is as effective as and less painful than conventional photodynamic therapy for actinic keratoses; a randomized, controlled, single-blinded study[END_REF]… have been proposed for dermatological PDT.

Among these alternative PDT protocols, several protocols involving exposure with daylight instead of irradiation with the Aktilite CL 128 lamp have been investigated over the last decade [START_REF] Wiegell | Continuous activation of PpIX by daylight is as effective as and less painful than conventional photodynamic therapy for actinic keratoses; a randomized, controlled, single-blinded study[END_REF][START_REF] Wiegell | Photodynamic therapy of actinic keratoses with 8% and 16% methyl aminolaevulinate and home-based daylight exposure: a double-blinded randomized clinical trial[END_REF][START_REF] Wiegell | A randomized, multicentre study of directed daylight exposure times of 1½ vs. 2½ h in daylight-mediated photodynamic therapy with methyl aminolaevulinate in patients with multiple thin actinic keratoses of the face and scalp[END_REF][START_REF] Wiegell | Daylight photodynamic therapy for actinic keratosis: an international consensus: International Society for Photodynamic Therapy in Dermatology[END_REF][START_REF] Wiegell | Daylight-mediated photodynamic therapy of moderate to thick actinic keratoses of the face and scalp: a randomized multicentre study[END_REF][START_REF] Rubel | Daylight photodynamic therapy with methyl aminolevulinate cream as a convenient, similarly effective, nearly painless alternative to conventional photodynamic therapy in actinic keratosis treatment: a randomized controlled trial[END_REF]. Most of these protocols involve an incubation with MAL for a maximum of 30 minutes followed by a daylight exposure for between 1.5 and 2.5 hours. From an European consensus [START_REF] Wiegell | Daylight photodynamic therapy for actinic keratosis: an international consensus: International Society for Photodynamic Therapy in Dermatology[END_REF], using a two hours daylight exposure within 30 minutes after MAL application leads to a protocol, hereinafter referred to as the European consensus protocol for daylight PDT, as effective as and better tolerated by patients than the conventional protocol for Aktilite CL 128 PDT. The European consensus protocol for daylight PDT, more manageable in clinical practice than the conventional protocol for Aktilite CL 128 PDT [START_REF] Wiegell | Continuous activation of PpIX by daylight is as effective as and less painful than conventional photodynamic therapy for actinic keratoses; a randomized, controlled, single-blinded study[END_REF], has been recently approved in Europe for the treatment of thin, non-hyperkeratotic AK [START_REF] Cantisani | MAL Daylight Photodynamic Therapy for Actinic Keratosis: Clinical and Imaging Evaluation by 3D Camera[END_REF]. Based on the two comparative clinical studies between a protocol for daylight PDT and the conventional protocol for Aktilite CL 128 PDT [START_REF] Wiegell | Continuous activation of PpIX by daylight is as effective as and less painful than conventional photodynamic therapy for actinic keratoses; a randomized, controlled, single-blinded study[END_REF][START_REF] Rubel | Daylight photodynamic therapy with methyl aminolevulinate cream as a convenient, similarly effective, nearly painless alternative to conventional photodynamic therapy in actinic keratosis treatment: a randomized controlled trial[END_REF], we set the incubation time for the European consensus protocol for daylight PDT to 30 minutes.

In this paper, we propose to compare the efficiency of the conventional protocol for Aktilite CL 128 PDT (light source: Aktilite CL 128, incubation time: three hours, light dose: 37 J/cm 2 ) to the one of the European consensus protocol for daylight PDT (light source: daylight, incubation time: 30 minutes, treatment duration: two hours) through a mathematical modeling already published in our previous works [START_REF] Vignion-Dewalle | Comparison of three light doses in the photodynamic treatment of actinic keratosis using mathematical modeling[END_REF][START_REF] Vignion-Dewalle | Red light photodynamic therapy for actinic keratosis using 37 J/cm 2 : fractionated irradiation with 12.3 mW/cm 2 after 30 minutes incubation time compared to standard continuous irradiation with 75 mW/cm 2 after three hours incubation time using a mathematical modeling[END_REF]. Two weather conditions for the European consensus protocol for daylight PDT have been considered for this comparison: a clear sunny day and an overcast day. The mathematical modeling that involves a logistic model for both the biological clearance of PpIX and the conversion of MAL into PpIX and an analytic model for the PpIX photobleaching enables the local damage induced by the therapy to be estimated [START_REF] Vignion-Dewalle | Red light photodynamic therapy for actinic keratosis using 37 J/cm 2 : fractionated irradiation with 12.3 mW/cm 2 after 30 minutes incubation time compared to standard continuous irradiation with 75 mW/cm 2 after three hours incubation time using a mathematical modeling[END_REF].

A detailed description of both the conventional protocol for Aktilite CL 128 PDT and the European consensus protocol for daylight PDT is presented in Section II while Section III describes the mathematical modeling enabling the quantification of the PDT local damage.

The PDT local damages achieved by the conventional protocol for Aktilite CL 128 PDT, the European consensus protocol for daylight PDT during a sunny day and the European consensus protocol for daylight PDT during an overcast day are compared in Section IV.

Based on the comparison results, some discussions and conclusions are drawn in Section V.

II. The two topical PDT protocols considered in this paper

Incubation time here refers to the time elapsed between the application of MAL cream and the beginning of the light irradiation.

A. The conventional protocol for Aktilite CL 128 PDT

The Aktilite CL 128 lamp (Galderma SA, Lausanne, Switzerland) is the most widely used device system for topical PDT in Europe. Equipped with 128 light emitting diodes (LEDs) arranged in a 8×16 array, this lamp emits red light with a fluence rate of 70-100 mW/cm 2 (at a distance from 5 cm to 8 cm) [20], a peak wavelength of 632 nm [START_REF] Moseley | Light distribution and calibration of commercial PDT LED arrays[END_REF] and a full width at half maximum (FWHM) of approximately 19 nm [START_REF] Moseley | Light distribution and calibration of commercial PDT LED arrays[END_REF]. The fluence can be adjusted at the control panel of the lamp and the irradiation time is calculated automatically accordingly: for the 37 J/cm 2 as recommended for dermatological MAL-PDT, the irradiation time varies between 6 and 10 minutes. Moreover, with an irradiation field size of 8×18 cm, the Aktilite CL 128 lamp enables quite large fields to be treated.

Sequentially involving gentle curettage of the lesions, MAL cream application with occlusive dressing to each lesion and 5 mm of surrounding tissue for an incubation time of three hours, removal of excess MAL cream, positioning of the Aktilite CL 128 lamp head 5-8 cm over the area to be treated and irradiation with a total light dose of 37 J/cm 2 , the so called conventional protocol for Aktilite CL 128 PDT has proven to be effective for the treatment of various skin malignancies [START_REF] Calzavara-Pinton | Methylaminolaevulinate-based photodynamic therapy of Bowen's disease and squamous cell carcinoma[END_REF][START_REF] Szeimies | A clinical study comparing methyl aminolevulinate photodynamic therapy and surgery in small superficial basal cell carcinoma (8-20 mm), with a 12-month follow-up[END_REF][START_REF] Lee | Photodynamic therapy: new treatment for recalcitrant Malassezia folliculitis[END_REF][START_REF] Kim | Photodynamic therapy with methyl-aminolaevulinic acid for mycosis fungoides[END_REF][START_REF] Kim | Photodynamic therapy with ablative carbon dioxide fractional laser for treating Bowen disease[END_REF][START_REF] Fernandez-Guarino | Six Years of Experience in Photodynamic Therapy for Basal Cell Carcinoma: Results and Fluorescence Diagnosis from 191 Lesions[END_REF] including actinic keratosis [START_REF] Morton | Intraindividual, right-left comparison of topical methyl aminolaevulinate-photodynamic therapy and cryotherapy in subjects with actinic keratoses: a multicentre, randomized controlled study[END_REF][START_REF] Szeimies | Photodynamic therapy using topical methyl 5-aminolevulinate compared with cryotherapy for actinic keratosis: A prospective, randomized study[END_REF][START_REF] Pariser | Topical methyl-aminolevulinate photodynamic therapy using red lightemitting diode light for treatment of multiple actinic keratoses: A randomized, doubleblind, placebo-controlled study[END_REF][START_REF] Szeimies | Topical methyl aminolevulinate photodynamic therapy using red light-emitting diode light for multiple actinic keratoses: a randomized study[END_REF]. With similar response rates compared with standard therapies in the treatment for AK, the conventional protocol for Aktilite CL 128 PDT has also demonstrated improved cosmetic outcomes [START_REF] Morton | Intraindividual, right-left comparison of topical methyl aminolaevulinate-photodynamic therapy and cryotherapy in subjects with actinic keratoses: a multicentre, randomized controlled study[END_REF][START_REF] Tyrrell | The effect of air cooling pain relief on protoporphyrin IX photobleaching and clinical efficacy during dermatological photodynamic therapy[END_REF].

B. The European consensus protocol for daylight PDT

Many studies on daylight PDT for the treatment of AK [START_REF] Wiegell | Continuous activation of PpIX by daylight is as effective as and less painful than conventional photodynamic therapy for actinic keratoses; a randomized, controlled, single-blinded study[END_REF][START_REF] Wiegell | Photodynamic therapy of actinic keratoses with 8% and 16% methyl aminolaevulinate and home-based daylight exposure: a double-blinded randomized clinical trial[END_REF][START_REF] Wiegell | A randomized, multicentre study of directed daylight exposure times of 1½ vs. 2½ h in daylight-mediated photodynamic therapy with methyl aminolaevulinate in patients with multiple thin actinic keratoses of the face and scalp[END_REF][START_REF] Wiegell | Daylight photodynamic therapy for actinic keratosis: an international consensus: International Society for Photodynamic Therapy in Dermatology[END_REF][START_REF] Wiegell | Daylight-mediated photodynamic therapy of moderate to thick actinic keratoses of the face and scalp: a randomized multicentre study[END_REF][START_REF] Rubel | Daylight photodynamic therapy with methyl aminolevulinate cream as a convenient, similarly effective, nearly painless alternative to conventional photodynamic therapy in actinic keratosis treatment: a randomized controlled trial[END_REF][START_REF] Wiegell | Weather conditions and daylight-mediated photodynamic therapy: protoporphyrin IX-weighted daylight doses measured in six geographical locations[END_REF][START_REF] Lane | Daylight photodynamic therapy: the Southern California experience[END_REF][START_REF] Perez-Perez | Daylight-mediated photodynamic therapy in Spain: advantages and disadvantages[END_REF][START_REF] See | Consensus recommendations on the use of daylight photodynamic therapy with methyl aminolevulinate cream for actinic keratoses in Australia[END_REF] have been published since the early work of Wiegell et al. [START_REF] Wiegell | Continuous activation of PpIX by daylight is as effective as and less painful than conventional photodynamic therapy for actinic keratoses; a randomized, controlled, single-blinded study[END_REF].

As stated above, the fundamental difference between the protocols for daylight PDT and the protocols for Aktilite CL 128 PDT is the exposure to daylight in place of the red light provided by the Aktilite CL128 lamp.

Furthermore, most of the studies on daylight PDT report an incubation with MAL cream for a maximum of 30 minutes and no excess cream removal before exposure to daylight. The combination of these two factors involves, during all the procedure, an accumulation of PpIX greatly reduced compared to that of the protocols for Aktilite CL 128 PDT. In fact, the incubation time prior to the daylight exposure that is shorter than the usual three hours used for the protocols for Aktilite CL 128 PDT, results in a low initial PpIX accumulation.

Thereafter, without excess cream removal and with a fluence rate lower than the one of the Aktilite CL 128 lamp, the protocols for daylight PDT allows for a balance between the development and the photodegradation of PpIX (the PpIX molecules are photoactivated / photodegraded as quickly they are formed) thus ensuring the maintain of a low PpIX accumulation.

Another difference between the protocols for daylight PDT and the protocols for Aktilite CL 128 PDT is the application of a chemical sunscreen to the treatment area to prevent sunburn.

This sunscreen allowing wavelengths activating PpIX to pass through is usually applied before the gentle curettage of the lesions [START_REF] Wiegell | A randomized, multicentre study of directed daylight exposure times of 1½ vs. 2½ h in daylight-mediated photodynamic therapy with methyl aminolaevulinate in patients with multiple thin actinic keratoses of the face and scalp[END_REF][START_REF] Wiegell | Daylight-mediated photodynamic therapy of moderate to thick actinic keratoses of the face and scalp: a randomized multicentre study[END_REF][START_REF] See | Consensus recommendations on the use of daylight photodynamic therapy with methyl aminolevulinate cream for actinic keratoses in Australia[END_REF].

Among the studies on daylight PDT, two randomized clinical trials each compared a protocol for daylight PDT with the conventional protocol for Aktilite CL 128 PDT [START_REF] Wiegell | Continuous activation of PpIX by daylight is as effective as and less painful than conventional photodynamic therapy for actinic keratoses; a randomized, controlled, single-blinded study[END_REF][START_REF] Rubel | Daylight photodynamic therapy with methyl aminolevulinate cream as a convenient, similarly effective, nearly painless alternative to conventional photodynamic therapy in actinic keratosis treatment: a randomized controlled trial[END_REF]. With an exposure to daylight PDT of 2.5 hours and 2 hours for the first clinical trial [START_REF] Wiegell | Continuous activation of PpIX by daylight is as effective as and less painful than conventional photodynamic therapy for actinic keratoses; a randomized, controlled, single-blinded study[END_REF] and the second one [START_REF] Rubel | Daylight photodynamic therapy with methyl aminolevulinate cream as a convenient, similarly effective, nearly painless alternative to conventional photodynamic therapy in actinic keratosis treatment: a randomized controlled trial[END_REF], respectively, the two involved protocols for daylight PDT (both with 30 minutes incubation time with MAL) have been demonstrated as effective as the conventional protocol for Aktilite CL 128 PDT. Furthermore, due to the above mentioned reduction in PpIX accumulation obtained with the protocols for daylight PDT, patients enrolled in these two trials reported less pain during the protocols for daylight PDT than during the conventional protocol for Aktilite CL 128 PDT [START_REF] Wiegell | Continuous activation of PpIX by daylight is as effective as and less painful than conventional photodynamic therapy for actinic keratoses; a randomized, controlled, single-blinded study[END_REF][START_REF] Rubel | Daylight photodynamic therapy with methyl aminolevulinate cream as a convenient, similarly effective, nearly painless alternative to conventional photodynamic therapy in actinic keratosis treatment: a randomized controlled trial[END_REF]. This experience of a nearly pain-free treatment for patients has also been reported in other studies on daylight PDT [START_REF] Wiegell | Photodynamic therapy of actinic keratoses with 8% and 16% methyl aminolaevulinate and home-based daylight exposure: a double-blinded randomized clinical trial[END_REF][START_REF] Wiegell | A randomized, multicentre study of directed daylight exposure times of 1½ vs. 2½ h in daylight-mediated photodynamic therapy with methyl aminolaevulinate in patients with multiple thin actinic keratoses of the face and scalp[END_REF][START_REF] Wiegell | Daylight photodynamic therapy for actinic keratosis: an international consensus: International Society for Photodynamic Therapy in Dermatology[END_REF][START_REF] Wiegell | Daylight-mediated photodynamic therapy of moderate to thick actinic keratoses of the face and scalp: a randomized multicentre study[END_REF][START_REF] Lane | Daylight photodynamic therapy: the Southern California experience[END_REF][START_REF] Perez-Perez | Daylight-mediated photodynamic therapy in Spain: advantages and disadvantages[END_REF][START_REF] See | Consensus recommendations on the use of daylight photodynamic therapy with methyl aminolevulinate cream for actinic keratoses in Australia[END_REF].

Finally, with a shorter time of clinic attendance than the conventional protocol for Aktilite CL 128 PDT, the protocols for daylight PDT are more convenient for both patients and clinicians.

With a maximum of 30 minutes of MAL incubation and no excess cream removal before daylight exposure for 2 hours, the so called European consensus protocol for daylight PDT has been adopted from a European consensus in 2012 [START_REF] Wiegell | Daylight photodynamic therapy for actinic keratosis: an international consensus: International Society for Photodynamic Therapy in Dermatology[END_REF]. Based on the above mentioned randomized clinical trials [START_REF] Wiegell | Continuous activation of PpIX by daylight is as effective as and less painful than conventional photodynamic therapy for actinic keratoses; a randomized, controlled, single-blinded study[END_REF][START_REF] Rubel | Daylight photodynamic therapy with methyl aminolevulinate cream as a convenient, similarly effective, nearly painless alternative to conventional photodynamic therapy in actinic keratosis treatment: a randomized controlled trial[END_REF], we set the incubation time for the European consensus protocol for daylight PDT to 30 minutes.

In this study, the outdoor temperature was assumed to be sufficient for the PDT process to occur [START_REF] Wiegell | Weather conditions and daylight-mediated photodynamic therapy: protoporphyrin IX-weighted daylight doses measured in six geographical locations[END_REF][START_REF] Mordon | A commentary on the role of skin temperature on the effectiveness of ALA-PDT in dermatology[END_REF] and thus no further consideration has been given to the temperature in this paper.

III. Mathematical modeling for topical PDT protocol

A. AK sample model

To account for both the confinement of AKs to the epidermis and the usual 100 μm thickness of epidermis [START_REF] Liu | A dynamic model for ALA-PDT of skin: simulation of temporal and spatial distributions of ground-state oxygen, photosensitizer and singlet oxygen[END_REF], the simplified AK sample model consists of a 10 µm wide and 100 μm thick parallelepiped (Figure 1). The AK tissue is assumed homogeneous and its optical properties are set to the values reported in Garcia-Uribe et al. [START_REF] Garcia-Uribe | In vivo diagnosis of melanoma and nonmelanoma skin cancer using oblique incidence diffuse reflectance spectrometry[END_REF].

A primary planar beam with fluence rate 0 S is assumed to perpendicularly irradiate the surface of the AK sample model as illustrated in Figure 1.

B. Models for the different fluence rates

As depending on the weather conditions [START_REF] Wiegell | Weather conditions and daylight-mediated photodynamic therapy: protoporphyrin IX-weighted daylight doses measured in six geographical locations[END_REF][START_REF] Spelman | Treatment of face and scalp solar (actinic) keratosis with daylight-mediated photodynamic therapy is possible throughout the year in Australia: Evidence from a clinical and meteorological study[END_REF][START_REF] O'gorman | Artificial White Light vs Daylight Photodynamic Therapy for Actinic Keratoses: A Randomized Clinical Trial[END_REF], the spectral fluence rate for the daylight is not unique. In this paper, two spectral fluence rates have been used: the first one corresponds to a clear sunny day (blue curve in Figure 2) while the second stands for an overcast day (cyan curve in Figure 2). These two spectral fluence rates were recorded in the study of O'Gorman et al. [START_REF] O'gorman | Artificial White Light vs Daylight Photodynamic Therapy for Actinic Keratoses: A Randomized Clinical Trial[END_REF], who kindly provided them to us. The use of these two spectral fluence rates for the daylight aims to quantify the effect of the weather conditions on the PDT local damage.

Among the various spectral fluence rates for the Aktilite CL 128 lamp we have at our disposal, the one measured at a distance of 8 cm from the lamp by O'Gorman et al. [START_REF] O'gorman | Artificial White Light vs Daylight Photodynamic Therapy for Actinic Keratoses: A Randomized Clinical Trial[END_REF] (red curve in Figure 2) was preferred. Indeed, recorded using the same measurement system as the two spectral fluence rates for daylight, this spectral fluence rate seems to be the most appropriate for a pertinent comparison.

C. Modeling of the PDT process

The modeling method is the same as in our previously validated work [START_REF] Vignion-Dewalle | Red light photodynamic therapy for actinic keratosis using 37 J/cm 2 : fractionated irradiation with 12.3 mW/cm 2 after 30 minutes incubation time compared to standard continuous irradiation with 75 mW/cm 2 after three hours incubation time using a mathematical modeling[END_REF] Based on our previous works [START_REF] Vignion-Dewalle | Comparison of three light doses in the photodynamic treatment of actinic keratosis using mathematical modeling[END_REF][START_REF] Vignion-Dewalle | Red light photodynamic therapy for actinic keratosis using 37 J/cm 2 : fractionated irradiation with 12.3 mW/cm 2 after 30 minutes incubation time compared to standard continuous irradiation with 75 mW/cm 2 after three hours incubation time using a mathematical modeling[END_REF], the modeling of the PDT process consists of two steps that are iteratively repeated: determination of the local fluence rate and updating of the PpIX absorption coefficient.

Determination of the fluence rate

The local total fluence rate at time t , depth z and wavelength  , denoted by  

  , , z t , is
given by equation 1 [START_REF] Vignion-Dewalle | Comparison of three light doses in the photodynamic treatment of actinic keratosis using mathematical modeling[END_REF][START_REF] Vignion-Dewalle | Red light photodynamic therapy for actinic keratosis using 37 J/cm 2 : fractionated irradiation with 12.3 mW/cm 2 after 30 minutes incubation time compared to standard continuous irradiation with 75 mW/cm 2 after three hours incubation time using a mathematical modeling[END_REF][START_REF] Farrell | Modeling of photosensitizer fluorescence emission and photobleaching for photodynamic therapy dosimetry[END_REF]:

                                                              end start z t z
         (1)
Where:

 The above defined 0 S is the spectral fluence rate of the primary planar broad beam,  Due to the accumulation of PpIX in AK lesions, the total absorption coefficient, P , depending on both the optical properties of the actinic keratosis and the boundary conditions at the actinic keratosis surface, are computed as described in [START_REF] Farrell | Modeling of photosensitizer fluorescence emission and photobleaching for photodynamic therapy dosimetry[END_REF].

Updating of the PpIX absorption coefficient

The updating formula for the PpIX absorption coefficient in an unit volume, U V , is expressed as follows (equation 2) [START_REF] Vignion-Dewalle | Comparison of three light doses in the photodynamic treatment of actinic keratosis using mathematical modeling[END_REF][START_REF] Vignion-Dewalle | Red light photodynamic therapy for actinic keratosis using 37 J/cm 2 : fractionated irradiation with 12.3 mW/cm 2 after 30 minutes incubation time compared to standard continuous irradiation with 75 mW/cm 2 after three hours incubation time using a mathematical modeling[END_REF]:

                                                                                    ~, , , , , , , , , , exp 1 1
exp 1 1 exp , , , , d z t V c z t dt z t V dt t k dt t k z L V z t z dt t PpIX a U PpIX a U U PpIX PpIX a PpIX a  (2)
Where:

 dt is the time increment,     
PpIX is the PpIX molar extinction coefficient for wavelength  ,   ,  and c are the Avogadro number (6.022×10 23 /mol), the Planck constant (6.626×10 34 J×s) and the speed of light (3×10 8 m/s), respectively,  L , k ,  and  are the parameters of the depth-dependent logistic growth related to both the biological clearance of PpIX and the conversion of 5-ALA into PpIX (2 nd term in the right hand side) [START_REF] Vignion-Dewalle | Red light photodynamic therapy for actinic keratosis using 37 J/cm 2 : fractionated irradiation with 12.3 mW/cm 2 after 30 minutes incubation time compared to standard continuous irradiation with 75 mW/cm 2 after three hours incubation time using a mathematical modeling[END_REF],

  and   ~ are the bimolecular rate constant for the reaction of singlet oxygen with PpIX and the singlet oxygen quantum yield, respectively. These two parameters allow the photobleaching process to be analytically modeled (3 rd term in the right hand side) [START_REF] Vignion-Dewalle | Red light photodynamic therapy for actinic keratosis using 37 J/cm 2 : fractionated irradiation with 12.3 mW/cm 2 after 30 minutes incubation time compared to standard continuous irradiation with 75 mW/cm 2 after three hours incubation time using a mathematical modeling[END_REF].

Iterative procedure

All the involved parameters were assigned to the values reported in our previous work and empirically determined in the literature [START_REF] Vignion-Dewalle | Comparison of three light doses in the photodynamic treatment of actinic keratosis using mathematical modeling[END_REF][START_REF] Vignion-Dewalle | Red light photodynamic therapy for actinic keratosis using 37 J/cm 2 : fractionated irradiation with 12.3 mW/cm 2 after 30 minutes incubation time compared to standard continuous irradiation with 75 mW/cm 2 after three hours incubation time using a mathematical modeling[END_REF] (Table 1). Assuming the initial distribution for the PpIX absorption coefficient reported in the last row in Table 1, .

              h 3 exp 1 k V U  5.3×10 9 l/mol/s   0.56     , , 0 , z PpIX a                 k z L V U PpIX exp 1 exp
From several studies on PDT [START_REF] Farrell | Modeling of photosensitizer fluorescence emission and photobleaching for photodynamic therapy dosimetry[END_REF][START_REF] Valentine | Monte Carlo modeling of in vivo protoporphyrin IX fluorescence and singlet oxygen production during photodynamic therapy for patients presenting with superficial basal cell carcinomas[END_REF], this cumulative quantity enables the quantification of the PDT local damage, D , over time (equation 3).

                                          1 0 ~, ~, , , , , i j start PpIX a U start start d z dt j t V c z dt j t dt z dt i t D           (3)

IV. Results

A. Comparison in terms of effective fluence

Given the spectral fluence rate for the Aktilite CL 128 lamp (Figure 2), several parameters have been computed (Table 2). First, from the integration of this spectrum over wavelength, the fluence rate has been estimated to be 85.39 mW/cm 2 that is consistent with the abovementioned 70-100 mW/cm 2 range provided by the manufacturer. By considering this fluence rate, the fluence of 37 J/cm 2 of the conventional protocol for Aktilite CL 128 PDT is achieved using an irradiation time of 433.3 seconds (Table 2). Moreover, weighting the spectral fluence rate with the normalized absorption spectrum for PpIX, which is derived from data measured by the Research Center for Automatic Control of Nancy (CRAN) (Figure 3.a), enables the effective or PpIX-weighted spectral fluence rate to be obtained (Figure 3.a) and the effective fluence rate to be deducted by integration over wavelength (approximately, 1.44 mW/cm 2 ) [START_REF] Wiegell | Continuous activation of PpIX by daylight is as effective as and less painful than conventional photodynamic therapy for actinic keratoses; a randomized, controlled, single-blinded study[END_REF][START_REF] O'gorman | Artificial White Light vs Daylight Photodynamic Therapy for Actinic Keratoses: A Randomized Clinical Trial[END_REF] (Table 2). Multiplying the effective fluence rate by the above determined irradiation time of 433.3 seconds leads to an effective fluence of 0.63 J/cm 2 for the conventional protocol for Aktilite CL 128 PDT.

Similar computations have been carried out for the sunny daylight spectral fluence rate and the overcast daylight spectral fluence rate (Figure 2). The integration over wavelength of these two spectral fluence rates yields fluence rates of 55.2 mW/cm 2 and 7.75 mW/cm 2 for the sunny daylight and the overcast daylight, respectively (Table 2). The effective or PpIXweighted spectral fluence rate for the sunny daylight (respectively, the overcast daylight), which was obtained by weighting the spectral for the sunny daylight fluence rate (respectively, the overcast daylight fluence rate) with the normalized absorption spectrum for PpIX, yields, when integrated over wavelength, an effective fluence rate of 5.11 mW/cm 2 (respectively, 0.80 mW/cm 2 ) for the sunny daylight (respectively, the overcast daylight) From these computations performed using a homemade software available online [START_REF] Vignion-Dewalle | A software for analyzing and comparing the light sources available for PDT in dermatology[END_REF], the effective fluence achieved using the European consensus protocol for daylight PDT during a sunny day is 58.40 and 6.40 times higher than those achieved using the conventional protocol for Aktilite CL 128 PDT and the European consensus protocol for daylight PDT during an overcast day, respectively. The effective fluence obtained using the European consensus protocol for daylight PDT during an overcast day is alsoalthough to a lesser extent than the European consensus protocol for daylight PDT during a sunny dayhigher (9.13 times higher) than that obtained using the conventional protocol for Aktilite CL 128. different scale is applied for (a) due to its much higher amplitude). All the three spectral fluence rates were provided by the authors of [START_REF] O'gorman | Artificial White Light vs Daylight Photodynamic Therapy for Actinic Keratoses: A Randomized Clinical Trial[END_REF]. computed from the spectral fluence rate provided by [START_REF] O'gorman | Artificial White Light vs Daylight Photodynamic Therapy for Actinic Keratoses: A Randomized Clinical Trial[END_REF]. The irradiation times used in this paper (row 3) and the corresponding standard and effective fluences are also indicated (rows 5 and 7).

(Table 2,

Light source Aktilite

B. Comparison in terms of photodynamic dose

All the computations were performed using a Matlab™ program on a standard personal computer (Intel Xeon CPU E3-1240 V2 3.40 GHz-8Go of RAM-Windows 7 64 bits).

The PDT local damages obtained using the conventional protocol for Aktilite CL 128 PDT, the European consensus protocol for daylight PDT during a sunny day and the European consensus protocol for daylight PDT during an overcast day are displayed as a function of depth in Figure 4 and as a function of time in Figure 5. From Figure 4, whatever the depth position in the AK sample model (Figure 1), the PDT local damage achieved at the end of the treatment using the European consensus protocol for daylight PDT during a sunny day is on average about 6.50 and 3.63 times higher than those achieved at the end of the treatment using the conventional protocol for Aktilite CL 128 PDT and using the European consensus protocol for daylight PDT during an overcast day, respectively. The ratio between the PDT local damage obtained using the European consensus protocol for daylight PDT during an overcast day and the one obtained using the conventional protocol for Aktilite CL 128 PDT is approximately 1.80 at 0 µm depth and 1.77 at 100 µm depth. From Figure 5, for the European consensus protocol for daylight PDT, the shape of the time courses of the PDT local damage tends to demonstrate an exponential trend (blue and cyan curves in Figure 5) while that for the conventional protocol for Aktilite CL 128 PDT suggests a very slight logarithmic trend (red curves). This results that, at least over the considered time intervals, the PDT local damage produced using the conventional protocol for Aktilite CL 128

PDT increases in time much faster than those produced using the European consensus protocol for daylight PDT (Figure 5).

C. Comparison in terms of number of PpIX molecules

The evolution of the number of PpIX molecules present in an unit volume, U V , at time t and depth z , deduced from equation 4, is illustrated for the conventional protocol for Aktilite CL 128 PDT, the European consensus protocol for daylight PDT during a sunny day and the European consensus protocol for daylight PDT during an overcast day in Figure 6. From Figure 6, regarding the conventional protocol for Aktilite CL 128 PDT, the irradiation (starting at time=10800 s) leads to a mean percent drop of 35.24% in the number of PpIX molecules while the number of PpIX molecules corresponding to the European consensus protocol for daylight PDT during an overcast day continues to increase even after the beginning of irradiation (occurring at time 1800 seconds). For the European consensus protocol for daylight PDT during a sunny day, the irradiation results, in terms of the number of PpIX molecules, in a transient slight decrease followed by a further controlled increase.

          PpIX U PpIX a PpIX V z t z t M    , , , , (4)
Furthermore, an important number of PpIX molecules is still present at the end of the irradiation for the conventional protocol for Aktilite CL 128 PDT and the European consensus protocol for daylight PDT during an overcast day (Figure 6). With regard to the European consensus protocol for daylight PDT during a sunny day, the number of PpIX molecules present at the end of the irradiation is slightly higher than that at the beginning of irradiation.

V. Discussion

In this paper, a comparison between the two most widely used MAL-PDT protocols in Europe for the treatment of actinic keratosis (AK) is performed using a mathematical modeling of the PDT process:

 The conventional protocol for Aktilite CL 128 PDT (light source: Aktilite CL 128, incubation time: three hours, light dose: 37 J/cm 2 ) [START_REF] Morton | Intraindividual, right-left comparison of topical methyl aminolaevulinate-photodynamic therapy and cryotherapy in subjects with actinic keratoses: a multicentre, randomized controlled study[END_REF][START_REF] Tyrrell | The effect of air cooling pain relief on protoporphyrin IX photobleaching and clinical efficacy during dermatological photodynamic therapy[END_REF][START_REF] Pariser | Topical methyl-aminolevulinate photodynamic therapy using red lightemitting diode light for treatment of multiple actinic keratoses: A randomized, doubleblind, placebo-controlled study[END_REF],

 The European consensus protocol for daylight PDT (light source: daylight, incubation time: 30 minutes, treatment duration: two hours) [START_REF] Wiegell | Daylight photodynamic therapy for actinic keratosis: an international consensus: International Society for Photodynamic Therapy in Dermatology[END_REF].

Two weather conditions for the European consensus protocol for daylight PDT have been considered for this comparison: a clear sunny day and an overcast day. The spectral fluence rates for the Aktilite CL 128, the clear sunny daylight and the overcast daylight have been provided by the authors of O'Gorman et al. [START_REF] O'gorman | Artificial White Light vs Daylight Photodynamic Therapy for Actinic Keratoses: A Randomized Clinical Trial[END_REF] (Figure 2).

The comparison is performed using an AK sample model consisting of a 100 μm thick parallelepiped (Figure 1) and a recently published modeling of the PDT process iteratively alternating determination of the fluence rate and updating of the optical properties [START_REF] Vignion-Dewalle | Comparison of three light doses in the photodynamic treatment of actinic keratosis using mathematical modeling[END_REF][START_REF] Vignion-Dewalle | Red light photodynamic therapy for actinic keratosis using 37 J/cm 2 : fractionated irradiation with 12.3 mW/cm 2 after 30 minutes incubation time compared to standard continuous irradiation with 75 mW/cm 2 after three hours incubation time using a mathematical modeling[END_REF].

The determination of the fluence rate involves solving the one-dimensional diffusion equation (equation 1, [START_REF] Vignion-Dewalle | Comparison of three light doses in the photodynamic treatment of actinic keratosis using mathematical modeling[END_REF][START_REF] Vignion-Dewalle | Red light photodynamic therapy for actinic keratosis using 37 J/cm 2 : fractionated irradiation with 12.3 mW/cm 2 after 30 minutes incubation time compared to standard continuous irradiation with 75 mW/cm 2 after three hours incubation time using a mathematical modeling[END_REF]) while the updating of the optical properties takes the biological clearance of PpIX, the conversion of MAL into PpIX and the PpIX photobleaching into account (equation 2, [START_REF] Vignion-Dewalle | Red light photodynamic therapy for actinic keratosis using 37 J/cm 2 : fractionated irradiation with 12.3 mW/cm 2 after 30 minutes incubation time compared to standard continuous irradiation with 75 mW/cm 2 after three hours incubation time using a mathematical modeling[END_REF]). All the parameters involved in equations 1 and 2 are set to published empirical values, which were obtained with PpIX and with either normal human epidermis or AK (Table 1, [START_REF] Vignion-Dewalle | Comparison of three light doses in the photodynamic treatment of actinic keratosis using mathematical modeling[END_REF]). The updating formula (equation 2) explicitly provides the number of singlet oxygen molecules produced at each iteration (integral in the right hand side).

Cumulating this number over the irradiation time enables the quantification of the PDT local damage (equation 3, [START_REF] Vignion-Dewalle | Comparison of three light doses in the photodynamic treatment of actinic keratosis using mathematical modeling[END_REF][START_REF] Vignion-Dewalle | Red light photodynamic therapy for actinic keratosis using 37 J/cm 2 : fractionated irradiation with 12.3 mW/cm 2 after 30 minutes incubation time compared to standard continuous irradiation with 75 mW/cm 2 after three hours incubation time using a mathematical modeling[END_REF]).

Applying to the conventional protocol for Aktilite CL 128 PDT, to the European consensus protocol for daylight PDT during a sunny day and to the European consensus protocol for daylight PDT during an overcast day, the model allows evaluation and comparison of their performance in terms of the PDT local damage. A comparison has also been performed in terms of effective or PpIX-weighted fluence.

From the results, regardless the two weather conditions considered in this study, the European consensus protocol for daylight PDT is more efficient in terms of both effective fluence and PDT local damage than the conventional protocol for Aktilite CL 128 PDT. First, this finding underlines the relevance of the European consensus protocol for daylight PDT, even when performed during an overcast day. Then this finding supports the non-inferiority in efficacy of the European consensus protocol for daylight PDT compared to the conventional protocol for Aktilite CL 128 PDT that has been demonstrated by a clinical trial, including all weather conditions except rain or cold [START_REF] Rubel | Daylight photodynamic therapy with methyl aminolevulinate cream as a convenient, similarly effective, nearly painless alternative to conventional photodynamic therapy in actinic keratosis treatment: a randomized controlled trial[END_REF]. This finding also suggests that, although a minimum effective fluence of 8 J/cm 2 is recommended with daylight PDT to result in an effective treatment of AK [START_REF] Wiegell | Photodynamic therapy of actinic keratoses with 8% and 16% methyl aminolaevulinate and home-based daylight exposure: a double-blinded randomized clinical trial[END_REF][START_REF] Wiegell | Weather conditions and daylight-mediated photodynamic therapy: protoporphyrin IX-weighted daylight doses measured in six geographical locations[END_REF], the effective fluence of 5.75 J/cm 2 related to the overcast day considered in this study (Table 2), is sufficient for the European consensus protocol for daylight PDT to perform better than the efficient conventional protocol for Aktilite CL 128 PDT. Finally, based on the well-known efficiency of the conventional protocol for Aktilite CL 128 PDT, this finding could reflect a potential over-treatment of the lesions when using the European consensus protocol for daylight PDT during a sunny day but also during an overcast day.

As already mentioned in many studies on daylight PDT [START_REF] Wiegell | Photodynamic therapy of actinic keratoses with 8% and 16% methyl aminolaevulinate and home-based daylight exposure: a double-blinded randomized clinical trial[END_REF][START_REF] Wiegell | A randomized, multicentre study of directed daylight exposure times of 1½ vs. 2½ h in daylight-mediated photodynamic therapy with methyl aminolaevulinate in patients with multiple thin actinic keratoses of the face and scalp[END_REF][START_REF] Wiegell | Daylight-mediated photodynamic therapy of moderate to thick actinic keratoses of the face and scalp: a randomized multicentre study[END_REF][START_REF] Wiegell | Weather conditions and daylight-mediated photodynamic therapy: protoporphyrin IX-weighted daylight doses measured in six geographical locations[END_REF][START_REF] O'gorman | Artificial White Light vs Daylight Photodynamic Therapy for Actinic Keratoses: A Randomized Clinical Trial[END_REF] and as evidenced by the above mentioned values, the effective fluence varies depending on weather conditions. This dependence on weather conditions is also evident here in terms of the PDT local damage:

the better the weather conditions, the more efficient is the European consensus protocol for daylight PDT. Nonetheless, based on the studies of Wiegell et al. [START_REF] Wiegell | Photodynamic therapy of actinic keratoses with 8% and 16% methyl aminolaevulinate and home-based daylight exposure: a double-blinded randomized clinical trial[END_REF][START_REF] Wiegell | Weather conditions and daylight-mediated photodynamic therapy: protoporphyrin IX-weighted daylight doses measured in six geographical locations[END_REF] and O'Gorman et al [START_REF] O'gorman | Artificial White Light vs Daylight Photodynamic Therapy for Actinic Keratoses: A Randomized Clinical Trial[END_REF] that have found no association between response rate and effective fluence in patients who received an effective fluence with daylight PDT higher than 8 J/cm 2 and 3.2 J/cm 2 , respectively, a sunny day is not necessarily required for the European consensus protocol for daylight PDT to be effective (an overcast daylight may be sufficient for the European consensus protocol for daylight PDT to destroy any cancer cells). This has also been highlighted by Rubel et al. [START_REF] Rubel | Daylight photodynamic therapy with methyl aminolevulinate cream as a convenient, similarly effective, nearly painless alternative to conventional photodynamic therapy in actinic keratosis treatment: a randomized controlled trial[END_REF], which have reported that, although patients received variable effective fluences during the European consensus protocol for daylight PDT, no correlation to efficacy was found. These clinical results with similar response rates whatever the weather conditions (except rainy or cold conditions) tend to support the above suggested hypothesis of overtreatment of the lesions by the European consensus protocol for daylight PDT.

outweigh the conversion of MAL into PpIX. Regarding the European consensus protocol for daylight PDT during a sunny day, despite the low initial PpIX accumulation, the 5.11 mW/cm 2 effective fluence rate appears to be appropriate for maintaining approximately constant the number of PpIX molecules.

Consideration also needs to be given to the important number of PpIX molecules still present at the end of the conventional protocol for Aktilite CL 128 PDT (Figure 6). This important number tends to demonstrate that the incubation time or the cream concentration in MAL for this protocol could be reduced [START_REF] Christiansen | 5-ALA for photodynamic photorejuvenation--optimization of treatment regime based on normal-skin fluorescence measurements[END_REF]. This observation can be extended to the European consensus protocol for daylight PDT during an overcast day…

The overall results emphasize the need to refine the parameters of both the conventional protocol for Aktilite CL 128 PDT and the European consensus protocol for daylight PDT.

This refinement (reduction of the incubation time, reduction of the irradiation time…) would allow for a similar efficiency but an improved tolerability and a more manageable clinical practice.

Only two weather conditions for the European consensus protocol for daylight PDT have been studied in this paper. These two weather conditions have been chosen because their spectral fluence rates discussed in the study of O'Gorman et al. [START_REF] O'gorman | Artificial White Light vs Daylight Photodynamic Therapy for Actinic Keratoses: A Randomized Clinical Trial[END_REF] have been kindly provided by the authors. Any other weather condition with an available spectral fluence rate could obviously be investigated using the same mathematical modeling… The subsequent results will either support the above suggested hypothesis of overtreatment of the lesions by the European consensus protocol for daylight PDT, or will allow to estimate the PDT local damage at which the treatment becomes effective.

In this study, the temperature was, as above mentioned, assumed to be at least 10°C so that the PDT process to occur [START_REF] Wiegell | Weather conditions and daylight-mediated photodynamic therapy: protoporphyrin IX-weighted daylight doses measured in six geographical locations[END_REF][START_REF] Mordon | A commentary on the role of skin temperature on the effectiveness of ALA-PDT in dermatology[END_REF]. If this is not the case, or if it's raining or windy, a greenhouse may be used [START_REF] Lerche | Alternatives to Outdoor Daylight Illumination for Photodynamic Therapy--Use of Greenhouses and Artificial Light Sources[END_REF]. However, the spectral fluence rate of the daylight is modified by the filtering effect of the greenhouse glass, and the mathematical modeling proposed in this paper could be used to quantify this filtering effect in terms of the PDT local damage. Another solution to get rid of this dependence on the weather conditions, geographical location, seasons…, is to use a light source with a spectral fluence rate close to that of daylight [START_REF] O'gorman | Artificial White Light vs Daylight Photodynamic Therapy for Actinic Keratoses: A Randomized Clinical Trial[END_REF][START_REF] Lerche | Alternatives to Outdoor Daylight Illumination for Photodynamic Therapy--Use of Greenhouses and Artificial Light Sources[END_REF]; a comparative study of these light sources for "indoor" daylight PDT could also be performed in terms of the PDT local damage as defined in this paper.

VI. Conclusion

In this paper, we have compared the conventional protocol for Aktilite CL 128 PDT (light source: Aktilite CL 128, incubation time: three hours, light dose: 37 J/cm 2 ) with the European consensus protocol for daylight PDT (light source: daylight, incubation time: 30 minutes, treatment duration: two hours). This comparison performed in terms of the effective fluence and the PDT local damage tends to demonstrate, whatever the two weather conditions considered for this comparison (i.e., a clear sunny day and an overcast day) the European consensus protocol for daylight PDT perform better than the conventional protocol for Aktilite CL 128 PDT.

Figure 1 :

 1 Figure 1: The AK sample model

Figure 2 :

 2 Figure 2: The spectral fluence rate for the Aktilite CL 128 lamp, the sunny daylight and the overcast daylight, respectively. All these three spectral fluence rates were provided by the authors of [39].

a

  , is the sum of the PpIX absorption coefficient, PpIX a,  , and the actinic keratosis absorption coefficient, The total transport coefficient, t  , is the sum of the total absorption coefficient, a  , and the actinic keratosis reduced scattering coefficient, The two parameters, b and

  calculated. Applying equation 2 then yields the distribution for the PpIX absorption coefficient at time dt , deduced. By iterating equations 1 and 2, all the PpIX absorption coefficients and local total fluence rates can be determined throughout the treatment.

Figure 3 .

 3 b and 3.c). Using an exposure of 2 hours as required with the European consensus protocol for daylight PDT, the sunny daylight fluence rate (respectively, the effective sunny daylight fluence rate) leads to a fluence of 397.44 J/cm 2 (respectively, an effective fluence of 36.79 J/cm 2 ) while the overcast daylight fluence rate (respectively, the effective overcast daylight fluence rate) provides a fluence of 55.8 J/cm 2 (respectively, an effective fluence of 5.76 J/cm 2 ).

Figure 3 :

 3 Figure 3: The spectral fluence rate (green curves) and the effective spectral fluence rate (red curve in a), blue curve in b) and cyan curve in c)) for a) the Aktilite CL 128 lamp at a distance of 8 cm from the lamp, b) the sunny daylight and c) the overcast daylight are scaled according to the left axis while the normalized absorption spectrum for PpIX (black curve) is plotted according to the right axis (Note: left axes in (b) and (c) are identically scaled whereas a

Figure 4 :

 4 Figure 4: Evolution in depth of the PDT local damage achieved at the end of the treatment using the conventional protocol for Aktilite CL 128 PDT (red curve), the European consensus protocol for daylight PDT during a sunny day (blue curve) and the European consensus protocol for daylight PDT during an overcast day (cyan curve).

Figure 5 :

 5 Figure 5: Evolution in time of the PDT local damage achieved using the conventional protocol for Aktilite CL 128 PDT (red curves), the European consensus protocol for daylight PDT during a sunny day (blue curves) and the European consensus protocol for daylight PDT during an overcast day (cyan curves). The bright solid curves (respectively, light dashed curves) represent the PDT local damages at 0 µm (respectively, at 100 µm) in depth in AK.

Figure 6 :

 6 Figure 6: Evolution in time of the number of PpIX molecules present at 50 μm in depth in AK using the conventional protocol for Aktilite CL 128 PDT (red curve), the European consensus protocol for daylight PDT during a sunny day (blue curve) and the European consensus protocol for daylight PDT during an overcast day (cyan curve).

  and therefore only outlines are referred to in this paper without further discussion.

	Let the incubation start at time	t 0 	s	and let the light irradiation be performed during the time
	interval  t , start t	end		:
		t	start		3	hours and	t	end		t	start		t 	(	t  is the irradiation time necessary to achieve the
		recommended 37 J/cm 2 using the spectral fluence rate for the Aktilite CL 128 lamp
		depicted in Figure 2) for the conventional protocol for Aktilite CL 128 PDT,
		t			5 . 0	hours and		t		t		7200	s	for the European consensus protocol for daylight
			start								end		start
		PDT.

Table 1 :

 1 Specification of the model parameters from[START_REF] Vignion-Dewalle | Red light photodynamic therapy for actinic keratosis using 37 J/cm 2 : fractionated irradiation with 12.3 mW/cm 2 after 30 minutes incubation time compared to standard continuous irradiation with 75 mW/cm 2 after three hours incubation time using a mathematical modeling[END_REF]

D. Quantification of the PDT local damage The

  integral in the last part of the right hand side of equation 2 represents the number of singlet oxygen molecules generated during the time interval   , located at depth z in the AK sample model when the PpIX molecules, excited by absorption of photons, return to the ground state. Summing this integral over the time intervals

																														t t  ;	dt	in an unit volume,
	 t	start	;	t	start		dt		,  t	start		dt	;	t	start		2		dt		, …,	 t	start		  i  1		dt	;	t	start		i		dt		provides the total
	cumulative singlet oxygen molecules produced during the time interval  t	start	;	t	start		i		dt	

U V

Table 2 :

 2 Standard (row 4) and effective (row 6) fluence rates for the Aktilite CL128 lamp (second column), the sunny daylight (third column) and the overcast daylight (fourth column)

	Sunny	Overcast
	daylight	daylight

Acknowledgements

The authors would like to thank M. Manley from the Department of Medical Physics and Clinical Engineering, Saint Vincent's University Hospital, Dublin, Ireland, UK, for providing the spectral fluence rates for 1) the Aktilite CL128 lamp (Galderma SA, Lausanne, Switzerland), 2) the daylight on a clear sunny day and 3) the daylight on an overcast day.

The information provided by the ratios between the effective fluences is consistent with that provided by the ratios between the PDT local damages. In fact, if a ratio between two considered effective fluences is higher than 1, then the ratio between the two corresponding PDT local damages is also higher than 1. Nonetheless, the ratios between the effective fluences are always higher than those between the PDT local damages. This is explained, to a large extent, by the fact that the effective fluence is computed using the normalized absorption spectrum for PpIX (Figure 3) whereas the PDT local damages involves the "actual" absorption spectrum for PpIX. With an incubation time longer than that for the European consensus protocol for daylight PDT, the conventional protocol for Aktilite CL 128 PDT leads to a higher initial PpIX accumulation (Figure 6) and subsequently to a higher initial absorption spectrum for PpIX. This in turn allows to partially offset the higher effective fluence rate for the European consensus protocol for daylight PDT in the iterative calculation of the PDT local damage (equations 1-3). In contrast, as the effective fluence is computed from the normalized PpIX absorption spectrum, no such offset is present and thus the information provided by the effective fluence might not be sufficient to predict the result of the MAL-PDT procedure. The PDT local damage, as defined in this paper, may therefore be a more appropriate predictor… Regarding the evolution in time of the PDT local damage (Figure 5), with an effective fluence rate about 3.55 times lower than that of the European consensus protocol for daylight PDT during a sunny day (Table 2), the conventional protocol for Aktilite CL 128 PDT, however, leads to a higher increase rate for the PDT local damage. This results from the above mentioned longer incubation time and subsequent higher initial PpIX absorption coefficient for the conventional protocol for Aktilite CL 128 PDT compared to the European consensus protocol for daylight PDT, leading to a higher initial photons absorption efficiency (equations 1-3).

From Figure 6, the beginning of the irradiation is clearly identifiable for the conventional protocol for Aktilite CL 128 PDT with a more than 30 percent drop in the number of PpIX molecules resulting from the above mentioned high initial photons absorption efficiency. On the contrary, the steady growth curve observed for the European consensus protocol for daylight PDT during an overcast day makes the identification of the beginning of the irradiation impossible: the 0.80 mW/cm 2 effective fluence rate combined with a low initial PpIX accumulation do not allow a photobleaching of the PpIX molecules important enough to


Abdeldjalil Aissa

El Bey

Amina Saleem

Karim Abed-Meraim

Azzedine Beghdadi

Abdeldjalil Aïssa-El-Bey

Watermark extraction using blind image separation and sparse representation

INTRODUCTION

Digital watermarking technology has evolved as an important technology in the recent years. The basic principle of most watermarking method involves the application of small, pseudo-random changes to the selected coefficients in the spatial or transform domain. Most of the watermark detection schemes use some kinds of correlating detector to verify the presence of the embedded watermark [START_REF] Hartung | Multimedia watermarking techniques[END_REF][START_REF] Stefan | Information Hiding Techniques for Steganography and Digital Watermarking[END_REF]. Digital image watermarking has its applications in copy rights protection, data tracking and monitoring [START_REF] Hartung | Multimedia watermarking techniques[END_REF]. Blind source separation (BSS) is an important area of research in signal and image processing. The BSS problem can be solved using sparse representations of the source signals. Solution for the blind separation of image sources using sparsity include the wavelet-transform domain method in [START_REF] Stefan | Information Hiding Techniques for Steganography and Digital Watermarking[END_REF] and the method in [START_REF] Bronstein | Separation of reflections via sparse ICA[END_REF] using projection onto sparse dictionaries and the iterative Blind source separation algorithm presented in [START_REF] Souidene | Blind image separation using sparse representation[END_REF]. This paper introduces a new image watermark extraction technique based on the iterative sparse blind separation algorithm (ISBS) and a new ISBS algorithm. The ISBS algorithm employs an p norm based contrast function for blind signal separation. When the images are sparse or sparsely representable, a smooth approximation of the absolute value function is a good choice for the cost function. The NISBS algorithm proposed in this paper employs the modeling of the distributions of sparse images using a family of convex smooth functions. The ISBS algorithm proposed in [START_REF] Souidene | Blind image separation using sparse representation[END_REF] and the NISBS algorithm presented in this paper are shown to be more efficient than other existing techniques in the literature and both lead to improved separation quality with lower computational cost. The performance of the proposed algorithms is compared to the performance of other ICA algorithms using the objective image quality measure inspired by the Human Visual System (HVS) proposed in [START_REF] Beghdadi | A new image distortion measure based on wavelet decomposition[END_REF]. It is shown that the ISBS and NISBS algorithms perform better in terms of PSNR-WAV. The new watermarking technique is presented with the principal assumptions of (i) image source sparsity, (ii) instantaneous mix-tures and (iii) the same number of mixtures and sources (three mixtures and three sources). In our proposed BSS based method, we do not have restrictions on the mixing process as well as the mixing coefficients. The paper is organized as follows. Section 2 presents the data model and assumptions of our system. The blind watermark extraction system using the sparsity based algorithms is described in Section 3. The simulations and the performance of the algorithms is discussed in Section 4. The conclusions are drawn in Section 5.

DATA MODEL AND ASSUMPTIONS

A generic watermark embedding system consists of the inputs which are the original data f1 , the watermark signal f2 and an optional public or secret key f3 each of size (m f , n f ). The key is used to enforce the security, that is, to prevent unauthorized party from recovering and manipulating the watermark. The proposed image watermarking system uses a watermark f2 and a secret key f3 for the purpose of conducting two levels of security, by using the special images as the watermark and the key, with the same size as the original image f1, to be embedded. Both the watermark f2 and the key f3 are inserted in the spatial domain of the original image f1. The watermarked image g1 is a linear mixture of the original image, key and watermark. That is,

g1(m, n) = f1(m, n) + a f2(m, n) + b f3(m, n) (1)
where a and b are the weighting coefficients [START_REF] Yu | Watermark detection and extraction using independent component analysis[END_REF]. To assure the identifiability of BSS model, it is required that the number of observed linear mixture inputs is at least equal to or larger than the number of independent sources. For the proposed watermark extraction scheme, at least three linear mixtures of the three independent sources are needed. Using the key image f3 and with the help of original image f1, two more mixed images are generated by adding them into the watermarked image

g2(m, n) = c g1(m, n) + d f3(m, n) (2) g3(m, n) = k g1(m, n) + l f1(m, n) (3)
where {c, d, k, l} are arbitrary real numbers. The latter mixtures can be modeled by the following linear system:

g(m, n) = Af (m, n) (4)
where,

f (m, n) = [f1(m, n), • • • , fN (m, n)] T is a N × 1 (with N = 3
) image source vector consisting of the stack of corresponding pixels of source images, A is the M × N full column rank mixing matrix (here,

M = N = 3), g(m, n) = [g1(m, n), • • • , gM (m, n)]
T is an M × 1 vector of mixture image pixels and the superscript T denotes the transpose operator. The purpose of blind image separation is to find a separating matrix, i.e. a N × M matrix B such that f (m, n) = Bg(m, n) is an estimate of original images.

BLIND WATERMARK EXTRACTION

As shown in [START_REF] Bronstein | Separation of reflections via sparse ICA[END_REF][START_REF] Zibulevsky | Blind source separation by sparse decomposition in signal dictionary[END_REF], exploiting the sparsity of some representations of the original images afford us to achieve the BSS problem. Indeed, the mixture destroys or 'reduces' the sparsity of the considered signals that is restored after source separation. Reversely, it is shown in [START_REF] Bronstein | Separation of reflections via sparse ICA[END_REF][START_REF] Zibulevsky | Blind source separation by sparse decomposition in signal dictionary[END_REF] that restoring (maximizing) the sparsity leads to the desired source separation. Based on this, we propose in the sequel a two-step BSS solution consisting in a linear pre-treatment that transforms the original sources into sparse signals followed by a BSS algorithm that minimizes the cost function of the transformed image mixtures using natural gradient technique.

Image pre-treatment

The algorithms proposed in this article are efficients for separating sparse sources. For some signals, one can assume that the spatial or temporal representation is naturally sparse, whereas for natural scenes, this assumptions falls down. We propose to make the image sparse by simply taking into account its Laplacian transform:

F = ∇f = ∂ 2 f ∂x 2 + ∂ 2 f ∂y 2 , (5)
or, in discrete form

F (m, n) = f (m + 1, n) + f (m -1, n) + f (m, n + 1) +f (m, n -1) -4f (m, n) .
Our motivation for choosing this transformation is two fold. First the Laplacian transform is a sparse representation of the image since it acts as an edge detector which provides a two-level image, the edges and the homogeneous background. Second, the Laplacian is a linear transformation. This latter property is 'critical' since the separation matrix estimated to separate the image mixtures is the same to separate the mixture of Laplacian images:

G = ∂ 2 Af ∂x 2 + ∂ 2 Af ∂y 2 = AF ( 6
)
where G is the Laplacian transform of the mixtures. In the literature, some other linear transformations were proposed in order to make the image sparse, including the projection into a sparse dictionary [START_REF] Zibulevsky | Blind source separation by sparse decomposition in signal dictionary[END_REF]. In Figure 1, the original cameraman image is displayed as well as its Laplacian transform and their respective histograms that clearly show the sparsity of the latter.

In the pre-treatment phase, we also propose an optional whitening step which aims to set the mixtures to the same energy level. Furthermore, this procedures reduces the number of parameters to be estimated. More precisely, the whitening step is applied to the Laplacian image mixtures before using our separation algorithm. The whitening is achieved by applying a N × M matrix Q to the Laplacian image mixtures in such a way Cov(QG) = I in the noiseless case, where Cov(•) stands for the covariance operator. As shown in [START_REF] Belouchrani | A blind source separation technique using second-order statistics[END_REF],

Q can be computed as the inverse square root of the noiseless covariance matrix of the Laplacian image mixtures (see [START_REF] Belouchrani | A blind source separation technique using second-order statistics[END_REF] for more details). In the following, we apply our separation algorithm on the whitened data:

G w (m, n) = QG(m, n).

Sparsity-based BSS algorithm

In this section, we propose an iterative algorithm for the separation of sparse signals, namely the ISBS for Iterative Sparse Blind Separation algorithm. It is well known that Laplacian image transform is characterized by its sparsity property in the spatial domain [START_REF] Cichocki | Adaptive Blind Signal and Image Processing[END_REF][START_REF] Zibulevsky | Sparse source separation with relative Newton method[END_REF]. This property can be measured by the p norm where 0 ≤ p < 2. More specifically, one can define the following sparsity based contrast function,

Gp(F ) = N i=1 [Jp(Fi)] 1 p (7)
where

Jp(Fi) = 1 m f n f m f m=1 n f n=1 |Fi(m, n)| p . (8)
The algorithm finds a separating matrix B such as,

B = arg min B {Gp(B)} (9)
where Gp(B) Gp(H) [START_REF] Zibulevsky | Sparse source separation with relative Newton method[END_REF] and H(m, n) BG w (m, n) represents the estimated image sources Laplacian. The approach we choose to solve (9) is inspired from [START_REF] Pham | Blind separation of mixture of independent sources through a quasi-maximum likelihood approach[END_REF]. It is a block technique based on the processing of m f n f observed image pixels and consists in searching the minimum of the sample version of [START_REF] Cichocki | Adaptive Blind Signal and Image Processing[END_REF]. Solutions are obtained iteratively in the form:

B (k+1) = (I + (k) )B (k) (11)
H (k+1) (m, n) = (I + (k) )H (k) (m, n) . ( 12
)
At iteration k, a matrix (k) is determined from a local linearization of Gp(BG w ). It is an approximate Newton technique with the benefit that (k) can be very simply computed (no Hessian inversion) under the additional assumption that B (k) is close to a separating matrix. This procedure is illustrated in the following. At the (k+1) th iteration, the proposed criterion ( 8) can be developed as follows:

Jp(H (k+1) i ) = 1 m f n f m f ,n f m,n=1 H (k) i (m, n) + N j=1 (k) ij H (k) j (m, n) p
Under the assumption that B (k) is close to a separating matrix, we have

| (k) ij | 1
and thus, a first order approximation of Jp(H

(k+1) i
) is given by:

Jp(H (k+1) i ) ≈ 1 m f n f m f ,n f m,n=1 |H (k) i (m, n)| p + p N j=1 (k) ij |H (k) i (m, n)| p-1 sign H (k) i (m, n) H (k) j (m, n) (13
) where sign(•) represents the sign value operator. Using equation (13), equation ( 7) can be rewritten in more compact form as:

Gp B (k+1) = Gp B (k) + T r (k) R (k)T D (k) ( 14
)
where T r(•) is the matrix trace operator, the ij th entry of matrix R (k) is given by:

R (k) ij = 1 m f n f m f ,n f m,n=1 |H (k) i (m, n)| p-1 sign H (k) i (m, n) H (k) j (m, n) .
and

D (k) = diag R (k) 11 , . . . , R (k) N N 1 p -1 . (15)
Using a gradient technique, (k) can be written as:

(k) = -µD (k) R (k) , (16)
where µ > 0 is the descent step. Replacing (16) into (14) leads to,

Gp B (k+1) = Gp B (k) -µ D (k) R (k) 2 , (17)
so µ controls the decrement of the criterion. Now, to avoid the algorithm's convergence to the trivial solution B = 0, one normalizes the outputs of the separating matrix to unit-power, i.e. ρ (k+1)

H i E |H (k+1) i (m, n)| 2 = 1, ∀ i.
Using first order approximation, this normalization leads to:

(k) ii = 1 -ρ (k) H i 2ρ (k) H i . ( 18
)
The final estimated separation matrix B = B (K) Q is applied to the image mixtures g to obtain an estimation of the original images. K denotes here the number of iterations that can be either chosen a priori or given by a stopping criterion the form B (k+1) -B (k) < δ where δ is a small threshold value.

New sparsity-based BSS algorithm

When the images are sparse or sparsely representable, a smooth approximation of the absolute value function might be a better choice for the cost function [START_REF] Bronstein | Blind deconvolution of images using optimal sparse representations[END_REF]. We, therefore, focus our attention on modeling distributions of sparse images using a family of convex smooth functions and propose the following cost function :

G λ (F ) = N i=1 J λ (Fi) (19)
with

J λ (F i ) = 1 m f n f m f ,n f m,n=1 |F i (m, n)| -λ log 1 + |F i (m, n)| λ ( 20
)
where λ is a positive smoothing parameter [START_REF] Bronstein | Blind deconvolution of images using optimal sparse representations[END_REF]. The algorithm finds a separating matrix B such as,

B = arg min B G λ (B) (21)
where

G λ (B) G λ (H) . (22)
In the same way as the ISBS algorithm presented in the previous section, the solutions are obtained iteratively as describe by the equations [START_REF] Pham | Blind separation of mixture of independent sources through a quasi-maximum likelihood approach[END_REF] and [START_REF] Bronstein | Blind deconvolution of images using optimal sparse representations[END_REF]. Therefore, at the (k+1) th iteration, the proposed criterion (20) can be developed as follows :

J λ (H (k+1) i ) = 1 m f n f m f ,nf m,n=1 H (k) i (m, n) + N j=1 (k) ij H (k) j (m, n) -λ log 1 + H (k) i (m,n)+ N j=1 (k) ij H (k) j (m,n) λ .
Under the assumption that B (k) is close to a separating matrix, we have

| (k) ij | 1
and thus, by using a first order approximation of J λ (H (k+1) i

), we can rewrite equation (19) in more compact form as:

G λ B (k+1) = G λ B (k) + T r (k) R (k)T (23)
where the ij th entry of matrix R (k) is given by: R

k) ij = 1 m f n f m f ,nf m,n=1 sign H (k) i (m, n) H (k) j (m, n) λ + |H (k) i (m, n)| . (
Using a gradient technique, (k) can be written as:

(k) = -µ R (k) , (24)
Replacing ( 24) into (23) leads to,

G λ B (k+1) = G λ B (k) -µ R (k) 2 . (25)
Now, to avoid the algorithm's convergence to the trivial solution B = 0, one normalizes the outputs of the separating matrix to unitpower (see equation ( 18)). 1 shows the PSNR-WAV between the original and extracted images for the example described in Fig. 2. We compare the performance of the proposed algorithms to the ICA algorithm. It is clearly shown that our algorithms (ISBS and New ISBS) perform better in terms of PSNR-WAV.

CONCLUSION

In this paper, a new watermarking technique based on BSS algorithm using sparsity property of images has been proposed. The proposed technique consists in a sparsification of the natural observed mixtures followed by a blind separation of the original images. The

Fig. 1 .

 1 Fig. 1. (a) Original image, (b) Laplacian transform, (c) Original image histogram, (d) Sparse Laplacian transform histogram

Fig. 2 .

 2 Fig. 2. (a) watermarked image, (b) and (c) generated mixture images, (d) extracted watermark, (e) extracted image, (f) extracted key.

Table 1 .

 1 Performance evaluation : Comparison of the PSNR-WAV for ICA, ISBS and NISBS algorithms sparsification is simply the Laplacian transform and has a low computational cost. The separation is performed using an iterative algorithm based on the minimizing of the sparsity cost function of the Laplacian image.

			PSNR-WAV (dB)
	BSS
	Algorithm	Parrot	Key	Watermark
	ICA	28.76	15.85	22.36
	ISBS	33.12	19.32	25.67
	NISBS	40.70	28.92	34.49


François-Xavier Socheleau
email: fx.socheleau@telecom-bretagne.eu

Sébastien Houcke
email: sebastien.houcke@telecom-bretagne.eu

Abdeldjalil Aissa

El Bey

Franc ¸ois-Xavier Socheleau

Abdeldjalil Aissa-El-Bey
email: abdeldjalil.aissaelbey@telecom-bretagne.eu

Philippe Ciblat
email: philippe.ciblat@telecom-paristech.fr

OFDM system identification based on m-sequence signatures in cognitive radio context

teaching and research institutions in France or abroad, or from public or private research centers.

I. INTRODUCTION

The increasing demand of wireless services faced to the limited spectrum resources constrains wireless systems to evolve towards more embedded intelligence. The Cognitive Radio (CR) concept [START_REF] Mitola | Cognitive radio: making software radios more personal[END_REF] appears as a key solution to make different systems coexist in the same frequency band. CR terminals have the ability to reconfigure themselves (i.e. to adapt the modulation parameters, carrier frequency, power, etc.) with regards to the surrounding radio environment and spectrum policy. Spectrum sensing and especially system identification is therefore a crucial step towards radio environment awareness. In this paper we focus on OFDM based systems as it becomes the physical layer for many wireless standards [START_REF]Air Interface for Fixed Broadband Wireless Access Systems[END_REF]- [5]. Identification of such systems has mainly been studied using various OFDM cyclostationary properties induced by the Cyclic Prefix (CP) [START_REF] Yucek | OFDM Signal Identification and Transmission Parameter Estimation for Cognitive Radio Applications[END_REF]- [START_REF] Jallon | An algorithm for detection of DVB-T signals based on their second order statistics[END_REF]. The performance of this approach fully depends on the cyclic prefix duration and on the length of the multipath propagation channel. Note also that such a method is totally inefficient in the presence of a zero-padded OFDM system which may be relevant in a cognitive context [START_REF] Khambekar | Utilizing OFDM Guard Interval for Spectrum Sensing[END_REF]. Moreover, considering the increasing interest in OFDM by the wireless designers, cyclostationary properties of such systems are likely to become closer and closer. For instance 3GPP/LTE [5] and Mobile WiMAX [START_REF]Air Interface for Fixed and Mobile Broadband Wireless Access Systems, Amendment 2: Physical and Medium Access Control Layers for Combined Fixed and Mobile Operations in License Bands and Corrigendum 1[END_REF] systems have already an intercarrier spacing only different from 4% which may prevent from getting an accurate system identification based on the intercarrier spacing estimation principle.

To overcome these limitations, methods that involve more particular signatures in OFDM systems are required. [START_REF] Maeda | Recognition Among OFDM-Based Systems Utilizing Cyclostationarity-Inducing Transmission[END_REF] suggested approaches using specific preambles or dedicated subcarriers with cyclostationary patterns. Preambles being usually intermittently transmitted and therefore difficult to intercept, the use of dedicated subcarriers is preferable. However, dedicating subcarriers to only embed signatures may have a cost as it adds overhead and thus reduces systems capacity. One way to address this issue is to jointly use pilot tones for synchronization and/or channel estimation (their initial usage) as well as for system identification.

In this contribution, we develop a solution relying on the use of m-sequence (MS) modulated pilot tones to embed signatures in OFDM signals. MSs show specific high order statistics relevant for system identification and avoid any additional overhead as they meet the requirements of usual training sequences (such m-sequences are indeed already used in existing standards for synchronization and/or channel estimation, see WiMAX [START_REF]Air Interface for Fixed Broadband Wireless Access Systems[END_REF] and Wifi [START_REF]IEEE Std. 802.11a, Part 11: Wireless LAN Medium Access Control (MAC) and Physical Layer (PHY) specifications High-speed Physical Layer in the 5 GHz Band[END_REF]). Thus, we here suggest to take advantage of signatures created as a side-effect of existing MS structures to identify standards such as [START_REF]Air Interface for Fixed Broadband Wireless Access Systems[END_REF] and [START_REF]IEEE Std. 802.11a, Part 11: Wireless LAN Medium Access Control (MAC) and Physical Layer (PHY) specifications High-speed Physical Layer in the 5 GHz Band[END_REF] and also advocate to generalize MSs use in a cognitive context.

The paper is organised as follows: Section II describes the MS useful properties. Section III introduces the OFDM identification scheme and especially the associated cost function based on MS signature characterization. In section IV the impact of synchronization impairments is analysed. Identification performance is assessed through simulations in Section V. Finally, conclusions are presented in Section VI.

II. MS PROPERTIES

A maximum length sequence, commonly called msequence, is a type of pseudorandom binary sequence generated using maximal linear feedback shift registers and modulo 2 addition. A necessary and sufficient condition that a sequence be of maximal length (i.e. sequence of length 2 p -1 for length-p registers) is that its corresponding generator polynomial, denoted by P M S , be primitive.

Any binary MS w k generated by length p registers of polynomial P M S = p i=0 α i X i verifies over GF(2)

w k = p i=1 α i w k-i ⇔ p i=0 α i w k-i = 0.
Moreover, let wk = 1 -2w k be the BPSK associated sequence. Thanks to [START_REF] Macwilliams | Pseudo-Random Sequences and Arrays[END_REF], one can see that

lim M →+∞ 1 M M -1 k=0 i∈B wk-i = 1, if B = A 1/(1 -2 p ) otherwise (1) where B is any subset of {0, • • • , 2 p -2} and A = {i ∈ {0, • • • , p}|α i = 1}.

III. OFDM SYSTEM IDENTIFICATION ALGORITHM

A. System model

In order to facilitate the identification of OFDM systems, we suggest to generalize the use of m-sequence modulated pilot tones. Assuming that a transmitted OFDM symbol consists of N subcarriers and N p comb-type pilot tones, the discrete-time baseband equivalent signal is given by

x(m) = E s N x d (m) + x r (m) , (2)
where

x d (m) = k∈Z N -1 n=0 n / ∈Ip a k (n)e 2iπ n N [m-D-k(N +D)] g[m-k(N +D)],
and

x r (m) = k∈Z n∈Ip wk (n)e 2iπ n N [m-D-k(N +D)] g[m-k(N +D)].
E s is the signal power, a k (n) are the transmit data symbols assumed to be independent and identically distributed (i.i.d), D is the CP length, m → g(m) is the pulse shaping filter. I p denotes the set of pilot subcarrier indexes. For each n ∈ I p , wk (n) is a BPSK pilot symbols sequence associated with one m-sequence obtained by the generator polynomial P M S (n). System signature is thus entirely characterized by I p and

{P M S (n)} n=0,••• ,Np-1 .
Notice that the number of primitive polynomials of degree p over GF(2) is given by φ(2 p -1)/p, where φ(.) is Euler's Totient function [START_REF] Morgan | Primitive Normal Polynomials Over Finite Fields[END_REF]. As an example, for N p = 1 and p ≤ 11, there are 335 different possible signatures which is much larger than the number of existing OFDM systems! Consequently each existing or future system may have its own system signature based on the knowledge of

I p and {P M S (n)} n=0,••• ,Np-1 .
Now, we consider that the signal propagates through a multipath channel. Let {h(l)} l=0,••• ,L be the base-band equivalent discrete-time channel impulse response of length L. The received samples of the OFDM signal are thus given by

y(m) = e -i(2πε m-τ N +θ) L-1 l=0 h(l)x(m -l -τ ) + η(m), ( 3
)
where ε is the carrier frequency offset (normalized by the intercarrier spacing), θ the initial arbitrary carrier phase, τ the timing offset and η(m) a zero mean circularly-symmetric complex-valued white Gaussian noise of variance σ 2 per complex dimension.

The k-th received symbol on subcarrier n is therefore written as

Y k (n) = 1 √ N N -1 m=0 y[k(N + D) + D + m]e -2iπ nm N .
In the case of perfect synchronization (i.e. ε = 0, τ = 0 and θ = 0) and for

n ∈ I p , Y k (n) simplifies to Y k (n) = H k (n) wk (n) E s + N k (n)
where H k (n) and N k (n) are respectively the channel frequency response and the noise at subcarrier n of the k-th received symbol.

B. Identification cost function

Thanks to Eq. ( 1), systems described in Eq. ( 2) can be discriminated by using the following criterion

J J = n∈Ip lim M →+∞ 1 M M -1 k=0   i∈A(n) wk-i (n)   2
where A(n) is the set of indexes associated with the nonnull components of P M S (n). For the sake of simplicity, we consider that the same MS generator polynomial is used for all pilot subcarriers. Consequently A(n) and P M S (n) are independent of n. Moreover, we limit P M S to trinomials of the form 1 + X l + X p (p > l). 1J can now be written as

J = n∈Ip lim M →+∞ 1 M M -1 k=0 wk (n) wk-l (n) wk-p (n) 2 .
In practice, J cannot be computed and the sequence wk (n) is only accessible via the observations Y k (n). Thus, the cost function is based on the following estimate

Ĵ = n∈Ip |Z(n)| 2 , (4)
where Z(n) is defined as

Z(n) = 1 M -p M -1 k=p Ỹk (n) Ỹ * k-l (n) Ỹk-p (n). (5)
M is the number of available OFDM symbols and the superscript "*" stands for complex conjugation which is added on the second term to mitigate the influence of the frequency offset ǫ on Ĵ (for more details, see Section IV). Moreover, in order to get the criterion Ĵ independent of the received signal gain, each term Y k (n) in Eq. ( 5) is normalized so that

Ỹk (n) = Y k (n) Var [Y (n)] , (6)
where Var[.] denotes the variance and

Var [Y (n)] = 1 M M -1 k=0 |Y k (n)| 2 . (7)

C. Decision statistics

In this section we assume perfect synchronization (i.e. ε = 0, τ = 0 and θ = 0). Synchronization impairments are studied in Section IV.

Our identification problem described in the previous subsection boils down to a standard detection problem for which we have to select the most likely hypothesis between the following two hypotheses

          
H 0 : y(m) writes as in Eq. (3) without MS structure or with MS structure for tones in I p associated with P ′ M S = 1 + X l ′ + X p ′ and (l ′ , p ′ ) = (l, p) H 1 : y(m) writes as in Eq. (3) with MS structure for tones in I p described by

P M S = 1 + X l + X p . ( 8
)
To decide the most likely hypothesis, we propose a detection test constrained by the asymptotic false alarm probability similar to what is suggested in [START_REF] Jallon | An algorithm for detection of DVB-T signals based on their second order statistics[END_REF]. The decision is made by comparing Ĵ to a positive threshold such that

ĴH1 > < H0 Λ,
with Λ defined as

F Ĵ|H0 (Λ) = 1 -P f a . (9)
F Ĵ|H0 is the cumulative distribution of Ĵ when H 0 holds and P f a is the tolerated false alarm probability.

As implied in ( 8), H 0 embodies two different subhypotheses respectively named H a 0 and H b 0 . H a 0 represents the case where y(m) does not have any MS structure that is to say that y(m) is any signal such that Y k (n), Y * k-l (n) and Y k-p (n) are mutually independent and E [Y k (n)] = 0. As for H b 0 , it corresponds to the case where the tones of y(m) belonging to the set I p follow a MS structure with generator polynomial P ′ M S = 1 + X l ′ + X p ′ ( = P M S ) where l ′ = l and/or p ′ = p. Thanks to the random, independent and centered nature of the vast majority of digital modulated signals, we can make the reasonable assumption that H 0 = H a 0 ∪ H b 0 . In order to find a relevant threshold Λ, we hereafter analyse the asymptotic statistical behavior of Ĵ under both hypotheses H a 0 and H b 0 . 1) Asymptotic probability density function of Ĵ under H a 0 : As shown in Eq.( 6), Ỹk (n) is expressed as a ratio of two random variables. The variance estimator introduced in Eq.( 7) being consistent, it converges almost surely to a constant denoted v n so that, thanks to the asymptotic theory developped in [START_REF] Brockwell | Time Series: theory and analysis[END_REF], Ỹk (n) converges in distribution to Y (n)/ √ v n . Moreover, Z(n) being a sum of i.i.d random variables when H a 0 holds, we deduce that Z(n)|H a 0 is asymptotically normal with

E [Z(n)|H a 0 ] = 0, Var [Z(n)|H a 0 ] = E |H k (n)| 2 ρ(n) + σ 2 /N 3 (M -p)v 3 n .
where ρ(n) is the signal power of subcarrier n. If we consider the multipath channel as static over the observation window (impacts of channel variation are discussed in section V) then

v n = |H(n)| 2 ρ(n) + σ 2 /N and Var [Z(n)|H a 0 ] simplifies to Var [Z(n)|H a 0 ] = 1 M -p .
Therefore, the asymptotic probability density function of Ĵ under H a 0 is given by 2(M -p) Ĵ ∼ χ 2 2Np , where χ 2 d denotes a chi-square distribution with d degrees of freedom.

2) Asymptotic probability density function of Ĵ under H b 0 : According to Eq. ( 1),

lim M →+∞ 1 M M -1 k=0 wk (n) wk-l (n) wk-p (n) = 1 1 -2 p ′
when H b 0 holds. Following the same approach described in Sec. III-C1, we then have

E Z(n)|H b 0 = |H(n)| 2 H(n)ρ(n) 3 2 (1 -2 p ′ )v 3 2 n , Var Z(n)|H b 0 = 1 (M -p) 2 v 3 n M -p-1 i,j=0
[C] i,j .

[C] i,j are the elements of the covariance matrix defined as

C = E (Υ -E{Υ})(Υ -E{Υ}) H
where the superscript H stands for transpose conjugate and

Υ = [Y k (n)Y * k-l (n)Y k-p (n), Y k+1 (n)Y * k-l+1 (n)Y k-p+1 (n), • • • , Y k+M -p (n)Y * k+M -p-l (n)Y k+M -2p (n)]
. By developing each product term of the covariance matrix and assuming that M ≤ 2 p ′ -1, we get

[C] i,j =          |H(n)| 2 ρ(n) + σ 2 /N 3 -|H(n)| 6 ρ(n) 3 (1-2 p ′ ) 2 , i = j σ 2 |H(n)| 4 ρ(n) 2 N (1-2 p ′ ) -|H(n)| 6 ρ(n) 3 2 p ′ (1-2 p ′ ) 2 , |i -j| = p -|H(n)| 6 ρ(n) 3 2 p ′ (1-2 p ′ ) 2
, otherwise.

In a realistic scenario, the probability density function of Ĵ under H b 0 cannot be easily estimated as it depends on H(n), ρ(n), σ 2 and p ′ which are unknown by the receiver. However, in practice MS degrees can be chosen large enough (e.g., p = 11 in [START_REF]Air Interface for Fixed Broadband Wireless Access Systems[END_REF] and [START_REF][END_REF]) to consider the covariance matrix C as diagonal. In that case,

Z(n)|H b 0 is asymptot- ically normal and Var Z(n)|H b 0 is well approximated by 1/(M -p). Furthermore, if we assume that M ≪ 2 p ′ -1 then E Z(n)|H b 0 ≪ Var Z(n)|H b
0 . Therefore, we can consider that Ĵ follows the same cumulative distribution under both hypotheses H a 0 and H b 0 , that is

F Ĵ|H0 (x) = γ (N p , (M -p) x) (N p -1)!
where γ(a, x) is the incomplete gamma function.

IV. EFFECT OF SYNCHRONIZATION IMPAIRMENTS

Timing missynchronization (τ = 0) and/or frequency offset (ε = 0) damage the observations Y k (n) as inter-symbol (ISI) and inter-carrier (ICI) interferences occur [START_REF] Sathananthan | Probability of error calculation of OFDM systems with frequency offset[END_REF]. In addition to interference, ε modifies the phase of

Y k (n)Y * k-l (n)Y k-p (n) and consequently makes E[Z(n)|H 1 ]
decrease (note that the complex conjugation in Eq.( 5) mitigates the phase variation speed). Therefore, as illustrated in Figure 1, the identification algorithm performance decreases dramatically in the case where ε = 0 and/or τ = 0. To overcome this issue, ǫ and τ can be estimated as

[ε, τ ] = argmax (ε,τ ) Ĵ. (10)
The direct use of the identification cost function Ĵ to estimate (ε, τ ) implies to change the detection threshold Λ. In the specific case where Eq. ( 10) is solved using a grid of K points and assuming that the different Ĵ on this grid are mutually independent, Λ is given by

F Ĵ|H0 (Λ) K = 1 -P f a .

V. SIMULATIONS

In the following, all the results are averaged over 1000 Monte Carlo runs. The system to be recognized is the Fixed WiMAX described in [START_REF]Air Interface for Fixed Broadband Wireless Access Systems[END_REF]. WiMAX embeds MS structures of polynomial P M S = 1+X 9 +X 11 . We recall that N = 256 and N p = 8. Unless otherwise stated, N/D = 32. The subcarriers are equipowered. The asymptotic false alarm probability P f a is fixed to 0.01. The Signal-to-Noise Ratio (SNR) is defined as SNR(dB) = 10log 10 E s /σ 2 .

In Figure 2, we plot the correct detection probability versus SNR in the context of AWGN channel. Different synchronization assumptions and various M are considered. We show that the performance of the MS criterion is significantly improved when the observation window increases. Moreover, we observe the impact of the synchronization method based on Eq. [START_REF] Khambekar | Utilizing OFDM Guard Interval for Spectrum Sensing[END_REF]. We see that the loss due to missynchronization decreases when M increases. For the simulation, uniformly distributed random ε and τ were generated with -0.5 ≤ ε ≤ 0.5 and -0.5(N + D) ≤ τ ≤ 0.5(N + D). ε and τ were estimated by maximizing Ĵ over a grid with a step of 4.10 -3 over ε and 0.1(N + D) over τ .

] = Ge -l/β for l = 0, • • • , L and G is chosen such that L l=0 E[|h(l)| 2 ] = 1).
Notice that β approximately corresponds to the root mean square (RMS) delay spread.

In Figure 3, we display the correct detection probability versus SNR for various RMS delay spread. We observe that the more frequency-selective the channel is, the better the performance is. This is due to the fact that Var[|h(l)| 2 ] decreases as the RMS delay spread increases. In Figure 4, we compare the correct detection probability versus SNR between the proposed MS criterion and the standard correlation based method for β = 0.5D and the various CP lengths handled by the WiMAX system.

To compare both methods, we consider that the correlation based detection is correct when

N -δ/2 ≤ argmax v∈[vmin,vmax] M ′ -v-1 m=0 y(m)y * (m + v) ≤ N +δ/2
where M ′ = M (N + D), v min = 32 and v max = 2048 which corresponds to searching systems from 32 to 2048 subcarriers, and where δ is the tolerated error on the subcarrier spacing. We choose δ to be conditioned by the P f a under the white gaussian noise hypothesis such that δ = (v max -v min )P f a . We observe that the proposed algorithm is not dependent on the CP length and outperforms the correlation based method as soon as N/D ≥ 8. Moreover, for a fair comparison, it is important to remind that as long as the MSs are different, our method can discriminate systems with the same intercarrier spacing whereas the correlation algorithm cannot.

In Figure 5, we plot the correct detection probability versus SNR when the frequency-selective channel becomes timevariant. Various values of maximum Doppler frequencies f d have been inspected. We see that our algorithm is quite robust to Doppler spread below 100Hz (at 3GHz, this corresponds to a relative velocity of 36kph) whilst above this frequency, performance degrades significantly.

VI. CONCLUSION

In this paper, we developed a new method based on msequence properties to embed signatures in OFDM systems without adding any overhead to standard pilot tones. We also studied the MS identification cost function and showed that this method exhibits excellent performance and is quite robust to channel impairements. Moreover, simulation results indicate that in addition to stronger discriminating properties, MS identification outperforms classical correlation based method in some relevant contexts.

Fig. 1 .

 1 Fig. 1. Effect of ε and τ on the correct detection probability (SN R = 0 dB, M = 50, P M S = 1 + X 9 + X 11 , Np = 8, P f a = 0.01).

Fig. 2 .

 2 Fig. 2. Effect of SNR, M , ε and τ on the correct detection probability.

Fig. 3 .

 3 Fig. 3. Effect of β on the correct detection probability (M = 50).

Fig. 4 .

 4 Fig. 4. Comparison between correlation based method and MS criterion (M = 50, β = 0.5D).

Fig. 5 .

 5 Fig. 5. Effect of Doppler spread on the correct detection probability (M = 50, β = 0.5D).

Limiting P M S to trinomials reduces the number of possible signatures but does not call the validity of the concept into question. For instance, there are still 19 different possible signatures for Np = 1 and p ≤ 11.