Jan Unkelbach
email: junkelbach@partners.org

Bjoern H Menze

Ender Konukoglu

Florian Dittmann

Nicholas Ayache

Helen Shi

Helen A Shih

Radiotherapy planning for glioblastoma based on a tumor growth model: implications for spatial dose redistribution

  L'archive ouverte pluridisciplinaire

Introduction

Gliomas differ from most other tumors in the sense that they do not form a solid tumor mass with a defined boundary. Gliomas infiltratively grow into the adjacent brain tissue, forming a diffuse spread of tumor cells. Most high grade gliomas, i.e. glioblastoma multiform (GBM), form a solid tumor mass which is locally destructive to the normal brain tissue and mainly consists of tumor cells. Reactive edema typically surrounds this core that histologically reveals largely normal brain tissue with an infiltration of tumor cells. It is well known that tumor cells infiltrate the adjacent brain tissue and can be found several centimeters beyond the enhancing tumor mass that is visible on MRI [START_REF] Kelly | Imaging-based stereotaxic serial biopsies in untreated intracranial glial neoplasms[END_REF][START_REF] Watanabe | Magnetic resonance imaging and histopathology of cerebral gliomas[END_REF][START_REF] Sj Price | Improved delineation of glioma margins and regions of infiltration with the use of diffusion tensor imaging: an image-guided biopsy study[END_REF]. Functional imaging modalities including amino acid Positron Emission Tomography (PET), including FET (Fluoro-Ethyl-Tyrosine) and MET (Methionine) [START_REF] Grosu | L-(methyl-11c) methionine positron emission tomography for target delineation in resected high-grade gliomas before radiotherapy[END_REF][START_REF] Miwa | Discrepancy between lesion distributions on methionine pet and mr images in patients with glioblastoma multiforme: insight from a pet and mr fusion image study[END_REF][START_REF] Niyazi | Fetpet for malignant glioma treatment planning[END_REF][START_REF] Götz | 18F]FET-PET imaging for treatment and response monitoring of radiation therapy in malignant glioma patients a review[END_REF], have potential to improve the definition of the gross tumor volume. However, these modalities also fail to identify areas of low tumor cell infiltration.

In this paper, we analyse a phenomenological tumor growth model based on the Fisher-Kolmogorov equation, which estimates the spatial distribution of tumor cells in the brain. The model is based on the assumption that tumor cells proliferate locally and migrate into neighboring tissue. This can be described mathematically by a reactiondiffusion equation with a logistic growth term. The model reflects two characteristics of glioma:

• An anisotropic spread of infiltrating tumor cells around the central core: The spatial distribution of tumor cells is modulated by anatomical boundaries like the ventricles and the falx, and by a reduced infiltration of gray matter compared to white matter [START_REF] Matsukado | The growth of glioblastoma multiforme (astrocytomas, grades 3 and 4) in neurosurgical practice[END_REF][START_REF] Coons | Anatomy and growth patters of diffuse gliomas[END_REF].

• A continuous fall-off of the tumor cell density with distance from the core and the lack of a defined tumor boundary [START_REF] Kelly | Imaging-based stereotaxic serial biopsies in untreated intracranial glial neoplasms[END_REF][START_REF] Watanabe | Magnetic resonance imaging and histopathology of cerebral gliomas[END_REF][START_REF] Sj Price | Improved delineation of glioma margins and regions of infiltration with the use of diffusion tensor imaging: an image-guided biopsy study[END_REF].

A comprehensive review of the growth characteristics of glioma can be found in [START_REF] Coons | Anatomy and growth patters of diffuse gliomas[END_REF]. The growth model can be used for radiotherapy planning in two stages:

• The shape of the target volume is modified as to match an isoline of the simulated tumor cell density.

• The prescription dose is redistributed based on the local tumor cell density in order to deliver less dose to regions of low tumor cell density, and (possibly) boost the regions with high tumor cell density.

This paper primarily addresses the second aspect, i.e. the implications of the tumor model for determining the optimal dose distribution prescribed to an inhomogeneous tumor. To that end, we couple the Fisher-Kolmogorov tumor growth model with the linear-quadratic cell survival model after irradiation. The first aspect, using the growth model for automatic delineation of a radiotherapy target volume, taking into account growth characteristics, is discussed in an accompanying paper [START_REF] Unkelbach | Radiotherapy planning for glioblastoma based on a tumor growth model: improving target volume deliniation[END_REF].

Clinical data and standard of care

Typical radiotherapy treatment protocols prescribe a homogeneous radiation dose of 60 Gy to the target volume defined as the abnormality on T2 weighted MRI plus approximately a 2 cm margin. Variations of this protocol prescribe a dose of 46 Gy to this volume, and a dose of 60 Gy to a smaller volume defined as the abnormality on T1 weighted MRI post contrast plus a margin of approximately 2 cm. Both strategies are considered standard of care for radiotherapy of glioblastoma [START_REF] Becker | Status Quo -Standard-of-care Medical and Radiation Therapy for Glioblastoma[END_REF]. The current standard of care is mainly driven by the following three observations:

1. Most tumors recur centrally: For the commonly delivered dose of approximately 60 Gy, most tumors recur at the location of the primary tumor [START_REF] Wallner | Patterns of failure following treatment for glioblastoma multiforme and anaplastic astrocytoma[END_REF][START_REF] Berger | The Gliomas[END_REF][START_REF] Chan | Survival and failure patterns of high-grade gliomas after three-dimensional conformal radiotherapy[END_REF][START_REF] Mcdonald | Pattern of failure after limited margin radiotherapy and temozolomide for glioblastoma[END_REF]. In rare cases, the tumor recurs distant from the primary tumor and outside of the treatment volume ‡. Whereas standard of care is based on a 2-3 cm margin around the visible tumor, substantially reduced margins are reported to not change the recurrence patterns [START_REF] Mcdonald | Pattern of failure after limited margin radiotherapy and temozolomide for glioblastoma[END_REF]. Some studies report changes of the recurrence patterns towards more distant recurrences for doses of 90 Gy [START_REF] Nakagawa | High-dose conformal radiotherapy influenced the pattern of failure but did not improve survival in glioblastoma multiforme[END_REF][START_REF] Fitzek | Accelerated fractionated proton/photon irradiation to 90 cobalt gray equivalent for glioblastoma multiforme: results of a phase ii prospective trial[END_REF]. However this was not found in all studies [START_REF] Chan | Survival and failure patterns of high-grade gliomas after three-dimensional conformal radiotherapy[END_REF] and not for doses in the range of 70-80 Gy [START_REF] Souhami | Randomized comparison of stereotactic radiosurgery followed by conventional radiotherapy with carmustine to conventional radiotherapy with carmustine for patients with glioblastoma multiforme: Report of radiation therapy oncology group 93-05 protocol[END_REF].

2. Dose escalation trials fail to demonstrate a benefit for doses above 60 Gy: An early trial comparing 60 Gy to 70 Gy showed no improvement in the high dose arm [START_REF] Df Nelson | Combined modality approach to treatment of malignant gliomas-re-evaluation of rtog 7401/ecog 1374 with long-term follow-up: a joint study of the radiation therapy oncology group and the eastern cooperative oncology group[END_REF]. A randomized clinical trial comparing standard fractionated radiotherapy to standard radiotherapy plus a stereotactic boost failed to demonstrate any improvement in survival [START_REF] Souhami | Randomized comparison of stereotactic radiosurgery followed by conventional radiotherapy with carmustine to conventional radiotherapy with carmustine for patients with glioblastoma multiforme: Report of radiation therapy oncology group 93-05 protocol[END_REF]. Similarly, trials comparing standard fractionated radiotherapy to standard radiotherapy plus brachytherapy failed to demonstrate a significant benefit for dose escalation [START_REF] Selker | The brain tumor cooperative group nih trial 87-01: a randomized comparison of surgery, external radiotherapy, and carmustine versus surgery, interstitial radiotherapy boost, external radiation therapy, and carmustine[END_REF]. Whether doses of 90 Gy yield improvements in survival is unclear. It has been reported in some studies [START_REF] Fitzek | Accelerated fractionated proton/photon irradiation to 90 cobalt gray equivalent for glioblastoma multiforme: results of a phase ii prospective trial[END_REF], but was not confirmed in others [START_REF] Chan | Survival and failure patterns of high-grade gliomas after three-dimensional conformal radiotherapy[END_REF][START_REF] Nakagawa | High-dose conformal radiotherapy influenced the pattern of failure but did not improve survival in glioblastoma multiforme[END_REF].

3. The risk of radiation necrosis increases for doses above 60 Gy: Main side effects of radiation include radiation necrosis as well as brain atrophy, which correlates with fatigue and cognitive decline. For standard fractionation, the QUANTEC review [START_REF] Richard Lawrence | Radiation dose-volume effects in the brain[END_REF] estimates a 5% risk of radiation necrosis at a dose of 72 Gy and a 10% risk at 90 Gy. This is consistent with the results reported in [START_REF] Ruben | Cerebral radiation necrosis: Incidence, outcomes, and risk factors with emphasis on radiation parameters and chemotherapy[END_REF]. For stereotactic radiosurgery treatments, the incidence of radiation necrosis increases with the treatment volume [START_REF] Richard Lawrence | Radiation dose-volume effects in the brain[END_REF]. The risk of radiation necrosis increases if chemotherapy is administered concurrently. It is commonly assumed that cognitive decline and fatigue worsen with increased integral dose, however, data to establish a doseresponse relation is weak [START_REF] Richard Lawrence | Radiation dose-volume effects in the brain[END_REF].

The prescribed dose of 60 Gy in current clinical practice is driven by three main factors: the established benefit of radiotherapy to improve median survival compared to surgery alone [START_REF] Michael D Walker | An analysis of dose-effect relationship in the radiotherapy of malignant gliomas[END_REF][START_REF] Michael D Walker | Evaluation of BCNU and/or radiotherapy in the treatment of anaplastic gliomas: a cooperative clinical trial[END_REF][START_REF] Kristiansen | Combined modality therapy of operated astrocytomas grade III and IV. Confirmation of the value of postoperative irradiation and lack of potentiation of bleomycin on survival time: a prospective multicenter trial of the Scandinavian Glioblastoma Study Group[END_REF]; the negative outcome of dose escalation trials beyond 60 Gy; and the increasing risk of radiation necrosis for larger doses. ‡ Most authors distinguish central, marginal, and distant recurrence. However, the definition of these terms is not always entirely clear and varies between authors.

Approach to optimizing spatial dose prescription

The observation that tumor cells can be found at almost any distance from the central core [START_REF] Kelly | Imaging-based stereotaxic serial biopsies in untreated intracranial glial neoplasms[END_REF][START_REF] Watanabe | Magnetic resonance imaging and histopathology of cerebral gliomas[END_REF][START_REF] Sj Price | Improved delineation of glioma margins and regions of infiltration with the use of diffusion tensor imaging: an image-guided biopsy study[END_REF] questions the concept of a target volume in general, and a homogeneous dose prescription in particular. A possible criterion for determining the ideal dose distribution prescribed to an inhomogeneous tumor is based on the idea of minimizing integral cell survival in the context of the linear-quadratic cell survival model [START_REF] Brahme | Optimal dose distribution for eradication of heterogeneous tumors[END_REF][START_REF] Webb | A model for calculating tumor control probability in radiotherapy including the effects of inhomogeneous distributions of dose and clonogenic cell density[END_REF][START_REF] Webb | A proof that uniform dose gives the greatest TCP for fixed integral dose in the planning target volume[END_REF].

