MD Lucie Peigné

MD Florence Godey

PhD Mathieu Le Gallo

François Le Gall

PhD Alain Fautrel

Jeff Morcet

MD Franck Jégoux
email: franck.jegoux@chu-rennes.fr

One-step nucleic acid amplification for detecting lymph node metastasis of head and neck squamous cell carcinoma

Keywords: Lymph node, OSNA, head and neck squamous cell carcinoma, pathological analysis, CK19, sentinel lymph node biopsy, micro metastasis

BACKGROUND:

In head and neck squamous cell carcinoma (HNSCC) 30% of cN0 patients have occult metastasis. LN invasion is a major prognostic factor. Sentinel lymph node (SLN) is an option for cN0 neck management. One-step nucleic acid amplification (OSNA) used to analyze SLN in breast cancer is also a candidate to get more reliable intraoperative HNSCC lymph node (LN) staging.

OBJECTIVE:

To compare OSNA analysis to pathological analysis in cN0 HNSCC.

MATERIALS AND METHODS:

157 LN from 26 cN0 HNSCC patients were prospectively analyzed (6.3LN/patient). Exclusion criteria were previous surgery or radiotherapy. Each node was cut into 4 equal pieces alternatively sent to pathological analysis and OSNA technique. IHC CK19 was performed on the primary tumor biopsy and RT-qPCR of CK19, PVA and EPCAM on the LN lysate of discordant cases.

RESULTS:

OSNA was able to provide intraoperative result in all patients. OSNA detected 21 metastases. There were 139 concordant LN (88.5%). There were 18 initial discordant LN (11.5%), 13 (8.3%) were OSNA positive/pathological analysis negative, 5 (3.2%) were OSNA negative/pathological analysis positive. After elimination of allocation bias, false negative rate was 1.3%, sensitivity and specificity were 90% and 95.6%, PPV and NPV were 75% and 98.5%.

CONCLUSION:

Our results suggest that OSNA should be considered to improve SNB analysis both for increasing micro metastasis diagnosis and offer extemporaneous results.

Study registered under clinicaltrials.gov database number NCT02852343.

INTRODUCTION

In head and neck squamous cell carcinoma (HNSCC) 30% of cN0 patients have occult lymph node (LN) metastasis [START_REF] Barrera | Detection of Occult Cervical Micrometastases in Patients With Head and Neck Squamous Cell Cancer[END_REF][START_REF] Gurney | Implications of a positive sentinel node in oral squamous cell carcinoma[END_REF]. Lymph node metastasis reduce the survival by up to 50% [START_REF] Tao | Detection of occult carcinomatous diffusion in lymph nodes from head and neck squamous cell carcinoma using real-time RT-PCR detection of cytokeratin 19 mRNA[END_REF][START_REF] Alkureishi | Joint Practice Guidelines for Radionuclide Lymphoscintigraphy for Sentinel Node Localization in Oral/Oropharyngeal Squamous Cell Carcinoma[END_REF]. The diagnosis of occult LN metastasis is relevant since metastatic patients may benefit from adjuvant therapies whereas observation may be sufficient in pN0. Actual treatment of LN is systematic neck dissection to treat all metastatic patients and to get accurate LN staging. Systematic neck dissection may increase morbidity in 70% of the patients [START_REF] Cruz | Elective versus Therapeutic Neck Dissection in Node-Negative Oral Cancer[END_REF]. To limit morbidity and to improve pathological analysis focusing it on a small number of nodes, Sentinel Lymph Node (SLN) analysis has been developed for cN0 neck management [START_REF] Stoeckli | Sentinel node biopsy for early oral carcinoma[END_REF][START_REF] Broglie | Long-Term Experience in Sentinel Node Biopsy for Early Oral and Oropharyngeal Squamous Cell Carcinoma[END_REF] . If the SLN biopsy is negative, neck dissection can be avoided, reducing the risk for post-operative complications such as nerve injuries, cervical hematomas and post-operative infections. However, intra operative SLN pathological analysis by frozen section or imprint cytology doesn't provide definitive result. After completion of pathological analysis, post operatively secondtime surgery is required in patient with metastatic SLN. A technique giving definitive and accurate intra-operative results is needed. One step nucleic acid amplification (OSNA) assay allows rapid amplification of CK19 mARN which is a biomarker of epithelial and adenocarcinoma cells absent in healthy LN. OSNA assay uses RT-LAMP (Loop isothermal Amplification) method. The main advantages of this technique are to give a definitive intra operative result within 30min thus allowing concurrent neck dissection and avoiding second time surgery. This technique has been approved for intraoperative SLN analysis in breast cancers [START_REF] Tsujimoto | One-step Nucleic Acid Amplification for Intraoperative Detection of Lymph Node Metastasis in Breast Cancer Patients[END_REF] and is widely used in SLN analysis of colon cancer in Europe and Japan [START_REF] Wild | Is it time for onestep nucleic acid amplification (OSNA) in colorectal cancer? A systematic review and meta-analysis[END_REF]. Three main studies investigated the feasibility of OSNA assay in HNSCC [START_REF] Goda | One-step nucleic acid amplification for detecting lymph node metastasis of head and neck squamous cell carcinoma[END_REF][START_REF] Suzuki | Diagnostic potential of 1-step nucleic acid amplification assay in patients with head and neck squamous cell carcinoma based on CK19 expression in a primary lesion[END_REF][START_REF] Matsuzuka | Intraoperative Molecular Assessment for Lymph Node Metastasis in Head and Neck Squamous Cell Carcinoma Using One-Step Nucleic Acid Amplification (OSNA) Assay[END_REF] whatever the neck LN invasion was. As OSNA assay is a candidate for intra operative SLN analysis in HNSCC the main objective of our study was to compare OSNA assay to standard pathological analysis in cN0 HNSCC patients.

MATERIALS AND METHODS

Study design

A monocentric, prospective non-interventional study of feasibility was conducted in our university tertiary care institution. This study of non-inferiority was comparative; each LN was its own control, two quarters of the LN was analyzed by pathological analysis and the two quarters left by OSNA assay.

Inclusion criteria were cN0 HNSCC patients over 18 years old requiring neck dissection or SLN. Noninclusion criteria were previous neck dissection or radiotherapy. Patients unable to consent weren't enrolled. Characteristics of the 26 patients included were summarized in Table 1. Mean age was 61.1 ±10.0 years old (range 29-79 years old) and male sex was predominant (80.8%). 158 LN were collected (6.3LN/patients), 1 LN wasn't analyzed due to technical issue. Institutional ethical committee approval and written informed consent from all patients was obtained. This study was recorded in clinicaltrials.gov database under the number NCT02852343.

Lymph nodes processing

From February 2015 to December 2016, lymph nodes were harvested during surgery. While patients underwent neck dissection, LN were harvested in a limited time of 15 minutes to prevent mRNA degradation. After removing extra nodal tissue and lipid, LN harvested from the neck dissection were cut into 4 equal pieces. The 4 pieces were alternatively sent to pathological examination and OSNA technique. The 2 pieces dedicated to OSNA technique were immediately snap-frozen in liquid nitrogen (Figure 1).

For pathological examination, LN samples were embedded in paraffin blocks and 5µm-sections were cut every 200µm. Each level was stained with HES. HES results were express according to Hermanek's classification [START_REF] Hermanek | Classification of isolated tumor cells and micrometastasis[END_REF] : micro metastasis ≤2mm and macro metastasis >2mm.

For OSNA assay LN samples were homogenized with Lynorhag Lysis buffer (Sysmex®, Kobe, Japan) and shortly centrifuged at 10.000xg at room temperature according to the manufacturer's instructions (Sysmex®, Kobe, Japan). OSNA analysis was carried out in duplicate with pure and diluted sample (1/10) of LN lysate without prior isolation and purification of mRNA. A 2µL sample of each lysate was subjected to a RT-LAMP reaction. Six CK19 primers were used for the reaction (5'-GGA GTT CTC AAT GGT GGC ACC AAC TAC TAC ACG ACC ATC CA-3' (CK19FA), 5'-GTC CTG CAG ATC GAC AAC GCC TCC GTC TC AAA CTT GGT TCG-3' (CK19RA), 5'-TGG TAC CAG AAG CAG GGG-3' (CK19F3), 5'-GTT GAT GTC GGC CTC CAC G-3' (CK19R3), 5'-AGA ATC TTG TCC CGC AGG-3' (CK19LPF) and 5'-CGT CTG GCT GCA GAT GA-3' (CK19LPR)) to avoid pseudo gene detection. After 16min amplification time using Lynoamp BC gene amplification reagent (Sysmex®, Kobe, Japan), the CK19 mRNA copy number per µL of lysate was determined by real-time monitoring of turbidity caused by increasing levels of magnesium pyrophosphate (a byproduct of amplification) with the RD-100i device (Sysmex®, Kobe, Japan). The node status was assessed as previously validated for breast cancer: copy number <250 = no metastasis, 250 ≤ copy number = metastasis [START_REF] Tsujimoto | One-step Nucleic Acid Amplification for Intraoperative Detection of Lymph Node Metastasis in Breast Cancer Patients[END_REF]. Leftover lysates are used for RNA extraction, and other molecular techniques.

Tumor biopsy analysis

To study CK19 protein expression all primary tumors were stained with IHC CK19 (Agilent Dako®, Santa Clara, USA). CK19 tumoral expression was considered as positive if there were more or equal to 1% of the cells stained [START_REF] Suzuki | Diagnostic potential of 1-step nucleic acid amplification assay in patients with head and neck squamous cell carcinoma based on CK19 expression in a primary lesion[END_REF]. Paraffin-embedded tissue was cut at 4 µm, mounted on positively charged slides and dried at 58°C for 60 minutes. Immunohistochemical staining was performed on the Discovery XT Automated IHC stainer using the Ventana® detection kit (Ventana® Medical Systems, Tucson, Arizona).

Following deparaffination with Discovery wash solution (Ventana® Medical Systems, Tucson, Arizona) at 75 °C for 8 minutes, antigen retrieval was performed using Tris EDTA-based buffer solution pH8, at 95°C to 100°C for 24 minutes. Endogen peroxidase was blocked 3% H2O2 for 10 minutes. After rinsing, slides were incubated at 37°C for 32minutes with primary antibody: mouse anti-CK19 (clone RCK108, dilution at 1/50). Signal enhancement was performed using a goat anti-mouse HRP. Secondary antibody was incubated at 37°C for 16 minutes using OmniMap Kit (Ventana® Medical Systems, Tucson, Arizona). Slides were then counterstained for 16 minutes with hematoxylin and rinsed. Slides were manually dehydrated and cover slipped.

Complementary analysis:

For all discordant samples complementary molecular analysis were done with OSNA lysate: RT-PCR for CK19 (target of OSNA Assay) and also two other makers EPCAM and PVA. PVA and EPCAM were chosen because Ferris et al. [START_REF] Ferris | Molecular Staging of Cervical Lymph Nodes in Squamous Cell Carcinoma of the Head and Neck[END_REF][START_REF] Ferris | Intraoperative qRT-PCR for Detection of Lymph Node Metastasis in Head and Neck Cancer[END_REF][START_REF] Trivedi | Correlation of tumor marker expression with nodal disease burden in metastatic head and neck cancer[END_REF] showed that PVA and EPCAM were present in HNSCC metastatic LN with a high accuracy.

Tissues and pathological examination

A second reading of all HES slides from all discordant LN in this study was realized by an experimented head and neck pathologist aware of the discordance of HES staining with OSNA assay. All discordant LN were stained with IHC CKAE1/3 with the same protocol described for CK19 above.

RNA isolation, cDNA synthesis and quantitative PCR

RNA was isolated from lymph nodes lysates used for OSNA assay (method described above) using NucleoSpin®RNA kit (Macherey-Nagel®, Düren, Germany) essentially as described by the manufacturer. Total RNA was quantified by spectrophotometry (260/280 nm ratio).

Reverse transcription (RT) was performed on 500ng of RNA, using random hexamer primers, using the Superscript III Reverse transcription kit (Thermo-fisher®, Massachusetts, USA) according to the manufacturer's instructions.

