M.D. Markus Stenner
email: markus.stenner@uk-koeln.de

M.dent.D Basima Yosef

Ph.D Christian U Huebbers

M.D Simon Florian Preuss

M.D Hans Peter Dienes

Ernst-Jan M Speel

Ph.D Margarete Odenthal

M.D Jens Peter Klussmann

Ph.D E J M Speel

Nuclear translocation of β-catenin and decreased expression of E-cadherin in HPV-positive tonsillar cancer: An early event in

Keywords: HPV, metastasis, nuclear beta-catenin, p16, tonsillar cancer HPV, metastasis, nuclear β-catenin, p16, tonsillar cancer

Dear Editor, Thank you very much for the work and the help with our manuscript. We were pleased to read the reviewers' comments, which were very helpful to revise our paper.

We would be grateful if you now consider this article to be published in

Histopathology.

Here is our list of changes according to the reviewers' comments (all changes were made using the track changes mode in MS Word): "This is an interesting manuscript that examines the ßcatenin and Ecadherin expression in HPV positive tonsillar cancer." (Response): -"The assessment of the immunohistochemistry is subjective, but uses two independent observers with a consensus agreement on discrepant scoring. The manuscript would be enhanced by an indication of the inter-observer variation, using for example Kappa agreement scores."

(Response): We now introduced Cohen's kappa coefficient as a measure of inter-observer agreement for all evaluated immunohistochemical markers (p16, E-cadherin, β-catenin, vimentin). We furthermore note that all cases of discrepancies between the two experienced pathologists were seen by Hans Peter Dienes, the head of the Institute for Pathology at the University Hospital Cologne to come to a decision.

"The photomicrographs in Figure 1 show examples of the ßcatenin staining. The description of panel B is not accurate.

The staining shows strong nuclear, cytoplasmic and membrane staining, not exclusive nuclear staining. In addition, the background staining is high when compared to the photomicrographs illustrated in A, C and D. Photomicrograph B raises the possibility that there are interpretive errors in the data." (Response): We are very grateful to this reviewer, as we find a misleading labelling error concerning the different localisations of β-catenin staining. The cut-off point for the group with prominent nuclear expression of β-catenin was >50% of positive nuclei irrespective of membranous or cytoplasmic staining. We now addressed this in the text and accordingly changed the name of this group into "nuclear expression" wherever necessary. We corrected this error in the legend as well. Because of the high background signal, we chose another picture for prominent nuclear staining in primary tumors (Figure 1 B).

"Figure 3B states that there are 'scattered vimentin positive carcinoma cells within the tumour tissue' it is possible that these cells represent transiting immune cells." (Response): All specimens were independently classified by two experienced pathologists under a microscope with magnification x400 (in contrast to mag. x100 in the figure), and positive cells within tumor tissue were rated as carcinoma cells. According to their common training the pathologists distinguished between immune cells and cancer cells during their evaluation. However, sections were additionally examined by Hans Peter Dienes, the head of the Institute for Pathology at the University Hospital Cologne, who finally decided whether a cell was counted as a vimentin-positive cancer cell or not. Nevertheless we agree with your critique. Therefore, we now did an additional immunohisochemistry for all specimens with an antibody to the Leucocyte Common Antigen (LCA) which is commonly known to detect about 95% of human leucocytes. All specimens were again seen by the two independent pathologists, and we did not find any leucocytes within areas of tumor tissue. We attached at the end of this document representative sections of the LCA immunohistochemistry of the same patient seen in figure 3 B and addressed this aspect within "methods".

"The results section states that 'For further analyses, the patients were categorized due to their HPV status determined by the above mentioned HPVspecific PCR. Categorizing patients into HPV-and p16-positive ones and all others did not lead to any different results.' The data should be analysed using the molecular classification described by Smeets et al 2007 (Ref 18). This would allow the reader to clearly see the association between 'biologically relevant' HPV infection (p16+/HPV DNA+ cases) and E-cadherin/ ßcatenin expression." PCR as gold standard. In their series the sensitivity and specificity for active HPV infection of this approach was 100%. As proposed, we now re-classified our complete data into p16+/HPV-16 DNA+ cases and all others (named HPV-related and HPV-unrelated). The results section as well as the table and figures 2 and 4 have been rewritten. As stated before, this did not lead to any different results.

Associate Editor:

"This is a potentially interesting study providing novel data on HPVrelated carcinoma of the oropharynx. I support the views of the main reviewer, including the comments on the photomicrographs which do not clearly demonstrate the features in the legends. The authors should also consider the following points:" (Response): (see statement concerning changes above) "The analysis of the results is based on dichotomous grouping of negative and positive cases, however, the methods indicates that scoring was based on 0, +, ++, +++. Was any degree of +/++/+++ regarded as positive for statistical analysis? or, was another threshold applied? Are 5% p16 positive carcinomas often HPV positive?" (Response): After primary scoring of the immunohistochemical markers, a grouping was done. "0" was grouped as "negative" whereas "+ / ++ / +++" was grouped as "positive". Less than 5% positivity was graded as negative in case of p16, as described before (Preuss et al. 2008 Br J Cancer). We found only four cases with less than 5% p16-positivity and a positive PCR for HPV-16 DNA in primary tumors and only two cases in metastases.

"Is the use of categorical data the best way to look at differences in expression between primary tumours and metastases? It would be interesting the the categorical analysis were supported by a comparison based on continuous data, since it is not apparent how much difference there was between primary and metastatic tumours." (Response): We made use of categorical data but used 4 categories as this is commonly done in literature and it gives some information on quantification. We admit that continuous data would be interesting, however we were not able re-evaluate all sections for continuous data.

"Table 1 is unnecessarily detailed -there are no other references to observations on individual cases, so grouped data could be substituted." (Response): There truly are no other references to observations on individual cases. As grouped data for tumor grading, HPV status, age and gender is already mentioned in the text in "methods" and in (former) Table 2, we now left away (former) Table 1. For a better understanding, we added the information on patients' pT and pN status out of (former) Table 1 in the text in "methods" as well.

"Table 3 Aims: High-risk human papillomaviruses (HPV) seem to be an important risk factor for tonsillar cancer. However, only few data exist concerning molecular changes during metastasis. Methods and Results: We examined 48 primary tonsillar carcinoma samples (25 HPV-16 DNA positive, 23 HPV-16 DNA negative) and their respective lymph node metastases for their HPV status and for the expression of p16, E-cadherin, β-catenin and vimentin. A positive HPV-specific PCR correlated significantly with p16 overexpression in both the primary tumor and metastasis (p<0.0001 both). In HPV-unrelated patients the expression of E-cadherin was significantly less abundant in metastases than in primary tumors (p<0.001). The expression of nuclear β-catenin was contrariwise significantly higher in metastases than in primary tumors (p=0.016). Compared to HPV-unrelated patients, in HPVrelated ones the nuclear localisation of β-catenin expression was already apparent in primary tumors (p=0.030). The expression of vimentin significantly correlated with the grading of the primary tumor (p=0.021). Conclusions: Our data newly indicate that the downregulation of E-cadherin together with an upregulation of nuclear β-catenin might be crucial steps during tumor progression of tonsillar cancer being already present in primary tumors in HPV driven cancers but become apparent in HPVunrelated tumors later during the process of metastasizing.

INTRODUCTION

Head and neck squamous cell carcinoma (HNSCC) is the fifth most common malignancy worldwide 1 . Among these, the oropharyngeal squamous cell carcinoma (OSCC) is the most common entity and the sixth most common cancer in the United States, with approximately 8,000-10,000 deaths each year 2 . Tonsillar cancer is the most frequent OSCC. It is characterized by local invasiveness and a marked propensity for dissemination to cervical lymph nodes. Despite advances in surgery, radiation, and chemotherapy, the five-year survival rate has not improved significantly over the past several decades and it remains at about 50 to 55% [START_REF] Ragin | The Epidemiology and Risk Factors of Head and Neck Cancer: a Focus on Human Papillomavirus[END_REF] .

The development of tonsillar cancer is closely related to both intensity and duration of alcohol and tobacco consumption, which represent the most significant risk factors [START_REF] Gillison | Current topics in the epidemiology of oral cavity and oropharyngeal cancers[END_REF] .

Since a few years, there is growing evidence, that high-risk human papillomaviruses (HPV) are an important risk factor for OSCC as well [START_REF] Gillison | Distinct risk factor profiles for human papillomavirus type 16-positive and human papillomavirus type 16-negative head and neck cancers[END_REF][START_REF] Klussmann | Prevalence, distribution, and viral load of human papillomavirus 16 DNA in tonsillar carcinomas[END_REF] . HPV can infect the epithelium of the upper aerodigestive tract in general, but the tonsil appears uniquely susceptible to transformation by the virus [START_REF] Venuti | Presence of HPV in head and neck tumors: high prevalence in tonsillar localization[END_REF][START_REF] Klussmann | Human papillomavirus-positive tonsillar carcinomas: a different tumor entity?[END_REF] . Although the genome of the prototypic high-risk types HPV-16 and -18 codes for several different proteins, two viral proteins, E6 and E7 are classified as oncogenes because they contribute to tumor progression by inactivating two tumor suppressor genes. E6 degrades p53 and E7 affects the retinoblastoma (Rb) protein, respectively [START_REF] Galloway | The disruption of cell cycle checkpoints by papillomavirus oncoproteins contributes to anogenital neoplasia[END_REF][START_REF] Scheffner | The HPV-16 E6 and E6-AP complex functions as a ubiquitin-protein ligase in the ubiquitination of p53[END_REF] . A recent review of the prognostic significance of HPV in survival of head and neck cancer patients has shown that the majority of studies have reported better survival in those patients with HPV-positive HNSCC when compared to patients with HPV-negative tumors [START_REF] Ragin | Survival of squamous cell carcinoma of the head and neck in relation to human papillomavirus infection: review and meta-analysis[END_REF][12][START_REF] Schlecht | Prognostic value of human papillomavirus in the survival of head and neck cancer patients: an overview of the evidence[END_REF] , even though HPV-related tumors are less differentiated than HPV-unrelated ones [START_REF] Badaracco | Molecular analyses and prognostic relevance of HPV in head and neck tumors[END_REF][START_REF] Wilczynski | Detection of human papillomavirus DNA and oncoprotein overexpression are associated with distinct morphological patterns of tonsillar squamous cell carcinoma[END_REF] and frequently appear with lymph node metastases. In OSCC the tumor suppressor protein p16 is so consistently overexpressed in HPV-positive tumors as opposed to HPV-negative ones 12,[START_REF] Li | The expression of key cell cycle markers and presence of HPV in squamous cell carcinoma of the tonsil[END_REF]17 that it has been proposed as a surrogate marker for highrisk HPV infection. Therefore p16 immunostaining and subsequent HPV-PCR has been recommended to detect HPV-related primary HNSCC [START_REF] Smeets | A novel algorithm for reliable detection of human papillomavirus in paraffin embedded head and neck cancer specimen[END_REF] . Only some studies exist describing the presence of HPV DNA in neck metastases of HNSCC [START_REF] Hoffmann | Human papillomaviruses in lymph node neck metastases of head and neck cancers[END_REF][START_REF] Goldenberg | Cystic lymph node metastasis in patients with head and neck cancer: An HPV-associated phenomenon[END_REF] and a correlation of p16 expression and HPV status has so far only been shown once in metastases of the oropharynx [START_REF] Begum | Detection of human papillomavirus in cervical lymph nodes: a highly effective strategy for localizing site of tumor origin[END_REF] .