For most solid tumors with a defined boundary, differences in the tumor cell density have often been considered irrelevant. This results from the exponential cell kill behavior: For a surviving fraction of 50% at a standard fraction dose of 2 Gy, a dose of 2 Gy is sufficient to compensate for a factor 2 difference in the initial cell density. Thus, heterogeneities in radiosensitivity (e.g. due to hypoxia) have been considered to be more important than differences in clonogen density. Gliomas play a special role in that respect: The cell density drops by orders of magnitude with distance from the core. For the Fisher-Kolmogorov tumor growth model, the tumor cell density drops approximately exponentially with distance from the solid tumor mass [START_REF] Murray | Mathematical biology: I. An introduction[END_REF][START_REF] Konukoglu | Extrapolating glioma invasion margin in brain magnetic resonance images: Suggesting new irradiation margins[END_REF]. This paper discusses the implications of the exponential fall-off for radiotherapy planning. In particular, this paper addresses the question, how the prescribed dose should be distributed spatially, given a limit on the integral dose delivered to the patient. Determining this integral dose constraint is not only determined by the tumor but depends also on the tolerance of the brain tissue; it it thus not the topic of this paper.

Contributions and organization of this paper

The remainder of this paper is organized as follows: In section 2 we summarize the essential features of the previously published Fisher-Kolmogorov tumor growth model, which are relevant for this work. Section 3 is the main contribution of this paper. Here, we discuss the implications of the tumor growth model regarding the optimal distribution of the prescription dose in the context of the linear-quadratic cell survival model. It is shown that an exponential fall-off of the tumor cell density suggests a linear fall-off of the prescription dose with distance. We introduce a dose fall-off rate parameter, which quantifies the steepness of the prescription dose fall-off, and replaces the concept of a target volume with a defined boundary and a uniform dose prescription.

In section 4, we discuss intensity modulated radiotherapy (IMRT) planning strategies and illustrate the differences in the corresponding dose distribution. In section 5, we discuss the results of this paper in the context of the available clinical data as well as limitations of the model and uncertainty in its parameters.

Tumor growth model

Tumor growth is modeled as local proliferation of tumor cells and migration into neighboring tissue. This is described mathematically via the Fisher-Kolmogorov equation. For a description of the model, the reader is referred to section 2 of the accompanying paper [START_REF] Unkelbach | Radiotherapy planning for glioblastoma based on a tumor growth model: improving target volume deliniation[END_REF] and the references therein. The property of the model that is most relevant to this paper is the exponential fall-off of the tumor cell density. Denoting the tumor cell density by c and the distance from the central tumor mass as |r| we qualitatively have the behaviour

c(|r|) ∝ exp - |r| λ w (1)
within a region of white matter. The parameter λ w is called the infiltration length and denotes the distance at which the cell density drops by a factor of e in white matter.

The infiltration length relates to the model parameters proliferation rate ρ and white matter diffusion coefficient D w via λ w = D w /ρ. Due to anatomical boundaries and a reduced infiltration of gray matter ( modeled through a shorter gray matter infiltration length λ g ≪ λ w ) the distribution of tumor cells is modulated spatially. Equation ( 1) is only valid within white matter, whereas the falloff is steeper in gray matter.

A detailed case study of spatial growth patterns can be found in the accompanying paper [START_REF] Unkelbach | Radiotherapy planning for glioblastoma based on a tumor growth model: improving target volume deliniation[END_REF]. To help the reader, we discuss the same patient in this paper and use the same slices for illustration. To make this paper self-contained, we include figure 1

Infiltration length

Ideally, the infiltration length is estimated for an individual patient based on the available imaging data. One approach to this problem is based on the difference in size between the contrast enhancing core on T1 gadolinium and the hyperintense region on T2-FLAIR [START_REF] Harpold | The evolution of mathematical modeling of glioma proliferation and invasion[END_REF][START_REF] Swanson | A mathematical modelling tool for predicting survival of individual patients following resection of glioblastoma: a proof of principle[END_REF]. To that end, it is assumed that the boundary of the enhancing core corresponds to a high tumor cell density close to the carrying capacity, whereas the hyperintense region on T2-FLAIR is associated with a lower tumor cell infiltration. However, a quantification of the tumor cell density at the T2-FLAIR boundary is problematic: For GBM patients, the hyperintense region on T2-FLAIR contains infiltrating parts of the tumor, but also large regions of vasogenic edema. Therefore, the boundary of the edema may correspond to a wide range of tumor cell densities. Quantification is ultimately based on histopathological analysis after biopsy. In [START_REF] Harpold | The evolution of mathematical modeling of glioma proliferation and invasion[END_REF] and [START_REF] Swanson | A mathematical modelling tool for predicting survival of individual patients following resection of glioblastoma: a proof of principle[END_REF], a tumor cell infiltration between 2% and 10% was assumed, which is consistent with biopsy data presented in [START_REF] Stadlbauer | Improved delineation of brain tumors: an automated method for segmentation based on pathologic changes of 1H-MRSI metabolites in gliomas[END_REF]. This allows for a coarse estimate of the infiltration length: Assuming that the cell density drops exponentially by a factor of exp(3) ≈ 20 between the edge of the enhancing core (where the cell density is close to 100%) and the boundary of the edema (where the cell density is close to 5%), this distance corresponds to three times the infiltration length. Thus, if (in white matter) the edema extends 1 cm beyond the enhancing core, the infiltration length is approximately 3 mm. For an edema that extends 2 cm beyond the core, it is approximately 6 mm. This concept of estimating the infiltration length is illustrated in figure 2 for the patient discussed in this paper. Figure 2a shows T1 post contrast image; in figure 2b, we overlaid the 10% and 2% isoline of the tumor cell density for λ w = 2.1 mm on the T2-FLAIR image. The blue contour is the 70% isoline which corresponds to the boundary of the contrast enhancing core. The value λ w = 2.1 for the infiltration length is appropriate to approximately match the 5% isoline of the tumor cell density to the extent of the edema in the left-anterior direction. Figure 2b also illustrates limitations due to the unspecific nature of the hyperintensity of the T2-FLAIR image: the edema does not extend equally far in all directions, even within white matter. Therefore, other imaging modalities including PET may prove valuable in improving the estimation of the infiltration length in the future (see section 5.2 for further discussion).

Redistribution of prescription dose

The desired radiotherapy dose distribution is determined by the conflicting objectives of minimizing side effects of radiation and maximizing time of survival. Major side effects of radiation include radiation necrosis and cognitive decline. The probability for developing radiation necrosis is a function of the maximum dose delivered to the target volume [START_REF] Richard Lawrence | Radiation dose-volume effects in the brain[END_REF]. Other side effects, like brain atrophy and subsequent cognitive decline, may depend on the integral dose delivered to the brain tissue. However, the published data still lacks a clear demonstration of cognitive decline in adults according to the recent Quantec review [START_REF] Richard Lawrence | Radiation dose-volume effects in the brain[END_REF]. For treatment planning, these objectives are formulated mathematically. For this section, we assume that in order to maximize time of patient survival, treatment planning should aim at minimizing integral tumor cell survival after the end of treatment. Radiation side effects are modeled via a maximum dose constraint and an integral dose constraint.

Cell survival after radiotherapy

Cell survival after radiation is commonly described by the linear-quadratic cell survival model with radiosensitivity parameters α and β. The density of surviving tumor cells s after irradiation with a cumulative dose d in N f fractions is given by

s = c exp (-αd(1 + κd)) (2)
where we have defined κ = β αN f . The α/β ratio for brain tumors is assumed to be in the order of 10 as for many other tumors. In this paper, we do not investigate effects of fractionation but consider absolute cell kill. In that case, the factor (1 + κd i ) can be considered as a small correction factor to the main radiosensitivity parameter α. We want to neglect the dose dependence of the correction factor and evaluate the correction factor at the mean prescription dose d p to obtain the effective radiosensitivity

ᾱ = α(1 + κ d p ) (3)
Assuming that the mean dose prescribed to the target volume is 60 Gy in 30 fractions, the correction factor is 1.2. If the dose in the target volume varies by ±50%, i.e. between 30 Gy and 90 Gy, the effective radiosensitivity varies by approximately 10%. This variation is relatively small compared to the overall uncertainty in the radiosensitivity parameter. Below, we therefore work with this approximation. Plausible values of the radiosensitivity parameter are discussed in sections 3.4 and 5.2.

Minimizing integral tumor cell survival

For treatment planning, we want to find the dose distribution d that minimizes integral cell survival. This can be formulated as the following optimization problem:

minimize d i∈T c i exp (-ᾱd i ) (4) subject to 1 N T i∈T d i ≤ d p (5)
d i ≥ 0 ∀i ∈ T (6)
where c i is the tumor cell density in voxel i obtained from the tumor growth model, d i is the dose delivered to voxel i, and T refers to the set of all voxels in the target volume. Equation 5 represents a constraint on the integral dose with the mean prescription dose d p . Without this constraint, the problem becomes unbounded, i.e. the optimal prescription dose is infinite.

Qualitative behavior

The optimal solution to the optimization problem (4-6) can be written in closed form. We note that the integral dose constraint (5) will be binding. We thus seek for stationary points of the Lagrange function

L = i∈T c i exp (-ᾱd i ) + µ i∈T d i -N T d p (7)
Setting the derivative of L with respect to d i to zero yields the optimal dose prescription

d i = max 0, 1 ᾱ ln ᾱc i µ (8)
The constant µ > 0 corresponds to the lagrange multiplier that arises from the integral constraint (5) and hence depends on the mean prescription dose d p . µ turns out to be given by

1 ᾱ ln ᾱ µ = d p - 1 ᾱN T i∈ T ln (c i ) (9)
where T denotes the set of target voxels with nonzero dose prescription, i.e. the set of voxels at which the positivity constraint ( 6) is not active. The optimal prescription dose is thus proportional to the logarithm of the tumor cell density. This implies that the optimal prescription corresponds to the dose distribution that yields homogeneous cell survival after therapy within the region of initially high enough cell density.

Within the tumor growth model, the tumor cell density drops approximately exponentially with distance from the visible tumor. Thus, the logarithm of the cell density falls off approximately linear. Due to the linear dependence of dose on logarithmic cell density, the optimal prescription dose falls off linearly with distance from the gross tumor volume. If we denote the distance from the visible tumor core by |r|, the normalized tumor cell density is approximately given by

c(|r|) = exp - |r| λ w (10)
with infiltration length λ w . Inserting this expression into the formula for the optimal prescription dose (8), we obtain

d(|r|) = d p - 1 ᾱN T i∈ T ln c i d p max - 1 λ w ᾱ |r| (11)
where d p max is the maximum prescribed dose to the center of the tumor with highest tumor cell density, which is larger than the mean prescribed dose d p . We conclude that the dose drops linearly with distance from the gross tumor, and that the proportionality constant is given by the inverse of the product of infiltration length and radiosensitivity. We therefore introduce the dose fall-off rate 1 ᾱλ w

which quantifies the steepness of the prescription dose fall-off in units of Gy per mm. The optimal prescription dose drops fast if the tumor cells are radioresistant or the tumor cell density drops fast. This result is intuitive: If ᾱ is large, i.e. tumor cells are radiosensitive, a relatively small dose is sufficient to compensate for initial differences in the tumor cell density and the dose fall-off is shallow. For small values of ᾱ, large doses are required to compensate for initial differences in cell density. Thus, the dose is more concentrated in the center of the tumor and the fall-off with distance is steeper.

Quantifying the dose fall-off rate

In order to quantify the dose fall-off rate, we require an estimate of the radiosensitivity.

The parameter ᾱ has been determined in in-vitro experiments for different human glioblastoma cell lines. The resulting values vary over a wide range between 0.05 and 0.4 for different cell lines under different conditions (see [START_REF] Dionysiou | Critical parameters determining standard radiotherapy treatment outcome for glioblastoma multiforme: a computer simulation[END_REF] and references therein).

An effective ᾱ = 0.35 corresponds to a surviving fraction of 50% at 2 Gy. Given the above range of radiosensitivity values and the infiltration length estimates in section 2.1, plausible parameters for the dose fall-off range from roughly 10 Gy per mm (assuming ᾱ = 0.05 1/Gy and λ w = 2 mm) to 5 Gy per cm (assuming ᾱ = 0.35 1/Gy and λ w = 6 mm).