Each 5 µL PCR reaction contained 2,5µL TaqMan® Fast Advance Master Mix (2X) (Applied Biosystem®, Foster City, USA), 0,25µL of TaqMan® Assay (20X), 1,75µL of nuclease free water and 2µL of 1/10 diluted cDNA. The amplification program consisted of one stage with an initial 50°C activation stage for 2 minutes followed by one 95°C polymerase activation stage for 2 minutes and 40 stages of 95°C denaturation for 1 second followed by one 60°C extend stage for 20 seconds.

All quantitative PCR was performed on the QuantStudio® 5 system (Thermofisher®, Massachusetts, USA). Relative expression of the marker genes was calculated using the ΔCT methods and with βglucuronidase as the endogenous control gene. All assays were designed for use with 5' nuclease hybridization probes.

The references of probes for β-glucuronidase, CK19, EPCAM and PVA (purchased from Thermofisher®, Massachusetts, USA) are listed in Supplementary table 1.

Data from the PCR melting curve were analyzed using the applied biosystems Relative Quantitation software on the Thermo Fisher Cloud. The first derivative of the melting curve was used to determine the product Tm and to establish the presence of the specific product in each sample. Samples were analyzed in triplicate PCR reactions and the average Ct value was used in the calculation of relative expression to obtain the highest value of background expression for the sample.

To validate the RT-qPCR method in our study, the correlation between PVA, EPCAM and CK19 was calculated (Figure 2).

Definition of final discordant rate

Final discordant rate was defined as the number of discordant nodes after exclusion of the allocation bias, meaning the node in which complementary analysis confirmed that metastases were strictly located in either the tissue used for OSNA or the tissue used for histology which renders a comparison of the 2 methods impossible.

Statistical analysis

Qualitative variables were expressed as numbers and percentages. Quantitative data were expressed as means ± standard deviations (SD) if normally distributed and medians (Q1-Q3) if not normally distributed and compared with matched Student t-tests and McNemar tests. HES was considered as the objective standard and those of the OSNA assay as the diagnostic test value. P values <0.05 were considered statistically significant. The sensitivity, specificity, positive and negative predictive values were calculated. Statistical analyses were performed using SAS® software version 9.4 (SAS Institute, Cary, NC, USA).

RESULTS

Comparison of OSNA assay with pathological analysis

OSNA was successful in 157 LN in a median time of 30'. With OSNA assay 136 LN were negative. OSNA analysis detected 21 metastases, representing 13.4% of the LN. OSNA number of copies range from <250 to 670 000 copies/µL of CK19 mRNA (Table 3). Average number of copies/µL was 48 702 ± 156 representing 8.3% of the LN. There was no significant difference regarding to the rate of metastasis detection between the two techniques (p=0.14). There were 139 (88.5%) concordant LN in which 131 (83.4%) were negative and 8 (5.1%) were positive for both OSNA and HES (Table 2). There were 18 initially discordant LN (11.5%): 13 OSNA+/HES-(8.3 %) and 5 OSNA-/HES+ (3.2%) (Table 2).

Analysis of primary tumor CK19 protein and mRNA expression

On the primary tumor, CK19 IHC was positive in 14/26 cases (53.8%). Four patients had an IHC CK19 staining at 100%, 7 had 75% staining, one had 50% staining and 2 had 25% staining (Table 4). The 12 remaining patients did not express CK19 at protein level on the primary tumor. Number of discordant cases was not influenced by IHC CK19 (p=0.12).

Identification of allocation bias and analysis of discordant cases

The second reading by an experimented pathologist did not change the results of the first reading. IHC CKAE1/3 was negative in all discordant LN except for patient 8.

Among the 18 nodes initially discordant, six were associated to concordant positive nodes (patients 8, 9, 24 and 26, Table 4) in adjacent levels leading to consider them as allocation bias (Figure 3). Among the 12 remaining LN, two were OSNA-/HES+ and 10 were OSNA+/HES-. The 2 OSNA-/HES+ were associated to a non-expression of CK19 on primary tumor biopsy, a negative CK19 RT-qPCR and a positive PVA and EPCAM RT-qPCR on the LN lysate confirming the metastatic status of the LN (patient 15, Table 4). These 2 LN were considered as false negative for OSNA analysis (Figure 3). Among the ten OSNA+/HES-LN, four LN had positive RT-qPCR detection of CK19, PVA or EPCAM, leading to consider them as true positive for OSNA analysis and allocation bias in favor to OSNA (patients 2, 5, 7, Table 4). RT-qPCR was negative for the remaining 6 OSNA+/HES-LN (patients 12, 13, 16, 21 and 22, Table 4) leading to consider them as OSNA false positive LN (Figure 3).

After identification of allocation bias the final discordant rate was 5.1% (n=8). The OSNA false negative rate was 1.3% (2 LN) and the sensitivity and specificity were respectively 90% and 95.6%. The positive predictive value and the negative predictive value were respectively 75% and 98.5%.

DISCUSSION

In our study, OSNA analysis was technically feasible in all specimen and allowed micro metastasis diagnosis with a low false negative rate of 1.3%.

Allocation bias are the main pitfall of histological tissue comparison studies. They are related to the necessity to split the node which lead to an isolation of micro metastasis that become strictly located in either the tissue used for OSNA or the tissue used for histology. This renders a comparison of the two methods impossible. In our study overall discordant rate before exclusion of tissue allocation bias was 11.5% with a majority of OSNA+/HES-discordant LN. Thus, real discordant samples were distinguished from allocation bias using a second reading of all HES slides, IHC CKAE1/3 staining and complementary molecular analysis. This last was based on EPCAM and PVA markers as Ferris et al. [START_REF] Ferris | Molecular Staging of Cervical Lymph Nodes in Squamous Cell Carcinoma of the Head and Neck[END_REF][START_REF] Ferris | Intraoperative qRT-PCR for Detection of Lymph Node Metastasis in Head and Neck Cancer[END_REF][START_REF] Trivedi | Correlation of tumor marker expression with nodal disease burden in metastatic head and neck cancer[END_REF] showed that among multiple markers these ones were present in HNSCC metastatic LN with a high accuracy. After these analyses, the final statistical analysis, only eight LN remain discordant. There were 6 LN with positive OSNA results despite negative HES analysis and negative RT-qPCR. OSNA assay is able to detect one tumoral cell among 10 6 non tumoral cells [START_REF] Ghossein | Molecular detection and characterisation of circulating tumour cells and micrometastases in solid tumours[END_REF] with more sensitive results when compared to HES or RT-qPCR [START_REF] Shaw | Intraoperative Sentinel Lymph Node Evaluation: Implications of Cytokeratin 19 Expression for the Adoption of OSNA in Oral Squamous Cell Carcinoma[END_REF] . False positive results for OSNA analysis could be due to the detection of heterotopic salivary, thyroid or epithelial inclusion in LN [START_REF] Shinohara | Heterotopic salivary gland tissue in lymph nodes of the cervical region[END_REF]. These inclusions account for only 1.6% of the LN [START_REF] Hamakawa | Keratin mRNA for detecting micrometastasis in cervical lymph nodes of oral cancer[END_REF] and are mainly located in the parotid gland or in level II and could express low level of CK19 potentially detected by OSNA. Actually, these inclusions are probably not detected as positive by OSNA regarding to the choosing threshold. Goda et al. studied OSNA results in salivary gland tissue and found that CK19 mRNA copies with OSNA ranged from 180 to 200 copies/µL [START_REF] Goda | One-step nucleic acid amplification for detecting lymph node metastasis of head and neck squamous cell carcinoma[END_REF] which is below the positive cut off (250copies/µL). Another cause for false positive is contamination by epithelial cells from skin. For this reason, as managed in the present study, LN should be carefully dissected from all surrounding non-lymphatic tissue with clean tools free from any prior use to avoid false positive results. Thus, it is likely that positive OSNA/negative HES LN suggest a superiority of OSNA on HES.

The overall false negative rate after identification of allocation bias was only 1.3% which is consistent with the literature. Goda et al. [START_REF] Goda | One-step nucleic acid amplification for detecting lymph node metastasis of head and neck squamous cell carcinoma[END_REF], Matsuzuka et al. [START_REF] Matsuzuka | Intraoperative Molecular Assessment for Lymph Node Metastasis in Head and Neck Squamous Cell Carcinoma Using One-Step Nucleic Acid Amplification (OSNA) Assay[END_REF] and Suzuki et al. [START_REF] Suzuki | Diagnostic potential of 1-step nucleic acid amplification assay in patients with head and neck squamous cell carcinoma based on CK19 expression in a primary lesion[END_REF] found respectively 2.6%, 3.4% and 3.7% of FN. We found a sensitivity and specificity similar to the literature, respectively of 90% and 95.6%. Goda et al. [START_REF] Goda | One-step nucleic acid amplification for detecting lymph node metastasis of head and neck squamous cell carcinoma[END_REF] found a sensitivity of 92% and a specificity of 95%. Matsuzuka et al. [START_REF] Matsuzuka | Intraoperative Molecular Assessment for Lymph Node Metastasis in Head and Neck Squamous Cell Carcinoma Using One-Step Nucleic Acid Amplification (OSNA) Assay[END_REF] found a sensitivity of 82.4% and a specificity of 99.3%. Suzuki et al. found a sensitivity and specificity of OSNA assay respectively at 86% and 100% in tumor expressing CK19 protein [START_REF] Suzuki | Diagnostic potential of 1-step nucleic acid amplification assay in patients with head and neck squamous cell carcinoma based on CK19 expression in a primary lesion[END_REF]. The studies of Goda [START_REF] Goda | One-step nucleic acid amplification for detecting lymph node metastasis of head and neck squamous cell carcinoma[END_REF], Matsuzuka [START_REF] Matsuzuka | Intraoperative Molecular Assessment for Lymph Node Metastasis in Head and Neck Squamous Cell Carcinoma Using One-Step Nucleic Acid Amplification (OSNA) Assay[END_REF] and Suzuki [START_REF] Suzuki | Diagnostic potential of 1-step nucleic acid amplification assay in patients with head and neck squamous cell carcinoma based on CK19 expression in a primary lesion[END_REF] were undertaken in a population with a variety of neck invasion and sites of HNSCC. None of these studies targeted a population fully eligible for SLN biopsy. As we included only cN0 nodes the low number histologically positive nodes could have impacted the sensitivity and specificity rates.

Most of the negative OSNA/positive HES were due to an allocation bias, only two of them were due to non-expression of CK19 mRNA. As found in our study, CK19 protein expression in HNSCC primary tumor range from 39% to 80% [START_REF] Suzuki | Diagnostic potential of 1-step nucleic acid amplification assay in patients with head and neck squamous cell carcinoma based on CK19 expression in a primary lesion[END_REF][START_REF] Shaw | Intraoperative Sentinel Lymph Node Evaluation: Implications of Cytokeratin 19 Expression for the Adoption of OSNA in Oral Squamous Cell Carcinoma[END_REF][START_REF] Noorlag | Cytokeratin 19 expression in early oral squamous cell carcinoma and their metastasis: Inadequate biomarker for one-step nucleic acid amplification implementation in sentinel lymph node biopsy procedure[END_REF]. Some authors suggested to restrain OSNA analysis to patient harboring a tumor positive for CK19 IHC [START_REF] Suzuki | Diagnostic potential of 1-step nucleic acid amplification assay in patients with head and neck squamous cell carcinoma based on CK19 expression in a primary lesion[END_REF][START_REF] Shaw | Intraoperative Sentinel Lymph Node Evaluation: Implications of Cytokeratin 19 Expression for the Adoption of OSNA in Oral Squamous Cell Carcinoma[END_REF] excluding tumor with no CK19 protein expression. However, we think that this selection criterion is not relevant because we found that CK19 mRNA can be detected both by OSNA and RT-qPCR even if the protein expression is not detected in the tumor burden as reported by previous studies [11][21] this was probably the case of patient 15.

