Andreia S Margarido

Ludovic Le Guen

Amandine Falco

Sandrine Faure

Norbert Chauvet

Pascal De

Santa Barbara

PROX1 is a specific and dynamic marker of sacral neural crest cells in the chicken intestine

Keywords: Enteric nervous system, sacral neural crest cell, nerve of Remak, colon, Prospero-related homeobox 1, chick embryo, pelvic plexus, RRID: AB_10063408, RRID: AB_63793, RRID: AB_221448, RRID: AB_10013821, RRID:AB_221568, RRID:AB_2762834, RRID:AB_141607, RRID:AB_2536183

The Enteric Nervous System (ENS) is a complex network constituted of neurons and glial cells that ensures the intrinsic innervation of the gastrointestinal tract. ENS cells originate from vagal and sacral neural crest cells that are initially located at the border of the neural tube. In birds, sacral neural crest cells (sNCCs) first give rise to an extramural ganglionated structure (the so-called Nerve of Remak (NoR)) and to the pelvic plexus. Later, sNCCs enter the colon mesenchyme to colonize and contribute to the intrinsic innervation of the caudal part of the gut. However, no specific sNCC marker has been described. Here, we report the expression pattern of Prospero-related homeobox 1 (PROX1) in the developing chick colon.

PROX1 is a homeobox domain transcription factor that plays a role in cell type specification in various tissues. Using in situ hybridization and immunofluorescence techniques, we showed that PROX1 is expressed in sNCCs localized in the NoR and in the pelvic plexus.

Then, using real-time quantitative PCR we found that PROX1 displays a strong and highly dynamic expression pattern during NoR development. Moreover, we demonstrated using in vivo cell tracing, that sNCCs are the source of the PROX1-positive cells within the NoR. Our results indicate that PROX1 is the first marker that specifically identifies sNCCs. This might help to better identify the role of the different neural crest cell populations in distal gut innervation, and consequently to improve the diagnosis of diseases linked to incomplete ENS formation, such as Hirschsprung's disease.

-INTRODUCTION