Illustration of dose redistribution

We illustrate the optimal dose distribution for the patient shown in figure 1. For that purpose, we first optimize an IMRT reference plan based on a target volume that is prescribed to a homogeneous dose of 60 Gy. The target volume is given by the isoline of the tumor cell density that yields the same total volume of the target as the contour that was drawn manually in the clinically applied plan. The treatment plan optimization problem is specified in section 4. For this plan, we calculate the mean dose delivered to brain tissue (i.e. any white matter, gray matter, or target voxel). This provides the mean prescribed dose d p in equation ( 5), considering all white and gray matter voxels as the set of target voxels T .

Figure 3a shows the optimal dose distribution in brain tissue that minimizes integral cell survival for λ w = 4.2 mm and ᾱ = 0.35 1/Gy. For these parameters, the dose fall-off rate according to (12) is 7 Gy per cm, which is approximately seen in figure 3a. The maximum prescribed dose in the center of the tumor is approximately 70 Gy, and the dose at the edge of the CTV is lowered from 60 Gy to approximately 40 Gy.

For smaller values of the infiltration length or the radiosensitivity parameter, the dose fall-off is steeper. For the same integral dose, this leads to higher doses in the center of the tumor. Figure 3b shows the optimal dose distribution for the parameters λ w = 2.1 mm and ᾱ = 0.35 1/Gy. This leads to a maximum dose of approximately 90 Gy in the center of the tumor and a dose fall-off rate of 14 Gy per cm. Almost the same result is obtained for an infiltration length of λ w = 4.2 mm and a radiosensitivity of ᾱ = 0.17.

For even smaller values of infiltration length and radiosensitivity parameter, the dose fall-off becomes steeper, suggesting even higher doses to core. In the limit, all the dose is concentrated to the region with highest cell density. In contrast, for very large values of infiltration length and radiosensitivity parameter, the optimal dose becomes more spread out over the brain. In the limit, the dose is spread out evenly over the brain tissue.

Including a maximum dose constraint

In current clinical practice, the maximum dose to the tumor is limited to 60 Gy for reasons discussed in section 1.1. We therefore discuss the optimal dose distribution that minimizes integral cell survival given a constraint on the maximum dose (in addition to the constraint on integral dose). For the parameter values used in figures 3a/b, this leads to the optimal dose distributions in figures 3c/d. For these parameters, the maximum dose constraint is binding in the center of the tumor. In the periphery, the prescribed dose is given by equation (8), i.e. the dose fall-off outside the region that is irradiated to 60 Gy is still given by 1/ᾱλ w . This has the following consequence: For very small values of the radiosensitivity parameter, the dose prescription approaches a binary dose prescription in which 60 Gy is delivered to a target volume, with a steep fall-off outside.

Treatment plan comparison

In sections 3.5 and 3.6, we compare the idealized prescription doses. In this section, we consider dose distributions that are physically realizable through intensity-modulated radiotherapy. We formulate the IMRT planning problem in section 4.1; in section 4.2 we illustrate the differences in the dose distributions for the conventional dose prescriptions and the model derived prescription doses. Below, we analyze the patient shown in sections 2 and 3; in the appendix we show results for a second patient with a multifocal GBM located in the temporal and frontal lobe.

IMRT planning formulation

We optimize an IMRT plan using 9 equally spaced coplanar 6MV photon beams and a piece-wise quadratic objective function. The optimization problem is solved using our own implementation of the L-BFGS quasi-newton method; dose-calculation is performed using the pencil-beam dose calculation algorithm, Quadratic Infinite Beam (QIB), embedded in CERR 3.0 Beta 3 (see [START_REF] Deasy | CERR: a computational environment for radiotherapy research[END_REF]). More specifically, we minimize the following objective function

f (d) = η w o η N η i∈Vη d i -d max η 2 + (13) + w u T N T i∈T (d pres i -d i ) 2 + + w o T N T i∈T (d i -d pres i ) 2 + ( 14
) + w u H N H i∈H (d pres i -d i ) 2 + ( 15
) + w o H N H i∈H (d i -d max i ) 2 + ( 16
)
The first term ( 13) denotes overdose objectives for the organs at risk (OAR). For the maximum doses d max η and weighting factors w o η , we use the values summarized in table 1 in the accompanying paper [START_REF] Unkelbach | Radiotherapy planning for glioblastoma based on a tumor growth model: improving target volume deliniation[END_REF]. The second term ( 14) denotes over and under dose objectives within a confined target volume T . For the comparison in this section, T is the model derived target, which has the same total volume as the CTV used in the clinical plan. The third term [START_REF] Mcdonald | Pattern of failure after limited margin radiotherapy and temozolomide for glioblastoma[END_REF] aims at delivering the prescribed dose to voxels outside of T ; and the forth term ( 16) represents a conformity objective that penalizes dose to unclassified voxels outside of T (including skull, brain tissue, ventricles). The conformality objective for an inhomogeneous prescription dose in defined in equation ( 9) in [START_REF] Unkelbach | Radiotherapy planning for glioblastoma based on a tumor growth model: improving target volume deliniation[END_REF].

Dose comparison

Based on the IMRT planning formulation (13-16), we optimize three treatment plans corresponding to three prescription dose distributions. Figure 4a shows the IMRT dose distribution for a homogeneous dose prescription of 60 Gy to the CTV. The prescribed dose outside of the CTV is set to zero. Figure 4c shows the IMRT plan for the prescription dose in figure 3c, which corresponds to a dose falloff rate of 7 Gy per cm and a maximum target dose of 60 Gy. Figure 4d shows the difference of the dose distributions in 4c and the homogeneous prescription in 4a. The dose distribution that uses a falloff rate of 7 Gy per cm leads to approximately the same dose distribution in the posterior region of the target volume. At the edge of the target volume (primarily on the contralateral side and anterior part on the ipsilateral side), the dose is reduced by 10-15 Gy. In regions outside of the target (but adjacent), the dose is increased by up to 15 Gy.

Figure 4e shows the IMRT plan for a the prescription dose in 3a, i.e. a dose falloff rate of 7 Gy per cm and no maximum dose constraint. The maximum prescription dose in the center of the target is approximately 70 Gy. This corresponds to a moderate dose escalation to the central region of the tumor, but keeping the integral dose to the brain constant. The dose difference compared to the homogeneous prescription case is shown in figure 4f. It illustrates dose increases in the center of the tumor and regions adjacent to the target; decreased doses are observed around the edge of the target volume.

In figure 4b we show the plan for the commonly used "cone down" approach in which 60 Gy is prescriped to the boost volume and 46 Gy to the CTV ¶. Comparing the dose distribution in figure 4c to the cone down approach in figure 4b shows that the dose distribution within the CTV is very similar, i.e. a dose fall-off rate of 7 Gy per cm lead to a similar result as the cone down approach. However, the IMRT plan in 4b aims at reducing the dose outside the target, whereas the plan in 4c has a non-zero prescription dose in regions adjacent to the target. Thus, the integral dose in 4c is higher compared to 4b.

Controlling conformity for complex prescription dose distributions

Anatomical barriers and gray matter surrounding sulci lead to complex spatial dose prescriptions. Figure 4 shows that IMRT can reproduce the main features of the idealized prescribed dose distribution. However, an appropriate formulation of the IMRT planning problem is necessary to control the trade-off between matching the prescribed dose distribution and sparing uninvolved healthy tissue. In the formulation in section 4.1, we aim at matching the prescribed dose relatively accurately in the region of the target volume, using a high weighting factor w u T . Outside of the target, the formulation also aims at delivering the prescribed dose, but a separate (lower) weight w u H is assigned. The use of a fall-off objective generalized to non-uniform prescription doses, together with separate weights w u T and w u H where found adequate for controlling the tradeoff between matching the prescription dose and reducing integral dose.

Discussion

Implications for dose prescription

Infiltrative glioblastoma lack a defined boundary. Individual tumor cells that infiltrate the brain parenchyma can be found at several centimeters distance from the contrast enhancing tumor mass. This questions the paradigm of a radiotherapy target volume that is treated to a homogeneous dose. Under the assumption of exponential cell kill and an exponential fall-off of the tumor cell density, minimizing integral tumor cell survival for a constrained integral dose suggest a linear fall-off of the prescription dose with distance from the tumor core. The proposed concept of a dose fall-off rate represents a ¶ This plan is optimized based on the objective detailed in section 4 in [START_REF] Unkelbach | Radiotherapy planning for glioblastoma based on a tumor growth model: improving target volume deliniation[END_REF]. quantitative approach to relate dose prescription to underlying parameters for radiation response (ᾱ) and the infiltrativeness of a tumor (λ). Dose prescription without maximum dose constraint: If no constraint on the maximum dose is imposed, this approach will typically suggest a redistribution of the dose towards higher doses in the center and lower doses at the periphery of the CTV as illustrated in figures 3a and 3b. For small values of the infiltration length and the radiosensitivity parameter, the dose is concentrated in the center of the tumor. The dose fall-off in the periphery is steep for the following reason: for a small infiltration length, the cell density drops fast with distance; and for a small radiosensitivity parameter, large doses are required to compensate for initial differences in the cell density. In contrast, for large values of the infiltration length or the radiosensitivity parameter, the dose fall-off becomes more shallow. This approach is intuitive, especially regarding the fact that the vast majority of tumors recur at the primary site. However, escalating the dose to the central tumor mass is a controversial topic due to the negative outcome of dose escalation trails [START_REF] Souhami | Randomized comparison of stereotactic radiosurgery followed by conventional radiotherapy with carmustine to conventional radiotherapy with carmustine for patients with glioblastoma multiforme: Report of radiation therapy oncology group 93-05 protocol[END_REF]. Therefore, we also discuss the application of dose fall-off rate concept given a maximum dose constraint.

Dose prescription with maximum dose constraint: Under the assumption that the maximum dose is kept at 60 Gy, this approach suggests the delivery of 60 Gy to an inner target, surrounded by a peripheral target in which the prescription dose falls off linearly. For radioresistant tumor cells, the dose fall-off is steep, resulting in an essentially uniform dose prescription (figure 3d). For example, assuming a radiosensitivity of ᾱ = 0.1 and an infiltration length of 3 mm, which appears realistic, the model suggests a dose falloff rate of approximately 30 Gy per cm. This roughly corresponds to the dose fall-off that is naturally achieved for an IMRT treatment plan for a uniform dose prescription -due to the physical characteristics of a photon beam. Thus, under these conditions, the optimal dose prescription is consistent with current clinical practice in which a homogeneous dose of 60 Gy is delivered to the CTV.

Comparison to cone-down approach: The alternative "cone down" approach commonly used, in which 46 Gy is prescribed to the CTV and 60 Gy to a boost volume, corresponds to larger values of either radiosensitivity or infiltration length. For example, assuming a radiosensitivity of α = 0.2 and an infiltration length of 3 mm, which appears realistic for some patients, the model suggests a dose fall-off of approximately 15 Gy per cm. This is consistent with the cone down approach under the assumption that the edema extends about 1 cm beyond the enhancing core and that the same margin is added to obtain the boost volume and the CTV from the enhancing core and the T2-FLAIR abnormality, respectively. In our approach, the linear fall-off of the prescription dose continues beyond the CTV contour prescribed to 46 Gy. Thus, our approach suggests to increase the dose in some regions outside the CTV contour, compared to the dose that is naturally delivered due to the depth dose curve of photon beams.

Model limitations and parameter uncertainty

An application of the dose fall-off rate concept requires an estimation of radiosensitivity and infiltration length. Furthermore, personalizing the prescribed dose distribution to the patient at hand requires methods to estimate these parameters based on patent data acquired prior to treatment. The estimation of both parameters is currently related to large uncertainties.