Furthermore, Norlaag et al reported that 36 to 43% of CK19 negative tumor have CK19 positive metastatic LN for [START_REF] Noorlag | Cytokeratin 19 expression in early oral squamous cell carcinoma and their metastasis: Inadequate biomarker for one-step nucleic acid amplification implementation in sentinel lymph node biopsy procedure[END_REF]. This underline that CK19 protein expression in the primary tumor is not well correlated to the one of metastatic LN. Thus, we do think that when dealing with OSNA analysis the question of under expression of CK19 has to be address but not only with the CK19 protein detection in the tumor. In a cohort of 87 primary tumor biopsies, Shaw et al. analyzed the detection of CK19 using OSNA assay, RT-qPCR, RNA-ISH and IHC. They found that OSNA was the most sensitive technique with 80% of tumor expressing CK19.

To solve the issue of the variability of CK19 expression in HNSCC, molecular tests targeting other markers in the node analysis could be discuss. Yamauchi et al. also showed that CK19 had an insufficient expression in tongue cancer and proposed to use p63 [START_REF] Yamauchi | Quantitative expression study of four cytokeratins and p63 in squamous cell carcinoma of the tongue: suitability for sentinel node navigation surgery using one-step nucleic acid amplification[END_REF] . In a pioneer study, Brennan et al. proposed p53 as a biomolecular target [START_REF] Brennan | Molecular assessment of histopathological staging in squamous-cell carcinoma of the head and neck[END_REF], afterwards Hamakawa et al. suggested to use SCCA [START_REF] Hamakawa | Genetic diagnosis of micrometastasis based on SCC antigen mRNA in cervical lymph nodes of head and neck cancer[END_REF].

Other molecular targets were also studied [START_REF] Elsheikh | Importance of molecular analysis in detecting cervical lymph node metastasis in head and neck squamous cell carcinoma[END_REF]: MUC1 [START_REF] Cai | Detection mucin gene in cervical lymph nodes in supraglottic cancer by reverse transcription polymerase chain reaction[END_REF], E48 transcripts [START_REF] Nieuwenhuis | Assessment and clinical significance of micrometastases in lymph nodes of head and neck cancer patients detected by E48 (Ly-6D) quantitative reverse transcriptionpolymerase chain reaction[END_REF] , CK13 [START_REF] Hamakawa | Keratin mRNA for detecting micrometastasis in cervical lymph nodes of oral cancer[END_REF][START_REF] Onishi | Basic and clinical studies on quantitative analysis of lymph node micrometastasis in oral cancer[END_REF], CK14 [START_REF] Shores | Clinical Evaluation of a New Molecular Method for Detection of Micrometastases in Head and Neck Squamous Cell Carcinoma[END_REF], CK17 [START_REF] Solassol | Pemphigus vulgaris antigen mRNA quantification for the staging of sentinel lymph nodes in head and neck cancer[END_REF][START_REF] Garrel | The Diagnostic Accuracy of Reverse Transcription-PCR Quantification of Cytokeratin mRNA in the Detection of Sentinel Lymph Node Invasion in Oral and Oropharyngeal Squamous Cell Carcinoma: A Comparison with Immunohistochemistry[END_REF] , CK20 [START_REF] Elsheikh | Importance of molecular analysis in detecting cervical lymph node metastasis in head and neck squamous cell carcinoma[END_REF], PTHrp [START_REF] Ferris | Intraoperative qRT-PCR for Detection of Lymph Node Metastasis in Head and Neck Cancer[END_REF]. Recently, Ferris et al. screened 40 potential markers using RT-qPCR for their ability to detect HNSCC LN metastasis [START_REF] Ferris | Molecular Staging of Cervical Lymph Nodes in Squamous Cell Carcinoma of the Head and Neck[END_REF][START_REF] Trivedi | Correlation of tumor marker expression with nodal disease burden in metastatic head and neck cancer[END_REF]. As they identified PVA [START_REF] Solassol | Pemphigus vulgaris antigen mRNA quantification for the staging of sentinel lymph nodes in head and neck cancer[END_REF] and EPCAM as the most specific markers they developed a rapid RT-qPCR targeting them using an automated RNA isolation and RT-qPCR instrument (Cepheid GeneXpert®, Sunnyvale, USA) with an excellent reproducibility and linearity [START_REF] Ferris | Rapid molecular detection of metastatic head and neck squamous cell carcinoma as an intraoperative adjunct to sentinel lymph node biopsy[END_REF]. At that time « this homemade » molecular technology could not be used for clinical practice instead of OSNA that is already implemented in clinical practice for breast, gastric and colon cancer. A 14% false negative rate have been reported using SNB analysis in T1T2N0 in the SENT European Trial [START_REF] Schilling | Sentinel European Node Trial (SENT): 3-year results of sentinel node biopsy in oral cancer[END_REF] which renders the 1.3% false negative rate of OSNA acceptable in clinical practice.

It could be suggested to analyze extemporaneously the LN with OSNA and to keep a central histological section that could be analyzed post-operatively with IHC staining targeting more ubiquitous cytokeratin like CKAE1/3. The great advantage with this strategy is to analyze extemporaneously SLN with good specificity and sensitivity in order to complete SLN dissection or not. The central section kept for histological analysis would allow to confirm that the tissue analyzed with OSNA assay is indeed lymphatic tissue and to diagnose the rare macro metastasis in cases without CK19 expression.

CONCLUSION

OSNA was able to provide intraoperative result in all specimen and to diagnose micro metastasis in cN0 HNSCC with a low false negative rate. This suggests that OSNA should be considered to improve SNB analysis both for increasing micro metastasis diagnosis and offer extemporaneous results. Further studies in SLN conditions with morbidity, prognosis, and survival analysis to address the benefit for patients are warranted.

ROLE OF THE FUNDING SOURCE

The University hospital of Rennes gave us a funding in order to finance the analysis and the technical support needed for the study. This funding was free from conflicts of interest.

Table 1 : Patients characteristics

 1

	Gender
	Male	21 (80.8%)
	Female	5 (19.2%)
	Age (year)	61.1 (29-79)
	Tumor location
	Oral Cavity	15 (57.7%)
	Oropharynx	1 (3.8%)
	Hypopharynx	2 (7.7%)
	Larynx	8 (30.8%)
	Differentiation
	Well	18 (69.2%)
	Medium	8 (30.8%)
	Low	0 (0.00%)
	Stage
	T1N0M0	7 (26.9%)
	T2N0M0	11 (42.3%)
	T3N0M0	2 (7.7%)
	T4N0M0	6 (23.1%)

Table 2 : Concordance between HES and OSNA.

 2


Emmanuel Gobet
email: emmanuel.gobet@imag.fr

C Éline

SOLVING BSDE WITH ADAPTIVE CONTROL VARIATE

Keywords: Backward stochastic differential equations, adaptive control variate, semilinear parabolic PDE AMS subject classifications. 60H10, 60H35, 65C05, 65G99

We present and analyze an algorithm to solve numerically BSDEs based on Picard's iterations and on a sequential control variate technique. Its convergence is geometric. Moreover, the solution provided by our algorithm is regular both w.r.t. time and space.

1. Introduction. Let (Ω, F , P) be a given probability space on which is defined a q-dimensional standard Brownian motion W , whose natural filtration, augmented with P-null sets, is denoted (F t ) 0≤t≤T (T is a fixed terminal time). We aim at numerically approximating the solution (Y, Z) of the following forward backward stochastic differential equation (FBSDE) with fixed terminal time T

-dY t = f (t, X t , Y t , Z t )dt -Z t dW t , Y T = Φ(X T ), (1.1)
where f : [0, T ] × R d × R × R q → R, Φ : R d → R and X is the R d -valued process solution of

X t = x + t 0 b(s, X s )ds + t 0 σ(s, X s )dW s , (1.2)
b : [0, T ] × R d → R d and σ : [0, T ] × R d → R d×q .
The main focus of this work is to provide and analyze an algorithm -based on Picard's iterations and an adaptive Monte Carlo method -to approximate the solution (Y, Z) of (1.1).

Several algorithms to solve BSDEs can be found in the literature. Ma, Protter and Yong [START_REF] Ma | Solving forward backward stochastic differential equations explicitly-a four step scheme[END_REF] present an algorithm to solve quasilinear PDEs (associated to forward BS-DEs) using a finite difference approximation. Concerning algorithms based on the dynamic programming equation, we refer to Bouchard and Touzi [START_REF] Bouchard | Discrete time approximation and Monte Carlo simulation of backward stochastic differential equations[END_REF], Gobet, Lemor and Warin [START_REF] Gobet | A regression-based Monte Carlo method to solve backward stochastic differential equations[END_REF], Bally and Pagès [START_REF] Bally | Error analysis of the optimal quantization algorithm for obstacle problems[END_REF] and Delarue and Menozzi [START_REF] Delarue | A forward-backward stochastic algorithm for quasi-linear PDEs[END_REF]. In [START_REF] Bouchard | Discrete time approximation and Monte Carlo simulation of backward stochastic differential equations[END_REF], the authors compute the conditional expectations appearing in the dynamic programming equation by using Malliavin calculus techniques, whereas [START_REF] Gobet | A regression-based Monte Carlo method to solve backward stochastic differential equations[END_REF] proposes a scheme based on iterative regression functions which are approximated by projections on a reduced set of functions, the coefficients of the projection being evaluated using Monte Carlo simulations. [START_REF] Bally | Error analysis of the optimal quantization algorithm for obstacle problems[END_REF] and [START_REF] Delarue | A forward-backward stochastic algorithm for quasi-linear PDEs[END_REF] use quantization techniques for solving reflected BSDEs and forward BSDEs respectively. Bender and Denk [START_REF] Bender | A forward scheme for backward SDEs[END_REF] propose a forward scheme which avoids the nesting of conditional expectations backwards through the time steps. Instead, it mimics Picard's type iterations for BSDEs and, consequently, has nested conditional expectations along the iterations. This work has some connections with our approach but it does not handle the error analysis for the conditional expectations.

Our algorithm works as follows. First, we use Picard's iterations to approximate the solution (Y, Z) of (1.1) by the solutions of a sequence of linear BSDEs converging geometrically fast to (Y, Z) (see El Karoui, Peng and Quenez [START_REF] Karoui | Backward Stochastic Differential Equations in Finance[END_REF] for more details).

Then, since we can link linear BSDEs and linear PDEs, we use the adaptive control variate method proposed by Gobet and Maire [START_REF] Gobet | Sequential control variates for functionals of Markov processes[END_REF]. This method approximates the solutions of linear PDEs, which can be written as expectations of functionals of Markov processes via Feynman-Kac's formula. The authors use a control variate, changing at each step of the algorithm, to reduce the variance of the simulations. The convergence of this technique is also geometric w.r.t. the iterations (see [START_REF] Gobet | Sequential control variates for functionals of Markov processes[END_REF] for more details). As a consequence, we can guess that combining these two methods in order to solve BSDEs will lead to a geometrically converging algorithm. As a difference with previous works, we provide an approximated solution to the semi-linear PDE that has the same smoothness as the exact solution, which is quite satisfactory and potentially useful. As another difference with other Monte Carlo approaches, the final accuracy does not depend much on the number of simulations, but rather on an operator P used for approximating functions. It means that our algorithm benefits of the ability of Monte Carlo methods to solve high dimensional problems, without suffering for their relatively low rate of convergence. The paper is organized as follows. In Section 1, we give some definitions and notations, and recall the link between BSDEs and semilinear PDEs. In section 2, we describe the two main ingredients (Picard's iterations and adaptive control variate) of our algorithm. In Section 3, we define the norm used to measure the convergence of the algorithm, and in Section 4, we present the operator P used in the algorithm to approximate functions, emphasizing the important properties that P should satisfy in order to make our algorithm converge. We give the main convergence result in Section 5, and its proof in Section 6. We present in Section 7 an example of an operator P based on kernel estimators. Finally, in Section 8, we expose some numerical results in the field of financial mathematics.

Definitions and Notations.

• Let C k,l b be the set of continuously differentiable functions φ : (t, x) ∈ [0, T ] × R d with continuous and uniformly bounded derivatives w.r.t. t (resp. w.r.t.

x) up to order k (resp. up to order l). The function φ is also bounded. • C k p denotes the set of C k-1 functions with piecewise continuous k-th derivative.