Lymph node metastasis seems to be the most important prognostic indicator for patients with OSCC [START_REF] Tanaka | pRb2/p130 protein expression is correlated with clinicopathologic findings in patients with oral squamous cell carcinoma[END_REF] . Some of the properties of metastatic cells have been linked to downregulation of the epithelial cadherin (E-cadherin) glycoprotein [START_REF] Wijnhoven | Genetic manipulation of E-cadherin expression by epithelial tumor cells reveals an invasion suppressor role[END_REF][24][25] , and in OSCC several studies revealed that abnormal expression of E-cadherin correlates with lymph node metastasis [START_REF] Kudo | Invasion and metastasis of oral cancer cells require methylation of E-cadherin and/or of membranous betacatenin[END_REF][START_REF] Ghosh | Loss of adhesion-regulated proteinase production is correlated with invasive activity in oral squamous cell carcinoma[END_REF] . Furthermore, E-cadherin was found to be a selective and strongly dominant prognostic factor in squamous cell carcinomas [START_REF] Bosch | E-cadherin is a selective and strongly dominant prognostic factor in squamous cell carcinoma: a comparison of Ecadherin with desmosomal components[END_REF] . E-cadherinmediated cell adhesion requires intracellular attachment to the actin cytoskeleton by interaction with catenins. Especially the β-catenin is not only an important protein linked to epithelial adhesion molecules but may also act as an oncogene, becoming one of the key downstream effectors in the Wnt and Akt signaling pathway. Dimers of β-catenin and one of the members of the TCF/LEF (T cell factor/lymphoid enhancer factor) family function as a transcription factor when located to the nucleus. Amongst other genes it transactivates cyclin D1. Some studies reported that reduced assembly of membranous β-catenin and nuclear translocation was found frequently in invasive areas and metastases of HNSCC [START_REF] Van Der Velden | Expression of Cytokeratin Subtypes and Vimentin in Squamous Cell Carcinoma of the Floor of the Mouth and the Mobile Tongue[END_REF] . Studies of epithelial malignancies have furthermore consistently demonstrated that loss of cell adhesion and the acquisition of mesenchymal features, a process termed epithelial-mesenchymal transition (EMT), precede their invasion and progression [START_REF] Huber | Molecular requirements for epithelial-mesenchymal transition during tumor progression[END_REF][START_REF] Thiery | Epithelial-mesenchymal transitions in tumor progression[END_REF] [START_REF] De Araujo | Vimentin in oral squamous cell carcinoma[END_REF] . In this context, a significant correlation between decreased E-cadherin and vimentin expression and aggressive tumor features including lymph node metastasis could recently be shown by immunophenotypic analysis of HNSCC [START_REF] Mandal | Epithelial to mesenchymal transition in head and neck squamous carcinoma: association of Src activation with E-21 cadherin down-regulation, vimentin expression, and aggressive tumor features[END_REF] . Underscoring the possible role of HPV for tumor progression in primary tumors, it remains unclear whether there is a link between HPV and tumor progression in metastases as well. In this context, it has recently been shown that HPV DNA could even be detected in histologically confirmed tumour-free neck lymph nodes of head and neck cancers [START_REF] Hoffmann | HPV16 DNA in histologically confirmed tumour-free neck lymph nodes of head and neck cancers[END_REF] .

Here we have studied the loss of cell adhesion and the acquisition of mesenchymal features in relation to the HPV status in primary tonsillar carcinomas and their respective metastases.

MATERIAL AND METHODS

Patient data

The study was conducted according to the Declaration of Helsinki on biomedical research involving human subjects. Written informed consent was obtained from each patient and the scientific protocol was approved by the Local Ethics Committee.

In a retrospective design, we selected from previous studies [START_REF] Klussmann | Genetic signatures of HPV-related and unrelated oropharyngeal carcinoma and their prognostic implications[END_REF]17 48 primary tumor samples (25 HPV-16 DNA positive and 23 HPV-16 DNA negative) obtained at the time of surgical resection as well as respective lymph node metastasis samples acquired during the neck dissection of these patients. A neck dissection was performed in patients with primary tumor stages T2-T4. Basis of selection was the availability of sufficient fresh frozen tumor tissue with at least 70% tumor cells, enough high quality tumor DNA and assessment of the presence or absence of HPV sequences. HPV DNA was detected using a highly sensitive group-specific nested PCR with degenerate primers A5/A10 and A6/A8 [START_REF] Wieland | Imiquimod treatment of anal intraepithelial neoplasia in HIV-positive men[END_REF] . HPV typing was performed as described previously [START_REF] Klussmann | Prevalence, distribution, and viral load of human papillomavirus 16 DNA in tonsillar carcinomas[END_REF] . All patients were treated at the Department of Otorhinolaryngology, Head and Neck Surgery, School of Medicine by JPK between 1997 and 2004. The treatment modality employed was definitive surgery with postoperative radiation or chemo radiation, and none of the samples examined were influenced by chemotherapeutic agents or preceding radiotherapy.

The patients' age ranged from 37 to 76 years (median 57 years) with a male to female ratio of 3 to 1. The classification of primary tumors and metastases as well as the histological grading was performed in a blinded manner following the World Health Organization 2005 criteria for squamous cell carcinomas of the oral mucosa.

The pathologic primary tumor and lymph node stages were as follows: pT1 13 ( 45.8%), pN2c 3 (6.3%), pN3 1 (2.1%). As tonsillar tumors tend to be non-keratinizing or basaloid, grading sometimes is difficult. Therefore, in 10 out of 48 cases, the grades 1-2 and 2-3 were introduced in which the assignment was unclear. All patients were cM0 as evaluated by means of a computed tomography scan of the thorax and an ultrasound examination of the abdomen.

Immunohistochemistry

Tissue specimens were fixed in 10% buffered formaldehyde and embedded in paraffin. Sections 4 µm thick were mounted on Superfrost Plus glass slides (Microm, Walldorf, Germany), deparaffinised and rehydrated using graded ethanols.

Immunostaining for p16, E-cadherin and vimentin was performed using the DAKO Tech-Mate 500 Plus Immunostainer (DAKO, Glostrup, Denmark). Antigen retrieval and antibody dilutions were optimized individually as follows. Using specific mouse monoclonal antibodies, a 1:100 dilution for p16 (clone DCS-50, NeoMarkers, Fremont, CA, USA), a 1:50 dilution for E-cadherin (clone NCH-38) and a 1:200 dilution for vimentin (clone Vim 3B4) (both DAKO) were applied for 25 minutes after microwave treatment twice for 7 minutes in citrate buffer (0.01 M, pH 6.0) at 600 W.

Developing was initiated by incubating the sections for 25 minutes with a biotinylated rabbit anti-mouse secondary antibody (DAKO). Slides were then immersed in Chem Mate Peroxidase-Blocking Solution (DAKO) three times for 2.5 minutes to block endogenous peroxidase activity. After 25 minutes of incubation with streptavidinhorseradish peroxidase (DAKO) and subsequent staining with chromogen AEC (DAKO) for 10 minutes, sections were counterstained with haematoxylin (DAKO), dehydrated and mounted.

Immunohistochemical staining for β-catenin was performed using the ABC-method.

Briefly, antigen retrieval was done using citrate buffer as mentioned above in a Deleted: (Table 1)

F o r P e e r R e v i e w 8
microwave oven (600 W) 3 times for 5 minutes. Endogenous peroxidase activity was blocked by incubating the sections in 0.3% H 2 O 2 in 0.001 M phosphate buffered saline for 30 minutes. This was followed by blocking avidin and biotin using the Avidin-Biotin Blocking Kit (Vector Laboratories, Burlingame, CA, USA). Primary monoclonal mouse antibody (BD Biosciences, Heidelberg, Germany) was applied at a dilution of 1:500 overnight at 4°C. A biotin-conjugated rabbit anti-mouse antibody (DAKO) was used as a secondary antibody at a dilution of 1:200 for 30 minutes.

Avidin-biotin complex-peroxidase (Vector Laboratories) was used for antibody detection and visualization with AEC was done as described above.

For quality control of staining, positive and negative controls were included in each staining run. The positive controls were normal tonsils for p16, normal breast tissue for E-cadherin, colorectal tumors for β-catenin and lymphoma tissue for vimentin. The vascular walls in the examined slides served as an internal positive control for vimentin as well. Staining in the absence of primary antibodies served as negative controls, in which no staining reaction was seen.

Scoring of immunohistochemical staining

All sections were independently examined by two experienced pathologists participating in the education of pathology at the International Academia of Pathology (IAP), and they were blinded to the patients' clinical information. First of all, both examiners determined the areas of carcinoma tissue within the samples. Then the staining pattern was evaluated semi quantitatively by comparing the intensity and cellular localisation of immunostaining with positive and negative controls. To determine percentage labelling indices, for each section usually more than 1000 cells were analysed in 5 randomly selected fields of tumor tissue using high-power (40x objective and 10x ocular) magnification. For each antibody, levels of immunoreaction were graded into 4 categories as follows: 0, <5%; 1+, 5-25%; 2+, 26-50%; 3+, >50% positive cells for E-cadherin, vimentin and p16, and 0, <1%; 1+, 1-10%; 2+, 11-25%; 3+, >25% positive cells for β-catenin, for practical and statistical purposes as based on the results. In tumors showing heterogeneity of staining, the immunostaining was judged according to the prominent pattern. Whenever there were some discrepancies in the evaluation between observers, a consensus was reached using a multiheaded microscope. In case of vimentin, we did an additional immunohistochemistry for the leucocyte common antigen (LCA) to rule out immune cells within tumor tissue.

Statistical analysis

All statistical analyses were performed using SPSS 17.0 (SPSS Inc., Chicago, IL, USA) for Microsoft Windows XP (Microsoft Corporation, Redmont, WA, USA). The association between experimental findings and clinicopathological variables was analysed using the χ 2 -test for categorical data. In case of individual datasets of primary tumor and metastasis of the patients, a paired Student's t-test was applied whenever suitable. A p-value < 0.05 was considered to be statistically significant.

RESULTS

Correlation between a positive HPV status and p16 expression was significant for primary tumors (p < 0.0001) as well as for metastases (p < 0.0001). In primary tumors, 84% of HPV-positive ones were p16-positive and 91% of HPV-negative ones were p16-negative. In metastases, 89% of HPV-positive ones were p16-positive and 97% of HPV-negative ones were p16-negative. There were 3 out of 48 patients with a mismatch of HPV-status and p16 expression in their metastasis. These 3 patients belonged to the group of 6 patients with a mismatch in their primary tumor.

Significant differences comparing the individual datasets were found for the expression of p16 being lower in metastases (35%) than in primary tumors (48%; p = 0.013) (Table 1) and for the patients HPV status being less frequently positive in metastases (38%) than in primary tumors (52%; p = 0.018). In other words, 7 out of 25 HPV-positive primary tumors had an HPV-negative metastasis and 5 out of these 7 were p16-positive in their primary tumor and p16-negative in their metastasis. These cases were of low T and N classification, of high grade and presented vimentin levels above the average and a prominent nuclear β-catenin staining. All 23 HPV-negative primary tumors had a negative HPV-specific PCR of their respective metastasis, and all 25 p16-negative primary tumors had p16-negative respective metastases.

For further analyses, according to Smeets et al. [START_REF] Smeets | A novel algorithm for reliable detection of human papillomavirus in paraffin embedded head and neck cancer specimen[END_REF] , the patients were categorized into p16-and HPV-positive ones (p16+/HPV+) and all others. For better understanding, the first group was classified as HPV-related tumors. The remaining cases were summarized as HPV-unrelated. Comparing individual datasets, this HPV-relation was less frequent in metastases (33%) than in primary tumors (44%; p = 0.024).

Deleted: (Table 2) Deleted: (Tables 2,3) All examined immunohistochemical markers could be detected in primary tumors and in metastases. The Cohen's kappa coefficients for the inter-observer agreements were as follows: 0.895 for p16, 0.690 for E-cadherin, 0.865 for β-catenin and 0.744 for vimentin. 13 out of 48 patients showed concurrent positive immunostaining for p16, β-catenin and vimentin both of the primary tumor and the metastasis respectively.

Deleted
Expression of E-cadherin was significantly less abundant in metastases (21%) than in primary tumors (56%; p < 0.001) (Table 1). This tendency became even more obvious in the subgroup of patients with an N2b, N2c or N3 status. Here, E-cadherin expression again was lower in metastases (23%) than in primary tumors (73%; p < 0.0001). E-cadherin expression in metastases of T3/4 primary tumors was significantly lower (0%) compared to metastases of T1/2 primary tumors (26%; p < 0.001).