Estimating radiosensitivity The radiosensitivity parameter ᾱ has been determined in in-vitro experiments for different human glioblastoma cell lines. The resulting values vary over a wide range between 0.05 and 0.4 for different cell lines under different conditions (see [START_REF] Dionysiou | Critical parameters determining standard radiotherapy treatment outcome for glioblastoma multiforme: a computer simulation[END_REF] and references therein). In addition, attempts to estimate a radiosensitivity parameter in vivo from imaging data have been proposed [START_REF] Rockne | Predicting the efficacy of radiotherapy in individual glioblastoma patients in vivo: a mathematical modeling approach[END_REF]. However, such approaches need further validation. The estimation of radiosensitivity based on patient data acquired prior to therapy is further complicated by the heterogeneity of the tumor. Glioblastoma consist of different cell populations with varying radiosensitivity. In particular, the infiltrating tumor cells may behave differently from cells in the tumor core, which itself is heterogeneous, typically showing regions of hypoxia and necrosis. Potential differences in radiosensitivity between core and periphery are relevant in the context of increasing the dose to the core, i.e. the scenario without maximum dose constraint. If, instead, the prescribed dose is limited to 60 Gy, we are only interested in the dose fall-off in the periphery. This is determined by the radiosensitivity of the infiltrating tumor cells rather than the radiosensitivity of the enhancing core.

Estimating the infiltration length Current approaches to estimate the infiltration length are based the difference of contrast enhancing tumor volume and the hyperintense region visible on T2-FLAIR [START_REF] Harpold | The evolution of mathematical modeling of glioma proliferation and invasion[END_REF][START_REF] Swanson | A mathematical modelling tool for predicting survival of individual patients following resection of glioblastoma: a proof of principle[END_REF]. The use of the T2-FLAIR image has inherent limitations due to its unspecific nature and a weak correlation between tumor infiltration and vasogenic edema. This results in large uncertainties in the infiltration length estimation. Modern imaging modalities, in particular PET and MR Spectroscopy have significant potential to improve the assessment of infiltrative disease and the estimation of the infiltration length. In PET imaging, tracers for amino acid metabolism (e.g. MET and FET) have been suggested as surrogares for tumor [START_REF] Grosu | L-(methyl-11c) methionine positron emission tomography for target delineation in resected high-grade gliomas before radiotherapy[END_REF][START_REF] Miwa | Discrepancy between lesion distributions on methionine pet and mr images in patients with glioblastoma multiforme: insight from a pet and mr fusion image study[END_REF][START_REF] Niyazi | Fetpet for malignant glioma treatment planning[END_REF][START_REF] Götz | 18F]FET-PET imaging for treatment and response monitoring of radiation therapy in malignant glioma patients a review[END_REF]; in MR Spectroscopy attempts to correlate the choline over NAA ratio (N-Acetylaspartat) have been made [START_REF] Stadlbauer | Improved delineation of brain tumors: an automated method for segmentation based on pathologic changes of 1H-MRSI metabolites in gliomas[END_REF]. However, PET and MR Spectroscopy both suffer from a low spatial resolution. In addition, more biopsy studies guided by these imaging modalities are required to establish a quantitative relationship between signal intensity and tumor cell density.

Further remarks

Margin expansion versus dose redistribution: It should be noted that the model discussed in this paper only predicts the steepness of the dose fall-off. The volume that is irradiated to the maximum dose, i.e. at what distance from the contrast enhancing core the fall-off is located, is determined by the integral dose constraint. This has to be determined based on a tradeoff between radiation side effects and the risk of marginal recurrences. McDonald [START_REF] Mcdonald | Pattern of failure after limited margin radiotherapy and temozolomide for glioblastoma[END_REF] analyzed outcomes for glioblastoma treatments with substantially reduced target volumes, and reported no deterioration in survival or increase in marginal recurrences compared to historical controls. This suggests that the integral dose, and therefore the volume that is irradiated to 60 Gy, can be reduced compared to current practice.

Low grade gliomas: In this paper, we discussed the implications of the proposed tumor growth model for glioblastoma patients. With minor modifications, the same model is applicable to low grade gliomas, which are assumed to show the same infiltrative growth patterns. Similar to high grade glioma, dose escalation trials have shown a point of diminishing return for doses above approximately 50 Gy. In contrast to glioblastomas, patients receiving radiation therapy for low grade gliomas survive many years to decades. Thus, determining the optimal tradeoff between tumor control and the minimization of integral dose (to reduce long term side effects) may be more important for low grade gliomas than it appears for glioblastoma. In order to apply the to low grade gliomas, a delineation of the tumor core (which is associated with an isoline of the cell density) is needed as input (see [10, section 2] or [START_REF] Konukoglu | Extrapolating glioma invasion margin in brain magnetic resonance images: Suggesting new irradiation margins[END_REF]). Unlike glioblastoma considered in this paper, low grade glioma typically do not show a contrast enhancing core. Consequently, a different imaging modality has to be used for that purpose. This can, for example, be the abnormality on T2-FLAIR.

Conclusions

The Fisher-Kolmogorov glioma growth model yields an exponential fall-off of the tumor cell density with distance from the tumor core. For an exponential cell kill model, minimizing integral tumor cell survival in radiotherapy planning suggests a linear falloff of the prescription dose with distance from the core. We introduce the dose fall-off rate in units of Gy per mm, which quantifies the steepness of the fall-off and is given by the inverse of radiosensitivity times infiltration length. The dose fall-off rate is a concept that reflects the absence of a defined boundary of infiltrating gliomas and could replace the paradigm of a confined target volume with constant dose prescription.

Without a constraint on the maximum dose, the approach suggests (for typical parameter values) a dose increase to the tumor core and a reduction of dose at the periphery where the tumor cell density is low. Although intuitive, dose increases to the core are controversial because clinical dose escalation trials for GBM failed to demonstrate a survival benefit. Therefore, we also discuss the application of the dose fall-off rate concept given a maximum dose constraint.

If the maximum dose is constrained to 60 Gy, the model suggests delivering 60 Gy to an inner target volume, which is surrounded by a linear dose fall-off region. The clinical approach, in which a homogeneous dose of 60 Gy is delivered to the CTV, corresponds to our approach under the assumption of weakly infiltrating tumors and radioresistant tumor cells. The "cone down" approach, in which 46 Gy is delivered to the CTV and 60 Gy to a boost volume, is consistent with our approach for more infiltrating tumors and radiosensitive tumor cells.

The dose fall-off rate represents a quantitative approach to relate dose prescription to underlying parameters for radiation response and the infiltrativeness of a tumor. Current uncertainties in radiosensitivity and infiltration length lead to a large interval of plausible dose fall-off rates. The dose distributions prescribed in current clinical practice fall within this range. However, if better methods to estimate radiosensitivity and infiltration length on a patient-by patient basis became available, the concept could be used to individualize dose prescription. target volume. In regions where the target volume is bounded by CSF or gray matter, the desired dose fall-off is much steeper and not determined by the dose fall-off rate. This is most apparent adjacent to the ventricles and the lateral sulcus.

  , which shows an example of the simulated tumor cell density for a patient with a glioblastoma located in the left parietal lobe §. Parameters for this simulation are λ w = 4.2 mm and λ w /λ g = 10. The figure illustrates the spatial modulation of tumor growth through the anatomical boundaries (falx cerebri and ventricles) and through the fast fall-off of the tumor cell density in gray matter. It also shows the approximately linear fall-off of the logarithmic tumor cell density along major white matter tracts.

Figure 1 .

 1 Figure 1. Simulated tumor cell density for parameters λ w = 4.2 mm and λ w /λ g = 10.The color bar refers to the base-10 logarithmic cell density, which is normalized to one in the center of the enhancing core.

  Figure 2. (a) T1 post gadolinium image of a GBM; (b) Isolines of the cell density (70%, 10%, 2%) overlaid on T2-FLAIR for λ w = 2.1 mm and λ w /λ g = 10.

Figure 3 .

 3 Figure 3. Dose distributions in brain tissue that minimize integral tumor cell survival for λ w = 4.2 mm (a,c) and λ w = 2.1 mm (b,d). For (c) and (d), the maximum dose is constrained to 60 Gy; for (a) and (b), no maximum dose constraint is imposed. ᾱ = 0.35 is assumed in all plots.

w o w u d grad d low Target 10

 10

  presciption dose in figure3c

  presciption dose in figure3a

Figure 4 .

 4 Figure 4. IMRT dose distributions in brain tissue for different prescription doses: (a) homogeneous prescription of 60 Gy; (b) 60/46 Gy prescribed to the boost/CTV; (c) for the prescription dose in figure 3c; (d) dose difference of (a) and (c); (e) for the prescription dose in figure 3a; (f) dose difference of (a) and (e).

  Figure A1. Illustration of dose redistribution for a patient with a multifocal GBM located in the left temporal and frontal lobe. (a) ideal dose prescription for α = 0.35 1/Gy, λ w = 3 mm and a maximum dose constraint of 70 Gy; (b) IMRT plan based on the dose prescription in (a); (c) IMRT treatment plan for a homogeneous 60 Gy dose prescription to the CTV; (d) dose difference of (b) and (c).

Table 1 .

 1 Objective function parameters used for IMRT optimization

Appendix A. Multifocal temporal/frontal lobe case

In this section, we demonstrate the concept of dose redistribution for an additional patient showing a multifocal GBM located in the left temporal and frontal lobe. The use of the tumor growth model for target delineation, and the differences compared to manual target delineation, are described in the accompanying paper ([10, Appendix A.1]). Here, we show an IMRT plan comparison using the same planning parameters as in section 4.1.

Figure A1b shows the ideal dose prescription for the parameters λ w /λ g = 10, ᾱ = 0.35 1/Gy, λ w = 3 mm, and a maximum dose constraint of 70 Gy. Figure A1a shows the IMRT plan for a uniform dose prescription of 60 Gy to the CTV; figure A1c shows the IMRT plan for the prescription dose in (b); and A1d displays the difference of the two dose distributions. The parameters for radiosensitivity and infiltration length yield a dose fall-off rate of approximately 10 Gy per cm. In figure A1, this becomes most apparent in the white matter regions in the contralateral frontal lobe and ipsilateral parietal lobe. Here, the dose falls from 70 Gy to approximately 50 Gy at the edge of the


Yogendra P Chaubey
email: chaubey@alcor.concordia.ca

Christophe Chesneau
email: chesneau@math.unicaen.fr

Hassan Doosti
email: h.doosti@tmu.ac.ir-hassan.doosti@unimelb.edu.au

Adaptive wavelet estimation of a density from mixtures under multiplicative censoring

Keywords: AMS 2000 subject classifications: Primary 62G07, 62G20; secondary 60K35 Density estimation, mixture, Besov balls, wavelets, hard thresholding

In this paper, a mixture model under multiplicative censoring is considered. We investigate the estimation of a component of the mixture (a density) from the observations. A new adaptive estimator based on wavelets and a hard thresholding rule is constructed for this problem. Under mild assumptions on the model, we study its asymptotic properties by determining an upper bound of the mean integrated squared error over a wide range of Besov balls. We prove that the obtained upper bound is sharp.

Introduction

First of all, let us present the standard density estimation problem under multiplicative censoring as described in [START_REF] Vardi | Multiplicative censoring, renewal processes, deconvolution and decreasing density: Nonparametric estimation[END_REF]. We observe n independent random variables Y 1 , . . . , Y n where

Y i = U i X i , i ∈ {1, . . . , n}, (1.1)
U 1 , . . . , U n are n unobserved i.i.d. random variables having the uniform distribution on [0, 1] and X 1 , . . . , X n are n unobserved i.i.d. random variables. We suppose that X 1 , . . . , X n and U 1 , . . . , U n are independent. We aim to estimate f, the unknown density of X 1 , from Y 1 , . . . , Y n .