• C k+α , α ∈]0, 1] is the set of C k functions whose k-th derivative is Hölder continuous of order α. • f v function. Let f v : [0, T ] × R d → R denote the following function f v (t, x) = f (t, x, v(t, x), (∂ x vσ)(t, x)),
where f denotes the driver of BSDE (1.1), σ denotes the diffusion coefficient of the SDE satisfied by X and v :

[0, T ] × R d → R is C 1 in space.
• Euler scheme. When it exists, we approximate the solution of (1.2) by its N -time-steps Euler scheme, denoted X N :

∀s ∈ [0, T ], dX N s = b(ϕ(s), X N ϕ(s) )ds + σ(ϕ(s), X N ϕ(s) )dW s , (1.3)
ϕ(s) := sup{t j : t j ≤ s} and

{0 = t 0 < t 1 < • • • < t N = T } is a regular subdivision of the interval [0, T ].
• Transition density function p(t, x; s, y). If σ is uniformly elliptic, the Markov process X admits a transition probability density p(t, x; s, y). Concerning X N (which is not Markovian except at times (t k ) k ), for any s > 0 X N s has a probability density p N (0, x; s, y) w.r.t. the Lebesgue measure. • Ψ(s, y, g 1 , g 2 , W ) and Ψ N (s, y, g 1 , g 2 , W ). We define the two following functions

Ψ(s, y, g 1 , g 2 , W ) = T s g 1 (r, X s,y r (W ))dr + g 2 (X s,y T (W )), Ψ N (s, y, g 1 , g 2 , W ) = T s g 1 (r, X N,s,y r (W ))dr + g 2 (X N,s,y T (W )),
where X s,y (resp. X N,s,y ) denotes the diffusion process solving (1.2) and starting from y at time s (resp. its approximation using an Euler scheme with N time steps), and W denotes the standard Brownian motion appearing in (1.2) and used to simulate X N , as given in (1.3).

• Constants c i,j (•). For any function φ in C i,j b , c i,j (φ) denotes i,j k,l=0 |∂ k x ∂ l t φ| ∞ . For i = j = 0, we set c 0 (φ) := c 0,0 (φ).
• Functions K(T ). K(•) denotes a generic function non decreasing in T which may depend on d, µ, β, on the coefficients b and σ (through

σ 0 , σ 1 , c 1,3 (σ), c 0,1 (∂ t σ), c 1,3 (b)
) and on other constants appearing in the Appendix A (i.e.

K A.1 (T ), α 1 , α 2 , c A.2 , K A.2 (T ) , c A.3 , K A.3 (T ), c A.4 , K A.4 (T ), c A.5 , K A.5 (T ) and K B.2 (T ))
. The parameter β is defined in Section 2.1, µ is defined in Section 3.2, σ 0 and σ 1 are defined in Hypothesis 1. • Functions K 0 (T ). K 0 (T ) are analogous to K(T ) except that they may also depend on the operator P (through c 1 (K t ) and c 2 (K x ), defined in Section 7). Hypothesis 1.

• The driver f is a bounded Lipschitz continuous function, i.e. for all

(t 1 , x 1 , y 1 , z 1 ), (t 2 , x 2 , y 2 , z 2 ) ∈ [0, T ] × R d × R × R d×q , |f (t 1 , x 1 , y 1 , z 1 ) -f (t 2 , x 2 , y 2 , z 2 )| ≤ L f (|t 1 -t 2 | + |x 1 -x 2 | + |y 1 -y 2 | + |z 1 -z 2 |).
• σ is uniformly elliptic on [0, T ] × R d : there exist two positive constants σ 0 , σ 1 s.t. for any vector ξ and any

(t, x) ∈ [0, T ] × R d σ 0 |ξ| 2 ≤ d i,j=1 [σσ * ] i,j (t, x)ξ i ξ j ≤ σ 1 |ξ| 2 . • Φ is bounded in C 2+α , α ∈]0, 1]. • b and σ are in C 1,3 b and ∂ t σ is in C 0,1 b . 1.2.
Link with semilinear PDE. According to [START_REF] Pardoux | Backward Stochastic Differential Equations and Quasilinear Parabolic Partial Differential Equations[END_REF]Theorem 3.1] (see also [START_REF] Karoui | Backward Stochastic Differential Equations in Finance[END_REF]Proposition 4.3] and [8, Theorem 2.1]), we can link (Y, Z)-the solution of the BSDE (1.1) -to u, the solution of the following PDE :

∂ t u(t, x) + Lu(t, x) + f (t, x, u(t, x), (∂ x uσ)(t, x)) = 0, u(T, x) = Φ(x), (1.4)
where L is defined by

Lu(t, x) = 1 2 i,j [σσ * ] ij (t, x)∂ 2 xixj u(t, x) + i b i (t, x)∂ xi u(t, x).
Theorem 1.1 ([8], Theorem 2.1). Under Hypothesis 1, the solution u of PDE (1.4)

belongs to C 1,2 b . Moreover, (Y t , Z t ) 0≤t≤T -solution of (1.1)-satisfies ∀t ∈ [0, T ], (Y t , Z t ) = (u(t, X t ), ∂ x u(t, X t )σ(t, X t )).
(1.5)

2. Description of the algorithm. As said in the introduction, the current algorithm is based on two ingredients: Picard's iterations and adaptive control variates. We present these ingredients in the following section, before describing the main algorithm in Section 2.3.

2.1. First ingredient: Picard's iterations. ¿From [10, Corollary 2.1], we know that under standard assumptions on Φ and f , the sequence ( Ŷ k , Ẑk ) k recursively defined by ( Ŷ0 = 0, Ẑ0 = 0) and

-d Ŷ k+1 t = f (t, X t , Ŷ k t , Ẑk t )dt -Ẑk+1 t dW t , Ŷ k+1 T = Φ(X T ) converges to (Y, Z), dP ⊗ dt a.s. (and in H 2 T,β (R) × H 2 T,β (R q )) as k goes to +∞, where H 2 T,β (R q ) := {φ ∈ H 2 T (R q ) such that E T 0 e βt |φ t | 2 dt < ∞} and β is such that 2(1 + T )L f < β.
∂ t ûk+1 + Lû k+1 + f (•, •, ûk , ∂ x ûk σ) = 0 and ûk+1 (T, •) = Φ(•).
It means that the sequence of solutions of linear PDEs (û k , ∂ x ûk ) k converges (in a L 2 norm) to (u, ∂ x u), solution of the semi-linear PDE (1.4).

Second ingredient:

Adaptive control variate. In their work [START_REF] Gobet | Sequential control variates for functionals of Markov processes[END_REF], Gobet and Maire present an adaptive algorithm to solve linear PDEs of type

∂ t v + Lv + g = 0 and v(T, •) = Φ(•).
Thanks to the Feynman-Kac formula, we know that the probabilistic solution of this PDE is v(t, x) = E t,x [Φ(X T ) + T t g(s, X s )ds] = E(Ψ(t, x, g, Φ, W )). Their idea is to compute a sequence of solutions (v k ) k by writing

v k+1 = v k + Monte Carlo evaluations of the error (v -v k ).

The probabilistic representation of the correction term c

k := v -v k is c k (t, x) = v(t, x) -v k (t, x) = E(Ψ(t, x, g + ∂ t v k + Lv k , Φ -v k (T, .), W )).
Their algorithm computes iterative approximations (v k ) k of the global solution v. These approximations rely on the computations of E(Ψ N (t, x, g, Φ, W )) (for data g and Φ possibly different from g +∂ t v k +Lv k and Φ-v k (T, .)) at some points (t i , x i ) 1≤i≤n ⊂ [0, T ] × R d . We briefly recall below their algorithm and the associated convergence result.

Algorithm 1. Begin with v 0 ≡ 0 and assume that an approximated solution

v k of class C 1,2 b is built at step k. • Evaluate c k (t i , x i ) using M independent simulations c M k (t i , x i ) = 1 M M m=1 Ψ N (t i , x i , g + ∂ t v k + L N v k , Φ -v k , W m ) where L N u(s, x) = 1 2 i,j [σσ * ] ij (ϕ(s), x)∂ 2 xixj u(s, x) + i b i (ϕ(s), x)∂ xi u(s, x).
• Build the global solution c M k (•) based on the values [c M k (t i , x i )] 1≤i≤n by using a linear approximation operator:

Pc M k (•) = n i=1 c M k (t i , x i )ω i (•) for some weight functions ω i . Then, get v k+1 = P(v k + c M k ).
(Examples of operator P: interpolation, projection, Kernel-based estimator...)

The main result of their paper is the following Theorem 2.1 (Theorems 3.1 and 3.2, [START_REF] Gobet | Sequential control variates for functionals of Markov processes[END_REF]). Define v-v k

2 2 := sup 1≤i≤n |E((v- v k )(t i , x i ))| 2 + sup 1≤i≤n Var(v k (t i , x i )). Then for any k ≥ 0, one has v -v k+1 2 2 ≤ ρ v -v k 2 2 + C v -Pv 2 2 ( 1 N + 1 M ),
where ρ < 1 (depending on M , N and P) and v -Pv 2 is a suitable norm related to the approximation error v -Pv (that we do not detail).

As for Picard's iterations, the algorithm converges at a geometric rate. Moreover, there is no need to take N and M large (in practice, in their experiments, they take M = 10 to get accurate approximations). The final error is strongly related to the ability of the operator P to approximate well the true solution v. In the following Section, we present an algorithm combining these two ingredients and leading to the same features.

2.3. Algorithm. We recall that we aim at numerically solving BSDE (1.1), which is equivalent to solving the semilinear PDE (1.4). The current algorithm provides an approximation of the solution of this PDE. Then, by simulating the diffusion X through an Euler scheme, we deduce from Equality (1.5) an approximation of the solution of BSDE (1.1). More precisely, let u k (resp. (Y k , Z k )) denote the approximation of u (resp. (Y, Z)) at step k, and let X N denote the approximation of X obtained with a N -time-steps Euler scheme. We write

(Y k t , Z k t ) = (u k (t, X N t ), ∂ x u k (t, X N t )σ(t, X N t )), for all t ∈ [0, T ]. (2.1)
where X N is described in Section 1.1. It remains to build u k+1 . Adaptive control variate. As in Algorithm 1, we write

u k+1 = u k + Monte Carlo evaluations of the error(u -u k ).
Combining Itô's formula applied to u(s, X s ) and to u k (s, X N s ) between t and T and the semilinear PDE (1.4) satisfied by u, we get that the correction term

c k := u -u k is c k (t, x) = E Ψ (t, x, f u , Φ, W ) -Ψ N t, x, -(∂ t + L N )u k , u k (T, .), W |G k .
Remark 2.2. As we will see later (see Remark 2.3), u k depends on several random variables. G k is the σ-algebra generated by the set of all random variables used to build u k . In the above equation, we compute the expectation w.r.t. the law of X and X N and not w.r.t. the law of u k , which is G k measurable. Picard's iteration. The correction term c k cannot be used directly: we have to replace u and ∂ x u (unknown terms) appearing in f by u k and ∂ x u k , as suggested by Picard's contraction principle:

ĉk (t, x) = E Ψ (t, x, f u k , Φ, W ) -Ψ N t, x, -(∂ t + L N )u k , u k , W |G k . (2.2)
We still have to replace the expectation by a Monte Carlo summation, and Ψ by Ψ N . Then, the algorithm computes iterative approximations of (u k ) k of the global solution u at some points (t k i , x k i ) 1≤i≤n ∈ [0, T ] × R d , which may change over the iterations.

Algorithm 2. We begin with u 0 ≡ 0. Assume that an approximated solution u k of class C 1,2 is built at step k.

• Evaluate ĉk (t k i , x k i ) using M independent simulations

ĉM k (t k i , x k i ) = 1 M M m=1 Ψ N t i , x i , f u k + (∂ t + L N )u k , Φ -u k , W m,k,i . • Build the global solution ĉM k (•) based on the values [ĉ M k (t k i , x k i )]
1≤i≤n by using a linear approximation operator:

P k c(•) = n i=1 c(t k i , x k i )ω i (•), (2.3)
where (ω k i ) i are some weight functions. Deduce the approximation of u at step k + 1