The gastrointestinal (GI) tract is a complex organ and its functioning relies on the coordinated work of many different tissue types. Indeed, GI smooth muscle, which is innervated by the enteric nervous system (ENS), allows the movement of the ingested food, while the intestinal epithelium and vascular system play a key role in nutrient absorption (Le [START_REF] Guen | Mesenchymal-epithelial interactions during digestive tract development and epithelial stem cell regeneration[END_REF].

Throughout development, the establishment and differentiation of all GI components are precisely regulated. The ENS, which is a complex network of neurons and glial cells, derives from Neural Crest Cells (NCCs) that colonize the immature GI tract from adjacent tissues [START_REF] Furness | The Enteric Nervous System[END_REF][START_REF] Uesaka | Development of the intrinsic and extrinsic innervation of the gut[END_REF]. Specifically, during development, NCCs located at the boundary of the neural plate and non-neural ectoderm undergo epithelial-mesenchymal transition, delaminate from the ectoderm, and migrate through specific pathways to reach their targets where they finally differentiate (Le [START_REF] Douarin | The Neural Crest (2nd edition)[END_REF]. Bird models have been used for decades to study NCC migratory routes (for review, [START_REF] Heanue | Enteric nervous system development in avian and zebrafish models[END_REF].

Experiments with chick-quail chimaeras have identified two NCC types that form the ENS of the GI tract: vagal NCCs (vNCCs) and sacral NCCs (sNCCs). vNCCs migrate ventrally from the neural tube at the level of somites 1-7 on embryonic day 3 (E3) [START_REF] Yntema | The origin of intrinsic ganglia of trunk viscera from vagal neural crest in the chick embryo[END_REF][START_REF] Douarin | The migration of neural crest cells to the wall of the digestive tract in avian embryo[END_REF][START_REF] Espinosa-Medina | Dual origin of enteric neurons in vagal Schwann cell precursors and the sympathetic neural crest[END_REF]. They first accumulate in caudal branchial arches before entering the foregut mesenchyme. In the GI tract, vNCCs migrate rostro-caudally, reaching the posterior end at E8-E8.5 in chicken embryos (Le [START_REF] Douarin | The migration of neural crest cells to the wall of the digestive tract in avian embryo[END_REF][START_REF] Burns | The sacral neural crest contributes neurons and glia to the post-umbilical gut: spatiotemporal analysis of the development of the enteric nervous system[END_REF]. Interestingly, their migratory behavior within the GI mesenchyme changes along the anterior-posterior axis. In the preumbilical gut, vNCCs are evenly distributed within the mesenchyme, while in the postumbilical gut they accumulate in the mesenchyme outermost layer. Finally, they reach the hindgut where they colonize the inner submucosa layer.

In chicken embryos, sNCCs originate in the neural tube from somite 28 onwards and give rise first to the so-called Nerve of Remak (NoR; an extramural ganglionated structure) and the pelvic plexus. Then, they migrate to the colon mesenchyme where they contribute to 20% of hindgut neurons (Le [START_REF] Douarin | The migration of neural crest cells to the wall of the digestive tract in avian embryo[END_REF][START_REF] Serbedzija | Vital dye labelling demonstrates a sacral neural crest contribution to the enteric nervous system of chick and mouse embryos[END_REF][START_REF] Catala | Organization and development of the tail bud analyzed with the quail-chick chimaera system[END_REF][START_REF] Burns | The sacral neural crest contributes neurons and glia to the post-umbilical gut: spatiotemporal analysis of the development of the enteric nervous system[END_REF][START_REF] Burns | Sacral neural crest cells colonise aganglionic hindgut in vivo but fail to compensate for lack of enteric ganglia[END_REF][START_REF] Nagy | Pelvic plexus contributes ganglion cells to the hindgut enteric nervous system[END_REF]. The pelvic plexus is found both in avian and mammalian embryos [START_REF] Nagy | Pelvic plexus contributes ganglion cells to the hindgut enteric nervous system[END_REF][START_REF] Serbedzija | Vital dye labelling demonstrates a sacral neural crest contribution to the enteric nervous system of chick and mouse embryos[END_REF]Kapur, 2010;Wang et al., 2011a). Conversely, the NoR is unique to birds, and runs in the mesentery, dorsal to the gut, from the cloaca to the duodenum [START_REF] Suzuki | Projections of neurons in the intestinal nerve of Remak to the chicken intestine[END_REF]. The NoR sends many nerve branches to the gut derivatives (intestine, cloaca, bursa of Fabricius), and also to other visceral organs, such as ovary and oviduct [START_REF] Liu | Distribution of vasoactive intestinal peptide and substance P messenger ribonucleic acid in intestinal nerve of Remak of chicken[END_REF].

Thus, vNCCs migrate in a rostro-caudal direction to the gut mesenchyme and colonize the entire chick gut by E8.5, whereas sNCCs enter the distal hindgut mesenchyme in a caudorostral direction [START_REF] Burns | The sacral neural crest contributes neurons and glia to the post-umbilical gut: spatiotemporal analysis of the development of the enteric nervous system[END_REF][START_REF] Burns | Enteric nervous system development: analysis of the selective developmental potentialities of vagal and sacral neural crest cells using quailchick chimeras[END_REF]. It was previously suggested that sNCC migration is dependent on vNCC colonization. However, in the absence of vNCCs after vagal neural tube ablation, sNCC migrate in an apparently normal manner [START_REF] Burns | Enteric nervous system development: analysis of the selective developmental potentialities of vagal and sacral neural crest cells using quailchick chimeras[END_REF]. Conversely, when the sacral neural tube is transplanted to the vagal area, sNCC migrate rostro-caudally along the gut and stop migrating at the level of the caecum before the colon [START_REF] Burns | In ovo transplantation of enteric nervous system precursors from vagal to sacral neural crest results in extensive hindgut colonisation[END_REF]. As the hindgut is the GI region most frequently affected by developmental disorders linked to NCC failure to colonize the gut, such as aganglionic megacolon in animals and Hirschsprung's disease in humans, the identification of intrinsic molecular differences between vNCC and sNCC is of great interest (Kapur, 2001;[START_REF] Burns | Enteric nervous system development: analysis of the selective developmental potentialities of vagal and sacral neural crest cells using quailchick chimeras[END_REF]. Key ENS genes are expressed in both NCC types without qualitative differences [START_REF] Delalande | The receptor tyrosine kinase RET regulates hindgut colonization by sacral neural crest cells[END_REF]. Only RET expression is 4-fold higher in vNCCs than in sNCCs. Moreover, RET overexpression in sNCCs leads to the earlier migration to the gut of a larger number of cells [START_REF] Delalande | The receptor tyrosine kinase RET regulates hindgut colonization by sacral neural crest cells[END_REF]. However, no specific marker for vNCCs or sNCCs has been identified so far.

Using a microarray approach to identify candidate genes involved in GI development [START_REF] Guen | Intermuscular tendons are essential for the development of vertebrate stomach[END_REF] that already allowed us to characterize regulators of GI mesenchyme development and differentiation [START_REF] Notarnicola | The RNA-binding protein RBPMS2 regulates development of gastrointestinal smooth muscle[END_REF][START_REF] Mckey | LIX1 regulates YAP1 activity and controls the proliferation and differentiation of stomach mesenchymal progenitors[END_REF], we identified Prospero-related homeobox 1 (PROX1), a homeodomain transcription factor with a key role in cell type specification in various organs and tissues [START_REF] Stergiopoulos | Prospero-related homeobox 1 (Prox1) at the crossroads of diverse pathways during adult neural fate specification[END_REF]. In mouse embryos, Prox1 plays an essential role in organ formation, especially brain, retina, lens, liver, pancreas, and the endothelial lymphatic system [START_REF] Oliver | Prox 1, a prospero-related homeobox gene expressed during mouse development[END_REF][START_REF] Tomarev | Chicken homeobox gene Prox 1 related to Drosophila prospero is expressed in the developing lens and retina[END_REF]Wigle and Oliver, 1999;Wigle et al., 1999;[START_REF] Sosa-Pineda | Hepatocyte migration during liver development requires Prox1[END_REF][START_REF] Dyer | Prox1 function controls progenitor cell proliferation and horizontal cell genesis in the mammalian retina[END_REF][START_REF] Wang | Prox1 activity controls pancreas morphogenesis and participates in the production of "secondary transition" pancreatic endocrine cells[END_REF][START_REF] Lavado | Prox1 expression patterns in the developing and adult murine brain[END_REF]. In the chick neural tube, PROX1 is strongly expressed in interneuron progenitors during development [START_REF] Misra | Prox1 regulates a transitory state for interneuron neurogenesis in the spinal cord[END_REF][START_REF] Kaltezioti | Prox1 regulates the notch1-mediated inhibition of neurogenesis[END_REF].

Moreover, PROX1 regulates cell cycle exit of neural progenitor cells and binary fate decision [START_REF] Kaltezioti | Prox1 regulates Olig2 expression to modulate binary fate decisions in spinal cord neurons[END_REF]. In the present study, we investigated PROX1 temporal and spatial expression in the developing chick colon. We found that PROX1 expression is strong and highly dynamic during the development of the NoR and the pelvic plexus. We also showed by in vivo cell tracing that sNCCs initially located in the neural tube from somite 28 onwards are at the origin of the PROX1-positive cells present in the NoR.

-MATERIALS AND METHODS

-Animals

Timed fertilized white Leghorn eggs (Haas Farm, France) were kept at 38°C in a humidified incubator (Coudelou, France) until use. Embryos were staged according to Hamburger and Hamilton (HH) [START_REF] Hamburger | A series of normal stages in the development of the chick embryo[END_REF]. GI tissues were dissected from chick embryos and staged according to the embryonic day (E) [START_REF] Southwell | Staging of intestinal development in the chick embryo[END_REF]. Dissected GI tissues were fixed in 4% paraformaldehyde at room temperature for 1h, washed in PBS, gradually dehydrated in methanol and stored at -20°C before processing for whole-mount in situ hybridization, as described below.

-Whole-mount in situ hybridization

Whole-mount in situ hybridization was carried out as previously described using antisense SOX10 [START_REF] Moniot | SOX9 specifies the pyloric sphincter epithelium through mesenchymal-epithelial signals[END_REF][START_REF] Bourret | Colonic mesenchyme differentiates into smooth muscle before its colonization by vagal enteric neural crest-derived cells in the chick embryo[END_REF] and PROX1 riboprobes. For this, E6 GI tract total RNA and primers specific for chick PROX1 (Gene Bank ID: NM_001005616.1; forward: 5' -TGT AAA GTT CAA CAG ATG CAT TAC C-3', reverse: 5' -ATG TTA AGG GTC TCG GGC AA-3') were used to amplify chick PROX1 cDNA. The resulting 616 base pair (bp) fragment was cloned in the pGEM-T Easy vector (Invitrogen) and sequenced to verify the identity. GI tissues were gradually rehydrated in PBS, washed in PBS/0.1% Tween (PBT) and incubated in 6% hydrogen peroxide for 1h (Sigma, France). Samples were then permeabilized by treatment with proteinase K (10 µg/ml) for 10 min, washed with glycine in PBT, and fixed in 4% paraformaldehyde/0.2% glutaraldehyde in PBT for 20 min. Then, samples were hybridized with digoxigenin-labeled (Roche) antisense SOX10 or PROX1 riboprobes at 70°C overnight. After washes at 70°C, samples were incubated in 10% sheep serum at room temperature for 2.5h, and finally mixed with a pre-adsorbed anti-digoxigenin antibody coupled with alkaline phosphatase (Roche) at 4°C overnight. Complexes were detected with BM Purple, a chromogenic substrate for alkaline phosphatase (Roche).

-Quantitative real-time RT-PCR

Total RNA was extracted from E5 to E12 chick colons with the HighPure RNA Isolation Kit (Roche). Reverse transcription (RT) was performed with the Verso cDNA Synthesis Kit (Thermo Scientific) and quantitative real time PCR was performed using the LightCycler technology (Roche Diagnostics). PCR primers (Supplementary Table 1) were designed using the LightCycler Probe Design software-2.0. For each developmental stage, each analyzed sample corresponded to total RNA pooled from three individual colons. Three samples were used at each stage to evaluate the expression of PROX1, SOX10 and SM22. Expression levels were determined with the LightCycler analysis software (version 3.5) relative to standard curves. Data were presented as the expression levels of PROX1, SOX10 and SM22 relative to the expression of the reference gene GAPDH. The relative mRNA expression was calculated using the 2 -ΔΔCT method [START_REF] Faure | Enteric neural crest cells regulate vertebrate stomach patterning and differentiation[END_REF][START_REF] Mckey | LIX1 regulates YAP1 activity and controls the proliferation and differentiation of stomach mesenchymal progenitors[END_REF]. Among the eight reference genes tested, GAPDH was selected based on its expression stability using the geNorm software. According to this program, GAPDH M value was 0.24, well below the accepted cutoff for expression stability (i.e., 0.5) [START_REF] Vandesompele | Accurate normalization of real-time quantitative RT-PCR data by geometric averaging of multiple internal control genes[END_REF].

-Antibodies

The anti-TUJ1 antibody (Covance, Princeton, NJ, USA; mouse monoclonal, cat.

#MMS-435P, RRID: AB_10063408) was raised against microtubules derived from rat brain and recognizes an epitope located within the last 15 C-terminal residues of neuron-specific class III β tubulin [START_REF] Lee | Posttranslational modification of class III betatubulin[END_REF]. It cross-reacts with avian class III β tubulin, but not with tubulin in glial cells. The anti-HNK-1 antibody (LabVision/Neo-Markers, CA, USA; mouse monoclonal, cat. #MS-1163-P1, RRID: AB_63793) recognizes a N-linked carbohydrate [START_REF] Abo | A differentiation antigen of human NK and K cells identified by a monoclonal antibody (HNK-1)[END_REF] and is a well-known avian neural crest marker [START_REF] Bronner-Fraser | Analysis of the early stages of trunk neural crest migration in avian embryos using monoclonal antibody HNK-1[END_REF][START_REF] Moniot | SOX9 specifies the pyloric sphincter epithelium through mesenchymal-epithelial signals[END_REF][START_REF] Nagy | Pelvic plexus contributes ganglion cells to the hindgut enteric nervous system[END_REF][START_REF] Faure | Enteric neural crest cells regulate vertebrate stomach patterning and differentiation[END_REF][START_REF] Bourret | Colonic mesenchyme differentiates into smooth muscle before its colonization by vagal enteric neural crest-derived cells in the chick embryo[END_REF].

The anti-HUC/D antibody (Invitrogen; mouse monoclonal, clone 16A11, cat. #A-21271, RRID: AB_221448) specifically recognizes neuronal perikarya [START_REF] Marusich | Hu neuronal proteins are expressed in proliferating neurogenic cells[END_REF][START_REF] Notarnicola | The RNA-binding protein RBPMS2 regulates development of gastrointestinal smooth muscle[END_REF]. The anti-PROX1 antibody (RELIA Tech, Germany; rabbit polyclonal, cat. #102-PA32, RRID: AB_10013821) used in this study was generated against the highly conserved C-terminus of recombinant human PROX1 produced in E. coli (amino acid 546 to 737, 100% identical to avian PROX1). In western blot experiments with E15 chick pancreas-duodenum samples, the anti-PROX1 antibody identified a single band of around 83 kDa that corresponds to the predicted size of chick PROX1 (Supplementary Figure 1a). Moreover, this antibody detected PROX1-positive cells in chick E15 pancreatic epithelium (Supplementary Figure 1b), consistent with a previous report [START_REF] Westmoreland | Pancreas-specific deletion of Prox1 affects development and disrupts homeostasis of the exocrine pancreas[END_REF]. Finally, PROX1 immunostaining was comparable to the results obtained by in situ hybridization with the antisense PROX1 probe.

-Immunofluorescence experiments

Chick gut samples were fixed in 4% paraformaldehyde, followed by gradual incubation in sucrose solution. Then, samples were embedded in OCT and stored at -80°C before cutting 10 µm cryosections collected on poly-d-lysine-coated slides (Thermo Fisher).

Immunofluorescence experiments were performed using antibodies against PROX1 (Reliatech; rabbit, dilution 1/400), TUJ1 (Covance; mouse, dilution 1/500), GFP (ThermoFisher; clone 3E6, mouse, cat. A-11120, RRID: RRID:AB_221568, dilution 1/400) and HNK-1 (NeoMarkers; mouse, dilution 1/200). Nuclei were labeled with Hoechst 33342 (Invitrogen, France). After washing, slides were incubated with Alexa 488 anti-mouse (Thermo Fischer Scientific; cat. #A-21202, RRID:AB_141607, 1:2000 dilution) and Alexa 555 anti-rabbit (Thermo Fischer Scientific; cat. #A32794, RRID:AB_2762834, 1:2000 dilution) or Alexa 647 anti-rabbit (Thermo Fischer Scientific; cat. A-31573, RRID:AB_2536183, 1:2000 dilution) for 30min. Sections were rinsed and mounted with Fluorescent Mounting Medium (Dako, France). Immunofluorescence control experiments were performed without primary antibodies (data not shown). Cells were imaged with a Zeiss AxioVision fluorescence microscope using standard filters, or a ZEISS LSM800 confocal laser scanning microscope. PROX1-positive cells in the developing NoR were quantified with the Fiji image analysis software. Image background was removed using on a "rolling ball" algorithm, and a median filter was applied followed by thresholding. Images were converted to a binary format and the areas of interest were selected and measured (TUJ1-positive area; µm2). PROX1 cells were then counted in such areas using the "analyze particle" plugin. The percentage of PROX1-positive cells per TUJ1-positive area (µm2) was calculated using the Prism7 software.

-In ovo electroporation to target sNCCs

Chick eggs were incubated for about 2.5 days until HH16 when integrative plasmids (pCAGGS-T2TP and pT2K-CAGGS-EGFP) (Wang et al., 2011b) were injected in the caudal neural tube (somite 25) using a micro-injector and the site was identified with Phenol Red Solution (Sigma, ref.: P0290). After microinjection, gold plated electrodes (Electroporation System ECM 830, BTX) were placed on both sides of the caudal neural tube and electroporation was performed bilaterally at 25V (3 pulses of 5ms) [START_REF] Fort | Activity of the RhoU/Wrch1 GTPase is critical for cranial neural crest cell migration[END_REF]. Then, eggs were sealed and incubated up to E6.5. After colon dissection, GFP signal was directly evaluated before immunofluorescence analysis.

-RESULTS

-PROX1 expression in the chick embryonic gastrointestinal tract

Among the genes with high expression level in the developing chick colon compared with the developing stomach identified by RNA microarray assays (Supplementary Table 2), PROX1 was strongly expressed in E5 and E6 chick colon, and then its expression level decreased in E9 colon. Conversely, it was faintly expressed in E6 and E9 stomach (Figure 1a). This pattern is unique compared with the expression of key NCC genes (SOX10 and RET) and markers of differentiated smooth muscle (SM22 and CNN1) (Figure 1a). Analysis of its expression pattern in the whole gastrointestinal tract at E6 by whole-mount in situ hybridization showed that PROX1 was expressed in the dorsal part of the colon and in pancreas (Figure 1b). As the NoR is situated in the dorsal part of the colon, the expression of SOX10, a marker of migrating NCCs, was also investigated by whole-mount in situ hybridization. At E6, SOX10 signal revealed the presence of vNCCs in the gizzard, their advancing front at the level of the small intestine, and the presence of sNCCs in the NoR (Figure 1b), where PROX1 was also expressed. Analysis of PROX1 expression in the developing chick colon using an antibody that recognizes also the Gallus gallus PROX1 protein (Supplementary Figure 1) showed that PROX1 was expressed in the NoR (Figure 1c). Moreover, in E6 colon sections, PROX1 colocalized with TUJ1-positive cells that define the NoR (Figure 1d) and also with HUC/Dpositive cells (a marker of neurons present in the NoR) (Figure 1e). Altogether, our data indicate that PROX1 is expressed by neurons of the NoR that are considered to originate from sNCCs (Le [START_REF] Douarin | The migration of neural crest cells to the wall of the digestive tract in avian embryo[END_REF][START_REF] Serbedzija | Vital dye labelling demonstrates a sacral neural crest contribution to the enteric nervous system of chick and mouse embryos[END_REF][START_REF] Catala | Organization and development of the tail bud analyzed with the quail-chick chimaera system[END_REF].

-PROX1 expression in the developing chick colon

As PROX1 was specifically expressed in the NoR, its temporal expression profile was thoroughly characterized and compared with that of SOX10 (NCC marker) in developing colon samples from E5 to E9 by whole-mount in situ hybridization (Figures 2a,2b). At E5 and E6, SOX10 expression in the colon was restricted to the NoR, thus identifying the sNCC population (Figure 2b). As vNCCs colonize the GI tract along an anterior-posterior migration axis, at E7 SOX10 expression in the small intestine, ceca and colon identified vNCCs. These cells then reached the most caudal part of the colon at E9 (Figure 2b). At E8, SOX10 expression highlighted the pelvic plexus surrounding the cloaca and connected to the NoR (Figure 2b and Supplementary Figure 2b, arrows). PROX1 was strongly expressed in the NoR from E5 (Figure 2a) and particularly at E6. In contrast to SOX10, PROX1 expression in the NoR quickly decreased from E7 onwards. At E8 PROX1 was faintly detected in the NoR, while its expression remained high in the pelvic plexus (Figure 2a and Supplementary Figure 2a, arrows). Moreover, PROX1 defined fibers that are present in the caudal part of the E8 colon and are connected to the NoR (Supplementary Figure 2a, arrowheads). At E9, PROX1 could be detected in the pelvic plexus (Figure 2a, arrow), but not in the NoR.

Quantification of PROX1, SOX10 (NCC marker) and SM22 (smooth muscle marker) expression levels by real-time quantitative RT-PCR (Figures 2c-e) using specific Gallus gallus primers (Supplementary Table 1) showed that SOX10 and SM22 transcript levels progressively increased during colon development from E5 to E12 (Figures 2d,2e). PROX1 was already strongly expressed at E5. Then, its expression increased at E6 and E7, and started to decrease from E8 (Figure 2c).

Immunofluorescence analysis of PROX1 and TUJ1 (NCC marker in the NoR) protein expression in E6 to E9.5 colon transverse sections (Figure 3, Supplemental Figure 3) showed that at E8, TUJ1-positive cells were in the inner layer of the colon mesenchyme. This stage corresponds to the time of colonization by vNCCs (Figure 3, Supplemental Figure 3). In addition, TUJ1-positive cells were detected also in the outer layer of the colon mesenchyme starting from E7.5, which corresponds to the presence of sNCC in the colon (Figure 3, Supplemental Figure 3). During this period, PROX1-positive cells were observed in the NoR, and at E9.5 they could not be detected in the colon mesenchyme and in sNCCs that have already migrated to the mesenchyme (Figure 3; Supplementary Figure 3). Although PROX1 expression was mainly restricted to the NoR, some PROX1-positive cells were located dorsally to the NoR and corresponded to developing lymphatic cells (white arrows, Figure 3). Quantification of PROX1-positive sNCCs in the NoR revealed a strong decrease in their number from E6 (10.7x10 -3 cells per µm 2 ) to E9 (0.86x10 -3 cells per µm 2 ) (Figure 4). These findings show that PROX1 is dynamically expressed in the developing chick colon and accurately defines the population of sNCCs that reside in the NoR and in the pelvic plexus from E6 to E9.

-Origin of the PROX1-positive cells detected in the nerve of Remak

As sNCCs originate from NCCs present in the neural tube adjacent to somite 28 caudalwards, PROX1 spatiotemporal expression pattern was examined also during early development from HH11 (E2) to HH24 (E4.5). Whole-mount in situ hybridization showed that PROX1 was expressed in the developing eye, cranial neural tube, heart, liver, and optic and otic placodes. At HH24, PROX1 was observed in somites and mid-neural tube (Figure 5a). However, at these stages, PROX1 expression could not be detected in the caudal neural tube where sNCCs originate (Figure 5a). Then, analysis of the expression of PROX1 and HNK-1 (a marker of migrating NCCs) by immunofluorescence in caudal colon transversal sections from HH18 to HH23 embryos (Figures 5b-d) showed that at HH18 and HH19, HNK-1-positive sNCCs delaminated from the neural tube and migrated to the embryo body to reach the dorsal part of the colon and form the future NoR (Figures 5b,5c). At HH18 and HH19, PROX1 was expressed in cells localized in the neural tube and in isolated cells that will constitute the lymphatic vessels (Figures 5b,5c), but not in sNCCs (Figures 5b,5c).

Conversely, at HH23, PROX1 was expressed in HNK-1-positive cells that constitute the developing NoR (Figures 5d,5e).

As PROX1 was not expressed in migrating sNCCs that have not reached the NoR yet, sNCCs located at the level of somite 28 and caudal-wards were labeled by bilateral electroporation of integrative GFP plasmids in the neural tube at HH16 before delamination (Figure 6a). At E6.5, most GFP-expressing cells were in the most dorsal part of the colon (Figure 6b). Immunofluorescence analysis of GFP and PROX1 expression in E6.5 colon samples after electroporation showed that GFP-positive cells in the NoR expressed also PROX1 (Figure 6c). Moreover, confocal microscopy analysis indicated that GFP-positive cells in the NoR expressed both PROX1 and HUC/D (Figure 6d), demonstrating that the PROX1-positive cells in the developing NoR are sNCCs from the caudal neural tube.

-DISCUSSION

The ENS derives from NCCs that originate from multiple positions to form a complex network of neurons and glial cells. Among them, sNCCs contribute to the intrinsic innervation of the posterior gut in vertebrates (Le [START_REF] Douarin | The migration of neural crest cells to the wall of the digestive tract in avian embryo[END_REF][START_REF] Serbedzija | Vital dye labelling demonstrates a sacral neural crest contribution to the enteric nervous system of chick and mouse embryos[END_REF][START_REF] Catala | Organization and development of the tail bud analyzed with the quail-chick chimaera system[END_REF][START_REF] Burns | The sacral neural crest contributes neurons and glia to the post-umbilical gut: spatiotemporal analysis of the development of the enteric nervous system[END_REF][START_REF] Burns | Sacral neural crest cells colonise aganglionic hindgut in vivo but fail to compensate for lack of enteric ganglia[END_REF][START_REF] Nagy | Pelvic plexus contributes ganglion cells to the hindgut enteric nervous system[END_REF][START_REF] Serbedzija | Vital dye labelling demonstrates a sacral neural crest contribution to the enteric nervous system of chick and mouse embryos[END_REF][START_REF] Kapur | Colonization of the murine hindgut by sacral crest-derived neural precursors: experimental support for an evolutionarily conserved model[END_REF]Wang et al., 2011a). A failure in NCC colonization of the gut leads to an incomplete formation of ENS (agangliogenesis) in the distal gut, a congenital obstructive disorder named the Hirschsprung's disease (de Santa Barbara et al., 2002;[START_REF] Amiel | Hirschsprung disease, associated syndromes and genetics: a review[END_REF][START_REF] Etchevers | The diverse neural crest: from embryology to human pathology[END_REF]. Unfortunately, the mechanisms leading to this pathology are not fully understood. A better knowledge about NCC development and the contribution of the different NCC populations to the distal gut innervation might improve their diagnosis. Although the origin and pathways of sNCC delamination and colonization are well described in avian and mouse models (Haenue et al., 2016;[START_REF] Kapur | Colonization of the murine hindgut by sacral crest-derived neural precursors: experimental support for an evolutionarily conserved model[END_REF]Wang et al., 2011a), no specific sNCC marker has been described so far (Haenue et al., 2016;[START_REF] Memic | Transcription and Signaling Regulators in Developing Neuronal Subtypes of Mouse and Human Enteric Nervous System[END_REF].

Using microarray, in situ hybridization and immunofluorescence approaches, we found that PROX1 is expressed in the NoR located dorsally to the colon and small intestine from E5 to E8. PROX1 is downregulated when sNCCs migrate to the colon/small intestine around E7.5-E8. As in situ hybridization and immunofluorescence approaches are not quantitative techniques, we quantified by RT-PCR the relative abundance of PROX1 transcripts during colon development. We confirmed PROX1 dynamic expression and its downregulation at E8-E9, when sNCCs migrate to the colon mesenchyme. Interestingly, at these stages, no expression of PROX1 was detected in SOX10-positive vNCCs that were migrating anteriorposteriorly. The integrative expression of GFP in sNCCs before their migration to the NoR confirmed that PROX1-expressing cells found in the NoR at E6 derive from sNCCs. In vertebrates (birds, humans and mice), sNCCs form the pelvic plexus that contribute to the hindgut ENS [START_REF] Serbedzija | Vital dye labelling demonstrates a sacral neural crest contribution to the enteric nervous system of chick and mouse embryos[END_REF][START_REF] Kapur | Colonization of the murine hindgut by sacral crest-derived neural precursors: experimental support for an evolutionarily conserved model[END_REF][START_REF] Nagy | Pelvic plexus contributes ganglion cells to the hindgut enteric nervous system[END_REF]Wang et al., 2011a).

PROX1 is also expressed in the pelvic plexus at E8 and E9. In vertebrate, PROX1 is implicated in cell type specification in various organs and tissues [START_REF] Stergiopoulos | Prospero-related homeobox 1 (Prox1) at the crossroads of diverse pathways during adult neural fate specification[END_REF], but its role in ENS organization has not been addressed yet. In heterozygous mice in which Prox1 was inactivated by in-frame insertion of the lacZ gene, lacZ expression was observed in the dorsal part of the caudal hindgut [START_REF] Burke | Prox1 is an early specific marker for the developing liver and pancreas in the mammalian foregut endoderm[END_REF], where the pelvic plexus is located (Wang et al., 2011a). Moreover, high resolution in situ hybridization in mouse embryo sections demonstrated Prox1 expression in the dorsal part of the caudal hindgut at E14.4 PROX1 expression in the NoR and also in the pelvic plexus suggests that PROX1 could be a specific marker of sNCCS.

PROX1 is a member of the homeobox transcription factor family and through its homeobox domain can regulate positively or negatively the transcription of target genes [START_REF] Stergiopoulos | Prospero-related homeobox 1 (Prox1) at the crossroads of diverse pathways during adult neural fate specification[END_REF]. For instance, Kaltezioti and colleagues showed that PROX1 represses NOTCH1 signaling in the chick and mouse neural tube [START_REF] Kaltezioti | Prox1 regulates the notch1-mediated inhibition of neurogenesis[END_REF]. As NOTCH1 signaling inhibits neurogenesis, PROX1 repression of NOTCH1 activity leads to neural progenitor cell differentiation. NOTCH1 and SOX10 maintain the ENS progenitor pool and later promote neurogenesis and/or gliogenesis through additional regulatory mechanisms at different developmental stages. Our results show that SOX10 is expressed in the NoR, in migrating sNCCs and in vNCCs from E5 to E9. Conversely, PROX1 is only expressed in the NoR between E5 and E8 before sNCC migration to the colon/small intestine. PROX1 could have a role in the mechanisms that regulate the dual function of SOX10 and NOTCH1. It could be hypothesized that in the NoR, SOX10 and NOTCH1 maintain the sNCC pool of progenitors before PROX1 expression. Upon PROX1 expression, NOTCH1 starts to be downregulated and this induces neuronal differentiation. Then, PROX1 downregulation might lead to upregulation of NOTCH1 that together with SOX10 can now induce gliogenesis.

However, although NOTCH1 can also contribute to NCC migration when PROX1 is still expressed in the NoR, some sNCCs are already migrating to the colon. Another pathway involved in ENS development is the HEDGEHOG signaling cascade that has a role in NCC migration and could inhibit the NOTCH pathway [START_REF] Liu | Hedgehog and Notch signaling in enteric nervous system development[END_REF][START_REF] Nagy | Sonic hedgehog controls enteric nervous system development by patterning the extracellular matrix[END_REF].

New studies are now needed to identify PROX1 downstream targets, which will also lead to the discovery of more sNCC-specific markers, providing new tools for the characterization of colon ENS-associated pathologies. Altogether, our results have identified a novel NCC marker that allows the distinction between vNCCs and sNCCs. Whole-mount in situ hybridization analysis of PROX1 (a) and SOX10 (b) in E5 to E9 chick colon samples. Scale Bars: 1 mm. PROX1 expression is detected in the NoR from E5 to E8, but not at E9. In addition, PROX1 expression is observed in the pelvic plexus (PP) at E8 and E9. SOX10 expression is restricted to the NoR at E5 and E6. Moreover, SOX10 is also observed in migrating vNCCs located in the small intestine and ceca at E7 and in the colon and pelvic plexus (PP) at E8 (arrowheads). At E9, the whole colon is completely colonized by NCCs and SOX10 is exclusively expressed by the glial lineage. Relative mRNA expression levels of PROX1 (c), SOX10, (d) and SM22 (e) in pooled E5 to E12 colon samples determined by quantitative RT-PCR. Data shown are the mean ± SD expression levels relative to GAPDH. Each measurement was done using total RNA pooled from three individual colons (three pools for each stage). PROX1 mRNA levels is highest at E6. SOX10 and SM22 expression levels progressively increase during colon development. (data from Figure 3). The number of cells in which PROX1 is detected in the nucleus was determined relative to the TUJ1-positive cell area in the NoR (n=2-4 colon/small intestine sections/stage). The number of PROX1-positive cells progressively decreased over time: 10.7x10 -3 cells/µm 2 at E6, 6.19x10 -3 cells/µm 2 at E7, 3.89x10 -3 cells/µm 2 at E8, and 0.86x10 -3 cells/µm 2 at E9. Data are the mean ± SD.

SUPPLEMENTAL FIGURE 3

Spatial and temporal PROX1 expression in the chick developing colon. Immunofluorescence analysis of PROX1 and TUJ1 expression in E6.5 to E9.5 chick colon transverse cryosections (focus on the NoR). TUJ1 is observed in sNCCs located in the NoR from E6.5 to E9.5. At E7.5, few TUJ1-positive sNCCs entered the outer part of the colon mesenchyme. At E6.5, PROX1 is observed in the majority of cells present in the NoR. PROX1 expression in these cells strongly decreased from E7.5 to E8.5 to finally be present only in some isolated cells at E9.5. PROX1 expression was also observed in lymphatic endothelial cells (white arrows) close to the NoR. Scale bars: 100 µm.

FIGURE 1

 1 FIGURE 1 Spatial and temporal expression of PROX1 mRNA during embryonic

FIGURE 2

 2 FIGURE 2 Spatial and temporal expression of PROX1 mRNA in the chick developing colon.

FIGURE 3

 3 FIGURE 3Spatial and temporal expression of PROX1 protein in the chick developing colon.

FIGURE 4

 4 FIGURE 4 Quantification of the number of PROX1-positive cells in the NoR at E6 to E9

FIGURE 5

 5 FIGURE 5 PROX1 expression in early developing chick colon. PROX1 mRNA expression

FIGURE 6

 6 FIGURE 6 Origin of PROX1-expressing cells present in the NoR. Electroporation procedure

ACKNOWLEDGMENTS

The authors thank all members of the INSERM Team "Development of visceral smooth muscle and Associated Pathologies" for technical help and comments. The authors thanks Alan Burns (UCL, London, UK) for fruitful discussion and advises.

Funding information LLG was supported by an ARC post-doctoral fellowship. ASM was supported by an Eramus fellowship (University of Lisbon -University of Montpellier). Research was supported by a Trampolin grant from the Association Française contre les Myopathies (AFM) to SF (N°15681), grants from AFM (N°18766 & N°20749) and the French Patients' Association POIC to PdSB. ORCID ://orcid.org/0000-0001-9040-2481

DATA AVAILABILITY STATEMENT

All plasmids developed during this study are available upon request.

CONFLICT OF INTEREST

The authors have nothing to disclose and no conflict of interest.

AUTHOR CONTRIBUTIONS

ASM and LLG performed in situ hybridization, immunofluorescence and quantitative RT-PCR experiments. LLG and PdSB performed electroporation experiments. AF analyzed the immunofluorescence data. SF and NC performed the microarray screen and analysis. All authors contributed to the writing and critical revision of the manuscript. stereomicroscope (b). Electroporated sNCCs are observed in the E6.5 NoR located dorsally to the colon (Col.). Immunofluorescence analysis of PROX1 and GFP expression in E6.5 colon transverse cryosections after electroporation at HH16 to visualize the NoR (c). Confocal immunofluorescence analysis of PROX1, HUC/D and GFP expression in E6.5 colon transverse cryosections after electroporation at HH16 to visualize the NoR (d). Cells that coexpress GFP, PROX1 and HUC/D are observed in the NoR. Epith: epithelium.

SUPPLEMENTAL TABLE 1

Gallus gallus primers used for RT-qPCR.

SUPPLEMENTAL TABLE 2

Microarray data used for producing the Heatmap (Figure 1a) that summarizes the expression of PROX1, SM22, CNN1, SOX10 and RET in chick colon (from E5, E6, and E9) and gizzard (E6 and E9) samples, with their corresponding p values.

Supplementary


L Galllimard

Adaptive reduced basis strategy for rare events simulations †

Keywords: Structural reliability, Finite element analysis, model reduction, reduced basis, error bounds, failure probability, Rare events simulations, cross-entropy

Monte-Carlo methods are well suited to characterize events of which associated probabilities are not too low with respect to the simulation budget. For very seldom observed events, these approaches do not lead to accurate results. Indeed, the number of samples are often insufficient to estimate such low probabilities (at least 10 +2 samples are needed to estimate a probability of 10 -with 10% relative deviation of the Monte-Carlo estimator). Even within the framework of reduced order methods, such as a reduced basis approach, it seems difficult to accurately predict low probability events. In this paper we propose to combine a cross-entropy method with a reduced basis algorithm to compute rare events (failure) probabilities.

INTRODUCTION

Many applications in structural analysis require taking into account stochastic properties of material, geometry or loads. Given a probabilistic description of the stochastic properties of the structure (i.e. a random vector Θ associated with a probability density function ), reliability analysis aims at computing the probability of failure of structures with respect to a prescribed limit state function . The failure is defined as an event = { (Θ) ≤ 0} and the failure probability is defined by = ({ (Θ) ≤ 0}) = ∫ ( )≤0 ( ) . The computation of this failure probability by a standard Monte-Carlo algorithm is too computationally expensive in practice (10 +2 finite element computations are need to estimate a probability of failure = 10 -). In order to reduce the number of simulations runs different alternatives have been proposed.

A first approach is the First Order Reliability Method (FORM) and the Second Order Reliability Method (SORM) [START_REF] Ditlevsen | Structural reliability methods[END_REF][START_REF] Lemaire | Structural reliability[END_REF] which consists in building a simple analytical approximation of the limit-state function around the so-called design point followed by a direct estimation of the failure probability [START_REF] Hasofer | An exact and invariant first order reliability format[END_REF][START_REF] Rackwitz | Structural reliability under random load sequences[END_REF][START_REF] Der Kiureghian | Efficient algorithm for second-order reliability analysis[END_REF][START_REF] Haukaas | Strategies for finding the design point in non-linear finite element reliability analysis[END_REF][START_REF] Gallimard | Error bounds for the reliability index in finite element reliability analysis[END_REF][START_REF] Gallimard | Coupling finite element and reliability analysis through proper generalized decomposition model reduction[END_REF] . A second approach consists in building a surface response as a surrogate model of the limit state function (Quadratic response surfaces, polynomial chaos expansions, kriging surrogates, ...) [START_REF] Colombi | Stochastic finite elements via response surface: fatigue crack growth problems[END_REF][START_REF] Pendola | Combination of finite element and reliability methods in nonlinear fracture mechanics[END_REF][START_REF] Sudret | Comparison of finite element reliability methods[END_REF][START_REF] Berveiller | Stochastic finite elements: a non intrusive approach by regression[END_REF][START_REF] Blatman | An adaptive algorithm to build up sparse polynomial chaos expansions for stochastic finite element analysis[END_REF] . The Monte-Carlo algorithm can then be applied on this surrogate model. A third approach consists in building a reduced order model (ROM) of the full order finite element model (FOM) [START_REF] Mohan | Multi-element stochastic reduced basis methods[END_REF][START_REF] Boyaval | Reduced basis approach for variational problems with stochastic parameters: application to heat conduction with variable robin coefficient[END_REF][START_REF] Boyaval | Reduced basis techniques for stochastic problems[END_REF][START_REF] Grigoriu | Reduced order models for random functions application to stochastic problems[END_REF][START_REF] Florentin | Adaptive reduced basis strategy based on goal oriented error assessment for stochastic problems[END_REF][START_REF] Haasdonk | Reduced basis methods for parametrized partial differential equations with stochastic influences using the Karhunen Loeve expansion[END_REF] , and using the Monte-Carlo method along with the reduced order model. This approach has been applied successfully in [START_REF] Chen | Accurate and efficient evaluation of failure probability for partial different equations with random input data[END_REF][START_REF] Gallimard | Towards error bounds of the failure probability of elastic structures using reduced basis models[END_REF] to compute the failure probability of structures in the case where is greater than 10 -4 , however in the case of rare events estimation the number of Monte-Carlo simulations needed to estimate accurately increases (10 8 simulations to estimate a failure probability of = 10 -6 ) and the computational cost becomes too important even in the context of the use of reduced order models. In order to reduce the number of samples different techniques have been developed [START_REF] Morio | a survey of rare event simulation methods for static input-output models[END_REF] , in particular the objective of the importance sampling (IS) method is to reduce the variance of the Monte-Carlo estimator. The main idea is to generate the samples Θ 1 , ..., Θ with an auxiliary probability density function ℎ that is able to generate more samples such that (Θ) ≤ 0 than the probability density function (PDF) and then to introduce a weight in the probability estimate to take into account the change in the PDF generating samples. The cross entropy (CE) optimization method [START_REF] Homem-De-Mello | A study on the cross-entropy method for rare-event probability estimation[END_REF] provides an efficient way to approximate this auxiliary probability density function. In [START_REF] Li | An efficient surrogate-based method for computing rare failure probability[END_REF] the authors propose to combine a CE method with a surrogate model of the limit state function and to use an exact computation of the limit state function only when the sample is near the failure limit. In [START_REF] Balesdent | Kriging-based adaptive importance sampling algorithms for rare event estimation[END_REF] the authors develop a CE method with an adaptive surrogate model based on a Kriging interpolation method for estimating the probability of rare events. In [START_REF] Peherstorfer | Multifidelity preconditioning of the cross-entropy method for rare event simulation and failure probability estimation[END_REF] the authors propose to use an a priori fixed hierarchy of surrogate models to accelerate the estimation of the auxiliary probability density function ℎ by the CE algorithm, then compute the IS probability estimate with the FOM.

In this work we develop an adaptive strategy, based on the computation of an error estimator of the limit state function, to simultaneously build the reduced basis and the auxiliary probability density function ℎ. The reduced basis vectors are chosen in order to control the quality of the estimation of the limit state function and the auxiliary PDF is built by minimizing the cross-entropy with respect to an ideal zero-variance distribution. An importance sampling Monte-Carlo simulation is then run with the surrogate model defined by the reduced basis and the auxiliary PDF. The control of the quality of the reduced basis is ensured by the computation of a lower bound and an upper bound of the failure probability. Compared to 24 , 26 , this paper introduces, during the CE optimization, an adaptive construction of the surrogate model based on an error estimator. Compared to [START_REF] Balesdent | Kriging-based adaptive importance sampling algorithms for rare event estimation[END_REF] it introduces surrogates models based on a reduced basis approximation which allows building bounds of the limit state function as well as of the computed Monte-Carlo failure probability.

The paper is organized as follows: Section 2 describes the problem to be solved. In Section 3 the cross-entropy optimization applied to importance sampling is presented. In Section 4, the reduced basis formulation used to solve the problem is presented as well as the technique used to compute bounds on quantities of interest. The reduced basis algorithm for rare events simulations is detailed in Section 5. Finally, Section 6 presents an application to multilayered structures in 2D plane elasticity.

PROBLEM FORMULATION

Linear elastic model

Let us consider a 2D elastic structure defined in a domain Ω bounded by Γ. The external actions on the structure are represented by a surface force density defined over a subset Γ of the boundary such that ∈ 2 (Γ ) [START_REF] Lemaire | Structural reliability[END_REF] and a body force density defined in Ω such that ∈ 2 (Ω) [START_REF] Lemaire | Structural reliability[END_REF] . We assume that a prescribed displacement = such that ∈ 1∕2 (Γ ) 2 is imposed on Γ = Γ -Γ . The material is assumed to be linear elastic, being the Hooke tensor. We consider that the problem is dependent of a vector = ( 1 , ..., ) ∈  ⊂ ℝ of uncertain independent parameters. These parameters are characterized in a probabilistic manner by means of a joint probability density function ( ) = ∏ =1 ( ), where is the probability function of . The problem can be formulated as:

1. For any ∈ , find a displacement field ∈  and a stress field defined in Ω which verify:

• the kinematic constraints:

( , ) = ( ) on Γ

• the equilibrium equations:

div ( , ) + ( , ) = in Ω and ( , ) = ( , ) on Γ (2)
• the constitutive equation:

( , ) = ( , ) ( ( , )) in Ω (3)
denotes the outer normal to Ω.  is the space in which the displacement field is being sought,  0 the space of the fields in  which are zero on Γ , and ( ) denotes the linearized deformation associated with the displacement: [ ( )] = 1∕2 ( , + , ).

Compute a quantity of interest:

( ) = ( ( , ))

where is a linear output of Remark: For the sake of simplicity assume that is a deterministic value (i.e. ( , ) = ( )). For an interested reader, the development of a reduced basis algorithm where is also described by a stochastic approach can be found in [START_REF] Kerfriden | Certification of projection-based reduced order modelling in computational homogenisation by the constitutive relation error[END_REF] .

The weak form formulation of the problem (Equations 1-3) is: find ∈  such that ( ( ), * ; ) = ( * ; ) ∀ * ∈  0

where

( ( ), * ; ) = ∫ Ω ( , ) ( ( , )) ∶ ( * ( )) Ω and ( * ; ) = ∫ Ω ( , ). * ( ) Ω + ∫ Γ ( , ). * ( ) Γ
To compute the solution ( ) of Equation ( 5), a finite element approximation ℎ of is introduced such that ℎ ∈  ℎ ⊂  . Let  ℎ be a partition of Ω into elements ( ∈ {1, ..., }). This partition formed by the union of all elements, is assumed to coincide exactly with the domain Ω and any two elements are either disjoint or share a common edge. We assume that can be represented by a displacement field in  ℎ . The discretized problem is: Find a displacement field ℎ ( ) ∈  ℎ and a stress field ℎ ( ) defined in Ω which verify:

• the kinematic constraints:

ℎ ( , ) = ( ) on Γ (6)
• the finite element equilibrium equations:

∫ Ω ℎ ( ) ∶ ( * ℎ ) Ω = ∫ Ω ( ). * ℎ Ω + ∫ Γ ( ). * ℎ Γ ∀ * ℎ ∈  0 ℎ ( 7
)
• the constitutive equation:

ℎ ( , ) = ( , ) ( ℎ ( , )) in Ω (8)
The classical weak form formulation is:

find ℎ ∈ { ∈  ℎ ; |Γ = } such that: ( ℎ ( ), * ℎ ; ) = ( * ℎ ; ) ∀ * ℎ ∈  0 ℎ ( 9
)
where

 0 ℎ = { ∈  ℎ ; | Γ = }.
Following [START_REF] Nguyen | Certified real-time solution of parametrized partial differential equations[END_REF] , we assume that the Hooke tensor , the body forces and the traction forces can be decomposed as the sum of functions of multiplied by deterministic functions.

( , ) = ∑ =1 Θ ( ) ̄ ( ), ( , ) = ∑ =1 Θ ( ) ̄ ( ) ( , ) = ∑ =1 Θ ( ) ̄ ( ) (10)
where Θ , Θ , Θ are known functions of , each ̄ ( ) is a fourth order tensor defined in Ω, ̄ ( ) are vector fields defined in Ω, and ̄ ( ) are vector fields defined on Γ .

Remark: An approximation of this decomposition can be obtained by the use of the Karhunen-Loeve expansion [START_REF] Ghanem | Stochastic Finite Elements: A spectral Approach[END_REF] as proposed in [START_REF] Chamoin | Robust goal-oriented error estimation based on the constitutive relation error for stochastic problems[END_REF] .

Furthermore, we assume that the inverse of the Hooke tensor can be decomposed as the sum of functions of multiplied by deterministic functions defined in Ω.

-1 ( ,

) = ∑ =1 Θ ( ) ̄ ( ) (11)
where each ̄ ( ) is a fourth order tensor defined on Ω. These hypotheses are required by the adaptive reduced basis method [START_REF] Gallimard | Towards error bounds of the failure probability of elastic structures using reduced basis models[END_REF] as shown in section 4.2.

Structural reliability model

In this section, we describe the problem of reliability assessment [START_REF] Ditlevsen | Structural reliability methods[END_REF][START_REF] Mohamed | Finite element and reliability: a method for compound variables: application on a cracked heating system[END_REF][START_REF] Pendola | Combination of finite element and reliability methods in nonlinear fracture mechanics[END_REF] . Let us denote by ( ) = ( ℎ ( , )) the computational output of interest and ( ) the threshold associated to this computational output. Classically the structure failure state is defined by a limit state function:

( ) = ( ) -( ) (12)
such that

• ( ) < 0 is a failure state for the structure,

• ( ) = 0 is the limit state,

• ( ) > 0 is a safe state for the structure.

According to this definition, the system fails when is lower or equal to zero. The failure probability is then given by:

= ∫  ={ ∈ℝ ; ( )≤0} ( ) (13)
where is the dimension of the vector and ( ) is a PDF. The evaluation of the integral defined by Equation ( 13) is not easy because it represents a very small quantity and because the integration domain is defined implicitly. Monte-Carlo algorithm is the main approach to solve the reliability problem. Recasting Equation ( 13) as:

= ∫ ℝ ( ( )≤0) ( ) ( ) = ( ( )≤0) (14)
where ( ( )≤0) is the failure indicator function being equal to one if ( ) ≤ 0 and zero otherwise. The probability of failure is equal to the expectation of ( ( )≤0) . The Monte-Carlo algorithm consists in generating a number of realizations of the random vector Θ using the actual probability density function, then in computing the Monte-Carlo estimator ̂ :

̂ = 1 ∑ =1 ( ( )≤0) ( ) (15)
According to the central limit theorem, this estimator is asymptotically unbiased and normally distributed with variance

̂ = ̂ (1 -̂ )
-1 When the failure probability is small the relative deviation of the estimator is

= √ ̂ ̂ ≈ 1 √ ̂ (16)
The relative deviation is consequently unbounded. From Equation ( 16), it can be seen that obtaining a relative deviation ≤ 10% for a probability of failure of 10 -requires about 10 +2 finite element simulations.

IMPORTANCE SAMPLING

Principles

The objective of importance sampling [START_REF] Melchers | Importance sampling in structural systems[END_REF][START_REF] Engelund | A benchmark study on importance sampling techniques in strucural reliability[END_REF][START_REF] Orsak | A class of optimum importance sampling strategies[END_REF] is to reduce the variance of the Monte-Carlo estimator ̂ . The idea is to generate the samples 1 , ..., with an auxiliary PDF ℎ such that the generated samples have a higher rate of falling in the failure region than with the PDF , because only these samples contribute to the evaluation of .

= ( ( )≤0) = ∫ ( ( )≤0) ( ) = ∫ ( ( )≤0) ( ) ℎ( ) ℎ( ) = ℎ ( ( )≤0) ( ) ℎ( ) (17)
Let us denote ( ) = ( ) ℎ( ) , then the failure probability is estimated by

̂ = 1 ∑ =1 ( ( )≤0) ( ) (18)
where ̂ is called the importance sampling (IS) estimator. The variance of this estimator is given by the following equation

( ̂ ) = ( ( )≤0) ( ) (19)
The optimal choice for ℎ is given by

ℎ * ( ) = ( ( )≤0) ( ) (20)
which leads to

( ( )≤0) ( ) = ( ( )≤0) ( ) ( ( )≤0) ( ) = (21)
and hence to ( ̂ ) = 0. The obvious difficulty of this approach is that ℎ * depends not only on the failure domain which is unknown, but also on the unknown probability . Nevertheless, a good sampling density ℎ will be close to the PDF ℎ * .

Cross-Entropy approach to determine a good Importance Sampling distribution

The idea of the Cross-Entropy (CE) [START_REF] Rubinstein | The cross-entropy method for combinatorial and continuous optimization[END_REF][START_REF] Boer | A tutorial on the cross-entropy method[END_REF][START_REF] Homem-De-Mello | A study on the cross-entropy method for rare-event probability estimation[END_REF] is to choose the IS density ℎ in a specified class of densities such that the crossentropy or Kullback-Leiber (K-L) divergence between the optimal importance sampling density ℎ * and ℎ is minimal. The K-L divergence between two PDF and is given by:

( , ) = ln ( ) ( ) = ∫ ( ) ln ( ) ( ) = ∫ ( ) ln ( ) -∫ ( ) ln ( ) (22)
Let us define ℎ a family of PDF indexed by the parameter ( could be the mean and the covariance matrix in case of Gaussian densities). The objective is then to find * which minimizes the Kullback-Leiber divergence (ℎ * , ℎ

) * = argmin (ℎ * , ℎ ) = argmin -∫ ℎ * ( ) ln ℎ ( ) (23)
substituting ℎ * ( ) from Equation (20) in Equation ( 23) we obtain

* = argmax ∫ ( ( )≤0) ( ) ln ℎ ( ) = argmax ( ( )≤0) ln ℎ ( ) (24)
Remark: In [START_REF] Homem-De-Mello | A study on the cross-entropy method for rare-event probability estimation[END_REF] , it is shown that solving the cross-entropy problem defined by Equation ( 23) yields an estimate with minimum variance. * can be estimated by solving the following stochastic program

* = argmax 1 ∑ =1 ( ( )≤0) ( ) ℎ ( ) ln ℎ ( ) (25)
where 1 , ..,. are random samples generated according to the probability density ℎ . In practice, ones does not solve directly Equation (25) since it requires the knowledge of several samples such that ( ) ≤ 0, which is not the case for rare event evaluation. To circumvent this difficulty, the optimization problem ( 25) is solved by a multilevel iterative method [START_REF] Rubinstein | The cross-entropy method for combinatorial and continuous optimization[END_REF][START_REF] Homem-De-Mello | A study on the cross-entropy method for rare-event probability estimation[END_REF] with a decreasing sequence of thresholds 0 > 1 > ... ≥ 0 chosen adaptively using a quantile definition. At each iteration, the value of -1 is available and ones determine in practice

= argmax 1 ∑ =1 ( ( )≤ ) ( ) ℎ -1 ( ) ln ℎ ( ) (26)
The algorithm is initialized by choosing a parameter (typically 10 -2 ≤ ≤ 10 -1 , see [START_REF] Rubinstein | The cross-entropy method for combinatorial and continuous optimization[END_REF][START_REF] Homem-De-Mello | A study on the cross-entropy method for rare-event probability estimation[END_REF] for more details), and by defining ℎ 0 = . At each iteration of the algorithm samples 1 , ..,. are generated according the probability density ℎ -1 , the values of the limit state function ( ) are computed and is set to the max(0, ) where is the value of the -quantile of ( 1 ), ... ( ) . Then, Equation ( 26) is solved and we set ℎ = ℎ . Let be the final iteration when = 0, samples 1 , ..,.

are generated according to the probability density ℎ and the importance sampling estimator is computed by

̂ = 1 ∑ =1 ( ( )≤0) ( ) ℎ ( ) (27)
The Cross-Entropy algorithm is described in Algorithm 1.

Remark: The maximization of Equation ( 26) can often be solved analytically, in particular when the probability density family is composed by independent Gaussian or log-normal random variables, as used in this paper (see section 6).

Simple 1D example

A simple 1D test problem proposed in [START_REF] Florentin | Adaptive reduced basis strategy based on goal oriented error assessment for stochastic problems[END_REF] is considered here to illustrate the CE algorithm. The unknown displacement ( , ) describes the deformation of a bar of section and length which is loaded by a traction on = and clamped at = 0 (see Fig 1 ). The random parameters are the Young modulus of the bar and the traction . They are functions of two independent Gaussian random variables ( 1 , 2 ) defined in Table 1 (

1 ) = 0 (1 + ( 1 )) with ( 1 ) = 2 arcsin(Erf( 1 √ 2
))

√ 2 -8 and ( 2 ) = 2 E( ) F( ) L FIGURE 1 Unidimensional bar Remark:
The nonlinear mapping is chosen such that the probability density function of ( 1 ) has a bounded support. Actually, the upper and lower limits for are

0 1 - √ 2 2 -8 ≤ ( 1 ) ≤ 0 1 + √ 2 2 -8
This definition allows precluding negative non-physical values of the Young modulus if

0 ≤ < √ 2 -8 2 .
Random variable distribution mean value standard deviation

1 Gaussian 1 1 2 Gaussian 2 2

TABLE 1 Statistical properties of the random variables

The limit state function is given by

( 1 , 2 ) = -( , 1 , 2 )
The analytical expression of the quantity of interest is easily available

( , 1 , 2 ) = ( 2 )
( 1 ) The PDF in this example is given as a function of the mean values and the standard deviations of 1 , 2 .

( 1 , 2 ; 1 , 1 , 2 , 2 ) = 2 ∏ =1 1 √ 2 exp - 1 2 - 2 (28)
The family of PDF to be optimized for the parameter = ( ̃ 1 , ̃ 1 , ̃ 2 , ̃ 2 ) will be defined by

ℎ( 1 , 2 ; ̃ 1 , ̃ 1 , ̃ 2 , ̃ 2 ) = 2 ∏ =1 1 √ 2 ̃ exp - 1 2 -̃ ̃ 2 (29)
As an example we consider the following values:

• Geometrical values: = 1, = 0.1.

• Young's modulus: 0 = 200, = 0.3, 1 = 0, 1 = 1.

• Applied traction: 2 = 1, 2 = 0.1.

• Critical displacement: = 0. Using the estimated optimal parameter vector 4 = (-3.10, 0.46, 0.73, 0.055), the final iteration with = 10 4 samples gives an estimate of ̂ = 1.51 × 10 -4 , with a relative deviation = 1.5%. In this case for an equivalent relative deviation we have reduced our simulation effort by a factor 10 3 . However, in practical mechanical computation it is mandatory to reduce the cost of generation of each sample. The aim of the section 4, is to introduce the so-called reduced basis method.

REDUCED BASIS METHOD

The purpose of a Reduced Basis Method is to provide a fast evaluation ( ) of the value of the displacement field ℎ ( ) and hence a fast evaluation ( ( ); ) of any quantity of interest ( ℎ ( ); ). In this section we recall a classical displacement approach which consists in performing a Galerkin projection onto a reduced basis space that is assumed to represent accurately the solutions of the problem to be solved. The construction of this reduced basis space, that will be detailed in section 5, is aimed to control the accuracy of the computed failure probability. To control this accuracy we define and compute an error estimator which is an upper bound of the error on the quantity of interest | ( ℎ ( ); ) -( ( ); )|. The computation of this upper bound necessitates the construction of a stress field that satisfies the finite element equilibrium equations (Equations 7). Following the approach developed in [START_REF] Gallimard | Towards error bounds of the failure probability of elastic structures using reduced basis models[END_REF] , in order to obtain a fast evaluation of this stress field, we propose in section 4.1 to build, from the displacement reduced basis, a stress reduced basis which satisfies the finite element equilibrium equations. In this paper, we use a reduced basis strategy which permits to control the quality of the reliability analysis which was proposed in [START_REF] Gallimard | Towards error bounds of the failure probability of elastic structures using reduced basis models[END_REF] . We present the main ingredients of this strategy in the following subsections.

Offline Phase

Let ∈  ℎ be a displacement field such that |Γ = . Let us introduce a set of samples in the parameter space  = { 1 , … , }, where ∈ , and for each compute a finite element solution 0 ℎ ( ) in  0 ℎ described by the corresponding vector of nodal values .

( 0 ℎ ( ), * ℎ ; ) = ( * ℎ ; ) -( , * ℎ ; ) ∀ * ℎ ∈  0 ℎ ( 30
) On this space, following [START_REF] Huynh | Reduced basis approximation and a posteriori error estimation for stress intensity factors[END_REF][START_REF] Rozza | Reduced basis approximation and a posteriori error estimation for affinely parametrized elliptic coercive partial differential equations[END_REF][START_REF] Martini | Reduced basis approximation and a-posteriori error estimation for the coupled stokes-darcy system[END_REF] , we perform a Gram-Schmidt orthonormalization process for the inner product ( , ; ̄ ) ( ̄ being a fixed value of the random vector which is selected a priori as shown in Section 5.1). The reduced basis space is then defined by

 0, = span { 1 , … , } ⊂  0 ℎ ( 31
)
The choice of the samples in the parameter space  and of the associated reduced basis  0, depends on the sampling strategy (see [START_REF] Rozza | Reduced basis approximation and a posteriori error estimation for affinely parametrized elliptic coercive partial differential equations[END_REF] for more details).

To construct the stress reduced basis, we follow the method presented in [START_REF] Gallimard | Towards error bounds of the failure probability of elastic structures using reduced basis models[END_REF] . The first step consists in building a stress ( ) which verifies the F.E. equilibrium (Equation 7) for all in  (We refer the reader to [START_REF] Gallimard | Towards error bounds of the failure probability of elastic structures using reduced basis models[END_REF] for more details about the construction of ( )). Let us consider the set of stress fields computed from the snapshot solutions (Equation 30)

= ( ) ( 0 ℎ ( ) + ) for ∈ {1, … , } (32)
and the set of stress fields defined by

Δ = - ( ) (33)
It follows that {Δ , for ∈ {1, … , }} is a set of stress fields equilibrated to zero in the FE sense. An orthonormal basis  0, = { 1 , … , } is built from Δ by a Gram-Schmidt process with respect to an internal product ∫ Ω 1 ∶ -1 ( ̄ ) 2 Ω.

Online Phase

The reduced basis approximation consists in solving Equation ( 5) in  0 + { }. A key point to justify the use of the reduced basis approximation is that is assumed to be much smaller than the number of degree of freedom of the F.E. model (i.e.

<<

). The reduced basis solution for the displacement field writes

( ) = + ∑ =1 (34)
The coefficients are computed by solving

( 0 ( ), * ; ) = ( * ; ) -( , * ; ) ∀ * ∈  0 (35)
which leads to an algebraic system

[ ( )] [ ] = [ ( )] (36
Thanks to the decomposition (Equation 10) of , and , and can be written as a linear combination of the functions Θ ( ), Θ ( ) and Θ ( )

( ) = ∑ =1 ̄ Θ ( ) and ( ) = ∑ =1 ̄ Θ ( ) + ∑ =1 ̄ Θ ( ) (38)
where

̄ = ∫ Ω ̄ ( ) ∶ ( ) Ω, ̄ = ∫ Ω ̄ . Ω, and ̄ = ∫ Γ ̄ . Γ
The reduced basis solution for the stress field writes

( ) = ( ) + ∑ =1 (39)
The coefficients are computed in order to minimize a distance between the stress field ( ) computed from the reduced basis solution ( ) through the constitutive equation ( ( ) = ( ) ( ( ))) and ( )

( ) = argmin ( 1 ,..., ) ( ) (40)
where

( ) = ∫ Ω ( ( ) - ( )) ∶ -1 ( )( ( ) - ( )) Ω (41)
The minimization of Equation ( 40) leads to the algebraic system

[ ( )] [ ] = [ ( )] (42)
The elements of [ ( )] and of [ ( )] are defined by

= ∫ Ω -1 ( ) ∶ Ω (43) and = ∫ Ω ∶ ( ( ) --1 ( ) ( )) Ω (44)
Thanks to the decompositions (Equation 10) and (Equation 11) and can be written as a linear combination of the functions Θ ( ), Θ ( ) and Θ ( )

= ∑ =1 ̄ Θ ( ) = ̄ 0 + ∑ =1 ̄ - ∑ =1 ∑ =1 ̄ Θ ( )Θ ( ) + ∑ =1 ̄ Θ ( )Θ ( ) where ̄ = ∫ Ω ∶ Ω, ̄ 0 = ∫ Ω ∶ ( ) Ω, ̄ = ∫ Ω ∶ ( ) Ω and ̄ = ∫ Ω ∶ ̄ ( ̄ ) ( ) Ω, ̄ = ∫ Γ ∶ ̄ ( ̄ ) ( ) Γ

Error estimator and bounds on the quantity of interest

The computation of error bounds on linear Quantities of Interest (QoI) has been developed by many authors [START_REF] Yang | A posteriori pointwise upper bound error estimates in the finite element method[END_REF][START_REF] Prudhomme | On goal-oriented error estimation for elliptic problems : application to the control of pointwise errors[END_REF][START_REF] Ohnimus | Local error estimates of fem for displacements and stresses in linear elasticity by solving local neumann problems[END_REF][START_REF] Prudhomme | Practical methods for a posteriori error estimation in engineering applications[END_REF][START_REF] Ladevèze | New bounding techniques for goal-oriented error estimation applied to linear problems[END_REF] for measuring the gap between an exact QoI and a QoI computed from a finite element analysis. The error bounds between a QoI computed from a finite element analysis and a QoI computed by a reduced modeling approach have been studied since early 2000 [START_REF] Machiels | Output bounds for reduced-order approximations of elliptic partial differential equations[END_REF][START_REF] Boyaval | Reduced basis techniques for stochastic problems[END_REF][START_REF] Haasdonk | Reduced basis methods for parametrized partial differential equations with stochastic influences using the Karhunen Loeve expansion[END_REF][START_REF] Ryckelynck | Estimation of the validity domain of hyper-reduction approximations in generalized standard elastoviscoplasticity[END_REF][START_REF] Hoang | A fast, certified and tuning free two-field reduced basis method for the metamodelling of affinely-parametrised elasticity problems[END_REF] . In this paper, we use the bounds proposed in [START_REF] Gallimard | A posteriori global error estimator based on the error in the constitutive relation for reduced basis approximation of parametrized linear elastic problems[END_REF][START_REF] Gallimard | Towards error bounds of the failure probability of elastic structures using reduced basis models[END_REF] within the framework of the application of the error in the constitutive relation to reduced basis computations. These bounds are obtained by using the parallelogram inequality as proposed in [START_REF] Prudhomme | On goal-oriented error estimation for elliptic problems : application to the control of pointwise errors[END_REF][START_REF] Machiels | Output bounds for reduced-order approximations of elliptic partial differential equations[END_REF] .

For any realization of the random vector the error introduced by the reduced approximation is given by

( ) = ℎ ( ) -( ) (45)
due to the linearity assumption, one has

( ( ); ) = ( ℎ ( ); ) -( ( ); ) (46)
Following [START_REF] Ladevèze | New bounding techniques for goal-oriented error estimation applied to linear problems[END_REF] , we consider the following auxiliary problem: find ℎ ∈  0 ℎ such that:

( * ℎ , ℎ ( ); ) = ( * ℎ ( ); ) ∀ * ℎ ∈  0 ℎ ( 47
)
and its solution in the reduced basis ( ( ), * ; ) = ( * ; ) ∀ * ∈  0 (48)

It can be shown that (we refer the reader to [START_REF] Prudhomme | Practical methods for a posteriori error estimation in engineering applications[END_REF][START_REF] Gallimard | A constitutive relation error estimator based on traction-free recovery of the equilibrated stress[END_REF][START_REF] Gallimard | A posteriori global error estimator based on the error in the constitutive relation for reduced basis approximation of parametrized linear elastic problems[END_REF][START_REF] Gallimard | Towards error bounds of the failure probability of elastic structures using reduced basis models[END_REF] for more details)

--( ) ≤ ( ( ); ) ≤ + ( ) (49)
where ( ) and ( ) are the stresses computed from the reduced basis solutions ( ) and ( ) through the constitutive equation, ( ) and , ( ) are stress fields equilibrated in the FE sense which are computed, for a given value of , as proposed in section (4.2). ‖ ‖ , is defined by

+ ( ) = 1 2 ‖ ( ) - ( )‖ , ‖ , ( ) - ( )‖ , + < ( ) - ( ), , ( ) - ( ) > , -( ) = 1 2 ‖ ( ) - ( )‖ , ‖ , ( ) - ( )‖ , -< ( ) - ( ), , ( ) - ( ) > ,
‖ ‖ 2 , =< , > , with < 1 , 2 > , = ∫ Ω 1 ∶ -1 ( ) 2 Ω
From Equations ( 46) and ( 49) it follows that the quantity of interest ( ℎ ( ); ) is bounded by

( ( ); ) --( ) ≤ ( ℎ ( ); ) ≤ ( ( ); ) + + ( ) (50)
Remark: If does not depend linearly on ℎ , bounding of the quantity interest is much more difficult and there are few works addressing the subject (see [START_REF] Huynh | Reduced basis approximation and a posteriori error estimation for stress intensity factors[END_REF][START_REF] Haasdonk | Reduced basis methods for parametrized partial differential equations with stochastic influences using the Karhunen Loeve expansion[END_REF] for quadratic outputs or [START_REF] Gallimard | Evaluation of the local quality of the von mises's stress and l2-norm of the stress[END_REF] for mean von Mises stress over a subdomain).

Bounds of the limit state and on the failure probability

The approximate limit state function computed from the reduced basis analysis is

( ) = ( ) -( ( ); ) (51)
The approximate failure probability computed by a Monte-Carlo method is given by

̂ = 1 ∑ =1 ( ( )≤0) ( ) (52)
and from Equation ( 50) the following computable bounds of the limit state function ( ) are obtained

-( ) ≤ ( ) ≤ + ( ) (53)
where

-( ) = ( ) -+ ( ) and + ( ) = ( ) + -( ) (54)
Let us introduce the failure domain associated with the lower bound and the upper bound of the limit-state function and the associated failure probabilities

 + = { ; ( ) -+ ( ) ≤ 0} and  -= { ; ( ) + -( ) ≤ 0} ̂ + = 1 ∑ =1  + ( ) and ̂ -= 1 ∑ =1  -( ) (55)
As  -⊂  = { ; ( ) ≤ 0} ⊂  + we obtain the following bounds for the failure probability

̂ -≤ ̂ ≤ ̂ + (56)
Remark: It must be noticed that the obtained bounds concern the failure probability ̂ computed from a Monte-Carlo algorithm.

The quality of the Monte-Carlo estimator will be classically assessed in the examples presented in section 6 via the central limit theorem and the computation of the relative deviation (Equation 16). Furthermore, to control the reduced basis algorithm an error estimator ( ), which is an upper bound on the error on the limit state function

( ) = | ( ) - ( )|, is introduced ( ) ≤ ( ) = max( -( ), + ( )) (57)
Properties:

1. If for all , ( ) = 0 then ̂ -= ̂ + = ̂ .

2. If ̂ -= ̂ + then ( ) = 0 for all and ̂ = ̂ -= ̂ + .

REDUCED BASIS ALGORITHM FOR RARE EVENTS SIMULATIONS

The Monte-Carlo method (see Section 2.2 and Equation 15) is the reference method to compute the failure probability. To evaluate a failure probability = 10 -with a relative deviation = 10% the Monte-Carlo method necessitates 10 +2 FOM computations. To reduce the computational effort it is interesting to introduce reduced order models based on the RB method [START_REF] Chen | Accurate and efficient evaluation of failure probability for partial different equations with random input data[END_REF][START_REF] Gallimard | Towards error bounds of the failure probability of elastic structures using reduced basis models[END_REF] . As shown in Section 4, this method involves the construction of a basis (the so-called reduced basis). The construction of this reduced basis can be performed offline in a training space as proposed in [START_REF] Chen | Accurate and efficient evaluation of failure probability for partial different equations with random input data[END_REF] or adaptively during the Monte-Carlo simulation as proposed in [START_REF] Gallimard | Towards error bounds of the failure probability of elastic structures using reduced basis models[END_REF] . In any case 10 +2 ROM computations are performed in order to evaluate the failure probability . For rare events simulations, it is mandatory to reduce not only the cost of the computations but also their number. The method developed here proposes to reduce the number of samples needed by using an IS estimator along with the reduced basis approximation. In Section 3, we presented the CE method for computing an IS optimized distribution. The idea of the method proposed in this paper is to build the reduced basis while the iterations of the CE method are performed. As, at each iteration the optimization process involves only the samples such that the limit-state function ( ) is below a given threshold (Equation 26), the reduced basis approximation is improved only when these samples are not computed with a sufficient accuracy in the actual reduced basis.

Compared to a classical MC-IS algorithm, the advantage of the proposed algorithm is that the computations of the IS failure probability is done with a ROM rather than a FOM. Compared to a classical MC-RB algorithm, the advantage of the proposed algorithm is that the number of computations of Monte-Carlo samples performed in the ROM is greatly reduced.

1. Offline stage : Construct the reduced bases while performing the iterations of the CE algorithm.

2. Online stage : Compute the importance sampling estimator, together with its inferior bound and superior bound, with the reduced bases defined in step 1.

Offline stage

The initial step of the algorithm consists in computing a displacement field satisfying the boundary conditions (Equation 6) and a stress field ( ) satisfying the equilibrium in the FE sense (Equation 7). The computation of and ( ) involves the choice of a particular value of denoted ̄ , a simple choice is to use the mean value of the random vector ̄ = [ ]. This value will also be used in the Gram-Schmidt orthonormalization algorithm. The idea is to construct the reduced bases along with the CE optimization process. For each sample generated during the CE optimization, an approximate solution ( ) of ℎ ( ) is computed in the reduced basis space  0, . The approximate limit state function ( ) is computed as well as its upper bound + ( ) and lower bound -( ). Two cases must be distinguished:

• If -( ) > -1 then ( ) > -1
and this sample will not contribute to the optimization problem (Eq 58).

• If -( ) ≤ -1 this sample may or may not contribute to the optimization problem (Eq 58).

In the first case the quality of value of ( ) provided by the reduced basis  0, is sufficient and the simulations are continued on the same reduced basis. In the second case, the reduced bases  0, and  0, are improved only if the error estimator on the limit state function is greater than a fixed tolerance ( ) ≥

; in this case the problem is solved in the finite element space and new reduced basis spaces  0, +1 and  0, +1 are constructed. When the values of ( 1 ), ..., ( ) have been computed, we compute their -quantile and we optimize the parameters in order to minimize the Kullback-Leibler divergence between the auxiliary PDF ℎ and the optimal auxiliary PDF for the limit-state function defined by the threshold : ( ) =

. The algorithm is given in Algorithm (2).

It must be noticed that the reduced bases are constructed iteratively for decreasing thresholds 1 > 2 > ... > > 0, whereas our actual limit-state function is defined by ( ) = 0. To obtain bases which accurately describe the limit-state function, a new iteration of the optimization process is added. In this iteration, we add vectors in the reduced bases when the RB accuracy is insufficient to certify the state of the structure, and we optimize the parameters in order to minimize the Kullback-Leibler divergence between the auxiliary PDF ℎ and the optimal auxiliary PDF for the limit-state function ( ) = 0. The algorithm is given in Algorithm 4. In this stage, the reduced bases are increased when the error estimator on the limit state function is greater than a fixed tolerance with ≤ , in order to describe accurately the limit-state defined by ( ) = 0.

Online stage and error assessement

During the online stage, the optimized values * of the parameters are used to generate the random variables such that a sufficient number of samples are in the failure domain. And the values of the limit-state function and its bounds are computed in the reduced bases computed along with the CE optimization which have been constructed to accurately represent the limit state function. The algorithm is given in Algorithm 5.

Then, the failure probability and the bounds are estimated by . At the end of the first step of optimization process (Algorithm 2), we obtain reduced bases of dimension 3, and the evolution of the PDF parameters presented in Table 3 . During the last step of optimization process (Algorithm 4) the dimension of the reduced bases increases to 5. The optimized PDF parameters are given in Table 4 . Using the estimated optimal parameter vector 4 = (115. 6, 11.5, 201.4, 42.6, 145.6, 22.6) and the size 5 reduced bases with = 5 × 10 3 samples gives an estimate of ̂ , = 2.48 × 10 -3 , with a relative deviation = 5.1%. We obtain as well a lower bound of the failure probability ̂ -, = 2.30 × 10 -3 and an upper bound ̂ +, = 2.52 × 10 -3 . We observe that, while 5 × 10 4 FOM resolutions are needed a standard MC method (SMC1) to obtain a relative deviation of = 8.6%, the proposed algorithm necessitates only 5 FOM resolutions to compute the RB and the IS optimized density during the offline stage and 5 × 10 3 ROM resolutions to compute the failure probability during the online stage. Compared to a standard MC method, the number of finite element simulations has been reduced by a factor 10 4 and the number of Monte-Carlo throws (during the online stage) has been reduced by a factor 10.

̂ , = 1 ∑ =1 ( ( )≤0) ( ) ℎ ( ) , ̂ +, = 1 ∑ =1 ( -( )≤0) ( ) ℎ ( ) and ̂ -, = 1 ∑ =1 ( + ( )≤0) ( ) ℎ ( ) (60)
̃ 1 ̃ 1 ̃ 2 ̃ 2 ̃ 3 ̃ 3 1 1.

Example 2

The second example concerns a composite structure composed of two layers linked by an interface as shown on Figure 3 . The structure is submitted, in plane strain, to a uniform normal traction applied along the lower and upper surfaces. The deterministic parameters are the value of the traction = 9.2 , the crack length = 0.15, the length of the structure = 0.6, the thickness of the layers = 0.14, the thickness of the interface = 0.02, the radius of the fibers = 0.05 and their positions. The random parameters, as defined in Table 5 , are the Young's moduli of the materials = ( 1 , 2 , 3 , 4 , , ) (where denotes the fibers, the matrix and the interface). The Poisson ratio's are set to 0.3. Figure 4 shows the 6-node triangular mesh used for the finite element analysis.

The Quantity of Interest is the mean value of the over the subdomain and the limit state function is given by

( ) = - 1 ( ) ∫ ( ) (61)
where is a deterministic value set to 180 . A standard Monte-Carlo method (SMC1) with 10 7 samples gives an estimate ̂ = 3.16 × 10 -4 with a relative deviation = 1.8%. For constructing the reduced bases and the CE optimization we use the following parameters: = 0.1, = 500, = 18 and = 1.8 . The optimization process converges in 4 iterations (Algorithms 2 and 4) and we obtain reduced bases of dimension 9. The optimized PDF parameters are given in Table 6 . Using optimized PDF parameters and the size 9 reduced bases with = 10 4 samples gives an estimate of ̂ , = 3.19 × 10 -4 , with a relative deviation = 2.9%. We obtain as well a lower bound of the failure probability ̂ -, = 2.81 × 10 -4 and an upper bound ̂ +, = 3.72 × 10 -4 . To compare with a Monte-Carlo simulation performed by a RB without CE optimization, we run the Monte-Carlo simulation with the reduced basis approximation computed during the offline step and the initial PDF , and we compute the failure probability ̂ with the Monte-Carlo estimator (Equation 15). The number of samples used is = 4 × 10 6 , and we obtain ̂ = 3.12 × 10 -4 , with a relative deviation = 2.8%. In this case, we have reduced the number of computations of the ROM during the online phase by a factor 4 × 10 for the two set of computations is shown in Figure 5 . It can be observed that the ratio is very sharp varying from 0.9 to 1.1 and that the estimate becomes more reliable when the tolerance decreases. For = 1.8 the size of the reduced bases varies from 9 to 10 for the different runs and for = 1.8 the size of the reduced bases varies from 15 to 17. Figure 6 displays the variation the failure's probability bounds ( ̂ +, and ̂ -, ) for 10 runs and shows that the bounds gaps become narrower as the tolerance decreases. , leading to = 28. It can be observed that in this case 10 3 samples are sufficient to obtain a relative deviation inferior to 10%.

Table 7 shows the evolution of the failure probability and the evolution of the number of the elements in the reduced bases when loading decreases. The algorithm is run with the following parameters: = 0.1, = 500, = 10 4 , = 1.8 , = 0.18 . The failure probability decreases rapidly with the force drop. However, it can be seen that there is

FIGURE 7

Evolution of the relative deviation as a function of no increase of the number of elements in the reduced bases and that less than 20 FE computations are sufficient to accurately estimate the failure probability. The relative deviation remains stable between 3% and 5%. If we consider that 10 9 finite element computations would be necessary to estimate the 7.4 × 10 -8 failure probability we observe that the number of finite element computations is reduced by a factor 5 × 10 7 and the number of computations in the reduced basis by a factor 10

CONCLUSIONS

In this paper, we have developed an adaptive computational strategy aiming to compute the failure probability of a structure in the context of rare events simulations. We use a cross-entropy method to construct an optimized probability density function for the importance sampling estimator. A surrogate model based on a reduced basis approach is constructed during the cross-entropy optimization and the reduced basis vectors are chosen in order to control the quality of the estimation of the limit state function.

In the numerical examples, we show that the algorithm permits estimating failure probabilities from 10 -4 to 10 -7 with less than 20 finite elements computations during the offline phase and only 10 4 calls to the surrogate model during the online phase. Finally, bounds of the Monte-Carlo estimator are computed which permit to assess the quality of the reduced basis. It is interesting to note that a possible improvement of the algorithm would be to use the bounds of the Monte-Carlo estimator to add vectors in the reduced basis during the online process in order to prescribe certified bounds for the estimate of the failure probability.

Algorithm 1

 1 Pseudocode for the IS based on a CE optimization input: the PDF , (number of samples used for the CE optimization), (number of samples used for the IS estimation), a parameter in [10 -2 , 10 -1 ] output: ̂ Define ℎ 0 = . Set 0 such that 0 > 0. Set = 0 (iteration counter) while > 0 do Set = + 1 Generate random samples 1 , .., according to the probability density ℎ -1 . Compute ( ) (for = 1,).Compute= max(0, ), where is the -quantile of ( 1 ), ..., ( ) Use 1 , .., and( 1 ), ..., ( ) to solve * = argmax 1 ∑ =1 ( ( )≤ ) ( ) ℎ -1 ( ) ln ℎ ( ) Set = * and ℎ = ℎ . end while Set ℎ = ℎ Generate samples according to the probability density ℎ and compute ̂

  )The elements of [ ( )] and of [ ( )] are defined by ( ) = ( , ; ) and ( ) = ( ; )

FIGURE 2

 2 FIGURE 2 Example 1: Thick plate studied

Fibers( 1

 1

FIGURE 3 FIGURE 4

 34 FIGURE 3 Example 2: Composite structure studied

FIGURE 5 Figure 7

 57 FIGURE 5 Evolution of =

FIGURE 6

 6 FIGURE 6Variation of the failure's probability bounds ( ̂ +, and ̂ -, ) for 10 runs

TABLE 2

 2 Convergence of the CE optimization optimization we use the following parameters: = 0.1, = 10 3 and = 10 4 . Table2displays the results of the CE optimization.

			̃ 1	̃ 1	̃ 2	̃ 2
	1	0.058	0	1	0.5	0.1
	2	0.022	-1.64 0.66 0.56 0.086
	3 0.00014 -2.59 0.37 0.64 0.060
	4	0	-3.10 0.46 0.73 0.055

1

A standard Monte-Carlo method with 10 6 samples gives an estimate ̂ = 1.48 × 10 -4 , with a relative deviation = 8.2%. With 10 7 samples we obtain an estimate ̂ = 1.54 × 10 -4 , with a relative deviation = 2.5%. For the IS method based on CE

TABLE 3

 3 Convergence of the CE optimization

			̃ 1	̃ 1	̃ 2	̃ 2	̃ 3	̃ 3
	4	0	115.6 11.5 201.4 42.6 145.6 22.6

TABLE 4

 4 Convergence of the CE optimization last iteration

TABLE 5

 5 Example 2: Statistical properties of the random variables

  [START_REF] Lemaire | Structural reliability[END_REF] .

	Random Variable	1	2	3	4
	Mean Value (GPa)	95.9 97.8 113.0 15.76 8.10 1.84
	Standard deviation (GPa) 19.6 19.4 27.2	2.88 2.19 0.25

TABLE 6

 6 Example 2: optimized PDF parameters Two series of tests are performed to check the robustness of the failure probability estimate. The algorithm is run 10 times with the following parameters = 0.1,

		= 500,		= 10 4 ,	= 18	,	= 1.8	, then with the same
							̂ ,
	parameters but with a smaller tolerance	= 1.8	,	= 0.18	. The evolution of the ratio =	̂

TABLE 7

 7 [START_REF] Der Kiureghian | Efficient algorithm for second-order reliability analysis[END_REF] . × 10 -4 8.4 × 10 -5 1.89 × 10 -5 3.4 × 10 -6 5.6 × 10 -7 7.4 × 10 -8 ̂ +, 3.7 × 10 -4 8.6 × 10 -5 1.94 × 10 -5 3.6 × 10 -6 6.3 × 10 -7 7.8 × 10 -8 ̂ -, 2.8 × 10 -4 8.1 × 10 -5 1.86 × 10 -5 3.3 × 10 -6 5.1 × 10 -7 7.1 × 10 -8 Computed failure probability and bounds

	(MPa)	9.2	9.0	8.8	8.6	8.4	8.2
	̂ , 3.2 Relative deviation (%)	2.9	2.9	4.0	3.6	4.7	3.0
	CE number of iterations	4	4	5	5	5	6
		15	16	17	18	16	19

† Adaptive reduced basis strategy for rare events simulations.

Algorithm 2 Pseudocode for the CE optimization algorithm input: the PDF , (number of samples used for the CE optimization), (Acceptable error on the estimation of the limit state function : ( ) ≤ ), a parameter in [10 -2 , 10 -1 ] output:  0, ,  0, and ℎ

Define ℎ 0 = . Set 0 such that 0 > 0. Set = 0 (iteration counter). Generate a random sample according to the probability density ℎ 0 . Compute ℎ ( ) on the finite element mesh. Compute 1 = ℎ ( ) -and 1 = ℎ ( ) -( ). Initialize the reduced bases  0,1 = span { 1 } and  0,1 = span { 1 }. = max(0, ), where is the -quantile of ( 1 ), ..., Compute -( ) and + ( ) from ( 49) Compute ( ) from (57) Compute ( ), -( ) and + ( ) from ( 51) and ( 54)

6 NUMERICAL EXAMPLES

Example 1

In this first example, we aim at validating the algorithm by a comparison with a standard Monte-Carlo algorithm. We consider a plate with two rectangular holes proposed in [START_REF] Peraire | Bounds for linear-functional outputs of coercive partial differential equations: local indicators and adaptive refinement[END_REF] . The structure is submitted, in plane strain, to a normal traction = 10 applied along the vertical edge and the plate is considered to be composed of three different materials. The symmetry of the problem allows to study only one fourth of the plate as shown in Figure 2 . The Poisson ratios are fixed 1 = 2 = 3 = 0. ( ( ), ) = 1 ( ) ∫ ( ).

The limit state function is given by ( ) = ̄ -( ( ), )