In HPV-unrelated primary tumors the protein expression of E-cadherin was higher (59%) than in HPV-related primary tumors (52%; p > 0.05) (Table 1), and the above observed significantly reduced expression of E-cadherin in metastases was seen only in HPV-unrelated patients (p = 0.002).

All primary tumors and all but one metastases were positively stained for β-catenin.

In both primary tumors and metastases, three different patterns of β-catenin immunostaining were found: predominantly in the nucleus (more than 50% positively stained nuclei, irrespective of membranous staining), exclusively in the membranes or a combination of nuclear and membranous localisation (Table 1, Figure 1). For better understanding, these groups were classified as nuclear, membranous and combined.

Deleted: There was no significant interobserver variability in the blinded analysis of the slides. 3). Comparing individual datasets, a significant difference was found for a nuclear βcatenin expression being higher in metastases (46%) than in primary tumors (27%, p = 0.002) (Table 1). This was valid only for HPV-unrelated patients. Here, the protein expression of nuclear β-catenin was 44% in metastases and 15% in primary tumors (p = 0.016). In HPV-related primary tumors the nuclear localisation of β-catenin expression was significantly higher than in HPV-unrelated ones (p = 0.030) (Table 1).

Deleted: exclusively

Looking at nuclear β-catenin-positive tumor samples, 69% of primary tumors and 36% of metastases were HPV-related whereas 31% of primary tumors and 64% of metastases were HPV-unrelated (p = 0.060).

The nuclear β-catenin expression significantly correlated contrary to the E-cadherin expression in the whole group (n = 96; p = 0.017) and in HPV-related primary tumors (n = 21; p = 0.015). In HPV-unrelated primary tumors one was found to be nuclear βcatenin-positive and E-cadherin-negative whereas in HPV-unrelated metastasis 13 were nuclear β-catenin-positive and E-cadherin-negative. In HPV-positive primary tumors 7 were found to be nuclear β-catenin-positive and E-cadherin-negative whereas in HPV-positive metastasis 6 were nuclear β-catenin-positive and Ecadherin-negative (p = 0.008) (Figure 2). A significant correlation could be found between the expression levels of vimentin and the differentiation grade of the primary tumors (p = 0.021) with higher vimentin levels in less differentiated tumors (Table 1, Figure 3). HPV status, our study revealed the following major findings. First, the process metastasis of tonsillar cancer is characterised by a downregulation of E-cadherin and an upregulation of nuclear β-catenin and seems to occur at different time points depending on the patients HPV status. Furthermore, the expression of the mesenchymal marker protein vimentin significantly correlates with the grading of tonsillar carcinomas (Figure 4). Second, metastasis of tonsillar cancer is accompanied by a loss of HPV and p16 in some Third, we could show a strong correlation of p16 expression and the presence of HPV DNA also for lymph node metastases.

P16 can be used as a surrogate marker for the patients' HPV status in OSCC 12 and in HNSCC metastases [START_REF] Begum | Detection of human papillomavirus-16 in fine-needle aspirates to determine tumor origin in patients with metastatic squamous cell carcinoma of the head and neck[END_REF] . In our study we could confirm the strong correlation of p16 and the patients' HPV status not only for the primary tumor but for the first time for the metastasis of tonsil cancer as well. The signaling pathways, namely the pRbpathway seems to work even in the partly disturbed molecular interior of the metastasis with the consequence of elevated p16 in HPV-positive metastases.

However we observed a loss of HPV-DNA in 7 cases during metastasis.

Interestingly, in most of these cases, the p16 expression was lost in metastases as well. It can be speculated that the loss of HPV-DNA during metastasis might be the consequence of loss of cellular DNA at sites with integrated HPV-DNA or more likely the consequence of loss of episomal HPV-sequences. Due to loss of E7 expression in these cases p16 might not be overexpressed anymore. For future analysis it will be Dysfunction of E-cadherin-mediated cell adhesion results in an increased metastatic potential of cancer cells [START_REF] Tanaka | Expression of E-cadherin, alphacatenin, and beta-catenin in the process of lymph node metastasis in oral squamous cell carcinoma[END_REF][START_REF] Rimm | Cloning and characterization of the human invasion suppressor gene E-cadherin (CDH1)[END_REF] . In our study E-cadherin expression was decreased in HPV-related primary tumors compared to HPV-unrelated ones. Furthermore, a decrease of E-cadherin during metastasis could be shown only for HPV-unrelated patients. This suggests, that parts of tumor progression during metastasis, characterized by downregulation of E-cadherin have already taken place in the primary tumor in HPV-related patients, whereas in HPV-unrelated patients these parts of tumor progression might happen later during the process of metastasis.

Interestingly in a recent study with primary tumors we found a loss of chromosome 16q as a specific genetic change in HPV-related OSCC which might explain the Ecadherin inactivation in this study [START_REF] Klussmann | Genetic signatures of HPV-related and unrelated oropharyngeal carcinoma and their prognostic implications[END_REF] since E-cadherin is located on chromosome 16q22.1 40 .

Our data for the first time show a significantly increased expression of nuclear βcatenin in HPV-related primary tumors compared to HPV-unrelated ones.

Furthermore, only within HPV-unrelated patients nuclear β-catenin is significantly increased in metastases compared to the corresponding primary tumors. This gives rise to the question, whether nuclear β-catenin plays a role in tumor progression at different time points in HPV-related and -unrelated patients as well. Acting as a cofactor for the LEF/TCF transcription factor, β-catenin activates target genes, including matrilysin, fibronectin, cyclin D1, c-myc and other oncogenes, resulting in deregulation of cell cycle progression and unscheduled proliferation [START_REF] Jung | The invasion front of human colorectal adenocarcinomas shows co-localization of nuclear beta-catenin, cyclin D1, and p16INK4A and is a region of low proliferation[END_REF] . In HPV-related cervical lesions, increased nuclear staining of β-catenin was detected in carcinoma samples, whereas in cervical dysplasia no nuclear staining was seen 42 . Uren et al. downregulates E-cadherin and membranous β-catenin expression [START_REF] Moustafa | E6/E7 proteins of HPV type 16 and ErbB-2 cooperate to induce neoplastic transformation of primary normal oral epithelial cells[END_REF] . This is accompanied by nuclear translocation of β-catenin, which is then acting as a transcriptional regulator with subsequent upregulation of oncoproteins being accountable for tumor progression. Analogical to the findings on E-cadherin our results suggest that this form of tumor progression characterized by elevated expression of nuclear β-catenin already takes place in the primary tumors of HPVrelated patients whereas in HPV-unrelated ones it happens later during metastasis.

Taken together, the combination of decreasing E-cadherin and increasing nuclear βcatenin might reflect an important event in HPV-related tumor progression.

Interestingly, recent studies on human colorectal carcinomas and endometrial carcinoma cells showed that p16 is under direct transcriptional control of nuclear βcatenin [START_REF] Jung | The invasion front of human colorectal adenocarcinomas shows co-localization of nuclear beta-catenin, cyclin D1, and p16INK4A and is a region of low proliferation[END_REF][START_REF] Wassermann | p16INK4a is a beta-catenin target gene and indicates low survival in human colorectal tumors[END_REF][START_REF] Saegusa | Induction of p16INK4A mediated by beta-catenin in a TCF4-independent manner: implications for alterations in p16INK4A and pRb expression during trans-differentiation of endometrial carcinoma cells[END_REF] . Therefore it could be speculated that the p16 overexpression in HPVrelated OSCC is not only depending on pRb inactivation by E7 but also a direct effect of nuclear β-catenin translocation. Here p16 overexpression does not lead to cell cycle suppression because of the disabled downstream effector pRb which is being blocked by E7, resulting in tumor progression. Elzagheid et al. observed a high nuclear expression of β-catenin in primary tumors and metastases being a significant predictor of a more favourable disease free survival in colorectal cancer [START_REF] Elzagheid | Nuclear beta-catenin expression as a prognostic factor in advanced colorectal carcinoma[END_REF] and recent studies have shown that p16 is a significant predictor for a favourable prognosis for OSCC which is in concordance with HPV DNA positivity in most cases [START_REF] Ragin | Survival of squamous cell carcinoma of the head and neck in relation to human papillomavirus infection: review and meta-analysis[END_REF]12 . Immunostaining revealed that the expression of vimentin in primary tumor cells increased significantly with a higher grading. This gives an indication that epithelial cells undergo the mesenchymal activation process known as EMT [START_REF] Mannic | Expression of Cytokeratin Subtypes and Vimentin in Squamous Cell Carcinoma of the Floor of the Mouth and the Mobile Tongue[END_REF] . Taken together, the vimentin expression corresponding to the grading as well as an upregulation of nuclear β-catenin and a downregulation of E-cadherin during tumor progression our observations give evidence that EMT might play an early event in dedifferentiation and progression of tonsillar cancer.

In conclusion, we have newly identified possible molecular changes involved in tumor progression in tonsillar primary tumors and their respective metastases at different timepoints during metastasis depending on the patients HPV status.

(

  Response): To assess a directly carcinogenesis-related HPV involvement in FFPE specimens of HNSCC Smeets et al. propose a test algorithm of p16 immunochemistry followed by a HPV-specific

Figure :

 : Figure: Expression of the marker protein LCA as seen in various representative sections of the same patient seen in original Figure 3 B, original magnification x400.In contrast to the carcinoma tissue, LCA-positive cells were only found within the connecting tissue.