u k+1 (t, x) = P k (u k + ĉM k )(t, x), (2.4)
where (P k ) k satisfies Hypothesis 2 (defined later in Section 4).

Remark 2.3. Since u k+1 is computed by using (2.4), u k+1 is a random function depending on the random variables needed to compute u k and W m,k,i , 1 ≤ m ≤ M, 1 ≤ i ≤ n, appearing in the computation of ĉM k . Moreover, P k may be random. In such a case, u k+1 also depends on the random variables used to build P k . Definition 2.4 (Definition of the σ-algebra G k ). Let G k+1 define the σ-algebra generated by the set of all random variables used to build u k+1 . Using (2.4) yields

G k+1 = G k ∨ σ(A k , S k ),
where A k is the set of random points used at step k to build the estimator P k , S k := {W m,k,i , 1 ≤ m ≤ M, 1 ≤ i ≤ n}, the set of independent Brownian motions used to simulate the paths X m,k,N (x k i ), and G k is the σ-algebra generated by the set of all random variables used to build u k .

3. Choice of the norm to measure the convergence. The choice of the norm to measure the convergence is not harmless. To prove the convergence of the algorithm we combine results on BSDEs stated in a norm leading to the integration w.r.t. e βs ds (see [START_REF] Karoui | Backward Stochastic Differential Equations in Finance[END_REF]), and results on the bounds for solutions of linear PDEs in weighted Sobolev spaces (leading to the integration w.r.t. e -µ|x| dx), coming from [START_REF] Bensoussan | Applications des inéquations variationnelles en contrôle stochastique[END_REF], and recalled in Theorem B.2. Although rather technical, this section is crucial in order to analyze the convergence of the algorithm and it is interesting for itself.

3.1. Norm of the convergence.

Definition 3.1 (Definition of H µ β (R q )).
For any β > 0, we define

(H µ β (R q ), • µ,β ) the set of processes V : Ω × [0, T ] × R d → R q that are Pr ⊗ B(R d )- measurable (where Pr is the σ-field of predictable subsets of Ω × [0, T ]) and such that V 2 µ,β := E T 0 R d e βs |V s (x)| 2 e -µ|x| dxds < ∞.
We use this norm to measure the error (Y -Y k , Z -Z k ), corresponding to the error we make at step k of the algorithm. Using (1.5) and (2.1), we get

Y -Y k 2 µ,β = E T 0 R d e βs |u(s, X x s ) -u k (s, X N,x s )| 2 e -µ|x| dxds , (3.1)
and

Z -Z k 2 µ,β = E T 0 R d e βs |(∂ x uσ)(s, X x s ) -(∂ x u k σ)(s, X N,x s )| 2 e -µ|x| dxds . Since u, u k , ∂ x u and ∂ x u k are bounded (see Theorem 1.1 and Hypothesis 2 defined later), Y -Y k 2 µ,β and Z -Z k 2 µ,β are finite.
Remark 3.2. We point out that the expectation appearing in the above definition of Y k -Y 2 µ,β and Z k -Z 2 µ,β is computed w.r.t. the law of X, X N and all the possible random variables used to compute u k .

3.2. Some other useful norms. The following Definitions introduce two norms strongly related to • µ,β . They will be useful in the proof of the main result.

Definition 3.3 (Space H m,µ β ). For any m ≤ 2, β > 0 and µ > 0, let H m,µ β define the space of functions v : [0, T ] × R d → R such that v 2 H m,µ β = T 0 e βs R d e -µ|x| k≤m |∂ k x v(s, x)| 2 dxds < ∞.
For m = 0, we set H 0,µ β = H µ β . Definition 3.4 (Function ν t µ ). For any s, t ∈ [0, T ] and any x, y ∈ R d such that t < s we define ν t µ (s, y) := R d e -µ|x| p(t, x; s, y)dx, where µ is a positive constant.

Definition 3.5 (Space H m,µ β, X ).
For any m ≤ 2, β > 0, µ > 0 and any diffusion process Xs , 0 ≤ s ≤ T starting from x at time 0, let H m,µ β,

X define the space of functions v : [0, T ] × R d → R such that v 2 H m,µ β, X = T 0 e βs R d e -µ|x| k≤m E|∂ k x v(s, Xx s )| 2 dxds < ∞.
For m = 0, we set H 0,µ β,

X = H µ β,
X . By using the definition of ν, we also get v 2

H µ β,X = T 0 e βs R d ν 0 µ (s, y)|v(s, y)| 2 dyds. Remark 3.6. We can interpret v 2 H µ β,X as T 0 e βs E|v(s, X µ s )| 2 ds
where (X µ s ) s stands for (X s ) s starting at time 0 with a random initial law having a density proportional to e -µ|x| .

Remark 3.7. In the above definition of v

T, d, µ, K A.1 (T ), α 1 , α 2 ) s.t. for every v ∈ L 2 ([0, T ] × R d , e βt dt ⊗ e -µ|x| dx) c v 2 H µ β ≤ v 2 H µ β,X ≤ c v 2 H µ β .
Proof. Proving Proposition 3.8 boils down to showing that there exist two con-

stants c i , i = 1, 2 depending on µ, d, α i s.t. ∀y ∈ R d 1 2 d K(T ) e -µ|y| e c1(s-t) ≤ ν t µ (s, y) ≤ 2 d K(T )e c2(s-t) e -µ|y| . (3.2)
The r.h.s. (resp. the l.h.s.) of the above inequality ensues from Proposition A.1 and Lemma C.1, with c = 1 2α2 (resp. c = 1 2α1 ). Moreover, we have

c i = dµ 2 αi 2 , i = 1, 2.
Proposition 3.9. Assume σ is uniformly elliptic and b and σ are in

C 0,2 b . For every v ∈ L 2 ([0, T ] × R d , e βt dt ⊗ e -µ|x| dx) it holds v 2 H µ β,X N ≤ K(T ) v 2 H µ β,X . Proof. Since v 2 H µ β,X N = T 0 e βs R d e -µ|x| R d v 2
(s, y)p N (0, x; s, y)dydxds, we apply successively Proposition A.5, Lemma C.1 and the l.h.s. of (3.2).

The following Proposition bounds the solution of a linear PDE in norm

• H µ β,X . Proposition 3.10. Assume σ is uniformly elliptic, σ ∈ C 1,1 b , b is C 1,1
and bounded, and g ∈ L 2 ([0, T ] × R d , e βt ν 0 µ (t, y)dt ⊗ dy). We also assume that v satisfies the linear PDE (∂ t + L)v(t, x) + g(t, x) = 0, with a terminal condition v(T, •) = 0. Then,

v 2 H 2,µ β,X + ∂ t v 2 H µ β,X ≤ K 3.10 (T ) g 2 H µ β,X
, where K 3.10 (T ) depends on K B.2 (T ), β, c, c and T .

Proof. Let v(t, x) := e βt/2 v(t, x). Using the PDE satisfied by v, we get

-∂ t v(t, x) -Lv(t, x) + β 2 v(t, x) = e βt 2 g(t, x), v(T, x) = 0.
Then, we apply Theorem B.2 to v to obtain v 2

H 2,µ β + ∂ t v 2 H µ β ≤ K(T ) g 2 H µ β
, where K(T ) denotes a constant depending on K B.2 (T ) and on β. Proposition 3.8 (applied to v and its derivatives) ends the proof.

4. Approximation operator P. P denotes the sequence of approximation operators (P k ) k satisfying the following Hypothesis. In Section 7, we give an example of such an operator P, based on kernel estimators.

Hypothesis 2 (Hypotheses on (P k ) k ). Let ǫ i (P), i = 1 • • • 4 denote some constants depending on P and tending to 0 when the parameters of P tend to infinity. For any k, the random operator P k satisfies the following properties.

1. Measurability and linearity. P k is linear and it writes (see (2.3))

P k c(•) = n i=1 c(t k i , x k i )ω i (•).
As said in Definition 2.4, the points (t k i , x k i ) i are σ(A k )-measurable and A k is independent of A j for k = j. In addition, the weights w k i (.

) : Ω×[0, T ]×R d → R are measurable w.r.t. σ(A k ) × B([0, T ] × R d ). 2. Regularity. For any function v : [0, T ] × R d → R, P k v ∈ C 1,2 (equivalently, the weights (w k i ) i are C 1,2 w.r.t. (t, x)).

Boundedness. For any bounded function

v : [0, T ] × R d → R, P k v, ∂ t (P k v), ∂ x (P k v) and ∂ 2
x (P k v) are bounded by a constant depending on c 0 (v) and on the parameters of the operator P k (the explicit form of the bounds is not needed). 4. Approximation. P k approximates well a function and its spatial derivatives.

Let v be a C 1,2 random function, independent of A k , from [0, T ]×R d in R. We also assume that v and ∂ x v are bounded (by c 0,1 (v)) and ∀t,

t ′ ∈ [0, T ], ∀x ∈ R d , |∂ x v(t, x) -∂ x v(t ′ , x)| ≤ c 1/2 (v) |t ′ -t|. Then, E v -P k v 2 H µ β,X + E ∂ x v -∂ x (P k v) 2 H µ β,X ≤ ǫ 1 (P)(E v 2 H 2,µ β,X + E ∂ t v 2 H µ β,X ) + ǫ 2 (P)(c 2 1/2 (v) + c 2 0,1 (v)). 5. Let v be a C 1,2 b function from [0, T ] × R d in R. Then E v -P k v 2 H µ β,X + E ∂ x v -∂ x (P k v) 2 H µ β,X ≤ ǫ 3 (P)(c 2 1/2 (v) + c 2 1,2 (v)
). 6. Stability and centering property for random functions. For any random func-

tion v from [0, T ] × R d in R, one has E P k v 2 H µ β,X + E ∂ x (P k v) 2 H µ β,X ≤ c 4 (P)E v 2 H µ β,X
, where c 4 (P) is a constant depending on P. If E(v(t, x)) = 0,

one has E P k v 2 H µ β,X + E ∂ x (P k v) 2 H µ β,X ≤ ǫ 4 (P)E v 2 H µ β,X .
5. Main result. The following Theorem states the convergence of (Y k , Z k ) k defined in (2.1) towards (Y, Z) in the • 2 µ,β norm. See Section 6 for a proof of it. Theorem 5.1. Assume that Hypotheses 1 and 2 hold. We also assume that f is a bounded Lipschitz function and Φ ∈ C 2+α b . Then, there exists a constant K(T ) such that

Y -Y k 2 µ,β + Z -Z k 2 µ,β ≤S k + K(T ) c 2 0,2 (u) N , where S k ≤ ηS k-1 + ǫ and η = K(T ) 4(1 + T )L 2 f β + ǫ 1 (P)K 3.10 (T )L 2 f + ǫ 4 (P) M , ǫ = K(T ) ǫ 3 (P)c 2 1,2 (u) + (ǫ 2 (P) + ǫ 3 (P))c 2 0 (f ) + ǫ 4 (P)c 2 0,2 (u) M N + c 4 (P) N 2 (c 2 0 (f ) + c 2 0 (Φ)) .
Corollary 5.2. Under the assumptions of Theorem 5.1, for β and P-parameters large enough so that η < 1, we have