It is shown in [START_REF] Vardi | Multiplicative censoring, renewal processes, deconvolution and decreasing density: Nonparametric estimation[END_REF] [see also [START_REF] Vardi | Large sample study of empirical distributions in a random multiplicative censoring model[END_REF]] how this model unifies several well-studied statistical problems, including non-para-metric inference for renewal processes, certain non-parametric deconvolution problems and estimation of decreasing densities. Among the methods used to solve this problem, there is the method of singular value decomposition (or orthonormal series based expansion) developed in [START_REF] Andersen | Multiplicative censoring: density estimation by a series expansion approach[END_REF], the method using the wavelet technique used in [START_REF] Abbaszadeh | Nonparametric estimation of a density under bias and multiplicative censoring via wavelet methods[END_REF] and the kernel method developed in [START_REF] Asgharian | Large-sample study of the kernel density estimators under multiplicative censoring[END_REF].

In this paper, we consider a more general model incorporating mixtures. Here, we observe n independent random variables Y 1 , . . . , Y n according to (1.1) where U 1 , . . . , U n are unobserved i.i.d. random variables having the uniform distribution on [0, 1] and X 1 , . . . , X n are unobserved independent random variables but not (necessarily) identically distributed: for any i ∈ {1, . . . , n}, X i depends on an unobserved random indicator I i taking its values in {1, . . . , m}. Applying the Bayes theorem, the density of X i is the following finite mixture:

h i (x) = m d=1 w d (i)f d (x), (1.2)
where w d (i) = P(I i = d) and f d is the conditional density of X i given {I i = d}.