Formatted

  : nearly Deleted: exclusive Deleted: due to their HPV status determined by the above mentioned HPV-specific PCR. Categorizing patients into Deleted: and p16-

  Deleted: (p < 0.001) (Table

39 Deleted: 13 DISCUSSION

 3913 There was no such correlation for nuclear β-cateninnegative tumor samples. Deleted: , Deleted: negative tumor samples Deleted: 53 Deleted: 20 Deleted: and in HPV-negative metastases (n = 30; p = 0.032).Examining the expression of mesenchymal and epithelial proteins as well as p16 in 48 tonsillar carcinomas and their respective metastases in relation to the patients'

  to collect more of these cases and examine them in more detail concerning the viral loss.

  high β-catenin levels together with activation of the Wnt pathway may serve as a marker to detect malignant progression in cervical cancer. Al Moustafa et al. could show, that overexpression of E6/E7 in combination with ErbB-2 activation

1 LEGENDSFigure 1 :

 11 Figure 1: Localisation of the marker protein β-catenin within carcinoma cells in primary tumors and metastases, original magnification x200. (A) exclusive membranous expression in primary tumor tissue; (B) nuclear expression in primary tumor tissue (short white arrows, exemplarily); (C) combined membranous (long black arrows, exemplarily) and nuclear (short white arrows, exemplarily) expression in primary tumor tissue; (D) combined membranous and nuclear expression in metastasis.

Figure 2 :

 2 Figure 2: Distribution of tumor samples being E-cadherin-negative and β-cateninpositive in relation to the patients' HPV status.

Figure 3 :

 3 Figure 3: Expression of the mesenchymal marker protein vimentin as seen in various differentiation grades of the primary tumor, original magnification x100. (A) G1, endothelial cells and connective tissue (long black arrows) are vimentin positive and carcinoma cells within the tumor tissue (short white arrows) are vimentin negative; (B) G2, scattered vimentin positive carcinoma cells within the tumor tissue (short white arrows); (C) G3, nearly all carcinoma cells are vimentin positive.

Figure 4 :

 4 Figure 4: Diagram summarizing the main results concerning molecular changes during tumor progression and metastasis of HPV-related and -unrelated tonsillar carcinomas.

  Figure 1 254x190mm (96 x 96 DPI)

Figure 3

 3 Figure 3 B 233x175mm (150 x 150 DPI)

Figure 3

 3 Figure 3 C 233x175mm (150 x 150 DPI)

  is probably not necessary as the data are either extracted from Table2or included in the text." (Response): We now left away Table2. All relevant information is stated in the text.

  . De Araujo et al.

	demonstrated that the expression of the mesenchymal molecule vimentin was most
	marked in OSCC with high tumor grade
	F o r
	P
	e
	e r
	R
	e
	v i e
	w
	5

Table 1 :

 1 Expression of p16, E-cadherin, β-catenin and vimentin in primary tumors and metastases.

	Clinicopathological	n	p16			E-cadherin		β β β β-catenin		β β β β-catenin localisation			Vimentin
	feature	-	+	p	-	+	p	-	+	p	mem	nuc	both	p	#	-	+	p
	Total		96 56 (58%) 40 (42%) >0.05 59 (61%) 37 (39%) >0.05 1 (1%) 95 (99%)	>0.05	28 (29%)	35 (36%)	33 (34%)	>0.05 54 (56%) 42 (44%) >0.05
	Primary tumors Metastases	48 25 (52%) 23 (48%) 48 31 (65%) 17 (35%)	0.013 †	21 (44%) 27 (56%) 38 (79%) 10 (21%)	<0.001 †	0 (0%) 1 (2%)	48 (100%) 47 (98%)	>0.05	15 (31%) (27%) 13	13 (27%) (46%) 22	20 (42%) (27%) 13	0.002 †	24 (50%) 24 (50%) 30 (63%) 18 (38%)	>0.05
	p16 + / HPV + all others	37 59			22 (60%) 15 (41%) 37 (63%) 22 (37%)	>0.05	0 (0%) 1 (2%)	37 (100%) 58 (98%)	>0.05	12 (32%) (27%) 16	17 (46%) (31%) 18	8 (22%) (42%) 25	>0.05	20 (54%) 17 (46%) 34 (58%) 25 (42%)	>0.05
	Primary tumors
		G1, G1-2, G2 24 13 (54%) 11 (46%) G2-3, G3 24 12 (50%) 12 (50%)	>0.05	13 (54%) 11 (46%) 8 (33%) 16 (67%)	>0.05	0 (0%) 0 (0%)	24 (100%) 24 (100%)	>0.05	6 (25%) (38%) 9	7 (29%) (25%) 6	11 (46%) (38%) 9	>0.05	16 (67%) 8 (33%) 8 (33%) 16 (67%)	0.021
		p16 + / HPV + 21 all others 27			10 (48%) 11 (52%) 11 (41%) 16 (59%)	>0.05	0 (0%) 0 (0%)	21 (100%) 27 (100%)	>0.05	7 (33%) (30%) 8	9 (43%) (15%) 4	5 (24%) (56%) 15	0.030	11 (52%) 10 (48%) 13 (48%) 14 (52%)	>0.05
	Metastases
		p16 + / HPV + 16 all others 32			12 (75%) 4 (25%) 26 (81%) 6 (19%)	>0.05	0 (0%) 1 (3%)	16 (100%) 31 (97%)	>0.05	5 (31%) (25%) 8	8 (50%) (44%) 14	3 (19%) (31%) 10	>0.05	9 (56%) 21 (66%) 11 (34%) 7 (44%)	>0.05
	mem: membranous expression
	nuc:	nuclear expression
	#:	for nuclear β-catenin expression
	†:	by means of a paired Student's t-test
	p-value named only if level of significance < 0.05
	Percentage values rounded

Published on behalf of the British Division of the International Academy of Pathology Histopathology

ACKNOWLEDGEMENTS

This work was supported by the program for Research & Education of the University Hospital of Cologne to H.P. Dienes and by funding from the Jean-Uhrmacher-Foundation.

We thank S. Sattler and K. Wendland for their excellent technical assistance.

Bassima Yosef holds a fellowship of the Syrian government.

There are no financial disclosures from any authors. There do not exist any conflicts of interest.


Aurélie Boisbunon

Stephane Canu

Dominique Fourdrinier

Criteria for variable selection with dependence

come

Criteria for variable selection with dependence Aurélie Boisbunon

Stéphane Canu Dominique Fourdrinier

Context

Most theoretical tools for model selection rely on one main assumption: the independence of noise components. However, in many real regression examples this assumption is too strong and does not fit well the reality (see for instance the discussion on this issue in [START_REF] Fourdrinier | Improved model selection method for a regression function with dependent noise[END_REF] and references therein).

A way to take dependence into account is to model noise as a multivariate spherically symmetric random variable. This general framework includes the well-known multivariate Student and Kotz distributions. Our work consists in integrating this idea into model selection problems.

The problem of model selection implies three steps, as discussed in [START_REF] Guyon | Model selection: beyond the bayesian/frequentist divide[END_REF]: (a) the definition of a way to explore models, (b) the estimation of parameters for each model, and (c) the evaluation of the models. If we consider the example of Lasso [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF], the L 1 -penalization solves both step (a), with the regularization path algorithm [START_REF] Efron | Least angle regression[END_REF], and step (b) simultaneously, while step (c) is usually performed with Mallows' C p [START_REF] Mallows | Some comments on cp[END_REF] or cross-validation [START_REF] Stone | Cross-validatory choice and assessment of statistical predictions[END_REF].

We propose a new procedure for step (c) based on loss estimation. Loss estimation is a data-driven approach closely related to Stein's Unbiased Risk Estimation (SURE) [START_REF] Stein | Estimation of the mean of a multivariate normal distribution[END_REF], and has been extended to spherically symmetric distributions (see for instance [START_REF] Fourdrinier | Estimation of a loss function for spherically symmetric distributions in the general linear model[END_REF]). Our estimator of loss is of the following form:

δ ρ (y) = δ 0 (y) -ρ(y), (1)
where y ∈ R n corresponds to the data, δ 0 (y) is an unbiased estimator of loss (close to Mallows C p and AIC [START_REF] Akaike | Information theory and an extension of the maximum likelihood principle[END_REF] in the Gaussian case) and ρ(y) is a correction function. We will give in the sequel some details on how to build ρ as well as an explicit formula in a more specific context.

Note that our procedure should not be considered alone. Indeed, as efficient as it can be, it may not give satisfactory results if step (a) and (b) are not well handled. For instance, methods such as Lasso give an interesting and low-computational solution for step (a) exploring models through a regularization path, but not all of them propose good estimates of the parameters (see [START_REF] Bruce | Understanding waveshrink: Variance and bias estimation[END_REF] for instance). We suggest to use Firm Shrinkage [START_REF] Bruce | Understanding waveshrink: Variance and bias estimation[END_REF], which leads to the same model exploration as Lasso for step (a) and whose estimates for step (b) are much less biased than those induced by the L 1 -penalty. The choice of a regularization path for step (a) is driven by the increasing size of datasets, but we will discuss the limitations induced by searching only on a small number of models.

Loss estimation as a variable selection criterion

We focus our study on the linear regression model

Y = Xβ + ε, (2)
where Y is a random vector in R n , X is a fixed design matrix containing p explanatory variables, β ∈ R p is the unknown regression coefficient, and ε is the error vector in R n with mean zero. We assume that the distribution of ε is spherically symmetric. In the case where it admits a density, this one is of the form t → f ( t 2 ) for a given function f mapping from R + into R + , with . the Euclidean norm. What characterizes our work is that we do not need to specify the form of f , thus our results are more robust in terms of distribution. Another important feature of this family of distributions is the dependence between the components of ε, with the exception of the Gaussian distribution. However, ε i and ε j are not correlated, but correlation could be handled by the more general elliptical case. The family of spherically symmetric distributions covers, among others, the well known Student and Kotz distributions.

Model selection includes the specific problem of variable selection. In order to address this problem, we assume that only a small subset of variables in X are relevant to explain Y. The literature on this issue, corresponding to step (a), covers three main ways to perform the selection: the first one is to consider all possible subsets, leading to an exhaustive exploration, like in Subset Selection; the second one is to assess the significativity of each variable through statistical tests in order to decide whether to keep it or drop it, as is done in Forward and Backward Selection; the last one, which has been widely studied in the last 20 years, is the construction of a regularization path adding one variable at a time with respect to some criterion. We can cite from this latter approach methods such as the LARS algorithm for Lasso [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF][START_REF] Efron | Least angle regression[END_REF] or Firm Shrinkage [START_REF] Bruce | Understanding waveshrink: Variance and bias estimation[END_REF]. These methods generally depend on some hyperparameter tuning the level of sparsity, reducing the problem of variable selection to the optimization of this hyperparameter. The huge number of comparisons of the first two ways discussed here has led researchers to focus their interest on regularization path methods, especially in a context of large datasets.

In the sequel, we consider the Firm Shrinkage estimator, a method developed by Bruce and Gao in 1997 that reduces the large bias of Lasso. This way, Firm Shrinkage achieves a good estimation of parameter β and shares a low computational cost with Lasso. It is an arguable choice since there exist other Lasso-type methods with low bias, but it turns out to perform quite well on our simulations. Nevertheless, a comparison of such methods in terms of selection and estimation could lead to even better results. Note that Firm Shrinkage requires X to be orthogonal. Its estimator of β is

βF S i =    0 | βLS i | ≤ λ a( βLS i -λsign( βLS i ))/(a -1) λ < | βLS i | ≤ aλ βLS i | βLS i | > aλ , (3)
where βLS is the least-squares solution, λ > 0 is the hyperparameter tuning the sparsity, and a > 1 a hyperparameter tuning the bias. We set a to 2 in order to avoid the estimation of a second parameter and because it yields good results (see [START_REF] At Walden | Spectrum estimation by wavelet thresholding of multitaper estimators[END_REF]).

We now define more precisely our procedure of selection. Loss estimation is a decision theory tool originally designed for the problem of estimating the location parameter µ of a distribution P µ , like the multivariate Gaussian (see for instance [START_REF] Stein | Estimation of the mean of a multivariate normal distribution[END_REF][START_REF] Johnstone | On inadmissibility of some unbiased estimates of loss[END_REF]). It aims at evaluating the quality of some estimator of µ by estimating its loss function L(μ, µ) through the observations only, the real loss L(μ, µ) being inaccessible because of its dependence on the unknown parameter µ. The idea of loss estimation is closely related to Stein's Unbiased Risk Estimation (SURE), but the main difference lies in the fact that it does not only consider unbiased estimates. Unbiasedness is indeed not always an interesting property, and potential improvement can be achieved in terms of quadratic risk by adding a bias. We formalize this idea as follows. Let δ 0 be an unbiased estimator of loss L(μ, µ), that is to say δ 0 verifies the property

E µ [δ 0 ] = E µ [L(μ, µ)]
, where E µ is the expectation with respect to spherically symmetric distributions with location parameter µ. We look for corrections ρ(y) in (1) such that

R µ (δ ρ ) = E µ [(δ ρ -L(μ, µ)) 2 ] = E µ [(δ 0 -ρ(y) -L(μ, µ)) 2 ] ≤ E µ [(δ 0 -L(μ, µ)) 2 ] = R µ (δ 0 ), (4)
and with strict inequality at least for one value of µ. A common loss for L is the quadratic loss L(μ, µ) = μ -µ 2 , since it allows easy computations. It is also interesting to use this loss function in combination with nearly unbiased estimators of µ, like Firm Shrinkage for µ = Xβ, as it has a high probability of reaching its minimum for the true model, as soon as this latter one belongs to the set of models explored.

From now on we will restrict the study to the estimation of the quadratic risk of the Firm Shrinkage estimator. The corrective estimator of loss is δ ρ (y) = δ 0 (y) -ρ(y)

where

δ 0 = 2df -n n -p y -X βLS 2 + y -X βF S 2 (6)
and ρ is a correction function. Note that the number of degrees of freedom for Firm Shrinkage is

df = k + #{i \ λ < |β LS i | ≤ aλ}/(a -1)
with k the number of selected variables. We consider correction functions of the form

ρ(y) = y -X βLS 4 γ(y) (7)
where γ(y) is a twice weakly differentiable function. In particular, we study the function

γ(y) = C k max i≤p {(q T i y) 2 \ |q T i y| ≤ λ} + j≤p (q T j y) 2 1 {|q T j y|≤λ} -1 , (8)
where q i corresponds to the i th column of matrix Q in the QR decomposition of X. The constant C leading to the lower approximate quadratic risk in ( 4) is

C = 2(p -2 -2k(k + 1)/p + 10 -4 λ 2 (n + 4)(n -2p + 2df )/(n -p)) (n -p + 4)(n -p + 6) (9)

Simulations

The following example is inspired by [START_REF] Fourdrinier | Comparaisons de procédures de sélection d'un modèle de régression: une approche décisionnelle[END_REF]: y is a vector of n = 40 observations from the random variable Y, and X contains p = 5 explanatory variables. The regression coefficient β is set to (2, 0, 0, 4, 0) T . The error vector ε is drawn from two different distributions with variance 1: the Gaussian distribution, which is the usual assumption for most criteria, and the Student distribution with ν = 4 degrees of freedom, corresponding to our assumption of spherical symmetry. We replicate this error vector 5000 times, this way we obtain 5000 regularization paths.

The estimators of loss δ 0 and δ ρ are those described in ( 6) and ( 5) with the choices of the corrective function in ( 7), ( 8) and ( 9). We compare their selection to the real quadratic loss L(μ, µ) as well as to the classical Akaike's Information Criterion (AIC) and Schwarz' Information Criterion (BIC) with a Gaussian assumption, and the distribution-free leave-one-out cross-validation (LOOCV).

Tables 1 and2 present for each method the empirical probabilities of selecting the subsets over the 5000 replicates. We iterated the experience ten times to estimate the means and standard deviations shown in the tables. Only the ten most voted subsets are displayed. Here, we can see that the corrective estimator δ ρ selects the right model {1, 4} with higher probability than the unbiased estimator δ 0 and the classical methods even for the usual i.i.d. Gaussian case, and is closer to the real loss results. Improvement might be even larger if we consider a more general form of correction, like δ * = α(δ 0 -y -X βF S 4 γ) for instance, as was done in [START_REF] Fourdrinier | Comparaisons de procédures de sélection d'un modèle de régression: une approche décisionnelle[END_REF]. But the most striking result is the low empirical probabilities of all the methods. The real loss L(μ, µ) manages to select the right subset only around half of the time, bounding from above the probability of selection with our criterion. On the contrary, the results in [START_REF] Fourdrinier | Comparaisons de procédures de sélection d'un modèle de régression: une approche décisionnelle[END_REF] showed that systematic exploration with subset selection leads to select {1, 4} with an empirical probability of 83% for their corrective estimator (which is slightly different from the one in ( 1)), that is to say twice the probability obtained with Firm Shrinkage. This important difference is a consequence of the regularization path, which sometimes introduces irrelevant variables first and fails to select the right subset. It also explains the selection by δ 0 and AIC of subsets of size 3 containing the two relevant variables, with empirical probabilities around 7%, or subsets of bigger size by δ ρ and LOOCV.

Discussion

We proposed a new data-driven procedure of model selection based on loss estimation and valid for the whole family of spherically symmetric distributions, allowing dependence between noise components. From this procedure, we derived a criterion for the Firm Shrinkage estimator, a regularization path method with low bias, evaluated by a quadratic loss. In the experience we drove, this criterion performs better than classical criterion like AIC, BIC and LOOCV, even under the usual i.i.d. Gaussian assumption for the noise.

Note that this procedure can be applied to other estimators and loss functions, and this way it could be adapted to problems such as classification or reinforcement learning.

Table 1 :

 1 Empirical probabilities (%) of selection with Firm Shrinkage (Gaussian case).

	Subset	δ0	δρ	AIC	BIC	LOOCV	L( μ, µ)
	{4}	20.18 (0.59) 26.12 (0.56) 20.18 (0.59) 40.05 (0.83) 14.42(16.18) 14.17 (0.43)
	{1,4}	39.02 (0.74) 44.41 (0.60) 39.02 (0.74) 39.37 (0.49) 32.71(12.27) 54.29 (0.56)
	{2,4}	2.09 (0.22)	3.08 (0.24)	2.09 (0.22)	1.51 (0.20)	3.66 (1.55)	0.00 (0.11)
	{3,4}	2.11 (0.21)	3.23 (0.32)	2.11 (0.21)	1.54 (0.16)	3.17 (1.09)	0.00 (0.15)
	{4,5}	2.05 (0.13)	2.74 (0.19)	2.05 (0.13)	1.47 (0.20)	3.60 (1.56)	0.00 (0.07)
	{1,2,4}	7.57 (0.34)	1.33 (0.15)	7.57 (0.34)	3.66 (0.26)	5.68 (3.06)	7.46 (0.33)
	{1,3,4}	7.83 (0.40)	1.30 (0.13)	7.83 (0.40)	3.73 (0.19)	6.93 (2.99)	7.63 (0.32)
	{1,4,5}	7.73 (0.40)	2.13 (0.20)	7.73 (0.40)	3.73 (0.36)	6.49 (3.73)	7.87 (0.27)
	{1,2,3,4}	2.62 (0.20)	5.03 (0.39)	2.62 (0.20)	0.00 (0.12)	2.56 (1.35)	1.96 (0.22)
	{1,2,3,4,5}	1.34 (0.13)	3.75 (0.23)	1.34 (0.13)	0.00 (0.04) 11.58 (6.84)	1.04 (0.18)
	Subset	δ0	δρ	AIC	BIC	LOOCV	L( μ, µ)
	∅	9.94 (0.65)	8.87 (0.43)	9.94 (0.65) 20.90 (0.74)	7.21 (3.12) 14.62 (0.45)
	{4}	15.77 (0.37) 19.11 (0.29) 15.77 (0.37) 24.33 (0.45) 12.63 (8.99) 14.88 (0.50)
	{1,4}	32.08 (0.74) 38.01 (0.62) 32.08 (0.74) 35.15 (0.82) 26.35(11.77) 46.08 (0.78)
	{1,2,4}	6.08 (0.21)	0.00 (0.14)	6.08 (0.21)	2.74 (0.16)	5.82 (2.93)	4.65 (0.21)
	{1,3,4}	5.97 (0.19)	0.00 (0.19)	5.97 (0.19)	2.75 (0.27)	5.62 (3.48)	4.72 (0.20)
	{1,4,5}	6.21 (0.36)	1.63 (0.22)	6.21 (0.36)	2.83 (0.20)	6.58 (3.39)	4.50 (0.16)
	{1,2,3,4}	2.12 (0.19)	4.03 (0.29)	2.12 (0.19)	0.00 (0.08)	2.84 (1.47)	1.34 (0.13)
	{1,2,4,5}	2.01 (0.25)	2.23 (0.23)	2.01 (0.25)	0.00 (0.06)	2.64 (1.14)	1.30 (0.15)
	{1,3,4,5}	2.09 (0.19)	2.15 (0.19)	2.09 (0.19)	0.00 (0.08)	2.67 (1.19)	1.37 (0.19)
	{1,2,3,4,5}	1.04 (0.13)	3.10 (0.16)	1.04 (0.13)	0.00 (0.05) 11.45 (3.47)	0.00 (0.11)

Table 2 :

 2

Empirical probabilities (%) of selection with Firm Shrinkage (Student case).


Hubert Stahn

Agnès Tomini

Rainwater Harvesting under Endogenous Capacity of Storage: a solution to aquifer preservation?

Keywords: Rainwater Harvesting, Conjunctive Use, Groundwater Management, Optimal Control JEL Classi…cation: Q25, C61, D61

This paper studies groundwater management in the presence of rainwater harvesting (RWH). We propose a two-state model that takes into account the standard dynamics of the aquifer and the dynamics of the storage capacity and we assume that the collection of rainwater reduces the natural recharge. We analyze the trade-o¤ between these two water harvesting techniques in an optimal control model. In particular, we show that when these techniques are pure substitutes, the development of RWH leads in the long run to a depletion of the water table even if pumping is reduced.

Introduction

Rainwater Harvesting (RWH) is a simple technique that has been used for thousands of years. Today, this practice is enjoying a revival in popularity and an international network, the International Water Association 1 , promotes and supports RWH initiatives worldwide as an important component in the sustainable provision of freshwater. However,

6

Corresponding author Email addresses: hubert.stahn@univmed.fr (Hubert Stahn ), agnes.tomini@univmed.fr (Agnès Tomini )

1 See the association's website http://www.iwahq.org/home/ the co-existence of this technique with groundwater withdrawals is not always a success story, as it can lead to depletion of the water table.

RWH was encouraged in the early 20th century in regions with no signi…cant river systems and remote and arid areas, sometimes on a large scale. For instance, the United Nations Environment Programme2 reports that in 2000, the Gansu province in China, which has a semi-arid to arid continental climate, reached a capacity of 2,183,000 rainwater tanks, generating 73.1 million cubic meters of drinking and supplementary irrigation water. [START_REF] Perrens | Design Strategy for Domestic Rainwater Systems in Australia[END_REF] estimates that in Australia approximately one million people rely on rainwater as their primary source of supply. Over and above these …gures, this alternative supply-side measure should contribute towards a number of Millennium Development Goals, especially that of halving the proportion of the population without sustainable access to safe drinking water by 2015 [START_REF] Lehmann | Low-Cost TEchnologies Toward Achieving the Millenium Development Goals : The Case of Rainwater Harvesting[END_REF]. The Rainwater Partnership was established in 2004 under the initiative of the UNEP not only to promote this strategy but also to integrate it into water resources management policies. This alternative water source can therefore co-exist with other water sources, providing decision-makers with a choice between withdrawing water from the ground and investing in development of a storage technology such as cisterns or rain tanks to harvest rainwater. Our paper thus aims to extend the canonical groundwater management model to take into account the development of RWH technology.

The literature on groundwater economics is quite well-developed. Apart from the strand dealing with quality aspects, studies focus mainly on pumping patterns and strategic behaviors, based on the dynamic model proposed by [START_REF] Gisser | Competition versus Optimal Control in Groundwater Pumping[END_REF] (e.g., [START_REF] Gisser | Groundwater : Focusing on the Real Issue[END_REF][START_REF] Provencher | A private property rights regime to replenish a groundwater aquifer[END_REF][START_REF] Provencher | The Externalities associated with the Common Property Exploitation of Groundwater[END_REF][START_REF] Rubio | Strategic Behavior and E¢ ciency in the Common Property Extraction of Groundwater[END_REF].

This literature often examines the conjunctive use of two sources of water, one of which is a stock and the other a ‡ow [START_REF] Provencher | Issues in the Conjunctive Use of Surface Water and Groundwater[END_REF]. The latter is considered as uncertain and groundwater then plays a role in protecting users against this uncertainty [START_REF] Tsur | The Stabilization Role of Groundwater when Surface Water Supplies are Uncertain : The Implications for Groundwater Development[END_REF][START_REF] Tsur | The Bu¤er Value of Groundwater with Stochastic Surface Water Supplies[END_REF][START_REF] Knapp | The Economics of Conjunctive Groundwater Management with Stochastic Surface Supplies[END_REF][START_REF] Provencher | The Externalities associated with the Common Property Exploitation of Groundwater[END_REF].

Our work follows this tradition, with an economic model of conjunctive use of two water sources. However, we seek to determine the optimal extraction of groundwater in addition to the optimal level of investment in RWH facilities without uncertainty.

This choice, as in the usual backstop literature 3 , will result from a trade-o¤ between groundwater pumping costs and RWH investment costs. However, we are dealing with two water sources that are typically interdependent: rainwater that is harvested cannot reach the aquifer. This implies that the quantity of RWH a¤ects the dynamics of the aquifer and even the marginal pumping cost, which depends on the aquifer head. The idea of this paper is to explore this complex dynamics and especially to show that in the long run, the introduction of RWH negatively a¤ects the head of the aquifer.

Our analysis is based on an optimal control framework. We propose a two-state model with pure state constraints in order to combine the standard dynamics of groundwater stock and the progressive development of RWH capacity. In fact, we extend the choices available to a social planner managing an aquifer by adding the opportunity to constitute RWH capacity. The groundwater extraction part is modeled in a standard and quite general way while the constitution of RWH capacity follows the literature on capital accumulation by identifying an investment cost and a progressive depreciation of the stock of this capacity. Within this context, we specify the optimality conditions and examine the long term e¤ects.

In a …rst step, we de…ne the necessary conditions for optimality and assume that a steady state exists. This gives us the opportunity to compare the long-term groundwater level to the level that would be obtained in the same model without the opportunity to invest in RWH capacity. We essentially show that the introduction of this storage capacity leads to depletion of the aquifer. However, this result relies on the existence of a steady state. In a second step, we display that these necessary conditions are also su¢ cient under a standard transversality condition and that all the trajectories that satisfy these conditions admit a unique steady state, which is even a local saddle point.

This paper is organized as follows. Section 2 introduces the setting of our model.

In section 3, we specify and discuss the necessary conditions for optimality. Section 4 studies the long-run e¤ects of the development of RWH capacity. In section 5, we focus 3 For the classical backstop literature the reader is referred to [START_REF] Heal | The Relationship between Price and Extraction Cost for a Resource with a Backstop Technology[END_REF], and [START_REF] Dasgupta | Economic Theory and Exhaustible Resources[END_REF]. Speci…c applications to groundwater resource can be found in [START_REF] Amigues | On the Optimal Order of Natural Resource Use when the Capacity of the Inhexhaustible Substitute is Limited[END_REF], [START_REF] Holland | Cadillac Desert Revisited : Property Rights, Public Policy, and Water-Resource Depletion[END_REF], [START_REF] Kim | Public Policies in Water-Resource Use : Their Effect on Groundwater Mining and Surface-Water Imports[END_REF], [START_REF] Koundouri | Dynamic adaptation to resource scarcity and backstop availability : Theory and Application to Groundwater[END_REF] or Kruce and al. (1997).

on the su¢ cient conditions in order to analyze the properties of the steady state. Section 6 provides some concluding comments.

The model

We start from a stylized dynamic and continuous time model of groundwater management. We assume that the aquifer is delineated by a ‡at bottom and perpendicular sides.

We denote by h(t) 2 0; h the water height where h > 0 stands for the maximum level of the aquifer. This water table simply rises because of recharge, which soaks through the soil to the aquifer, and falls because of water pumping4 . In this stylized world, a constant quantity of RW R > 0 is available at each t. Ordinarily, this quantity goes to the aquifer, but it can be stored in tanks if this storage capacity is developed. So let us denote by w s (t) the capacity accumulated up to time t. The rest of the water needs are covered by groundwater (GW) pumping denoted by w g (t). In this basic case, the dynamics of the aquifer is given by:

_ h(t) = (R w s (t)) w g (t) (1)
The total amount of water available at time t is given by w(t) = (w g (t) + w s (t)). The instantaneous bene…t induced by water use is estimated by a function F (w) that behaves like a standard production function.

Assumption 1 (A1). The social bene…t of the use of water is measured by a C 1 function F : R + ! R + which satis…es (i)8w > 0, F 0 (w) > 0, (ii) F " (w) < 0, (iii) F (0) = 0, (iv)

lim w!0 + F 0 (w) = +1 and lim w!+1 F 0 (w) = 0
Groundwater withdrawals induce instantaneous pumping costs C(w g ; h) depending on the amount that is pumped and on the height of the water table. As usually in this literature [START_REF] Koundouri | Current Issues in the Economics of Groundwater Resource Management[END_REF][START_REF] Fisher | Adjusting to Climate Change : Implications of Increased Variability and Asymmetric Adjustment Costs for Investment in Water Reserves[END_REF][START_REF] Rubio | Strategic Behavior and E¢ ciency in the Common Property Extraction of Groundwater[END_REF], we assume that @ 2 wg;h C(w g ; h) < 0 in order to capture two basic principles related to water extraction. The …rst stems from the fact that when the water table is lower, it is costlier to extract water because the resource must be pumped further. In other words, the marginal extraction cost is decreasing with the water height. The second principle is related to the dynamics of the model. Using an additional unit of water, at a given period of time, lowers the water table and so raises future extraction costs. Therefore, we simply say that the marginal e¤ect of a lowering of the water table on the future extraction cost increases with the level of pumping; this helps to slow down pumping.

Apart from this speci…c assumption, we also introduce usual assumptions on this cost function. This function is increasing with the level of extraction, decreasing with the height of water table and strictly convex:

Assumption 2 (A2). The groundwater extraction costs are given by a C 1 function

C : R + 0; h ! R + (i) 8(w g ; h), @ 2
wg;h C(w g ; h) < 0, (ii) 8(w g ; h), @ wg C(w g ; h) > 0; @ h C(w g ; h) < 0 and @ 2 C(w g ; h) is positive de…nite

The alternative water source requires an investment I(t) in order to build and to maintain a storage capacity w s (t). This investment is made ineach period at some cost (I) and the dynamics of the storage capacity takes into account an instantaneous depreciation which is measured by the function (w s ). The dynamics of the capital stock across time is therefore given by the relation:

_ w s (t) = I(t) (w s (t)) (2)
The depreciation function (w s ) and the cost function (I) satisfy the following conditions:

Assumption 3 (A3). The water harvesting technique is characterized by a C 1 investment cost : R ! R + and a C 1 depreciation function : R + ! R + which respectively satisfy:

(i) 8I; 0 (I) > 0, 00 (I) > 0 and lim I!0 (I) = lim I!0 0 (I) = 0

(ii) 8w s , 0 (w s ) > 0, " (w s ) > 0 and lim ws!0 (w s ) = 0

Within this stylized framework, we consider the problem of a social planner who must choose the optimal GW extraction and investment path to maximize the discounted social net bene…t. The discounting rate is given by > 0. We say that:

De…nition 1. An optimal GW extraction and investment path (w g (t); I (t)) solves the following control problem:

max w(t);I(t) Z 1 0 [F (w g (t) + w s (t)) C (w g (t); h(t)) (I(t))] e t dt (3)
subject to :

8 > > > > > > < > > > > > > : _ h(t) = R (w g (t) + w s (t)) ; with h(0) = h _ w s (t) = I(t) (w s (t)) with w r (0) = 0 w g (t) 0, h(t) 0, h h(t) 0 w s (t) 0, R w s (t) 0
We …nally introduce some restrictions on our parameters R, h and to maintain the consistency of the model. We say that:

If the aquifer is full (h = h) and the natural recharge is extracted (w g = R), there is always an incentive, at least marginally, to pump an additional quantity of water even taking into account the increase in future extraction costs induced by a change in the water table (i.e. 1 @ h C(R; h))

Conversely, the marginal cost of extracting one unit of water when the aquifer is empty (h ! 0) is very high, at least higher than the marginal productivity of the recharge even for very small extraction levels (w g ! 0). This quantity is nevertheless bounded by the long-term marginal investment cost for a capacity …xed at the recharge R. Apart from the mathematical convenience, this last restriction simply ensures that both sources of water are always used, since it is too costly to develop a capacity that captures the whole recharge.

These assumptions can be written as:

Assumption 4 (A4). Let us assume that:

(i) F 0 (R) > @ wg C(R; h) 1 @ h C(R; h) (ii) 0 ( (R)) ( + 0 (R)) > @ wg C(0; 0) > F 0 (R)

The Hamiltonian necessary conditions

Let us …rst study the social planner problem. We are dealing with an autonomous optimal control problem with mixed and pure state constraints and we decide to adopt the direct adjoining approach (see [START_REF] Hartl | A Survey of the Maximum Principles for Optimal Control Problems with State Constraints[END_REF]. If p 1 (t) and p 2 (t)

stand for the two co-state variables related respectively to the dynamics of the aquifer and of the harvesting capacity in de…nition 1, the present value Hamiltonian associated with this program is given by5 :

H = [F (w g + w s ) C(w g ; h) (I)] + p 1 [R (w g + w s )] + p 2 [I (w s )]
It remains for us to associate the Lagrangian multipliers with each constraint. Since we have one mixed and three pure state constraints, we denote these multipliers by (q i (t)) 5 i=1 and de…ne the following Lagrangian:

L = H + q 1 w g + q 2 h + q 3 h h + q 4 w s + q 5 (R w s )
Before going further, we …rst verify that::

Lemma 1. The problem given by de…nition 1 satis…es the constraint quali…cation property. Under A1(ii) and A2(ii), the Hamiltonian is also strictly concave with respect to the controls Since the quali…cation constraints are veri…ed, we can now begin to specify the necessary optimality conditions 6 . As usual, we know that the controls maximize the Hamiltonian along the optimal path of the state variables and verify the mixed constraints associated with the program. Since the Hamiltonian is concave, these necessary conditions become:

8 > > < > > :
@ wg L =F 0 (w g + w s ) @ wg C(w g ; h) p 1 + q 1 = 0

@ I L = 0 (I) + p 2 = 0
q 1 w g = 0 and q 1 0 (CTR)

Moreover, we know that the dynamics of the state and co-state variables satisfy the following conditions:

8 > > > > < > > > > : _ h = R (w g + w s ) ; with h(0) = h _ w s = I (w s ) with w r (0) = 0 _ p 1 = p 1 @ h L = p 1 + @ h C (w g ; h) q 2 + q 3 _ p 2 = p 2 @ ws L = p 2 ( + 0 (w s )) F 0 (w g + w s ) + p 1 q 4 + q 5 (DYN)
In this problem, we also have to manage pure state constraints. We therefore have the following slackness conditions:

8 > > < > > :
q 2 h = 0; and q 3 h h = 0, q 4 w s = 0 and q 5 (R w s ) = 0

(q i ) 5 i=2 0 (SLC)
Finally, for any entry, exit or contact time where the pure state constraints are binding, the co-state variables may have a discontinuity. At this junction point, there exists

(k i; ) 4 i=1
0 with the property that:

( p 1 ( + ) = p 1 ( ) k 1; + k 2; p 2 ( + ) = p 2 ( ) k 3; + k 4; and ( k 1; h( ) = 0, k 2; h h( ) = 0 k 3; w r ( ) = 0, k 4; (R w r ( )) = 0 (JMP)

The long-term e¤ect of RWH on the aquifer

In this section, we highlight the fact that the introduction of RWH lowers the longterm level of the water table. This result is obtained by considering the previous necessary conditions and by assuming that a steady state exists. Under this assumption, we actually show that the water table and the level of groundwater extraction must verify certain relations over the long run. In a second step, we brie ‡y recall the relations that must hold between these variables when there is no RWH. By comparing these two sets of relations, we …nally conclude that the implementation of RWH depletes the aquifer.

Some useful properties of the steady state

We begin with the conditions that characterize our steady state. By de…nition, if such a point exists, the controls must be chosen in a optimal way. It follows by (CT R) that:

p 1 = F 0 w s + w g @ wg C(w g ; h ) + q 1 and p 2 = 0 ( (w s )) (4)
Moreover, the stock variable are invariant, thus by (DY N ), we have:

8 > > < > > : R = w s + w g , I = (w s ) p 1 = 1 @ h C w g ; h q 2 + q 3 p 2 ( + 0 (w s )) = F 0 (w s + w g ) p 1 + q 4 q 5 (5)
Finally, we also know that the slackness conditions given by (SLC) are met.

Under our conjunctive use assumption, we immediately observe that the total water use w s + w g must be equal to the recharge R. One of these variables can therefore be forgotten, for instance the RWH level. Moreover, if we bear in mind that system (4) de…nes the co-state variables, we can reduce the previous system to a set of conditions by using the two last equations of (5), which only involves the water table, the level of groundwater extraction and the Lagrangian multipliers. By adding the slackness conditions, we can therefore say that the water table and the groundwater extraction level satisfy7 :

8 > > > > < > > > > : F 0 (R) @ wg C(w g ; h ) + 1 @ h C w g ; h q 2 + q 3 + q 1 = 0 0 ( (R w g )) + 0 (R w g
@ wg C(w g ; h ) + q 1 q 4 + q 5 = 0 and ( q 1 w g = 0, q 2 h = 0; q 3 h h = 0 q 4 R w g = 0; q 5 w g = 0; (q i ) 5 i=1 0

In the spirit of conjunctive use, we should even go a step further and ensure that both sources of water are really used at the steady state. For instance, if no rainwater is used we are back to the standard model and this comparison makes no sense. This was exactly the motivation behind the introduction of assumption A4(ii). Since w s + w g = R, we should be able to say that:

Lemma 2. Under A2(i), A3 and A4(ii), we observe that w g 2 ]0; R[. It follows that we can set q 1 = q 4 = q 5 = 0 in system (6)

Proof :Let us …rst assume that w g = 0. By the previous slackness conditions, we can say that q 4 = 0 and the second equation of ( 6) becomes q 1 + q 5 = @ wg C(0; h ) 0 ( (R)) ( + 0 (R). But under A2(i)

and A4(ii), we observe that:

q 1 + q 5 < @ wg C(0; 0) 0 ( (R)) ( + 0 (R) < 0
which contradicts the non-negativity of the multipliers. Now suppose that w g = R > 0, it follows that q 1 = q 5 = 0 and,under A3, the second equation of ( 6) becomes q 4 = @ wg C(R; h ) < 0, which is a contradiction.

From all these observations, we can say that:

Remark 1. At a steady state, the water table h 2 0; h and the groundwater extraction level w g 2 ]0; R[ satisfy the following set of relations

8 > > > < > > > : F 0 (R) @ wg C(w g ; h ) + 1 @ h C w g ; h q 2 + q 3 = 0 0 ( (R w g )) + 0 (R w g @ wg C(w g ; h ) = 0
q 2 h = 0; q 3 h h = 0 and q 2 ; q 3 0 (7)

The steady state without RHW

This steady state follows from an optimal control problem which is in some sense a subcase of the problem stated in de…nition 1 and is given by: max

$(t) Z 1 0 [F ($(t)) C ($(t); (t))] e t dt (8)
subject to :

( _ (t) = R $(t) with (0) = h $(t) 0, (t) 0, h (t) 0
To avoid confusion, we denote the water table by (t) and the level of GW extraction by $(t). If we denote by the co-state variable and by ( i ) 3 i=1 the Lagrangian multipliers associated with the three non-negativity constraints, it is immediate that the necessary conditions for optimality are simply a subset of the more general conditions in which RWH is allowed. A quick inspection of these conditions leads to the following de…nition of the steady state. The unique control must be optimal, the water consumption must be equal to the recharge, and _ = 0.

> > <

> > :

= F 0 ($ ) @ wg C($ ; ) + 1 $ = R = 1 ( @ C ($ ; ) + 2 3 )
and the following slackness conditions have to hold:

1 $ = 0, 2 = 0, 3 h = 0 with ( i ) 3 i=1 0
Since the extraction of groundwater w g must be equal to the total recharge R, we can say that:

Remark 2. The water table without RWH satis…es:

8 < : F 0 (R) @ wg C(R; ) + 1 (@ C (R; ) 2 + 3 ) = 0 2 = 0, 3 h = 0 with ( i ) 3 i=2 0 (9)
Now we can compare the water table that solves (9) with h , which solves (7).

The e¤ect of RHW

Let us now verify that the introduction of RWH leads to a depletion of the aquifer.

This result is trivial in two following cases.

Remark 3. By using the slackness conditions, we can observe that:

(i) if q 2 > 0, the aquifer with RWH is empty, i.e. h = 0 solves system (7). It follows trivially that h < unless = 0.

(ii) if 3 > 0; the aquifer without RWH is full, i.e. = h veri…es system (9) and trivially

h < unless h = h
We can, in the rest of the discussion, set q 2 = 0 and 3 = 0 (7) and ( 9) respectively

These two remarks simplify the argument since the …rst equations of ( 7) and ( 9)

become respectively:

( F 0 (R) = @ wg C(w g ; h ) 1 @ h C w g ; h q 3 F 0 (R) = @ wg C(R; ) 1 @ C (R; ) + 2
Let us now recall by lemma 2 that w g < R. If we also bear in mind that the total pumping cost, by A2(ii), is increasing in w g , the …rst equation becomes:

F 0 (R) < @ wg C(R; h ) 1 @ h C (R; h ) q 3
It follows by using the second equation that:

@ wg C(R; ) 1 @ C (R; ) + 2 < @ wg C(R; h ) 1 @ h C (R; h ) q 3
Since 2 ; q 3 0, this implies that:

@ wg C(R; ) 1 @ C (R; ) < @ wg C(R; h ) 1 @ h C (R; h )
Finally, we know by A2(i)(ii) that the total pumping cost is decreasing in the water table.

We can therefore say that:

Proposition 1. Under our assumptions, the introduction of RW harvesting, even if it is welfare improving, leads to environmental damage to the aquifer since h < unless the aquifer is either empty or full.

The existence of a steady state

At this point, we know that any optimal path which solves our water management problem and admits a steady state has the property that the long-term water table is lower than that obtained without RWH. To strengthen this result, we now show that path that satis…es the su¢ cient optimality conditions admits a unique steady state which is locally a saddle point. To simplify this discussion, we restrict our analysis to interior solutions. Up to now we have not used assumption A4(i), which, as we will see, ensures that the aquifer is never full or empty.

According to the Mangasarian su¢ cient conditions (see Hartl, Sethi and Vickson, 1995, theorem 8.4 p.203), we simply make sure that the Hamiltonian is strictly concave in the controls and the state variable while the constraints are quasi-concave in these variables.

By imposing a transversality condition, we can state that:

Lemma 3. A path wg ; Ĩ; wr ; h that satis…es CRT, DYN, SLC is an optimal path provided that:

lim t!+1 p 1 (t) h(t) h(t) + p 2 (t) (w r (t) wr (t)) = 0
Let us now show that a unique steady state exists. From the discussion that we have conducted in sub-section 4.1, we know that the steady state veri…es system (7). Moreover if we now use assumption A4, we can verify that:

Lemma 4. Using now A4(i), the water table is never full nor empty. This means that the steady state solves:

8 < : F 0 (R) @ wg C(w g ; h ) + 1 @ h C w g ; h = 0 0 ( (R w g )) + 0 (R w g @ wg C(w g ; h ) = 0 (10)
Proof Let us assume that the water table that solves system (7) of subsection 4.1 is respectively full or empty, and let us exhibit a contradiction on the sign of Lagrangian multipliers. So, if h = h, the slackness conditions entail that q 2 = 0 and the …rst equation of system (7) becomes:

1 q 3 = F 0 (R) + @ wg C(w g ; h) 1 @ d C w g ; h
But under A2(ii), the right-hand side of this equality is increasing in w g . Moreover, by lemma 2 w g < R and A4(i), we obtain the following contradiction:

1 q 3 < F 0 (R) + @ wg C(R; h) 1 @ d C R; h < 0
Now assume that h = 0, then q 3 = 0 and the …rst equation of ( 7) can be written:

1 q 2 = F 0 (R) @ wg C(w g ; 0) + 1 @ h C w g ; 0
Since @ h C w g ; 0 0, and @ 2 wgwg C(w g ; h) > 0 , q 2 is bounded from above by F 0 (R) @ wg C(0; 0)

which is negative under A4(ii)
The question of the existence and the uniqueness of the steady state can therefore be reduced to the study of system (10). Given our general setting, the proof of these results relies on di¤erential topology and uses a degree theory argument (see for instance Mas-Colell, 1985, p.207-08).

Proposition 2. There exists a unique interior steady state.

It remains for us to verify the stability properties of this steady state by linearizing the di¤erential system given by (DYN). Since we have an interior equilibrium, we can set all the multipliers to zero and obtain the following system:

0 B B B B @ _ h _ w s _ p 1 _ p 2 1 C C C C A = 2 6 6 6 6 6 6 6 4 @ 2 wg ;h C F " @ 2 wg ;wg C @ 2 wg ;wg C F " @ 2 wg ;wg C 1 F " @ 2 wg ;wg C 0 0 0 0 1 " @ 2 wg ;h C 2 F " @ 2 wg ;wg C + @ 2 h;h C F " @ 2 wg ;h C F " @ 2 wg ;wg C + @ 2 wg ;h C F " @ 2 wg ;wg C 0 F " @ 2 wg ;h C F " @ 2 wg ;wg C p 2 " + F " @ 2 wg ;wg C F " @ 2 wg ;wg C @ 2 wg ;wg C F " @ 2 wg ;wg C + 0 3 7 7 7 7 7 7 7 5 0 B B B B @ h h w s w s p 1 p 1 p 2 p 2 1 C C C C A
Following [START_REF] Dockner | On the Optimality of Limit Cycles in Dynamic Economic Systems[END_REF], we can show that:

Proposition 3. This steady state is locally a saddle point

Concluding remarks

This paper extends the standard groundwater extraction model by introducing the opportunity to access a new source of water (rainwater) through investment in capital.

Rainwater harvesting requires investment to build and maintain a storage and irrigation capacity that can be used jointly with groundwater. The main conclusion of this model is that the level of the aquifer, at the steady state, will be lower in the presence of this RWH capacity. We derive this result from the Hamiltonian necessary conditions of an intertemporal water management program by assuming, in a …rst step, the existence of a steady state. In a second step, we prove that these conditions are also su¢ cient, provided that a standard tranversality condition holds, and we verify that all these trajectories have a unique steady state which is locally a saddle point.

However, to isolate the e¤ect of rainwater harvesting on groundwater extraction as well as on the level of the aquifer depth, we consider the simplest possible dynamic setting with (i) a simple "bathtub"aquifer, i.e. a ‡at bottom with parallel sides, (ii) the social planner approach (iii) complete information on hydrological characteristics (iv) no uncertainty on capital.

These simpli…cations call for future extensions. Namely, in line with the literature relaxing some Gisser-Sanchez assumptions, it could be interesting to incorporate a more accurate depiction of groundwater hydrology and rainwater variability. For instance, Brozovic et al. (2003) or [START_REF] Saak | Groundwater Use under Incomplete Information[END_REF] integrate spatially variable features such as the speed of lateral ‡ow or di¤erences in the elevation of the bottom of the aquifer.

Thus, in our setting, one might expect the introduction of a two-cell aquifer where the elevation of bottom di¤ers across location to impact on our result through the trade-o¤ based on marginal costs.

Along the same lines, one could also consider more complex dynamics including, for instance, return ‡ows from irrigation to the aquifer, a more complex in…ltration process and water evaporation during storage. This would surely mitigate our clear-cut result that the water table is always depleted. For instance, [START_REF] Stahn | A drop of Rainwater against a drop of Groundwater : Does Rainwater Harvesting allows us really to spare groundwater? GREQAM Working Paper[END_REF] show that there exist certain hydrological characteristics under which the impact of RWH is not negative.

But these results are obtained in a context where RWH is a backstop technology requiring no capital accumulation.

Another more detailed characterization could lead us to incorporate uncertainties about rainfall variability, following [START_REF] Fisher | Adjusting to Climate Change : Implications of Increased Variability and Asymmetric Adjustment Costs for Investment in Water Reserves[END_REF], for instance, who model water resources as a stochastic process and focus on the determination of long-term water storage capacity. Failure to take uncertainty into account can lead to costly errors. In other words, by reckoning random capital in order to capture ‡uctuations in precipitations, one might expect the water level of the aquifer in the steady state to depend on risk behavior as well as the level of uncertainty.

Thus, various re…nements can be applied to this model to achieve a more detailed approach. Nevertheless, the results should be fairly similar to our main …nding about the impact on the aquifer water level. It is therefore interesting to wonder about the signi…cance of this result with respect to the principle of sustainable development. Groundwater also maintains the health of the ecosystem, which gives it a conservation value. In other words, the question must be addressed of whether the implementation of this technology might not jeopardize the sustainable level of groundwater for all its di¤erent functions.

APPENDIX

Appendix A. Proof of Lemma 1

(i) The quali…cation constraint.

Following [START_REF] Grass | Optimal Control of Non linear Processes with Applications to Drugs, Corruption and Terror[END_REF] th. 3.60, we have to verify that:

Q 1 = @ wg;I (w g ); w g and Q 2 = 2 6 6 6 6 6 6 4 @ !g;I @ h (h) _ h + @ wr (h) _ w r h 0 0 0 @ !g;I @ h h h _ h + @ wr h h _ w r 0 h h 0 0 @ !g;I @ h (w r ) _ h + @ wr (w r ) _ w r 0 0 w r 0 @ !g;I @ h (R w r ) _ h + @ wr (R w r ) _ w r 0 0 0 R w r 3 7 7 7 7 7 7 5
are matrices of full rank. Moreover, it is immediate that these matrices are given by:

Q 1 = h 1 0 w g i , Q 2 = 2 6 6 6 6 4 1 0 d 0 0 0 1 0 0 h h 0 0 0 1 0 0 w r 0 0 1 0 0 0 R w r 3 7 7 7 7 5
Q 1 is obviously of full rank. Concerning Q 2 , let us remember that h; R > 0. This means that we can always choose a non-zero vector when we respectively consider the columns 3; 4 and 5; 6. If we add to this choice the 2 …rst columns we can conclude that Q 2 is of rank 4.

(ii) The concavity of the Hamiltonian with respect to the controls.

For de…nition 1, a simple computation shows that: @ 2 (wg;I) H = " F "(w g + w s ) @ 2 wg;wg C(w g ; h) 0 0 0 (I)

#

Moreover, under A1(ii) and A2(ii), @ 2 (wg;I) H is negative de…nite:

Appendix B. Proof of lemma 3

It remains for us to verify that the Hamiltonian H w g ; I; h; w r ; (p i ) 2 i=1 is strictly concave in (w g ; I; h; w r ) and that the di¤erent constraints are quasi-concave in these variables. The last condition is always sat-is…ed since our constraints are linear. So let us now compute the Hessian of this Hamiltonian for each

(p i ) 2 i=1
. We obtain by taking the following order of the variables (w g ; w r ; h; I)

@ 2 H = 2 6 6 6 6 4 F " @ 2 wg;wg C F " @ 2 wg;h C 0 F " F " p 2 " 0 0 @ 2 wg;d C 0 @ 2 d;d C 0 0 0 0 " 3 7 7 7 7 5
Bearing in mind that 0 (I) = p 2 , we observe, under assumption A1(ii), A2(i)(ii) and A3(i)(ii), that:

D 1 = F " @ 2 wg;wg C < 0 D 2 = F " @ 2 wg;wg C F "
F " F " 0 " = F " @ 2 wg;wg C + 0 " + 0 "@ 2 wg;wg C > 0

D 3 = F " @ 2 wg;wg C F " @ 2 wg;d C F " F " 0 " 0 @ 2 wg;h C 0 @ 2 h;h C = F " 0 "@ 2 h;h C + (F " 0 ") @ 2 wg;wg C@ 2 h;h C @ 2 wg;h C 2 < 0
and

D 4 = @ 2 H 2 = "D 3 > 0
It follows that the Hamiltonian is strictly concave for all (p i )

2 i=1 :
Appendix C. Proof of proposition 2

Let : [0; R] 0; h ! R 2 be de…ned by:

(w g ; h) = 1 (w g ; h) 2 (w g ; h) ! = F 0 (R) @ wg C(w g ; h) + 1 @ h C(w g ; h) ( (R ! g )) ( + 0 (R ! g )) @ wg C(w g ; h)
! and let us apply a standard degree theory argument (see for instance Mas-Colell, 1985, p.207-08). For that purpose, let us introduce G (w g

; h) = h A B w g ! with (A; B) 2 0; h ]0; R[. Obviously, G (w g ; h) = 0
admits a unique solution given by (w g ; h) = (B; A) and det @ (w;h) G (B;A) = 0 1

1 0 = 1, hence deg(G) = 1.
Let us now introduce the following homotopy H (w g ; h; ) = G (w g ; h) + (1 ) (w g ; h) ; with 2 [0; 1]. A simple computation shows that:

8 > 0; rank @ (A;B) H (w g ; h; ; (A; B)) = rank " 0 0 #! = 2
It follows by the generic transversality theorem (see Mas-colell, 1985, I22, p.45) that for almost all (A; B) 2 0; h ]0; R[, and > 0, @ (wg;h; ) H (w g ; h; ; A; B) is of full rank. So let us …x one of them. If we want to make sure that H is a regular homotopy, we must verify that for = 0, @H (w g ; h; ) is also of full rank. This is a simple exercise of computation, since:

det @ wg;h H (w g ; h; 0) = det (@ ) = " 1 @ 2 h;wg C @ 2 wg;wg C 1 @ 2 h;h C @ 2 wg;h C dD dwg @ 2 wg;wg C @ 2 wg;h C # = 1 @ 2 wg;wg C @ 2 h;h C @ 2 h;wg C 2 dD dw g 1 @ 2 h;h C @ 2 wg;h C with D(w g ) = ( (R ! g )) ( + 0 (R ! g )).
We deduce from our assumptions on the cost function (A2(i)(ii)) and the investment function (A3) that det (@ ) > 0: it follows that @H (w g ; h; 0) is of full rank and more generally that the H is a regular homotopy.

It follows from the classi…cation theorem of 1-manifold that H 1 (0) is di¤eomorphic to segments or circles. Since "a segment starts at the regular solution of G (w g ; h) = 0", it "reaches a solution of (w g ; h) = 0" if we verify that there exists no sequence w n g ; h n ; n 2 H 1 (0) and w n g ; h n ; n ! w 0 g ; h 0 ; 0 2 f0; Rg 0; h

[0; 1] (C.1)
If this condition is satis…ed, the existence part of our proposition is proved. But we can go a step further.

Since det (@ ) > 0, i.e is a regular function, 1 (0) only contains isolated points. It is moreover a closed subset of a compact set given by [0; R] 0; h , hence a …nite set. By using a standard index theory argument (Mas-Colell 1985, J2 p46), we can say that P (wg;h)2 1 (0) ind (det (@ )) = 1 where ind(x) = 1 if x > 0 and ind(x) = 1 if x < 0. But the sign of det (@ ) is invariant. We can therefore also conclude that the solution to (w g ; h) not only exists but is also unique.

To conclude the proof it therefore remains for us to verify that the boundary condition given in equation (C.1) is true. To do so, let us …rst observe that: Assumption A4(i) and (ii) and the fact that @ wg 1 (w g ; h) < 0 bring us to the conclusion that:

8 > > > > < > > > > : 8w g 2 [0; R] ,
1 (w g ; 0) = F 0 (R) @ wg C(w g ; 0) + 1 @ h C(w g ; 0) < F 0 (R) @ wg C(0; 0) + 1 @ h C(0; 0) < F 0 (R) @ wg C(0; 0) < 0 8w g 2 [0; R] , 1 (w g ; h) = F 0 (R) @ wg C(w g ; h) + 1 @ h C(w g ; h) > F 0 (R) @ wg C(R; h) + 1 @ h C(R; h) > 0 By A4(ii) and the fact that @ d 2 (w g ; d) < 0, we can say that: Hence @ w n g ; d n ; n 2 H 1 (0) with w n g ! w 0 g 2 f0; Rg. We can therefore say that the boundary condition given in equation (C.1) holds.

Appendix D. Proof of proposition 3

We know from [START_REF] Dockner | On the Optimality of Limit Cycles in Dynamic Economic Systems[END_REF] that this system exhibits saddle-point stability with either real or conjugate complex root if (i) the determinant D of this system is positive and (ii) K = M 2 r 2 < 0 where M 2 stands for the sum of the second order principal minors of this system Let us …rst verify that K > 0. By computation, we have: K = @ 2 wg ;h C F " @ 2 wg ;wg C 1 F " @ 2 wg ;wg C @ 2 wg ;h C 2 F " @ 2 wg ;wg C + @ 2 h;h C + @ 2 wg ;h C F " @ 2 wg ;wg C + 0 1 " p 2 " + F " @ 2 wg ;wg C F " @ 2 wg ;wg C + 0 + 2 @ 2 wg ;wg C F " @ 2 wg ;wg C 0 F " @ 2 wg ;h CF " F " @ 2 wg ;wg C 0 = @ 2 wg;h C F " @ 2 wg;wg C + @ 2 h;h C F " @ 2 wg;wg C 0 ( + 0 ) p 2 F " @ 2 wg;wg C " F " @ 2 wg;wg C < 0

It remains fur us to verify that the determinant D of this linearized system is positive. If we develop this expression with respect to the second line we can say that D = D 1 + D 2 with:

D 1 = 0 ( + 0 )
@ 2 wg ;h C F " @ 2 wg ;wg C 1 F " @ 2 wg ;wg C @ 2 wg ;h C 2 F " @ 2 wg ;wg C + @ 2 h;h C + @ 2 wg ;h C F " @ 2 wg ;wg C = 0 ( + 0 )

F " @ 2 wg;wg C h @ 2 h;h C @ 2 wg;h C i > 0

As the second term D 2 is more tedious to compute, we skip this computation and simply gives the result:

D 2 = 1 " 2 6 6 6 6 6 4
@ 2 wg ;h C F " @ 2 wg ;wg C @ 2 wg ;wg C F " @ 2 wg ;wg C 1 F " @ 2 wg ;wg C @ 2 wg ;h C 2 F " @ 2 wg ;wg C + @ 2 h;h C F " @ 2 wg ;h C F " @ 2 wg ;wg C + @ 2 wg ;h C F " @ 2 wg ;wg C F " @ 2 wg ;h C F " @ 2 wg ;wg C p 2 " + F " @ 2 wg ;wg C F " @ 2 wg ;wg C @ 2 wg ;wg C F " @ 2 wg ;wg C 3 7 7 7 7 7 5 = p 2 @ 2 wg;h C @ 2 h;h C F " @ 2 wg;wg C + @ 2 wg;h C 2 @ 2 wg;wg C @ 2 h;h C " F " @ 2 wg;wg C > 0

  8d 2 [0; R] , 2 (0; d) = ( (R)) ( + 0 (R)) @ wg C(0; d) > 0 Finally, A3 and A2(ii) lead to:8d 2 [0; R] , 2 (R; d) = @ wg C(R; d) < 0 It follows, from the …rst observation, that 8 (w g ; ) 2 [0; R] [0; 1] ( H 1 (w g ; 0; ) = d 1 + (1 ) 1 (w g ; 0) < 0 H 1 (w g ; d; ) = d d 1 + (1 ) 1 (w g ; d) > 0It is therefore impossible that there exists a sequence w n g ; d n ; n 2 H 1 (0) with d n ! d 0 2 0; d . Let us now move to the the second and the third observations. We respectively conclude that 8 (d; ) 2 d; ) = w 1 + (1) 2 (0; d) > 0 H 2 (R; d; ) = (w 1 R) + (1) 2 (R; d) < 0

See http://www.unep.or.jp/ietc/publications/urban/urbanenv-2/9.asp)

Obviously, the amount of water reaching the ground depends on various factors. One usually identi…es some geological parameters like soil moisture and porosity as well as some climatic factors like evaporation or the intensity of precipitation. All these parameters de…ne a speci…c hydraulic system characterized by runo¤s, return ‡ows and evaporation. However, we ignore these characteristics in order to focus on the main argument. These extensions are discussed in our concluding remarks.

In order to simplify the notation, we omit the fact that all the arguments are functions of t whenever possible..

For these necessary conditions the reader is referred to[START_REF] Seierstad | Chap.6 : Mixed and Pure State Constraints[END_REF] theorem 9 p.381 and note 6 p.374, Grass et al. (2008) theorem 3.60 p.149 or Hartl et al. (1995) Theorem 4.1 p186.

The reader may be surprised by the …rst and the last slackness conditions. But the …rst follows from the non-negativity of w g while the last is related to the fact that w s R (remember that w g + w s = R).