lim sup k→∞ Y -Y k 2 µ,β + Z -Z k 2 µ,β ≤ ǫ 1 -η + K(T ) c 2 0,2 (u) N .
Remark 5.3. Using a sequential Monte Carlo method leads to an error of 1 M appearing in the contraction term η. If we had implemented a non adaptive method (i.e. using only Picard's iterations), M -1 would have appeared in ǫ and this would have led us to choose a much larger M , while practically M = 100 does the trick. A quick analysis of the form of ǫ shows that it is sufficient to take M and N of the same magnitude (taking N = M = 100 gives usually small related errors). Then, the final algorithm accuracy relies heavily on the quality of the operator P, through the convergence rates ǫ 2 (P) and ǫ 3 (P). These rates are related to those obtained in the approximation of a function v and its first partial derivative ∂ x v, when roughly speaking the function v is C 1,2 b and ∂ x v is 1 2 -Hölder in time. 6. Proof of Theorem 5.1. In this section, we assume that Hypotheses 1 and 2 hold.

We want to measure the error

Y -Y k 2 µ,β + Z -Z k 2 µ,β .
To do so, we split the different sources of errors: Euler scheme approximation, Picard's approximation, approximation operator P and Monte Carlo simulations. We prove Theorem 5.1 by combining Proposition 6.1, Equation (6.2), Propositions 6.3, 6.4 and 6.8.

6.1. First source of error: Euler scheme approximation. Proposition 6.1. It holds

Y -Y k 2 µ,β + Z -Z k 2 µ,β ≤c 2 0,2 (u)
K(T ) N + K(T )( Y -Ỹ k 2 µ,β + Z -Zk 2 µ,β ),
where

( Ỹ k s , Zk s ) := (u k (s, X s ), (∂ x u k σ)(s, X s )).
Proof. We start from (3.1) and we introduce ±u(s,

X N,x s ) in |u(s, X x s ) -u k (s, X N,x s )|. We get |u(s, X x s ) -u k (s, X N,x s )| 2 ≤ 2|u(s, X x s ) -u(s, X N,x s )| 2 + 2|u(s, X N,x s ) -u k (s, X N,x s )| 2 . Since u is a C 1
-Y k 2 µ,β ≤ c 2 0,1 (u) K(T ) N + 2E u -u k 2 H µ β,X N
. The same technique enables to bound Z -Z k 2 µ,β . Finally, we apply Proposition 3.9 and use the equality

E u -u k 2 H µ β,X = Y -Ỹ k 2
µ,β (see Remark 6.2) to end the proof.

Remark 6.2. We recall that the expectation appearing in the definition of Y -Ỹ k 2 µ,β is computed w.r.t. the law of X and all the variables used to compute u k , whereas the one appearing in the definition of

• 2 H µ β,X is computed w.r.t.
the law of X, that's why we have

Y -Ỹ k 2 µ,β = E u -u k 2 H µ β,X . u -u k 2 H µ β,X is G k -measurable.

Let us now work on the remaining term

Y -Ỹ k 2 µ,β + Z -Zk 2 µ,
Y k t = Φ(X T ) + T t f (s, X s , Ỹ k-1 s , Zk-1 s )ds - T t Z k s dW s . (6.1) We split Y -Ỹ k 2 µ,β and Z -Zk 2 µ,β by introducing Y k and Z k : Y -Ỹ k 2 µ,β + Z -Zk 2 µ,β ≤2( Y -Y k 2 µ,β + Z -Z k 2 µ,β ) + 2( Y k -Ỹ k 2 µ,β + Z k -Zk 2 µ,β ). (6.2) Let us study Y -Y k 2 µ,β + Z -Z k 2 µ,β . The study of Y k -Ỹ k 2 µ,β and Z k -Zk
Y -Y k 2 µ,β + Z -Z k 2 µ,β ≤ 2(1 + T )L 2 f β Y -Ỹ k-1 2 µ,β + Z -Zk-1 2 µ,β ,
where (Y k , Z k ) is the solution of (6.1).

Proof.

Let
T 0 e βs E[|Y s -Y k s | 2 G k-1 ]ds ≤ T β T 0 e βs E[|f (s, X s , Y s , Z s ) - f (s, X s , Ỹ k-1 s , Zk-1 s )| 2 G k-1 ]ds. Since f is Lipschitz, T 0 e βs E |Y s -Y k s | 2 G k-1 ds ≤ 2T L 2 f β T 0 e βs E |Y s -Ỹ k-1 s | 2 + |Z s -Zk-1 s | 2 G k-
Y k -Ỹ k 2 µ,β + Z k -Zk 2 µ,β ≤ η( Y -Ỹ k-1 2 µ,β + Z -Zk-1 2 µ,β ) + 3E P k-1 (ĉ k-1 -ĉM k-1 ) 2 H µ β,X +3c 2 0 (σ)E ∂ x P k-1 (ĉ k-1 -ĉM k-1 ) 2 H µ β,X
+ǫ where η = 3ǫ 1 (P)K 3.10 (T )L 2 f and ǫ = K(T )(ǫ 3 (P)c 2 1,2 (u) + (ǫ 2 (P) + ǫ 3 (P))c 2 0 (f )). Before proving Proposition 6.4, we introduce two sequences (u k ) k and (H k ) k . Definition 6.5. ∀k ≥ 1, u k denotes the solution of the linear PDE

(∂ t + L)u k (t, x) + f u k-1 (t, x) = 0, u k (T, x) = Φ(x). Then, ∀t ∈ [0, T ], (Y k t , Z k t ) = (u k (t, X t ), (∂ x u k σ)(t, X t )), where (Y k t , Z k t ) k
is the solution of BSDE (6.1). Definition 6.6. ∀k ≥ 0, H k denotes the solution of the linear PDE

(∂ t + L)H k (t, x) + f u k (t, x) -f u (t, x) = 0, H k (T, x) = 0. (6.3)
Clearly, one has u k = u + H k-1 for any k ≥ 1.

Proof. [Proof of Proposition6.4] Let us work on Y

k -Ỹ k 2 µ,β = E u k -u k 2 H µ β,X . Using (2.4) gives u k -u k = u k -P k-1 (u k-1 + ĉM k-1 ) = u k -P k-1 (u k-1 + ĉk-1 ) + P k-1 (ĉ k-1 -ĉM k-1 ). Since u k-1 + ĉk-1 = u k = u + H k-1 (see (2.
2)), we get

u k -u k = (u -P k-1 u) + (H k-1 -P k-1 H k-1 ) + P k-1 (ĉ k-1 -ĉM k-1 ). (6.4)
Then, since u and its derivatives are bounded and since H k-1 satisfies PDE (6.3), we combine Lemma 6.7 and the features of P k-1 (see Hypothesis 2) to get

E u k -u k 2 H µ β,X ≤ 3ǫ 3 (P)(c 2 1/2 (u) + c 2 1,2 (u)) + 3ǫ 1 (P)(E H k-1 2 H 2,µ β,X +E ∂ t H k-1 2 H µ β,X ) + 3ǫ 2 (P)(c 2 1/2 (H k-1 ) + c 2 0,1 (H k-1 )) + 3E P k-1 (ĉ k-1 -ĉM k-1 ) 2 H µ β,X . (6.5)
We bound the second term on the r.h.s. by using Proposition 3.10:

E H k-1 2 H 2,µ β,X +E ∂ t H k-1 2 H µ β,X ≤ K 3.10 (T )E f u k-1 -f u 2 H µ β,X .
Since f is Lipschitz, the result follows. The proof is similar for

Z k -Zk 2 µ,β , except that (6.5) contains c 2 0 (σ)E ∂ x P k-1 (ĉ k-1 -ĉM k-1 ) 2 H µ β,X
(and not

E P k-1 (ĉ k-1 -ĉM k-1 ) 2 H µ β,X
).

Lemma 6.7 (resp. Proposition B.3) enables to replace constants

c 2 1/2 (H k-1 ) + c 2 0,1 (H k-1 ) (resp. c 2 1/2 (u)) by K(T )c 2 0 (f ). (resp. K(T )(c 2 0 (f ) + c 2 2 (Φ))). Since Φ(•) = u(T, •), we get c 2 2 (Φ) ≤ c 2 0,2 (
u) and the result follows.

Lemma 6.7. For all k, H k and ∂ x H k are bounded by a constant of the form K(T )c 0 (f ).

Moreover, for all t, t ′ ∈ [0, T ],

|∂ x H k (t ′ , x) -∂ x H k (t, x)| ≤ c 1/2 (H k ) |t ′ -t| with c 1/2 (H k ) = K(T )c 0 (f )
Proof. [Proof of the lemma] First, we prove that H k and ∂ x H k are bounded. Using Feynman-Kac's formula yields

H k (t, x) = E T t (f u k (s, X t,x s ) -f u (s, X t,x s ))ds|G k . Since f is bounded, we get that H is bounded by 2T c 0 (f ). To prove that ∂ x H k is bounded, we write H k (t, x) = T t
R d fk (s, y)p(t, x; s, y)dyds, where p denotes the transition density of X and fk := f u kf u . We differentiate H k w.r.t. x and we use

Proposition A.2 to get |∂ x H k (t, x)| ≤ 2K(T )c 0 (f ) T t ds √ s-t R d e -c|x-y| 2 s-t dy (s-t) d 2
. Our statement on ∂ x H k readily follows. The second assertion ensues from Proposition B.3.

Fourth source of error: Monte Carlo simulations.

In this section, we study

E P k (ĉ k -ĉM k ) 2 H µ β,X and E ∂ x P k (ĉ k -ĉM k ) 2 H µ β,X . Proposition 6.8. Let ∆c k denote ĉk -ĉM k . It holds E P k (∆c k ) 2 H µ β,X +E ∂ x P k (∆c k ) 2 H µ β,X ≤ǫ 4 (P) K(T ) M ( Y -Ỹ k 2 µ,β + Z -Zk 2 µ,β ) +ǫ 4 (P) K(T ) M N c 2 0,2 (u) + c 4 (P) K(T ) N 2 (c 2 0 (f ) + c 2 0 (Φ)).
. Proof. We split ∆c k into two terms: the bias and the noise

P k (∆c k ) = P k (E(∆c k |G k )) + P k (∆c k -E(∆c k |G k )) := P k (E(∆c k |G k )) + P k (ε k ),
where

ε k := ∆c k -E(∆c k |G k ) = ĉM k -E(ĉ M k |G k ).
The same decomposition holds for ∂ x P k (∆c k ). It remains to apply Propositions 6.9 and 6.10 to end the proof.

6.4.1. Bias terms P k (E(∆c k |G k )) and ∂ x P k (E(∆c k |G k )). Proposition 6.9. It holds E P k (E(∆c k |G k )) 2 H µ β,X + E ∂ x P k (E(∆c k |G k )) 2 H µ β,X ≤ c 4 (P) K(T ) N 2 .
Proof. First, we use the properties of P (see Point 6, Hypothesis 2) to get

E P k (E(∆c k |G k )) 2 H µ β,X +E ∂ x P k (E(∆c k |G k )) 2 H µ β,X ≤ c 4 (P) E(∆c k |G k ) 2 H µ β,X
. Then, we define ∆(s, y)

:= Φ(X s,y T ) -Φ(X N,s,y T ) + T s (f u k (r, X s,y r ) -f u k (r, X N,s,y r
))dr. We have

E(∆c k |G k ) = ĉk -E(ĉ M k |G k ) = E(∆|G k ).
Let us bound E(∆|G k ). First, we work on the first term:

E[Φ(X s,y T ) - Φ(X N,s,y T )] = R d Φ(y ′ )(p(s, y; T, y ′ )-p N (s, y; T, y ′ ))dy ′ . By using Proposition A.4, we get E[Φ(X s,y T ) -Φ(X N,s,y T )] ≤ K(T ) N √ T -s R d |Φ(y ′ )| 1 (T -s) d 2 exp -cA.4|y ′ -y| 2 T -s dy ′ . Then, |E[Φ(X s,y T ) -Φ(X N,s,y T )]| ≤ c 0 (Φ) K(T ) N √ T -s. The same proof holds for T s (f u k (r, X s,y r ) -f u k (r, X N,s,y r ))dr: E|f u k (r, X s,y r ) -f u k (r, X N,s,y r )| ≤ K(T ) N T -s √ r-s c 0 (f ).
Integrating both sides w.r.t. r between s and T leads to

|E[∆(s, y)|G k ]| ≤ K(T ) N √ T -s(c 0 (Φ) + (T -s)c 0 (f ))
, and the result follows.

6.4.2. Noise terms P k (ε k ) and ∂ x P k (ε k ). Proposition 6.10.

It holds

E P k (ε k ) 2 H µ β,X + E ∂ x P k (ε k ) 2 H µ β,X ≤ ǫ 4 (P) K(T ) M ( c 2 0,2 (u) N + Y -Ỹ k 2 µ,β + Z -Zk 2 µ,β ).
Proof. Since E(ε k |G k ) = 0, the properties of P (Point 6, Hypothesis 2) give

E( P k ε k 2 H µ β,X |G k ) + +E( ∂ x P k (ε k ) 2 H µ β,X |G k ) ≤ ǫ 4 (P)E( ε k 2 H µ β,X |G k ). Since E( ε k 2 H µ β,X |G k ) = T 0 e βs R d ν 0 µ (s, y)E(|ε k (s, y)| 2 |G k )dsdy = T 0 e βs R d ν 0 µ (s, y)Var(ε k (s, y)|G k )dsdy, it remains to bound Var(ε k (s, y)|G k ) = 1 M Var(Ψ N (s, y, f u k + (∂ t + L N )u k , Φ -u k , W )|G k ).
By applying Itô's formula to u(r, X N r ) between s and T and to u k (r, X N r ) between s and T , we get Ψ N (s, y, f

u k +(∂ t +L N )u k , Φ-u k , W ) = (u-u k )(s, y)+ T s f u k (r, X N r )+ (∂ t + L N )u(r, X N r )dr + T s (∂ x u(r, X N r ) -∂ x u k (r, X N r ))σ(ϕ(r), X N ϕ(r) )dW r . Then, Var(ε k (s, y)|G k ) ≤ 2 M E   T s f u k (r, X N r ) + (∂ t + L N )u(r, X N r )dr 2 |G k   + 2 M E   T s (∂ x u(r, X N r ) -∂ x u k (r, X N r ))σ(ϕ(r), X N ϕ(r) )dW r 2 |G k   .
To bound the first term, we introduce ±f u (r, X N r ). Since f is Lipschitz continuous with constant L f and u satisfies (

∂ t + L)u + f = 0, we get |f u k (r, X N r ) + (∂ t + L N )u(r, X N r )| ≤ L f (|u k -u|(r, X N r ) + |∂ x u k σ -∂ x uσ|(r, X N r )) + d i=1 |b i (r, X N r ) -b i (ϕ(r), X N ϕ(r) )||∂ xi u(r, X N r )| + 1 2 d i,j=1 |[σσ * ] ij (r, X N r ) - [σσ * ] ij (ϕ(r), X N ϕ(r) )||∂ 2 xixj u(r, X N r )|.
Since b, σ and u are bounded with bounded derivatives, we get |f u k (r,

X N r ) + (∂ t + L N )u(r, X N r )| ≤ L f (|u -u k |(r, X N r ) + |∂ x u k σ -∂ x uσ|(r, X N r )) + K(T )c 0,2 (u)( T N + |X N r -X N ϕ(r) |). Thus, Var(ε k (s, y)|G k ) ≤ K(T ) M c 2 0,2 (u) N + T s E(Θ(r, X N r )|G k )dr , where Θ(r, z) = |u -u k | 2 (r, z) + |∂ x u -∂ x u k | 2 (r, z). Finally, it remains to compute T 0 e βs R d ν 0 µ (s, y) T s E(Θ(r, X N r )|G k )drdyds.
To do so, we write E(Θ(r, X N r )|G k ) = R d dzΘ(r, z)p N (s, y; r, z) and we use successively the r.h.s. of (3.2), Proposition A.5, Lemma C.1 and the l.h.s. of (3.2).