We suppose that all these densities are unknown, the weight w d (i) is known, and X 1 , . . . , X n and U 1 , . . . , U n are independent. For a fixed ν ∈ {1, . . . , m}, we aim to estimate f ν i.e. the ν-th component of the mixture, from Y 1 , . . . , Y n . Note that, when m = 1, w 1 (1) = . . . = w 1 (n) = 1 and f ν = f 1 , we arrive at the density estimation problem under multiplicative censoring described in the first paragraph. The problem of estimating f ν from X 1 , . . . , X n (not Y 1 , . . . , Y n ) has been considered in several papers [see e.g. [START_REF] Maiboroda | Estimators of components of a mixture with varying concentrations[END_REF], [START_REF] Hall | Nonparametric estimation of component distributions in a multivariate mixture[END_REF], [START_REF] Pokhyl'ko | Wavelet estimators of a density constructed from observations of a mixture[END_REF] and Prakasa [START_REF] Rao | Wavelet linear estimation for derivatives of a density from observations of mixtures with varying mixing proportions[END_REF]]. However, to the best of our knowledge, the estimation of f ν from Y 1 , . . . , Y n is a new statistical problem that has potential applications in signal processing, biology, industry and telecommunications. For instance, consider the example of a sensor network (i.e. a collection of spatially distributed autonomous sensors intended to measure physical phenomena at diverse locations). Suppose that it contains m nodes, each sensor records only one physical value of the phenomena and only one sensor information is collected at a time. Now, given n multiplicative noisy versions Y 1 , . . . , Y n of the physical phenomena from non-necessarily identified sensors (i.e. unobserved I 1 , . . . , I n ), the goal is to estimate the density function of the physical phenomena at any sensor ν ∈ {1, . . . , m}. An overview of sensor networks can be found in [START_REF] Akyildiz | Wireless sensor networks: a survey[END_REF]. The reader may be referred to [START_REF] Chabert | Optimal wavelet for abrupt change detection in multiplicative noise[END_REF] and references therein for applications in multiplicative noise set-up in signal processing.

Among the numerous existing methods (such as Kernel, Spline, Wavelets and others) [see e.g. Prakasa [START_REF] Rao | Nonparametric functional estimation[END_REF][START_REF] Rao | Nonparametric functional estimation: an overview[END_REF] and [START_REF] Tsybakov | Introduction à l'estimation non-paramétrique[END_REF]], we investigate the estimation of f ν via wavelet methods. They are attractive for nonparametric function estimation because of their spatial adaptivity, computational efficiency and asymptotic optimality properties. They can achieve near optimal convergence rates over a wide range of function classes (such as Besov balls) and enjoy excellent mean integrated squared error (MISE) properties when used to esti- [START_REF] Antoniadis | Wavelets in statistics: a review (with discussion)[END_REF] and [START_REF] Härdle | Wavelet, Approximation and Statistical Applications[END_REF]].

In this paper, we develop an adaptive wavelet estimator for f ν based on the hard thresholding rule [introduced in Donoho andJohnstone (1994, 1995)].

In its construction, we use the new methodology developed in [START_REF] Abbaszadeh | Nonparametric estimation of a density under bias and multiplicative censoring via wavelet methods[END_REF] for the density estimation under bias and multiplicative censoring, the "observations thresholding technique" introduced by [START_REF] Delyon | On minimax wavelet estimators[END_REF], and some tools on mixtures described in [START_REF] Maiboroda | Estimators of components of a mixture with varying concentrations[END_REF]. We evaluate the performance of the proposed density estimator by considering the MISE and assuming that f ν belongs to a Besov ball. Under mild assumptions on the weights of the mixtures, we prove that our estimator attains a sharp rate of convergence: it is the one attained by the "best" non-adaptive linear wavelet estimator up to a logarithmic term, where "best" is used in the sense of minimum MISE.

The outline of the paper is as follows. We provide detailed assumptions on the mixture model considered here in Section 2. Section 3 describes the wavelet basis on [0, 1] and the Besov balls. The estimators are presented in Section 4 and the main results are provided in Section 5. The proofs are gathered in Section 6.

Assumptions

Additional assumptions on our model are presented below.

• Without loss of generality, we assume that the support of f d is [0, 1] and f d ∈ L 2 ([0, 1]) for any d ∈ {1, . . . , m}, where, for any p ≥ 1,

L p ([0, 1]) = u : [0, 1] → R; u p = 1 0 |u(x)| p dx 1/p < ∞ .
Note that, for any i ∈ {1, . . . , n}, the density of Y i is

g i (x) = 1 x h i (y) y dy, x ∈ [0, 1],
where h i is as given in (1.2). We suppose that there exists a known constant

C * > 0 such that max i∈{1,...,n} sup x∈[0,1] g i (x) ≤ C * . (2.1)
Such a boundedness assumption is standard in the density estimation problems from mixtures [see e.g. [START_REF] Maiboroda | Estimators of components of a mixture with varying concentrations[END_REF][START_REF] Pokhyl'ko | Wavelet estimators of a density constructed from observations of a mixture[END_REF] and Prakasa [START_REF] Rao | Wavelet linear estimation for derivatives of a density from observations of mixtures with varying mixing proportions[END_REF]].

• We now formulate some assumptions on the weights of the mixture. Consider the matrix We assume that det(Γ n ) > 0.

Γ n = 1 n n i=1 w k (i)w (i) (k, )∈{1,...,
Then the solution of the following quadratic objective with linear constraints,

min b∈R n 1 n n i=1 b 2 i such that 1 n n i=1 b i w d (i) = δ ν,d , for d ∈ {1, • • • , m}, (2.3)
where δ ν,d denotes the Kronecker delta, is unique and is denoted by (a ν (1), . . . , a ν (n)). Specifically, we have

a ν (i) = 1 det(Γ n ) m k=1 (-1) k+ν ζ n ν,k w k (i), (2.4)
where ζ n ν,k denotes the determinant of the minor (ν, k) of the matrix Γ n . Details concerning the above results can be found in [START_REF] Maiboroda | Estimators of components of a mixture with varying concentrations[END_REF] and [START_REF] Pokhyl'ko | Wavelet estimators of a density constructed from observations of a mixture[END_REF].

Remark 2.1. The definition of (a ν (1), . . . , a ν (n)) (2.3) allows us to construct unbiased wavelet coefficients estimators well adapted to the structure of our mixture model (see Proposition 4.1), obtain sharp rate of convergence (the quantity (1/n) n i=1 a 2 ν (i) will play a central role); see Theorems 5.1 and 5.2.

We set

ρ n = 1 n n i=1 a 2 ν (i) (2.5)
and, for technical reasons, we suppose that lim n→∞ n/ρ n = ∞.

Wavelets and Besov balls

Wavelets

For the purposes of this paper, we use the compactly supported wavelet basis on [0, 1] briefly described below.

Let N ≥ 10 and φ and ψ be the initial wavelet functions of the Daubechies wavelets db2N . In particular, these functions are compactly supported with supp(φ) ⊆ [1 -N, N ] and supp(ψ) ⊆ [1 -N, N ] (for the sake of simplicity), and belong to C 1 .

We define the functions φ j,k and ψ j,k by dilation and translation of φ and ψ respectively, and normalization:

φ j,k (x) = 2 j/2 φ(2 j x -k), ψ j,k (x) = 2 j/2 ψ(2 j x -k).
Then, with an appropriate treatment at the boundaries, there exists an integer τ such that, for any integer ≥ τ , is an orthonormal basis of L 2 ([0, 1]).

B = {φ ,k , k ∈ {0, . . . , 2 -1}; ψ j,k ; j ∈ N -{0, . . . , -1}, k ∈ {0, . . . ,
For any integer ≥ τ , a function u ∈ L 2 ([0, 1]) can be expanded on B as

u(x) = 2 -1 k=0 α ,k φ ,k (x) + ∞ j= 2 j -1 k=0 β j,k ψ j,k (x), x ∈ [0, 1],
where α j,k and β j,k are the wavelet coefficients of u defined by

α j,k = 1 0 u(x)φ j,k (x)dx, β j,k = 1 0 u(x)ψ j,k (x)dx. (3.1)
We refer to [START_REF] Cohen | Wavelets on the interval and fast wavelet transforms[END_REF] and [START_REF] Mallat | A wavelet tour of signal processing[END_REF] for more details.

Besov balls

As is traditional in the wavelet estimation literature, we will investigate the performances of our estimators by assuming that the unknown function f ν belongs to a Besov ball defined below. Let M > 0, s > 0, p ≥ 1 and r ≥ 1. Set, for every measurable function u on [0, 1] and ≥ 0,

∆ (u)(x) = u(x + ) -u(x), ∆ 2 (u)(x) = ∆ (∆ (u))(x) and identically, ∆ N (u)(x) = ∆ N -1 (∆ (u))(x) for any positive integer N . Let ρ N (t, u, p) = sup ∈[-t,t] ||∆ N (u)|| p .
Then, for s ∈ [0, N ), we define the Besov ball B s p,r (M ) of radius M > 0 by

B s p,r (M ) = u ∈ L p ([0, 1]); 1 0 ρ N (t, u, p) t s r dt t 1/r ≤ M , with the usual modifications if p = ∞ or r = ∞.
The parameter s is related to the smoothness of the function u, whereas p and r are norm parameters.

We have the following equivalence (see [START_REF] Härdle | Wavelet, Approximation and Statistical Applications[END_REF], Corollary 9.1): u ∈ B s p,r (M ) with s ∈ (0, N ) if and only if there exists a constant M * > 0 (depending on M ) such that the associated wavelet coefficients given by (3.1) satisfy

2 τ (1/2-1/p) 2 τ -1 k=0 |α τ,k | p 1/p +    ∞ j=τ   2 j(s+1/2-1/p)   2 j -1 k=0 |β j,k | p   1/p    r    1/r ≤ M * ,
with the usual modifications if p = ∞ or r = ∞. For a particular choice of s, p and r, B s p,r (M ) contains the Hölder and Sobolev balls [see [START_REF] Meyer | Ondelettes et Opérateurs[END_REF]].

imsart-generic ver. 2009/12/15 file: bj4_yogen_rev.tex date: December 12, 2013 4. Estimators

Wavelet coefficients estimators

The first step in estimating f ν consists of expanding f ν on B and estimating its unknown wavelet coefficients. Here we adapt the methodology of [START_REF] Abbaszadeh | Nonparametric estimation of a density under bias and multiplicative censoring via wavelet methods[END_REF]. We consider the operator T :

C 1 ([0, 1]) → C 0 ([0, 1]) defined by T (u)(x) = (xu(x)) = u(x) + xu (x), u ∈ C 1 ([0, 1]), x ∈ [0, 1]. (4.1)
For any integer j ≥ τ and any k ∈ {0, . . . , 2 j -1}, we estimate the wavelet coefficients of f ν , i.e.

α j,k = 1 0 f ν (x)φ j,k (x)dx and β j,k = 1 0 f ν (x)ψ j,k (x)dx, by respectively α j,k = 1 n n i=1 a ν (i)T (φ j,k )(Y i ) (4.2)
and

β j,k = 1 n n i=1 a ν (i)T (ψ j,k )(Y i )1 {|aν (i)T (ψ j,k )(Yi)|≤ηj } , (4.3)
where Y i is defined by (1.1), a ν (i) by (2.4), T by (4.1), 1 A is the indicator function on A for any random event A,

η j = θ2 j nρ n ln(n/ρ n ) , θ = 2C * 1 + N 1-N (ψ (y)) 2 dy , (4.4)
C * is the constant in (2.1) and ρ n is defined by (2.5).

The wavelet coefficient estimators α j,k and β j,k are motivated by the properties of the constants a ν (i) defined in (2.3) and (2.4). We can write

f ν (x) = 1 n n i=1 a ν (i)h i (x)
and use the fact that (see the details in the proof of Proposition 4.1)

1 0 h i (x)φ j,k (x)dx = E(T (φ j,k )(Y i )).
Thus α j,k provides an unbiased estimator of α j,k where as β j,k is based on the "observations thresholding technique" that had been firstly introduced in [START_REF] Delyon | On minimax wavelet estimators[END_REF] to estimate the regression function in a regression imsart-generic ver. 2009/12/15 file: bj4_yogen_rev.tex date: December 12, 2013 model with general errors. In our study, the role of its selection is to provide a good estimator of β j,k under mild assumptions on the distributions of I 1 , . . . , I n . To be more specific, we have no boundedness assumptions on the eigenvalues of the matrix Γ n (as in [START_REF] Autin | Test on the components of mixture densities[END_REF]) and, more generally, on a ν (i). This point is developed in Remarks 5.3 and 5.4. Our statistical approach has been used for other mixture models in Chesneau (2011a,b) and [START_REF] Chesneau | Adaptive wavelet estimation of a function in an indirect regression model[END_REF].

Proposition 4.1. For any integer j ≥ τ and any k ∈ {0, . . . , 2 j -1}, let α j,k and β j,k be the wavelet coefficients of f ν as given in (3.1). Then

• α j,k defined by (4.2) is an unbiased estimator of α j,k , • we have E 1 n n i=1 a ν (i)T (ψ j,k )(Y i ) = β j,k .
Remark 4.1. In addition to the unbiased properties described in Proposition 4.1, α j,k and β j,k satisfy powerful moments and concentration inequalities, i.e., Moments inequality: there exists a constant C > 0 such that

E ( α j,k -α j,k ) 2 ≤ C2 2j ρ n n and E β j,k -β j,k 4 ≤ C2 4j (ρ n ln(n/ρ n )) 2 n 2 .
Concentration inequality: for a large enough κ, we have

P | β j,k -β j,k | ≥ (κ/2)θ2 j ρ n ln(n/ρ n ) n ≤ 2 ρ n n 2 ,
where θ is defined by (4.4).

The technical details can be found in the proof of Proposition 6.2.

The next subsections are devoted to the description of our main estimators for f ν .

Wavelet linear estimator

We define the linear estimator f l ν by

f l ν (x) = 2 j 0 -1 k=0 α j0,k φ j0,k (x), x ∈ [0, 1], (4.5)
where α j,k is defined by (4.2) and j 0 is an integer which will be chosen later (see Theorem 5.1). In the case m = 1, w 1 (1) = . . . = w 1 (n) = 1 and f ν = f 1 , f l ν becomes the linear wavelet estimator developed in [START_REF] Abbaszadeh | Nonparametric estimation of a density under bias and multiplicative censoring via wavelet methods[END_REF].

For a survey on linear wavelet estimators for various density models (including the density model under multiplicative censoring), we refer to [START_REF] Chaubey | Wavelet Linear Density Estimation: A review[END_REF].

Wavelet hard thresholding estimator

We define the hard thresholding estimator f h ν by

f h ν (x) = 2 τ -1 k=0 α τ,k φ τ,k (x) + j1 j=τ 2 j -1 k=0 β j,k 1 {| β j,k |≥κλj } ψ j,k (x), x ∈ [0, 1], (4.6)
where α j,k is defined by (4.2), β j,k by (4.3), τ is the primary resolution level of the wavelet basis B, j 1 is the integer satisfying

1 2 n ρ n 1/3 < 2 j1 ≤ n ρ n 1/3 , κ ≥ 2 + 8/3 + 2 4 + 16/9, λ j = θ2 j ρ n ln(n/ρ n ) n .
and θ is defined by (4.4).

The definitions of the thresholds η j and λ j are based on theoretical considerations; they are chosen in such a way that the MISE of f h ν is as small as possible, while being independent of f ν (and its degree of smoothness). This independence makes f h ν adaptive. Technical details can be found in the proof of Theorem 6.1.

The idea of the considered thresholding is to set to zero all wavelet coefficients that are less than a particular threshold and to estimate the others. The resulting estimators often possess optimal (or near-optimal) rates of convergence in the minimax sense. Moreover, they are typically implemented via fast algorithms which makes them appealing in practice. For the constructions of hard thresholding wavelet estimators in the standard density model, see e.g. [START_REF] Donoho | Density estimation by wavelet thresholding[END_REF], [START_REF] Delyon | On minimax wavelet estimators[END_REF] and [START_REF] Härdle | Wavelet, Approximation and Statistical Applications[END_REF], and for the standard density mixture model, see [START_REF] Pokhyl'ko | Wavelet estimators of a density constructed from observations of a mixture[END_REF]. For the multiplicative censoring density model from biased data, see [START_REF] Abbaszadeh | Nonparametric estimation of a density under bias and multiplicative censoring via wavelet methods[END_REF]. To the best of our knowledge, f h ν is the first adaptive wavelet estimator proposed for the density estimation from mixtures under multiplicative censoring.

imsart-generic ver. 2009/12/15 file: bj4_yogen_rev.tex date: December 12, 2013

Results

Theorem 5.1 below explores the asymptotic performance of f l ν under the MISE over Besov balls B s p,r (M ). We distinguish the "homogeneous zone" corresponding to p ≥ 2 and the "inhomogeneous zone" corresponding to p ∈ [1, 2), following the classification of [START_REF] Härdle | Wavelet, Approximation and Statistical Applications[END_REF].

Theorem 5.1. Consider (1.1) under the assumptions of Section 2.

(i) Suppose that f ν ∈ B s p,r (M ) with s ∈ (0, N ), p ≥ 2 and r ≥ 1. Let f l ν be given by (4.5) with j 0 satisfying

1 2 n ρ n 1/(2s+3) < 2 j0 ≤ n ρ n 1/(2s+3)
.

(5.1)

Then there exists a constant C > 0 such that

E f l ν -f ν 2 2 ≤ C ρ n n 2s/(2s+3) . (ii) Suppose that f ν ∈ B s p,r (M ) with s ∈ (1/p, N ), p ∈ [1, 2) and r ≥ 1. Let f l ν
be given by (4.5) with j 0 satisfying

1 2 n ρ n 1/(2s * +3) < 2 j0 ≤ n ρ n 1/(2s * +3) , (5.2)
where s * = s + 1/2 -1/p. Then there exists a constant C > 0 such that

E f l ν -f ν 2 2 ≤ C ρ n n 2s * /(2s * +3)
.

The proof of Theorem 5.1 (see §6) uses moment inequalities on (4.2) and (4.3), and a suitable decomposition of the MISE.

Let us mention that, in each zone of the Besov balls, the integer j 0 is chosen to minimize the MISE of f l ν in a nonadaptive way. Remark 5.1. The variability of the performance of linear estimators according to the homogeneous nature of the zones of the Besov balls is a well known phenomenon. In particular, in the context of the standard nonparametric smoothing (such as for density and regression), the rates are often non-optimal (in the minimax sense) over the inhomogeneous zone. We refer to [START_REF] Härdle | Wavelet, Approximation and Statistical Applications[END_REF] Theorem 5.2. Consider (1.1) under the assumptions of Section 2. Let f h ν be given by (4.6). Suppose that f ν ∈ B s p,r (M ) with r ≥ 1, {p ≥ 2 and s ∈ (0, N )} or {p ∈ [1, 2) and s ∈ (3/p, N )}. Then there exists a constant C > 0 such that

E f h ν -f ν 2 2 ≤ C ρ n ln(n/ρ n ) n 2s/(2s+3)
.

The proof of Theorem 5.2 uses a general result derived to Chesneau (2011b) (see the proof of his Theorem 2).

Theorem 5.2 shows that, in addition to being adaptive, f h ν attains a sharp rate of convergence. More specifically, adapting the terminology of Härdle et al.

(1998),

• for the homogeneous zone (p ≥ 2): f h ν attains the same rate of convergence to f l ν , i.e. v n = (ρ n /n) 2s/(2s+3) , up to the logarithmic term (ln(n/ρ n )) 2s/(2s+3) .

• for the inhomogeneous zone (p ∈ [1, 2)): f h ν attains a better rate of convergence to f l ν . Remark 5.2. Let us mention that, for the standard nonparametric models, other types of adaptive procedures are able to remove the extra logarithmic term in the rate of convergence, for example,

• the Lepski procedure [see e.g. [START_REF] Lepski | On a problem of adaptive estimation in Gaussian white noise[END_REF][START_REF] Lepski | Asymptotically minimax adaptive estimation I : upper bounds. Optimally adaptive estimates[END_REF][START_REF] Lepski | Asymptotically minimax adaptive estimation II : statistical models without optimal adaptation[END_REF] and [START_REF] Lepski | Local adaptation to inhomogeneous smoothness : resolution level[END_REF]], and • the BlockJS procedure [see e.g. [START_REF] Cai | Adaptive wavelet estimation: a block thresholding and oracle inequality approach[END_REF][START_REF] Cai | On block thresholding in wavelet regression: adaptivity, block size and threshold level[END_REF]].

Although the adaptation of these procedures to our statistical problem is possible, their asymptotic properties are not established; new important theoretical difficulties appear and it is not immediately clear how to solve them.

Remark 5.3. Suppose that the smallest eigenvalue of the matrix Γ n (given in (2.2)) is larger than or equal to K with K ∈ (0, 1). Then Lemma 6.1 of [START_REF] Autin | Test on the components of mixture densities[END_REF] gives

1 n m d=1 n i=1 a 2 d (i) ≤ m K .
Hence ρ n ≤ m/K and the rate of convergence in Theorem 5.2 becomes (ln n/n) 2s/(2s+3) . It corresponds to the optimal one (up to a logarithmic term) for the standard multiplicative density model [see [START_REF] Andersen | Multiplicative censoring: density estimation by a series expansion approach[END_REF]].

Remark 5.4. Suppose that there exists a constant C > 0 satisfying max d∈{1,...,m} max i∈{1,...,n}

|a d (i)| ≤ C.
Then, with j 1 such that for a suitable constant c * , the indicator function in (4.3) is equal to one and Theorem 5.2 can be proved with the rate of convergence (ln n/n) 2s/(2s+3) .

c * 2 n ln n 1/3 ≤ 2 j1 ≤ c * n ln n 1/3 imsart-

Conclusion and perspectives.

We construct a new adaptive wavelet estimator to estimate the function f ν from (1.2). Under mild assumptions, we prove that it attains a sharp rate of convergence for a wide class of functions. Possible perspectives are to

• determine the optimal lower bound of MISE for (1.1), perhaps using the Fano Lemma presented in [START_REF] Tsybakov | Introduction à l'estimation non-paramétrique[END_REF]. However, it is not immediately clear how to extend this to the mixture models. • improve the estimation of f ν by removing the extra logarithmic term (see Remark 5.2), • investigate the case where the distributions of I 1 , . . . , I n are unknown,

• consider other estimation problems as the adaptive estimation of the q-th derivatives of f ν (as in Prakasa [START_REF] Rao | Nonparametric estimation of the derivatives of a density by the method of wavelets[END_REF] for the standard density model and in Prakasa [START_REF] Rao | Wavelet linear estimation for derivatives of a density from observations of mixtures with varying mixing proportions[END_REF] for the standard mixture model).

All these aspects need further investigations that we leave for a future work.

Proofs

In this section, we consider (1.1) under the assumptions of Section 2. Moreover, C denotes any constant that does not depend on j, k and n. Its value may change from one term to another and may depend on φ or ψ.

Proof of Proposition 4.1

Proof of Proposition 4.1. Recall that the density of Y i can be expressed in terms of h i in the following manner:

g i (x) = 1 x h i (y) y dy, x ∈ [0, 1]. So h i (x) = -xg i (x),
x ∈ [0, 1] a.s.. It follows from integration by parts, g i (1) = 0, and φ j,k (0) and φ j,k (1) exist, that

1 0 h i (x)φ j,k (x)dx = - 1 0 xg i (x)φ j,k (x)dx = -[xg i (x)φ j,k (x)] 1 0 + 1 0 g i (x) (φ j,k (x) + x(φ j,k ) (x)) dx = 1 0 g i (x)T (φ j,k )(x)dx = E(T (φ j,k )(Y i )), (6.1)
where T denotes the operator in (4.1). Using (6.1), (1.2) and ( 2.3), we obtain

E( α j,k ) = 1 n n i=1 a ν (i)E(T (φ j,k )(Y i )) = 1 n n i=1 a ν (i) 1 0 h i (x)φ j,k (x)dx = m d=1 1 0 f d (x)φ j,k (x)dx 1 n n i=1 a ν (i)w d (i) = 1 0 f ν (x)φ j,k (x)dx = α j,k .
Similarly, taking ψ instead of φ, we prove the second point that completes the proof of Proposition 4.1.

6.2. Proofs of Theorems 5.1 and 5.2

Proof of Theorem 5.1.

(i) We expand the function f ν on B at the level j 0 given in (5.1) :

f ν (x) = 2 j 0 -1 k=0 α j0,k φ j0,k (x) + ∞ j=j0 2 j -1 k=0 β j,k ψ j,k (x),
where

α j0,k = 1 0 f ν (x)φ j0,k (x)dx and β j,k = 1 0 f ν (x)ψ j,k (x)dx. We can write E f l ν -f ν 2 2 = 2 j 0 -1 k=0 E ( α j0,k -α j0,k ) 2 + ∞ j=j0 2 j -1 k=0 β 2 j,k . (6.2)
Let us now bound the first term. Since, by Proposition 4.1, α j,k is an unbiased estimator of α j,k , we have

E ( α j0,k -α j0,k ) 2 = 1 n 2 n i=1 a 2 ν (i)V(T (φ j0,k )(Y i )) ≤ 1 n 2 n i=1 a 2 ν (i)E (T (φ j0,k )(Y i )) 2 .
By the definition of T in (4.1), for any i ∈ {1, . . . , n}, we have Since, by (2.1), max i∈{1,...,n} sup x∈[0,1] g i (x) ≤ C * , we have

E (T (φ j0,k )(Y i )) 2 ≤ 2(A + B), where A = E φ 2 j0,k (Y i ) , B = E Y 2 i ((φ j0,k ) (Y i ))
A = 1 0 φ 2 j0,k (x)g i (x)dx ≤ C * 1 0 φ 2 j0,k (x)dx = C * .
On the other hand, using again (2.1), (φ j0,k ) (x) = 2 3j0/2 φ (2 j0 x -k), making the change of variables y = 2 j0 x -k and owing to supp(φ) ⊆ [1 -N, N ], we obtain

B = 1 0 x 2 ((φ j0,k ) (x)) 2 g i (x)dx ≤ C * 1 0 x 2 ((φ j0,k ) (x)) 2 dx ≤ C * 1 0 ((φ j0,k ) (x)) 2 dx = C * 2 3j 1 0 (φ (2 j0 x -k)) 2 dx ≤ C * 2 2j0 N 1-N (φ (y)) 2 dy.
Hence (6.3) where

E (T (φ j0,k )(Y i )) 2 ≤ C o 2 2j0 ,
C o = 2C * 1 + N 1-N (φ (y)) 2 dy . Therefore E ( α j0,k -α j0,k ) 2 ≤ C2 2j0 1 n 1 n n i=1 a 2 ν (i) = C2 2j0 ρ n n .
For j 0 satisfying (5.1) it holds

2 j 0 -1 k=0 E ( α j0,k -α j0,k ) 2 ≤ C2 j0 2 2j0 ρ n n ≤ C ρ n n 2s/(2s+3)
. (6.4)

On the other hand, for p ≥ 2, we have

f ν ∈ B s p,r (M ) ⊆ B s 2,∞ (M ) [see Härdle et al. (1998), Corollary 9.2], which implies ∞ j=j0 2 j -1 k=0 β 2 j,k ≤ C2 -2j0s ≤ C ρ n n 2s/(2s+3)
. (6.5)

It follows from (6.2), (6.4) and (6.5) that

E f l ν -f ν 2 2 ≤ C ρ n n 2s/(2s+3)
.

The part (i) of Theorem 5.1 is proved. (ii) Let p ∈ [1, 2). The proof is similar to (i) with s * instead of s: we use the expansion (6.2) at the level j 0 given in (5.2), and the inclusion This completes the proof of Theorem 5.1.

B s p,r (M ) ⊆ B s * 2,∞ (
Proof of Theorem 5.2. First of all, we present a general result (see Theorem 6.1) derived from Theorem 2 of [START_REF] Chesneau | Adaptive wavelet estimator for a function and its derivatives in an indirect convolution model[END_REF]. Its proof is given in §6.2. In this theorem, we consider a more general form of the hard thresholding estimator f h ν denoted by fH for estimating an unknown function f ∈ L 2 ([0, 1]) from n independent random variables W 1 , . . . , W n defined on a probability space (Ω, A, P) (We use the notations of Section 3):

f H (x) = 2 τ -1 k=0 α τ,k φ τ,k (x) + j1 j=τ 2 j -1 k=0 β j,k 1 {| β j,k |≥κϑj } ψ j,k (x), (6.6)
where

α j,k = 1 υ n n i=1
q i (φ j,k , W i ), (6.7)

β j,k = 1 υ n n i=1 q i (ψ j,k , W i )1 {|qi(ψ j,k ,Wi)|≤ςj } , (6.8)
ς j = θ ψ 2 ωj υ n √
µ n ln µ n (6.9)

ϑ j = θ ψ 2 ωj ln µ n µ n , (6.10)
κ ≥ 2 + 8/3 + 2 4 + 16/9 and j 1 is the integer satisfying

1 2 µ 1/(2ω+1) n < 2 j1 ≤ µ 1/(2ω+1) n . (6.11)
Here, we suppose that there exist

• n functions q 1 , . . . , q n with q i : L 2 ([0, 1]) × R → R for any i ∈ {1, . . . , n},

• two sequences of real numbers (υ n ) n∈N and (µ n ) n∈N satisfying lim n→∞ υ n = ∞ and lim n→∞ µ n = ∞ such that, for γ ∈ {φ, ψ}, (A1). any integer j ≥ τ and any k ∈ {0, . . . , 2 j -1},

E 1 υ n n i=1 q i (γ j,k , W i ) = 1 0 f (x)γ j,k (x)dx.
(A2). there exist two constants, θ γ > 0 and ω ≥ 0, such that, for any integer j ≥ τ and any k ∈ {0, . . . , 2 j -1},

n i=1 E (q i (γ j,k , W i )) 2 ≤ θ 2 γ 2 2ωj υ 2 n µ n .
imsart-generic ver. 2009/12/15 file: bj4_yogen_rev.tex date: December 12, 2013

Remark 6.1. The estimators α j,k and β j,k defined in the above theorem coincide with those in (4.2) and (4.3), respectively by choosing v n = n. The particular form of the estimator chosen will be clear from the context. Now we are ready to state the theorem required for proving Theorem 5.2.

Theorem 6.1. Let fH be (6.6) under (A1) and (A2). Suppose that f ∈ B s p,r (M ) with r ≥ 1, {p ≥ 2 and s ∈ (0, N )} or {p ∈ [1, 2) and s ∈ ((2ω + 1)/p, N )}. Then there exists a constant C > 0 such that

E f H -f 2 2 ≤ C ln µ n µ n 2s/(2s+2ω+1)
.

Remark 6.2. Theorem 6.1 is still valid

• with j 1 such that c * 2 µ n ln µ n 1/(2ω+1) < 2 j1 ≤ c * µ n ln µ n 1/(2ω+1)
, c * > 0, instead of (6.11).

• without θ ψ in the definitions of (6.9) and (6.10), but with κ ≥ 2θ 2 ψ + 8/3 + 2 4θ 2 ψ + 16/9.

Let us now investigate the assumptions (A1) and (A2) of Theorem 6.1 for the model (1.1). Set, for any γ ∈ {φ, ψ}, q i (γ j,k , y) = a ν (i)T (γ j,k )(y).

By Proposition 4.1, for any γ ∈ {φ, ψ}, we have

E 1 n n i=1 q i (γ j,k , Y i ) = 1 0 f ν (x)γ j,k (x)dx.
Proceeding as in (6.3), we obtain

1 n 2 n i=1 E (q i (γ j,k , Y i )) 2 ≤ θ 2 γ 2 2j 1 n 2 n i=1 a 2 ν (i) = θ 2 γ 2 2j ρ n n , with θ γ = 2C * 1 + N 1-N
(γ (y)) 2 dy . (6.12)

Thus the assumptions (A1) and (A2) in Theorem 6.1 are satisfied with, for any i ∈ {1, . . . , n}, υ n = n, ω = 1, µ n = n/ρ n and θ γ is (6.12). This implies the following result: Let f ν ∈ B s p,r (M ) with r ≥ 1, {p ≥ 2 and s ∈ (0, N )} or {p ∈ [1, 2) and s ∈ (3/p, N )}, then the hard thresholding estimator f h ν (4.6) satisfies

W i = Y i = U i X i , q i (γ j,k , y) = a ν (i)T (γ j,k )(
E f h ν -f ν 2 2 ≤ C ρ n ln(n/ρ n ) n 2s/(2s+3)
.

This completes the proof of Theorem 5.2.

6.3. Proof of Theorem 6.1

The proof of the theorem is broken into two major parts. In the first part, we present two auxiliary results on the statistical properties of (6.7) and (6.8).

Then, following the idea of [START_REF] Härdle | Wavelet, Approximation and Statistical Applications[END_REF], Proof of Proposition 10.3, these results are used in the second part on a suitable decomposition of the MISE of f H .

Part I: Auxiliary results

Proposition 6.1. Consider β j,k as defined in (6.8). We have

| β j,k -β j,k | ≤ 1 υ n n i=1 U i + ϑ j ,
where

U i = q i (ψ j,k , W i )1 {|qi(ψ j,k ,Wi)|≤ςj } -E(q i (ψ j,k , W i )1 {|qi(ψ j,k ,Wi)|≤ςj } ). (6.13)
Proof of Proposition 6.1. By (A1), we have where U 1 , . . . , U n are given by (6.13).

β j,k = E 1 υ n n i=1 q i (ψ j,k , W i ) = 1 υ n n i=1 E(q i (ψ j,k , W i )1 {|qi(ψ j,k ,Wi)|≤ςj } ) + 1 υ n n i=1 E(q i (ψ j,k , W i )1 {|qi(ψ j,k ,Wi)|>ςj } ). Hence | β j,k -β j,k | ≤ 1 υ n n i=1 U i + 1 υ n n i=1 E |q i (ψ j,k , W i )|1
The Markov inequality and (A2) imply that

1 υ n n i=1 E |q i (ψ j,k , W i )|1 {|qi(ψ j,k ,Wi)|>ςj } ≤ 1 ς j υ n n i=1 E (q i (ψ j,k , W i )) 2 ≤ 1 ς j θ 2 ψ 2 2ωj υ n µ n
Using the definitions of ς j and ϑ j from (6.9) and (6.10) respectively, the right hand side of the above inequality becomes

1 θ ψ 2 ωj υn √ µn ln µn θ 2 ψ 2 2ωj υ n µ n = θ ψ 2 ωj ln µ n µ n = ϑ j
that completes the proof of Proposition 6.1.

Proposition 6.2. Consider α j,k and β j,k as defined in (6.7) and (6.8), respectively.

(i) There exists a constant C > 0 such that

E ( α j,k -α j,k ) 2 ≤ C2 2ωj 1 µ n .
(ii) There exists a constant C > 0 such that

E β j,k -β j,k 4 ≤ C2 4ωj ln µ n µ n 2 .
(iii) For any κ ≥ 2 + 8/3 + 2 4 + 16/9, we have

P | β j,k -β j,k | ≥ κϑ j /2 ≤ 2 1 µ 2 n . (6.14)
Proof of Proposition 6.2.

(i) Noting that α j,k is an unbiased estimator of α j,k , by (A1), the independence of W 1 , . . . , W n and (A2) yield (ii) Using Proposition 6.1 we have (6.15) where U 1 , . . . , U n are given by (6.13).

E ( α j,k -α j,k ) 2 = V ( α j,k ) = 1 υ 2 n n i=1 V(q i (φ j,k , W i )) ≤ 1 υ 2 n n i=1 E (q i (φ j,k , W i )) 2 ≤ θ 2 φ 2 2ωj 1 µ n . imsart-
E β j,k -β j,k 4 ≤ C   1 υ 4 n E   n i=1 U i 4   + ϑ 4 j   .
Note that E(U i ) = 0 and

E (U i ) 4 ≤ 16ς 2 j E (q i (ψ j,k , W i )) 2 .
Hence using the Rosenthal inequality (see the Appendix), (A2) and the definition of ς j from (6.9), we have

E   n i=1 U i 4   ≤ C   n i=1 E U 4 i + n i=1 E U 2 i 2   ≤ C   ς 2 j n i=1 E (q i (ψ j,k , W i )) 2 + n i=1 E (q i (ψ j,k , W i )) 2 2   ≤ C ς 2 j θ 2 ψ 2 2ωj υ 2 n µ n + θ 4 ψ 2 4ωj υ 4 n µ 2 n ≤ C 2 4ωj υ 4 n µ 2 n ln µ n + 2 4ωj υ 4 n µ 2 n ≤ C2 4ωj υ 4 n µ 2 n . (6.16)
It follows from (6.15), (6.16) and the definition of ϑ j from (6.10) that

E β j,k -β j,k 4 ≤ C2 4ωj ln µ n µ n 2 .
(iii) Proposition 6.1 yields (6.17) where U 1 , . . . , U n are given by (6.13) and κ * = κ/2 -1.

P | β j,k -β j,k | ≥ κϑ j /2 ≤ P n i=1 U i ≥ κ * υ n ϑ j ,
Observe that E(U 1 ) = 0, |U 1 | ≤ 2ς j and, by (A2),

n i=1 E U 2 i ≤ n i=1 E (q i (ψ j,k , W i )) 2 ≤ θ 2 ψ 2 2ωj υ 2 n µ n .
Hence, using the Bernstein inequality (see the Appendix) and the defini- tions of ϑ j and ς j from (6.9) (6.10), respectively, we have

P n i=1 U i ≥ κ * υ n ϑ j ≤ 2 exp - κ 2 * υ 2 n ϑ 2 j 2 n i=1 E (U 2 i ) + 2 3 κ * υ n ϑ j ς j ≤ 2 exp   - κ 2 * υ 2 n θ 2 ψ 2 2ωj ln µn µn 2 θ 2 ψ 2 2ωj υ 2 n µn + 2 3 κ * θ 2 ψ 2 2ωj υ 2 n µn   = 2µ - κ 2 * 2 ( 1+ 2 3 κ * ) n .
Taking κ such that κ 2 * / (2 (1 + (2/3)κ * )) ≥ 2, i.e. κ ≥ 2+8/3+2 4 + 16/9, we obtain

P | β j,k -β j,k | ≥ κϑ j /2 ≤ 2 1 µ 2 n .
This completes the proof of Proposition 6.2.

Part II: Proof of the theorem

We decompose f on the wavelet basis B as

f (x) = 2 τ -1 k=0 α τ,k φ τ,k (x) + ∞ j=τ 2 j -1 k=0 β j,k ψ j,k (x),
where α τ,k = 1 0 f (x)φ τ,k (x)dx and β j,k = 1 0 f (x)ψ j,k (x)dx. Using the orthonormality of the wavelet basis B, we can decompose the MISE of f H in three terms as: (6.18) where

E f H -f 2 2 = E + F + G,
E = 2 τ -1 k=0 E ( α τ,k -α τ,k ) 2 , F = j1 j=τ 2 j -1 k=0 E β j,k 1 {| β j,k |≥κϑj } -β j,k 2 and G = ∞ j=j1+1 2 j -1 k=0 β 2 j,k .
We now exhibit the bounds of each of these terms.

Bound for E : Bound for F : For the sake of simplicity, set Bj,k = β j,k -β j,k . The term F can be decomposed as (6.20) where Owing to the previous inequalities, the Cauchy-Schwarz inequality and Proposition 6.2, we obtain Let j * be the integer (6.22). Then F 4 can be bound as Observe that, using the definition of ϑ j , we have

F = F 1 + F 2 + F 3 + F 4 ,
F 1 = j1 j=τ 2 j -1 k=0 E B2 j,k 1 {| β j,
F 1 + F 3 ≤ C j1 j=τ 2 j -1 k=0 E B2 j,k 1 {| Bj,k |>κϑj /2}

  M ) (thanks to p ∈ [1, 2) and the definition of s * , i.e. s * = s + 1/2 -1/p, [see Härdle et al. (1998), Corollary 9.2]), we obtain the desired result. imsart-generic ver. 2009/12/15 file: bj4_yogen_rev.tex date: December 12, 2013

  k |≥κϑj , |β j,k |<κϑj /2} , 1 {| β j,k |≥κϑj , |β j,k |≥κϑj /2} , 1 {| β j,k |<κϑj , |β j,k |≥2κϑj } and 1 {| β j,k |<κϑj , |β j,k |<2κϑj } . Bound for F 1 + F 3 : Observe that 1 {| β j,k |<κϑj , |β j,k |≥2κϑj } ≤ 1 {| Bj,k |>κϑj /2} , 1 {| β j,k |≥κϑj , |β j,k |<κϑj /2} ≤ 1 {| Bj,k |>κϑj /2} and 1 {| β j,k |<κϑj , |β j,k |≥2κϑj } ≤ 1 {|βj,k|≤2| Bj,k |} .

F

  1 {|β j,k |<2κϑj } .

FF.

  r ≥ 1 and p ≥ 2, since f ∈ B s p,r (M ) ⊆ B s 2,∞ (M ) (see Härdle et al. (1998), Corollary 9.2), we have For r ≥ 1, p ∈ [1, 2) and s > (2ω + 1)/p, owing to the Markov inequality, the definition of ϑ j (from (6.10)), f ∈ B s p,r (M ) and (2s + 2ω + 1)(2 -p)/2 + (s + imsart-generic ver. 2009/12/15 file: bj4_yogen_rev.tex date: December 12, 2013

  generic ver. 2009/12/15 file: bj4_yogen_rev.tex date: December 12, 2013

  {|qi(ψ j,k ,Wi)|>ςj } ,

	imsart-generic ver. 2009/12/15 file: bj4_yogen_rev.tex date: December 12, 2013

  So, for r ≥ 1, {p ≥ 2 and s > 0} or {p ∈ [1, 2) and s > (2ω + 1)/p}, we have

	imply that
	F 2,2 ≤ C	ln µ n µ n	j1 j=j * +1	2 2ωj 1 j ϑ p	2 j -1 k=0	|β j,k | p
	≤ C	ln µ n µ n		(2-p)/2	∞ j=j * +1	2 jω(2-p) 2 -j(s+1/2-1/p)p
	≤ C	ln µ n µ n		(2-p)/2	2 -j * (s+1/2-1/p+ω-2ω/p)p ≤ C	ln µ n µ n	2s/(2s+2ω+1)	.
						F 2 ≤ C	ln µ n µ n	2s/(2s+2ω+1)	.	(6.23)
	Bound for F 4 :
	We have								j1	2 j -1
						F 4 ≤	β 2 j,k 1 {|β j,k |<2κϑj } .
									j=τ	k=0
	≤ C	j1	2 j -1		E B4 j,k	1/2	P | Bj,k | > κϑ j /2	1/2
			j=τ	k=0
	≤ C	ln µ n µ 2 n	j1 j=τ	2 j(1+2ω) ≤ C	ln µ n µ n	≤ C	ln µ n µ n	2s/(2s+2ω+1)	. (6.21)

imsart-generic ver. 2009/12/15 file: bj4_yogen_rev.tex date: December 12, 2013

imsart-generic ver. 2009/12/15 file: bj4_yogen_rev.tex date: December 12, 2013

Acknowledgements

The authors are thankful to the reviewers for their comments which have helped in improving the presentation.

* Supported in part by a Discovery Grant from NSERC of Canada. † Supported in part by ANR grant NatImages, ANR-08-EMER-009.

Bound for F 2 : Using the Cauchy-Schwarz inequality and Proposition 6.2, we obtain

Let j * be the integer defined by 1 2

.

(6.22)

Then F 2 can be bounded as

where

and

We have

.

For r ≥ 1 and p ≥ 2, the Markov inequality, definition of ϑ j (from (6.10)) and [START_REF] Härdle | Wavelet, Approximation and Statistical Applications[END_REF], Corollary 9.2) yield

.

For r ≥ 1, p ∈ [1, 2) and s > (2ω + 1)/p, the Markov inequality, the definition of 1/2 -1/p + ω -2ω/p)p = 2s, we obtain

.

Thus, for r ≥ 1, {p ≥ 2 and s > 0} or {p ∈ [1, 2) and s > (2ω + 1)/p}, we have

. (6.24)

Therefore it follows from (6.20), (6.21), (6.23) and (6.24) that

. (6.25)

Bound for G: For r ≥ 1 and p ≥ 2, we have f ∈ B s p,r (M ) ⊆ B s 2,∞ (M ). Since 2s/(2s + 2ω + 1) < 2s/(2ω + 1), we have

.

For r ≥ 1 and p ∈ [1, 2), we have [START_REF] Härdle | Wavelet, Approximation and Statistical Applications[END_REF], Corollary 9.2). Since s > (2ω+1)/p, we have (s+1/2-1/p)/(2ω+1) > s/(2s + 2ω + 1). Hence

.

Hence, for r ≥ 1, {p ≥ 2 and s > 0} or {p ∈ [1, 2) and s > (2ω + 1)/p}, we have

. (6.26)

Combining (6.18), (6.19), (6.25) and (6.26), we have, for any r ≥ 1, {p ≥ 2 and s ∈ (0, N )} or {p ∈ [1, 2) and s ∈ ((2ω + 1)/p, N )},

.

This completes the proof of Theorem 6.1.

See [START_REF] Rosenthal | On the subspaces of L p (p ≥ 2) spanned by sequences of independent random variables[END_REF] for details.

Lemma 6.2 (Bernstein's inequality). Let n ∈ N * and (A i ) i∈{1,...,n} be n zero mean independent random variables such that there exists a constant M > 0 satisfying, for any i ∈ {1, . . . , n}, |A i | ≤ M < ∞. Then, for any y > 0, we have

See [START_REF] Petrov | Limit Theorems of Probability Theory: Sequences of Independent Random Variables[END_REF] for details.