We get

E( ε k 2 H µ β,X |G k ) ≤ K(T ) M c 2 0,2 (u) N + E( u -u k 2 H µ β,X |G k ) + E( ∂ x u -∂ x u k 2 H µ β,X
|G k ) . It remains to use the ellipticity condition on σ to end the proof.

An example of operator P based on kernel estimators.

In this Section we present an operator satisfying Hypothesis 2. It is based on a non parametric regression technique called local averaging. We refer to Györfi et al [16, Chapter 2] and Härdle [START_REF] Härdle | Applied nonparametric regression[END_REF] for more details on non parametric regression.

Definition 7.1. We approximate a function v(t, x) by

P k v(t, x) = r k n (t, x) f k n (t, x) g(2 d+1 T λ(B)f k n (t, x)), (7.1)
where • r k n (t, x) = 1 nhth d x n i=1 K t ( t-T k i ht )K x ( x-X k i hx )v(T k i , X k i ), • f k n (t, x) = 1 nhth d x n i=1 K t ( t-T k i ht )K x ( x-X k i hx ), • the points (T k i , X k i ) 1≤i≤n are uniformly distributed on [0, T ] × B where B := [-a, a] d , • λ(B) = (2a) d , • g is such that g(y) =    0 if y < 0, 1 if y > 1, -y 4 + 2y 2 if y ∈ [0, 1]. (7.2) • The kernel function K t is defined on the compact support [-1, 1], bounded, even, non-negative, C 2 p and R K t (u)du = 1. • The kernel function K x is defined on the compact support [-1, 1] d , bounded,
and

such that ∀ y = (y 1 , • • • , y d ) ∈ R d , K x (y) = Π d j=1 K j x (y j ), where for j = 1, • • • , d K j x : R → R is an even non-negative C 2 p function and R K j x (y)dy = 1. • δ n denotes 1 nhth d x
, and T λ(B)δ n << 1.

• h x << a and h t << T 2 . Since we study the convergence when h t and h x tend to 0, we assume in the following that h t ≤ 1 and h x ≤ 1.

Remark 7.2. Since we want to solve the PDE on [0, T ] × R d , we have to choose the interval [-a, a] large enough. Besides, g(2 d+1 T λ(B)f k n (t, x)) is nothing more than a regularizing term. If we had defined P k v(t, x) as r k n (t,x) f k n (t,x) , P k v(t, x) would not have satisfied the regularity property of Hypothesis 2. Now, let us check that P k satisfies the properties of Hypothesis 2. Since we pick new random points at each iteration, the required measurability property is easily satisfied (with

A k := (T k i , X k i ) 1≤i≤n
). The linearity property ensues from the definition of r n and the regularity property comes from the definition of g (see Remark 7.2). Concerning the boundedness, the following Proposition holds Lemma 7.3. For any bounded function v, P k v, ∂ x P k v, ∂ 2 x P k v and ∂ t P k v are bounded by some constants depending on h t , h x , c 0 (v), c 1 (K t ), c 2 (K x ), T and λ(B).

Proof. It is sufficient to note that x -→ g(x)

x is bounded with bounded derivatives up to order 2. To get more details on the bounds, we refer the reader to [START_REF] Labart | BSDE: analysis of the discretisation error and simulation using adaptive Monte Carlo methods; Domain perturbations for American options[END_REF]Proposition 11.8].

The convergence rates are proved in the companion paper [START_REF] Gobet | Solving BSDE with adaptive control variate[END_REF] and are stated in the next theorem.

Theorem 7.4. Assume that the coefficients σ, b are bounded measurable functions on [0, T ] × R d , Lipschitz w.r.t. x and that σ satisfies the ellipticity condition. The sequence (P k ) k satisfy the Hypothesis 2 with:

ǫ 1 (P) = K 0 (T )(h 2 t + h 2 x + T λ(B)δ n h 2 x ), ǫ 2 (P) = K 0 (T )(h t + e -µa a d-1 h x + e -µa √ d + T λ(B)δ n h 2 x ), ǫ 3 (P) = ǫ 1 (P) + ǫ 2 (P), ǫ 4 (P) = K 0 (T ) T λ(B)δ n h 2 x , c 4 (P) = K 0 (T ) h 2 x .
8. Numerical results.

8.1. Choice of the parameters a, n, h x , h t , M and N and complexity of the algorithm. By using the operator P described in Section 7, the constants η and ǫ defined in Theorem 5.1 become

η = K(T ) 4(1+T )L 2 f β + K 0 (T )(1 + K 3.10 (T )L 2 f )(h 2 t + h 2 x + T λ(B)δn h 2 x (1 + M -1 ) and ǫ = K 0 (T )(c 2 1,2 (u) + c 2 0 (f ) + c 2 0 (Φ)) h t + h 2 x + T λ(B)δn h 2 x (1 + (M N ) -1 ) + e -µa a d-1 hx + e -µa √ d + 1 N 2 h 2 x
. From this, we deduce the assumptions to be imposed on the coefficients to get a converging algorithm:

h 2 t ≪ 1, h 2 x ≪ 1, T λ(B)δ n ≪ h 2 x , N -2 ≪ h 2
x , a large enough and M as small as we want. In particular, it means n ≫ T λ(B)h

-1 t h -d-2
x . The complexity of Algorithm 2 is of order K max nM (n + N ), where K max denotes the number of iteration steps. Numerically, one can choose N ≪ n, then the complexity is equivalent to K max M n 2 . We don't need to choose M large, however getting an accurate approximating operator requires a quite large value of n. The limiting factor of the complexity is undoubtedly n.

Application to financial problems.

In this Section we present two applications of Algorithm 2 to financial problems. As presented in [START_REF] Karoui | Backward Stochastic Differential Equations in Finance[END_REF], BSDEs appear in numerous financial problems, like the pricing and the hedging of European and American options.

First, in order to check the convergence of the algorithm, we consider a problem for which the exact solution is known. Let us consider the BSDE (1.1)-(1.2) in dimension d = 1 and with the following parameters: b(t, x) = µ 0 -σ 2 2 , σ(t, x) = σ, f (t, x, y, z) = -ry -µ0-r σ z and Φ(x) = (e x -K) + . Although the driver and the terminal condition are not bounded and Φ is not C 2+α , it seems that the convergence of our algorithm still holds.

We explicitly know the value of the solution (Y, Z) of the BSDE: Y t = F (t, e Xt ) and Z t = ∂ x F (t, e Xt )σe Xt , where F is the price function of a standard call option in the Black and Scholes model.

To bring out the result of Theorem 5.1, we compute in Table 8. divided by 100 (resp. by 4) between the first and the fifth iterations. The huge difference between the two reduction coefficients is due to the fact Z is linked to the "derivative" of Y : it is well known that a function is always better approximated than its derivatives. Moreover, the errors drastically decrease between the first and the second iterations. From iteration 2, the algorithm does not improve so much the result anymore. Figure 8.2 represents the level-sets of the pointwise error on Y at iteration 10. Time goes from 0 to 1 and space varies between 4.5 and 4.81 (which means that the starting point for e X belongs to [90,120]). We notice that the error is quite small, except around the point t = 1 and e x = 100. This corresponds to the fact that at maturity time, the solution Y T is equal to Φ(X T ), which is continuous

(but not C 1 ) in (1, K).
Second, let us consider the pricing of a basket call option in dimension 3 in the Black and Scholes model with the following parameters: b(t, x) = rx, σ(t, x) = σx, f is kept the same as above, and Φ(x) = ( x1+x2+x3 3 -K) + . We aim at computing the price at time 0 and for X 0 = (100, 100, 100). The reference price is given by the approximated formula given in [START_REF] Carmona | Generalizing the Black-Scholes formula to multivariate contingent claims[END_REF]. The option parameters are µ 0 = r = 0.02, σ = 0.2, T = 1, β = 0 and µ = 1. The algorithm parameters are N = M = 100, K max = 10 and 2a = 127 for each space component, i.e. we solve the BSDE on [0, T ] × [55, 182] 3 . As above, we use a truncated Gaussian kernel. Figure 8 on the wealth or portfolio processes. For example, we can use it to hedge with higher interest rate for borrowing, which boils down to solving a non linear BSDE. We can also deal with the pricing and hedging of American options by using a penalization method. Since pricing American options is equivalent to solving a reflected BSDE, we approximate the solutions of the RBSDEs by a sequence of standard BSDEs with penalizations (see [START_REF] Karoui | Reflected solutions of backward SDE's, and related obstacle problems for PDE's[END_REF] for more details). We refer to [START_REF] Labart | BSDE: analysis of the discretisation error and simulation using adaptive Monte Carlo methods; Domain perturbations for American options[END_REF]Chapter 15] for more precisions on these applications, and for other financial applications.

Appendix A. Properties of the transition density of a diffusion process. Proposition A.1 (Aronson [START_REF] Aronson | Bounds for the fundamental solution of a parabolic equation[END_REF]). Assume that the coefficients σ and b are bounded measurable functions and that σ is elliptic. There exist positive constants K A.1 (T ), α 1 , α 2 s.t. for any s, t ∈ [0, T ](s > t) and any x, y ∈ R

d K -1 A.1 (T ) (2πα 1 (s -t)) d 2 e -|x-y| 2 2α 1 (s-t) ≤ p(t, x; s, y) ≤ K A.1 (T ) (2πα 2 (s -t)) d 2 e -|x-y| 2 2α 2 (s-t) .
The constant K A.1 (T ) depends on d and on the suprema of the coefficients σ, b, and α 1 and α 2 depend on σ 0 and σ 1 .

Proposition A. < ∞.

For m = 0, H µ := H 0,µ .

• Space L 2 (0, T ; F ). Let (F, • F ) be a Banach space. L 2 (0, T ; F ) defines the space of functions φ from [0, T ] into F s.t. φ 2 L2(0,T ;F ) = T 0 φ(t) 2 F dt < ∞. The following result is a direct application of Bensoussan and Lions [5, Theorem 6.12 page 130].

Theorem B.2 (Bensoussan and Lions [START_REF] Bensoussan | Applications des inéquations variationnelles en contrôle stochastique[END_REF]). Assume σ is uniformly elliptic, σ is in C 1,1 b , b is in C 1,1 , b is bounded and k is bounded from below. We also assume that g ∈ L 2 (0, T ; H µ ). Then, the solution v of (∂ t + L)v(t, x) + k(t, x)v(t, x) + g(t, x) = 0, with terminal condition v(T, •) = 0, is in L 2 (0, T ; H 2,µ ) and ∂ t v ∈ L 2 (0, T ; H µ ). Furthermore, we have Proof. First, we do the proof in the special case Π ≡ 0. Assume without loss of generality that t < t

(λt ′ + (1 -λ)t)) √ s -t ′ ) -1 dλ ≤ K(T ) √ t ′ -t -1 . T t ′ ds (s-(λt ′ +(1-λ)t)) √ s-t ′ = 2 √ T -t ′ 0 du u 2 +(1-λ)(t ′ -t) = 2 √ (1-λ)(t ′ -t) q T -t ′ (1-λ)(t ′ -t) 0 dv 1+v 2 . Thus, T t ′ ds 1 0 1 (s-(λt ′ +(1-λ)t)) √ s-t ′ dλ ≤ π √ t ′ -t 1 0
dλ √ 1-λ , and the result follows. In the general case, we introduce v 0 (t, x) := v(t, x) -Π(x). v 0 satisfies the PDE (∂ t + L)v 0 (t, x) + LΠ(x) + g(t, x) = 0, with terminal condition v 0 (T, •) = 0. Then, the first part of the Proposition yields that for all t, t ′ ∈ [0, T ], it holds |∂ x v 0 (t ′ , x) -∂ x v 0 (t, x)| ≤ c 1/2 (v 0 ) |t ′ -t|, where c 1/2 (v 0 ) is of the form K(T )(c 0 (g) + c 2 (Π)). Since ∂ x v 0 (t ′ , x) -∂ x v 0 (t, x) = ∂ x v(t ′ , x) -∂ x v(t, x), we get the result.

, 2 b| 2

 22 scheme to get Y -Y k 2 µ,β ≤ c 2 0,1 (u) K(T ) N + 2E uu k

2 µ,β is postponed to the next step. Proposition 6 . 3 .

 263 For any β > 0 and any µ > 0, it holds

Ỹ k 2 µ,β and Z k -Zk 2 µ,β . 6 . 3 .Ỹ k 2 µ,β and Z k -Zk 2 µ

 226322 1 ds. By integrating w.r.t. e -µ|x| dx and taking the expectation, the result follows.It remains to measure the error Yk -Third source of error: approximation operator P. In this part, we split the remaining errors Y k -,β in three error terms (see(6.4) for the decomposition of Y ). The first two are approximation errors connected to P and the last one is studied in the next section. Proposition 6.4. It holds

1 √Fig. 8 . 1 .

 181 Fig. 8.1. Evolution of e(Y k -Y ) and e(Z k -Z) w.r.t. k

|∂ 2

 2 tx p(t, x; s, y)| ≤ K A.3 (T )(st) -d+3 2 e -cA.3 |x-y| 2 s-t , where c A.3 depends on σ 0 , σ 1 and K A.3 (T ) depends on d and on the suprema of σ, b . Proposition A.4 (Gobet and Labart[START_REF] Gobet | Sharp estimates for the convergence of the density of the Euler scheme in small time[END_REF]). Assume σ is uniformly elliptic, b and σ are inC 1,3 b and ∂ t σ is in C 0,1 b . Then, ∀(s, x, y) ∈ [0, T ] × R d × R d ,there exist two positive constants c A.4 and K A.4 (T ) s.t. |p(0, x; s, y)p N (0, x; s, y)| ≤ K A.4 (c A.4 depends on σ 0 , σ 1 and K A.4 (T ) depends on the dimension d, on the upper bounds of σ, b and their derivatives. Proposition A.5 (Gobet and Labart [12]). Assume σ is uniformly elliptic and b and σ are in C 0,2 b . Then, for any x, y ∈ R d , s ∈ [0, T ], there exist two positive constants c A.5 and K A.5 (T ) s.t. p N (0, x; s, y) ≤ K A.5 (T ) s d/2 e -cA.5 |x-y| 2 s , where c A.5 depends on σ 0 , σ 1 and K A.5 (T ) depends on d, on the upper bounds of b, σ and their derivatives. Appendix B. Bounds for linear PDEs in weighted Sobolev spaces. Definition B.1 (Space H m,µ and L 2 (0, T ; F )). • For any m ≤ 2 and µ > 0, H m,µ defines the space of functions v : R d → R s.t. v H m,µ =   k≤m R d e -µ|x| |∂ k v(x)| 2 dx

∂

  t v L2(0,T ;H µ ) + v L2(0,T ;H 2,µ ) ≤ K B.2 (T ) g L2(0,T ;H µ ) . Proposition B.3. Let v be the solution of (∂ t + L)v(t, x) + g(t, x) = 0, with terminal condition v(T, •) = Π(•). Assume g ∈ C 0,0 b , Π ∈ C 2 b and σ and b are in C 1,2 b . Then, for all t, t ′ ∈ [0, T ], it holds |∂ x v(t ′ , x) -∂ x v(t, x)| ≤ c 1/2 (v) |t ′ -t|, where c 1/2 (v) is of the form K(T )(c 0 (g) + c 2 (Π)).

  This sequence of linear BSDEs can be linked to a sequence of linear PDEs: by writing Ŷ k t = ûk (t, X t ) and Ẑk t = ∂ x ûk (t, X t )σ(t, X t ), one has

  Norm equivalence results. The following proposition gives a norm equivalence result between • H µ β,X and • H µ β . It can be compared to Bally and Matoussi[2, Proposition 5.1]. Proposition 3.8 (Norm equivalence). Assume that the coefficients σ, b are bounded measurable functions on [0, T ] × R d , Lipschitz w.r.t. x and that σ satisfies the ellipticity condition. There exist two constants c > 0 and c > 0 (depending on

	3.3.
	2 H m,µ β, X , we compute the expectation w.r.t. the law of X. If v is a random function, v 2 H m,µ β, X is random too.

  us do the proof for the error on Y .

										Using
	a	priori	estimates	(see	[10,	Proposition	2.1])	and	the
	standard	BSDE	(1.1)	leads	to

  .[START_REF] Bally | Error analysis of the optimal quantization algorithm for obstacle problems[END_REF] shows the evolution of Y 0 Evolution of the pointwise error on Y (i.e the price of the option at time 0) w.r.t. the number of iterations k, when n = 250, n = 500, n = 1000 and n = 2000. The values of h t and h x depend on the value of n. They are computed with the following formula: h t = ( T n ) 1/3 and h x = 2a( T n ) 1/3 . We notice that the larger is n, the better is the approximation and the faster is the convergence.This algorithm also enables to price and hedge contingent claims with constraints

			-0.06 0.43			-0.54
			-0.54							0.43
		4.78 4.78
		4.74 4.74
					-1.03
			-0.06	-0.54							-0.54
		4.70 4.70	-1.03
										-0.06
	space	4.66 4.66	0.43			-0.54	0.43
		4.62 4.62				-0.06		-0.06				0.43 0.91 1.40 1.88 2.37 2.86 3.34
		4.58 4.58
		4.54 4.54
		0 0	0.1 0.1	0.2 0.2	0.3 0.3	0.4 0.4	0.5 0.5	0.6 0.6	0.7 0.7	0.8 0.8	0.9 0.9	1.0 1.0
							time
		1 Fig. 8.2. ref price 2 5.5 6.0 6.5 7.0 7.5		3	4	5	6	app. price n=1000 7 8		9	10
			app price n=250				app price n=2000
			app price n=500

Fig. 8.3. Evolution of Y 0 w.r.t. the iterations

  2. [Ladyzenskaja et al.[START_REF] Ladyzenskaja | Linear and quasi-linear equations of parabolic type[END_REF], pages 376-377] Assume that σ is uniformly elliptic and that the coefficients σ, b are bounded, Hölder continuous of order α in x and α 2 in t. There exist two positive constants c A.2 (depending on σ 0 , σ 1 ) and K A.2 (T ) (depending on d, α and on the suprema of σ, b), s.t. for any s, t ∈ [0, T ] (s > t) and any x, y ∈ R d ≤ |a| + |b| ≤ 2, m = 0, 1, stated by Friedman (see[11, page 261]). We refer to [18, Proposition 6.9] for a detailed proof.Corollary A.3. Assume that σ and b are in C 1,2 b and that σ is uniformly elliptic. It holds for any s, t ∈ [0, T ] (s > t) and any x, y ∈ R d

	|∂ r t ∂ q x p(t, x; s, y)| ≤ K A.2 (T )(s -t) -(d+2r+q) 2	e -cA.2 |x-y| 2 s-t , where 2r + q ≤ 2.
	The following Corollary bounds ∂ 2 tx p(t, x; s, y). It ensues from a result bounding ∂ m+a x ∂ b y p(t, x; s, y), 0

  ′ . Feynman-Kac's formula gives v(t, x) = E[ T t g(s, X s )ds]. Using the transition density function of X s leads to ∂ xi v(t, x) = R d dy T t g(s, y)∂ xi p(t, x; s, y)ds. Then,∂ xi v(t ′ , x) -∂ xi v(t, x) = (s, y)[∂ xi p(t ′ , x; s, y) -∂ xi p(t, x; s, y)]ds , y)∂ xi p(t, x; s, y)ds := A 1 -A 2 .Let us establish upper bounds for A 1 and A 2 .Upper bound for A 2 . Applying Proposition A.2 yields|A 2 | ≤ K(T )c 0 (g) R d dyUpper bound for A 1 . Applying Taylor's formula gives∂ xi p(t ′ , x; s, y) -∂ xi p(t, x; s, y) = (t ′t) |A 1 | ≤ K(T )c 0 (g)(t ′ -t) R d dy -(λt ′ + (1λ)t)) -3 2 dλ. Since λt ′ + (1λ)t ∈ [t, t ′ ], we can write (s -(λt ′ + (1λ)t)) -3 2 ≤ ((s -(λt ′ + (1λ)t)) √ st ′ ) -1 . It remains to prove that

				T
			dy
		R d	t ′
			t ′
		-
	t ′ t	ds (s-t) (d+1) 2	e -c |x-y| 2 s-t	≤ K(T )c 0 (g) √	t ′ -t.
	Then, K(T )c 0 (g)(t ′ -t) 0 (s T we T t ′ ds 1 t ′ ds 1 0 ((s -	T t ′ ds	1 0 ∂ 2 tx p(λt ′ + (1 -λ)t, x; s, y)dλ. Corollary A.3: 0 (s-(λt ′ +(1-λ)t)) -d+3 use 1 2 e -c |x-y| 2 (s-(λt ′ +(1-λ)t)) dλ ≤

gt g(s

Appendix C. Technical results. Lemma C.1.

, where c > 0. For any s, t ∈ [0, T ] and any x, y ∈ R d such that t < s, the following assertion holds

Proof. Using a change of variables in

. The last term being bounded, we get the upper bound for I. Since E[e

, the same proof gives the lower bound.