Audrey Lejart

Siham Zentout

Catherine Chapuis

Ostiane D'augustin

Rebecca Smith

Gilles Salbert
email: gilles.salbert@univ-rennes1.fr

Sébastien Huet
email: sebastien.huet@univ-rennes1.fr

The N-terminal domain of TET1 promotes the formation of dense chromatin regions refractory to transcription

Keywords: Ten eleven translocation enzyme, DNA methylation, DNA hydroxymethylation, chromatin, transcription, fluorescence microscopy

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

In the nucleus of eukaryotic cells, DNA is wrapped around histones to form the chromatin fibre which, itself, displays a hierarchy of folding steps. The complex organization of chromatin is essential to ensure that the entire genetic material of the cell can be contained within the nucleus, where it occupies 20 to 50% of the nuclear volume [START_REF] Ou | ChromEMT: Visualizing 3D chromatin structure and compaction in interphase and mitotic cells[END_REF]. Despite this high level of packaging, the chromatin structure is not fixed and constantly evolves in relation to multiple cellular processes which require access to DNA such as transcription, replication or repair [START_REF] Fierz | Biophysics of Chromatin Dynamics[END_REF]. Key actors in the regulation of chromatin folding are epigenetic marks such as post-transcriptional modifications found on histones and base modifications found in DNA. These dynamic marks attract or repel effectors of DNA transactions and control chromatin organisation either directly, by modulating for example nucleosome stability, or indirectly, by recruiting remodelers or chromatin scaffolding proteins [START_REF] Talbert | Histone variants on the move: substrates for chromatin dynamics[END_REF].

While epigenetic marks found on histones are extremely diverse, only a handful of modifications are found at the DNA level, the main one being 5-methylcytosine (5-mC) (Raush et al., 2019). This mark, mainly found at CpG dinucleotides, has a repressive impact on transcription, in particular through the recruitment of repressor complexes and by regulating the binding of transcription factors to DNA at regulatory regions [START_REF] Watt | Cytosine methylation prevents binding to DNA of a HeLa cell transcription factor required for optimal expression of the adenovirus major late promoter[END_REF][START_REF] Iguchi-Ariga | CpG methylation of the cAMP-responsive enhancer/promoter sequence TGACGTCA abolishes specific factor binding as well as transcriptional activation[END_REF]. 5-mC is written by several DNA methyl-transferase enzymes and its removal arises either by passive dilution due to DNA replication or by active DNA demethylation [START_REF] Bhutani | DNA demethylation dynamics[END_REF][START_REF] Mulholland | Recent evolution of a TET-controlled and DPPA3/STELLA-driven pathway of passive DNA demethylation in mammals[END_REF]. The main active demethylation pathway is based on the successive oxidation of 5-mC into 5-hmC (5-hydroxymethylcytosine), 5-fC (5-formylcytosine) and finally 5-caC (5-carboxylcytosine), which are carried out by members of the TET (ten-eleven translocation) enzyme family [START_REF] Ito | Tet proteins can convert 5-methylcytosine to 5-formylcytosine and 5-carboxylcytosine[END_REF]. While the 5-fC and 5-caC forms are rapidly excised by the DNA repair machinery to restore an unmethylated cytosine [START_REF] Maiti | Thymine DNA glycosylase can rapidly excise 5-formylcytosine and 5carboxylcytosine: potential implications for active demethylation of CpG sites[END_REF], 5-hmC has a relatively long lifetime and could thus constitute a fully-fledged epigenetic mark [START_REF] Münzel | 5-Hydroxymethylcytosine, the sixth base of the genome[END_REF].

The TET family is composed of three members showing distinct expression patterns hinting at different cellular functions [START_REF] Lorsbach | TET1, a member of a novel protein family, is fused to MLL in acute myeloid leukemia containing the t(10;11)(q22;q23)[END_REF]. All members share a catalytic domain responsible for 5-mC oxidation that is localized at the C-terminus of the proteins [START_REF] Tahiliani | Conversion of 5-methylcytosine to 5hydroxymethylcytosine in mammalian DNA by MLL partner TET1[END_REF][START_REF] Hu | Crystal structure of TET2-DNA complex: insight into TETmediated 5mC oxidation[END_REF]. Furthermore, TET1 and TET3, but not TET2, bear a CXXC domain that binds to unmethylated CpGs and is thought to be responsible for the localization of these proteins at CpG-rich genomic sites like CpG islands, protecting them from aberrant methylation [START_REF] Tahiliani | Conversion of 5-methylcytosine to 5hydroxymethylcytosine in mammalian DNA by MLL partner TET1[END_REF]. By controlling the balance of 5-mC/5-hmC levels within the genome, the TET enzymes are key regulators of transcriptional processes associated with cell differentiation, tumorigenesis and neurodegeneration [START_REF] Rasmussen | Role of TET enzymes in DNA methylation, development, and cancer[END_REF][START_REF] Zhang | Tet3 Reads 5-Carboxylcytosine through Its CXXC Domain and Is a Potential Guardian against Neurodegeneration[END_REF]. However, the role of TET ACCEPTED MANUSCRIPT -CLEAN COPY enzymes is not limited to their catalytic activity. Indeed, they have been shown to serve as a binding platform for effectors that regulate transcription such as the histone methyltransferase complex PRC2 [START_REF] Wu | Dual functions of Tet1 in transcriptional regulation in mouse embryonic stem cells[END_REF], the histone acetyltransferase hMOF [START_REF] Zhong | TET1 modulates H4K16 acetylation by controlling autoacetylation of hMOF to affect gene regulation and DNA repair function[END_REF] or the O-linked Nacetylglucosamine transferase OGT [START_REF] Shi | Ten-eleven translocation 1 (Tet1) is regulated by O-linked Nacetylglucosamine transferase (Ogt) for target gene repression in mouse embryonic stem cells[END_REF][START_REF] Bauer | Phosphorylation of TET Proteins Is Regulated via O-GlcNAcylation by the O-Linked N-Acetylglucosamine Transferase (OGT)[END_REF]. Importantly, the existence of a shorter isoform of TET1 (TET1short) has been recently reported. While the long isoform of TET1 (TET1long) is expressed mainly in undifferentiated cells, TET1short levels seem to increase upon differentiation as well as during tumorigenesis [START_REF] Zhang | Isoform Switch of TET1 Regulates DNA Demethylation and Mouse Development[END_REF][START_REF] Good | A novel isoform of TET1 that lacks a CXXC domain is overexpressed in cancer[END_REF]. Despite lacking the N-terminal domain of TET1long, including the CXXC domain, TET1short is still able to demethylate CpG islands [START_REF] Zhang | Isoform Switch of TET1 Regulates DNA Demethylation and Mouse Development[END_REF]. Nevertheless, TET1short cannot fully complement the longer isoform, suggesting that the N-terminal domain of the full-length protein fulfils specific cellular functions that remain to be clarified.

In the present work, we demonstrate that the N-terminal region of the TET1 protein, which is lacking in TET1short, is able to induce a global reorganisation of chromatin in three different human cell lines (U2OS, HeLa and MCF-7). This restructuring process is seen at the larger chromatin folding scales but also at the molecular level since histone exchange rates are impacted by the TET1 N-terminal domain. We also show that the chromatin reorganisation is associated with a global shut-down of transcription. Altogether, our findings uncover a previously unidentified function of TET1 in the regulation of the chromatin structure with potential consequences regarding its roles in the regulation of transcription.

RESULTS

Overexpression of TET1 leads to a global reorganisation of the chromatin structure

By using confocal fluorescence imaging, we observed that overexpression of GFP-TET1 and GFP-TET2, but not GFP-TET3, led to a global reorganisation of the chromatin within the nucleus of human osteosarcoma U2OS cells as seen both by Hoechst and anti-DNA staining (Fig. 1a andS1a).

While the chromatin in control cells occupied the entire nuclear space relatively evenly, the expression of GFP-TET1 and GFP-TET2 led to chromatin condensation into large bundles, leaving some nuclear areas largely devoid of staining. Since this phenotype was more pronounced for TET1, we focused on the impact of TET1 expression on chromatin structure. To better understand how chromatin reorganisation arises upon TET1 overexpression, we transiently transfected GFPtagged TET1 or free GFP in U2OS cells stably expressing H2B-mCherry and monitored the progressive increase in the GFP fluorescence levels in parallel with changes in the fluorescence contrast of the chromatin pattern, providing a readout of the chromatin compaction state in living cells. While no significant change in the chromatin pattern was observed in cells expressing GFP (Fig. S1b,c), a rapid increase of the fluorescence contrast of H2B-mCherry was observed upon ACCEPTED MANUSCRIPT -CLEAN COPY expression of GFP-TET1 (Fig. 1b,c). The tight correlation between the timing of the chromatin reorganisation and the rising of TET1 levels within the nucleus suggests that the change in the chromatin pattern is directly linked to the presence of TET1 rather than a secondary consequence of cellular mechanisms controlled by this enzyme.

TET1 is responsible for 5-mC oxidation to allow its later replacement by unmethylated cytosine [START_REF] Tahiliani | Conversion of 5-methylcytosine to 5hydroxymethylcytosine in mammalian DNA by MLL partner TET1[END_REF]. However, when expressing a catalytically inactive form of TET1 we observed the same chromatin pattern as with the active enzyme (Fig. 1d). Furthermore, this reorganisation phenotype also did not rely on the CXXC domain of TET1 (Fig. 1d). Therefore, these data show that the change in chromatin structure promoted by TET1 is not the consequence of 5-mC oxidation processes or further DNA demethylation steps.

The N-terminal domain of TET1 is sufficient to induce chromatin reorganisation

To better delineate the domains of TET1 responsible for chromatin reorganisation, we analyzed the impact of the expression of two truncated forms of the full-length protein: a C-terminal fragment We also assessed the affinity of these two truncated TET1 constructs for chromatin by monitoring the kinetics of fluorescence recovery after photobleaching (FRAP) for each GFP tagged construct.

The NTER fragment displayed a much slower recovery than the CTER construct (Fig 2d-f), in agreement with previous findings showing that it is principally the NTER region that controls TET1 binding to chromatin [START_REF] Zhang | Isoform Switch of TET1 Regulates DNA Demethylation and Mouse Development[END_REF]. Nevertheless, the NTER fragment appeared less bound to chromatin than the full length TET1, which was also confirmed by a reduced colocalization of the former with Hoechst (Fig S2d). Altogether, these results demonstrate that the dynamic association of the N-terminal domain of TET1 with chromatin, which displays characteristic exchange times of ACCEPTED MANUSCRIPT -CLEAN COPY about 10 seconds, alters higher-order chromatin states, as shown by light microscopy, ultimately leading to a spatial redistribution of the genomic material within the nucleus.

Large-scale chromatin reorganisation promoted by the N-terminal domain of TET1 is associated with changes in histone association with chromatin

We next assessed whether the changes in chromatin folding induced by the N-terminal domain of TET1 were also associated with altered dynamics of its molecular constituents. Histones are key chromatin scaffolding proteins regulating the spatial organization of the genome within the nucleus, and their exchange rates correlate with the chromatin compaction state [START_REF] Misteli | Dynamic binding of histone H1 to chromatin in living cells[END_REF][START_REF] Bancaud | Molecular crowding affects diffusion and binding of nuclear proteins in heterochromatin and reveals the fractal organization of chromatin[END_REF]. Therefore, we assessed the impact of the N-terminal domain of TET1 on the exchange rates of two different histones: the core histone H2B as well as the linker histone H1.

Using FRAP, we observed that expression of the N-terminal domain of TET1 decreased H2B exchange rate (Fig 3a,b), in line with the fact that this domain seems to promote the local condensation of chromatin. Conversely, the linker histone H1 displayed reduced association with chromatin in cells expressing the N-terminal domain of TET1 compared to controls (Fig. 3c,d).

Therefore, the dense chromatin domains generated by the expression of the N-terminal domain of TET1 differ from transcriptionally silent heterochromatic areas in which H1 displays increased binding [START_REF] Misteli | Dynamic binding of histone H1 to chromatin in living cells[END_REF][START_REF] Bancaud | Molecular crowding affects diffusion and binding of nuclear proteins in heterochromatin and reveals the fractal organization of chromatin[END_REF].

Expression of the N-terminal domain of TET1 leads to a global shutdown of transcription

The results obtained so far demonstrate that the N-terminal domain of TET1 regulates the architecture of chromatin at multiple scales. Since this enzyme was shown to play key roles in the regulation of transcription, we next tested whether the changes in chromatin organisation induced by the N-terminal domain of TET1 were also associated with a modulation of the transcriptional state. By quantifying 5-ethynyl uridine (EU) staining both by confocal imaging and flow cytometry to assess the level of nascent RNA in cells expressing the N-terminal domain of TET1, we found that the presence of this domain of TET1 led to a partial inhibition of transcription, similarly to the overexpression of the full-length protein (Fig. 4a-c andS3b). Furthermore, the negative colocalization between the EU and Hoechst staining in cells expressing the N-terminal domain of TET1 suggests that transcription is confined at the interface between the dense chromatin domains and the nucleoplasm (Fig. S3a). Concomitant to the reduction in RNA synthesis rates, immunofluorescence staining also showed that the expression of the N-terminal domain of TET1 increased the levels of trimethylation of the lysine 27 of histone H3 (H3K27me3) (Fig 4d-e), an epigenetic mark written by the PRC2 complex and found within facultative heterochromatin [START_REF] Cao | Role of histone H3 lysine 27 methylation in Polycomb-group silencing[END_REF]. These results were confirmed by flow cytometry, which also revealed a positive ACCEPTED MANUSCRIPT -CLEAN COPY correlation between the amount the N-terminal domain of TET1 expressed in the cells and the levels of H3K27me3 mark (Fig. 4f andS3b). Altogether, these findings demonstrate that the Nterminal domain of TET1 is able to shut-down transcription by establishing a chromatin environment displaying some of the characteristics found in heterochromatin regions.

DISCUSSION

Our data show that overexpression of TET1 promotes chromatin remodeling in three different human cell lines (U2OS, MCF-7 and HeLa), which leads to the establishment of a condensed chromatin architecture refractory to transcription. These dense chromatin domains observed in the presence of TET1 share some features found in heterochromatin, such as increased levels of H3K27me3 [START_REF] Wiles | H3K27 methylation: a promiscuous repressive chromatin mark[END_REF], but nevertheless differ from these compact regions. Indeed, while the histone linker H1 was shown to be crucial for the condensation of the chromatin fiber [START_REF] Willcockson | H1 histones control the epigenetic landscape by local chromatin compaction[END_REF], we observed reduced association of this histone in the compact chromatin conformation promoted by TET1. Interestingly, senescence associated heterochromatin foci, which arise from heterochromatic areas detaching from the nuclear lamina, also show reduced H1 binding due to its replacement by the high mobility group protein A2 (HMGA2) [START_REF] Funayama | Loss of linker histone H1 in cellular senescence[END_REF]. It would be useful to assess whether TET1 overexpression could promote cell senescence.

The TET1-mediated chromatin condensation process reported here is independent of both the catalytic activity of this enzyme and its targeting to CpG islands via the CXXC domain. The Nterminal region of the protein, which triggers this condensation, displays a strong affinity for chromatin. Therefore, it could serve as a platform to recruit factors modulating chromatin conformation, as it was recently proposed for OGT [START_REF] Zhang | Isoform Switch of TET1 Regulates DNA Demethylation and Mouse Development[END_REF]. TET1 has been shown to interact with the histone-deacetylases Sin3A and HDAC1 as well as the chromatin remodeler CHD4, three actors known to promote chromatin compaction [START_REF] Shi | Ten-eleven translocation 1 (Tet1) is regulated by O-linked Nacetylglucosamine transferase (Ogt) for target gene repression in mouse embryonic stem cells[END_REF]. Future experiments would be needed to assess whether the recruitment of these proteins to chromatin through their binding to TET1 could contribute to the establishment of the compact genome organisation observed in cells expressing the N-terminal domain of TET1. Furthermore, the increase in H3K27me3 staining in these dense domains may be promoted by the association between TET1 and the PRC2 complex [START_REF] Neri | TET1 is controlled by pluripotency-associated factors in ESCs and downmodulated by PRC2 in differentiated cells and tissues[END_REF], contributing to a global silencing of transcription.

In this work, we uncovered a new role of the N-terminal region of TET1 in the regulation of the chromatin organisation and transcription inhibition. The inhibitory impact of this region of the protein might be needed to balance the positive influence of the catalytic domain on transcription, where unrestrained activity of the recently identified short isoform of TET1 is associated with ACCEPTED MANUSCRIPT -CLEAN COPY tumorigenesis [START_REF] Good | A novel isoform of TET1 that lacks a CXXC domain is overexpressed in cancer[END_REF]. Further work will be needed to better delineate the characteristics of this coordination mechanism between the two domains of TET1.

METHODS

Plasmids

Plasmids coding for the full length TET1, TET2 and TET3 tagged with GFP on their N-terminal end are gifts from Heinrich Leonhardt and Cristina Cardoso [START_REF] Frauer | Different binding properties and function of CXXC zinc finger domains in Dnmt1 and Tet1[END_REF]. The plasmids coding for the C-terminal and N-terminal fragments of TET1, both tagged with GFP on their Cterminal end, as well as the TET1 mutant lacking the CXXC domain, tagged with GFP on its Nterminal end, were provided by Wei Xie [START_REF] Zhang | Isoform Switch of TET1 Regulates DNA Demethylation and Mouse Development[END_REF]. The catalytically-inactive mutant of full length TET1 fused to GFP on its N-terminal end was obtain by site-directed mutagenesis on GFP tagged TET1 with primers 5'-gttgtgcatgttgtgaatggccttgtaagaatgggcacaaaaatcc-3' and 5'ggatttttgtgcccattcttacaaggccattcacaacatgcacaac-3' using Phusion polymerase (Thermo Scientific).

Haico Van Attikum [START_REF] Rother | CHD7 and 53BP1 regulate distinct pathways for the re-ligation of DNA double-strand breaks[END_REF] provided the plasmid coding for GFP fused to a NLS (nuclear localisation sequence), which promotes GFP accumulation in the nucleus.The H2B-mCherry was a gift from Jan Ellenberg (Euroscarf P30632, [START_REF] Neumann | Phenotypic profiling of the human genome by time-lapse microscopy reveals cell division genes[END_REF]). mEGFP was cut from H1-mEGFP (gift from Gyula Timinszky) and replaced with mCherry using the restriction sites BshTI and Bsp1407I (ThermoScientific), to generate the H1-mCherry construct.

Cell culture

The experiments were performed on three different human cell lines: osteosarcoma U2OS cells, cervical carcinoma HeLa cells and breast cancer MCF-7 cells. All cell lines were cultured in DMEM (4.5 g/L glucose, Sigma) supplemented with 10% fetal bovine serum (Life Technologies), 2 mM glutamine (Sigma), 100 µg/mL penicillin and 100 U/mL streptomycin (Sigma) in a humidified incubator at 37 °C with 5% CO 2 . For transient expression of fluorescent constructs, cells were transfected 5 hours after seeding into glass-bottom chambers (Thermo Scientific and Zell Contact) with XtremeGENE HP (Roche) or Xfect (Takara Bio) according to the manufacturer's instructions.

U2OS stably expressing H2B-mCherry were generated by transfection of the H2B-mCherry plasmid followed by a two-week selection in growing medium supplemented with puromycin (0.25μg/ml, Invivogen). Immediately prior to imaging live cells, growth medium was replaced with CO 2 -independent imaging medium (Phenol Red-free Leibovitz's L-15 medium (Life Technologies) supplemented with 20% fetal bovine serum, 2 mM glutamine, 100 µg/mL penicillin and 100 U/mL streptomycin). All experiments were completed with unsynchronized cells.

Immunofluorescence staining

ACCEPTED MANUSCRIPT -CLEAN COPY

Cells were fixed using 4% paraformaldehyde and permeabilized with 0.1% Triton X-100 at room temperature. For 5-hmC staining, chromatin was denaturated by adding hydrochloric acid at 2 M for 30min. Cells were blocked with blocking buffer (BSA 3% Tween 0.05% in PBS) before incubation overnight at 4°C with primary antibodies recognizing double-strand DNA (Abcam, ab27156), 5-hmC (Active Motif,ref 39769) or the trimethylation of histone H3 at lysine 27 (Upstate, ref 07-449), diluted in the blocking buffer. The secondary antibodies labeled with Alexa Fluor 555 (ThermoScientific, ref A21428) diluted in blocking buffer was added for one hour at room temperature. DNA was stained for 15 min with Hoechst 33342 at 1µg/mL (ThermoScientific).

Nascent RNA staining by 5-ethynyl uridine

For nascent RNA staining, cells were incubated with 5-ethynyl uridine (EU) at 1 mM in culture medium for 1 h before being fixed with 4% paraformaldehyde. Cells were permeabilised with 0.1% Triton X-100 for 15 min before RNA was labeled using the Click-iT® RNA Alexa Fluor 594 Imaging Kit according to the manufacturer's instructions (ThermoScientific).

Confocal imaging and fluorescence recovery after photobleaching

Fluorescence imaging of live and fixed cells was on a LSM880 point-scanner confocal microscope (Zeiss) using a 60x/1.4 N.A. oil-immersion objective lens. The pinhole was set to one Airy unit and fluorescence detection was performed on a GaAsP detector array. The fluorescence of the Hoechst was excited at 405 nm and detected within a window ranging from 415 to 470 nm. The fluorescence of the GFP was excited at 488 nm and detected within a window ranging from 500 to 550 nm. The fluorescence of the mCherry and Alexa Fluor 555 was excited at 561 nm and detected within a window ranging from 580 to 650 nm. The fluorescence of the Alexa Fluor 594 was excited at 594 nm and detected within a window ranging from 605 to 695 nm. The pixel size was set to 70 nm.

Laser intensities and detector gains were chosen to avoid detector saturation and the acquisition settings were kept constant within the course of a given experiment. FRAP of H2B-mCherry was performed on the same setup. Photobleaching was achieved by 30 iterations of scanning half the nucleus at 100% power of the 561 nm laser. Then cells were imaged for 11 h at 1 image per hour to monitor the progressive recovery within the bleached area. FRAP of H1-mCherry as well as the full-length or truncated forms of TET1 tagged with GFP was completed on a Ti-E inverted microscope from Nikon equipped with a CSU-X1 spinning-disk head from Yokogawa, a Plan APO 60x/1.4 N.A. oil-immersion objective lens and a sCMOS ORCA Flash 4.0 camera. A 4 µm diameter circle within the nucleus was bleached by 40 consecutive scans at 100% of 561 or 488 nm and subsequent fluorescence recovery was monitored by imaging the cells for 25 to 50 seconds at 4 frames/second. For all live-cell experiments, cells were maintained at 37°C with a heating chamber.

ACCEPTED MANUSCRIPT -CLEAN COPY

Image analysis

Image analysis was performed using the ImageJ software. For the quantification of fluorescence signals, the mean intensity within manually defined region delineating the nucleus was measured and subtracted from background. Colocalization between two channels was performed by calculating the Pearson coefficient using the JACoP plugin [START_REF] Bolte | A guided tour into subcellular colocalization analysis in light microscopy[END_REF]. To assess chromatin compaction, we estimated the fluorescence contrast on images of nuclei expressing H2B-mCherry or stained with Hoechst, using the plugin GLCM_Texture written by Julio E. Cabrera. The contrast is one of the Haralick features [START_REF] Haralick | Textural Features for Image Classification[END_REF] and measures the mean squared intensity difference between pixels separated by a given distance d. This distance d was set to 10 pixels, which corresponded to the characteristic size of the dense chromatin domains generated by TET1 overexpression. Since we found that the measured contrast was independent of the cardinal directions, this parameter was only calculated along the east direction. The contrast estimated for each nucleus was the average of the values measured for three regions of 20*20 pixels randomly chosen within this nucleus, taking care that these regions do not overlap with the nucleus border or the nucleoli. For the timelapse acquisitions shown on figures 1b and S1b the changes in fluorescence contrast were measured from the images of the H2B-mCherry patterns and the levels of expression of unfused GFP or GFP-TET1 was assessed by the mean intensity in the GFP channel. The expression level was then rescaled to range between 0 and 1. Before averaging, the individual curves for the chromatin contrast and the expression level were aligned in time using as a reference the timepoint at which the normalized expression level equalled 0.5. The time was then estimated relative to this reference timepoint.

To analyze the FRAP acquisitions, regions of interest were defined manually to measure the fluorescence intensities within the bleached area (Ibleach), the whole nucleus (Inuc) and outside the cell (Ibg), along the entire image sequence. Based on these intensity measurements, the intensity within the bleached area was corrected from background and from the loss of total nuclear intensity due to both local photobleaching and timelapse imaging [START_REF] Bancaud | Fluorescence perturbation techniques to study mobility and molecular dynamics of proteins in live cells: FRAP, photoactivation, photoconversion, and FLIP[END_REF], using the following expression:

I frap = I bleach -I bg I nuc -I bg
After subtracting the signal estimated immediately after the bleach, the intensity was normalized to the pre-bleach signal. To estimate the characteristic recovery time, the post-bleach curves were fitted with the following single exponential recovery model:

I rec (t )=1-exp(-t /τ ) ,
with t the time after bleaching and τ the fitted characteristic recovery time.

(

  CTER) bearing the catalytic domain shared between the different TET proteins, and a N-terminal fragment (NTER) which remains poorly characterized except for the CXXC domain (Fig 2a). Importantly, the CTER is reminiscent of the recently identified short isoform of TET1 found expressed in differentiated cells. While the CTER is still able to oxidize 5-mC into 5-hmC (Fig S2a,b), in line with what is reported for TET1short (Good et al., 2017), its expression in U2OS cells did not lead to visible changes in the chromatin organisation (Fig 2b). In contrast, expression of the NTER fragment, while not affecting 5-hmC levels within the nucleus (Fig S2a,b), led to a global reorganisation of the chromatin resembling the one observed when overexpressing the full-length protein (Fig 2b). The amount of chromatin reorganisation, assessed by the fluorescence contrast of the Hoechst staining, correlated with the levels of expression of the NTER fragment of TET1, similar to what was observed for the full-length protein (Fig S2c). Importantly, this chromatin restructuring phenotype upon expression of the NTER fragment was also observed in human HeLa and MCF-7 cells, demonstrating that it is not specific to U2OS cells (Fig 2c).

FIGURE LEGENDS Figure 1 .

 LEGENDS1 FIGURE LEGENDS Figure 1. The overexpression of GFP-TET1 leads to chromatin reorganisation. a. Images of Hoechst-stained UO2S cells expressing either unfused GFP, or GFP-tagged TET1, TET2 or TET3. Cells were imaged 24h after transfection. More than 30 cells were imaged per condition and a typical image is shown for each of them. b. Timelapse image sequence of live U2OS cells stably expressing H2B-mCherry and transiently expressing GFP-TET1. The images were pseudocolored according to the look-up table displayed below the images. The timelapse acquisition started 6 h after transfection c. Images from b. were quantified to monitor the changes in GFP-TET1 expression levels (black curve), assessed by GFP fluorescence intensity, in parallel to the chromatin spatial distribution, assessed by the fluorescence contrast estimated from the images of the H2B-mCherry signal (red curve). The curves show a mean ± SD of 17 cells. Before averaging, the individual curves for the chromatin contrast and the expression level were aligned in time using as a reference the timepoint at which the normalized expression level equalled 0.5. The time was then estimated relative to this reference timepoint. d. Images of Hoechst-stained U2OS cells expressing either unfused GFP, or full-lenght TET1 (GFP-FL), a catalytically-inactive mutant of TET1 (GFP-FLmut) or a TET1 construct lacking the CXXC domain (GFP-ΔCXXC), all tagged with GFP. Cells were imaged 24h after transfection. For panels a., b. and d., scale bars = 10 µm.

Figure 2 .

 2 Figure 2. The N-terminal domain of TET1 promotes chromatin reorganisation. a. Schematic representation of the GFP-tagged TET1 protein domains. Cys-rich corresponds to a cysteine-rich region and DSBH stands for double-stranded β-helix. b. Images of Hoechst-labeled U2OS cells expressing either unfused GFP, full-length (GFP-FL), N-terminal domain (NTER-GFP) or Cterminal domain (CTER-GFP) of TET1, all tagged with GFP. More than 100 cells were imaged per condition and a typical image is shown for each. c. Images of Hoechst-labeled HeLa and MCF-7 cells expressing unfused GFP or the N-terminal of TET1 tagged with GFP. More than 30 cells were imaged per condition and a typical image is shown for each. d. Image sequences of the fluorescence redistribution after photobleaching a circular area of the nucleus (dotted red circles) of U2OS cells

Acknowledgements

ACCEPTED MANUSCRIPT -CLEAN COPY

We thank the engineers of the Microscopy Rennes Imaging Center (BIOSIT, Université Rennes 1), Stéphanie Dutertre and Xavier Pinson for technical help on the microscopes. We are also grateful to the engineers of the flow cytometry facility (BIOSIT, Université Rennes 1), Laurent Deleurme and Alexis Aimé for technical assistance. This work would not have been possible without the generous donation of plasmids from Heinrich Leonhardt, Cristina Cardoso, Wei Xie, Haico Van Attikum, Gyula Timinszky and Jan Ellenberg.

Funding

This work was funded by the Ligue contre le Cancer du Grand-Ouest, the Conseil Régional de Bretagne and the Institut Universitaire de France. The Microscopy-Rennes Imaging Center (BIOSIT, Université Rennes 1), a member of the national infrastructure France-BioImaging supported by the French National Research Agency (ANR-10-INBS-04).

Availability of data and materials

The raw datasets used during the current study are available from the corresponding authors on reasonable request.

Code availability

Not applicable

Flow cytometry

For staining with the Click-iT® RNA reagent (ThermoScientific), U2OS cells were seeded into 6well plates and transfected with the corresponding GFP tagged TET1 expression plasmid. After 16h, cells were treated for 1h at 37°C with 1mM of 5-ethynyl uridine diluted in DMEM medium.

Then, cells were harvested by trypsinization and fixed with 4% paraformaldehyde for 20min at room temperature. After permeabilization for 15 min using the Click-iT® saponin-based permeabilization reagent, the cells were stained with the click-chemistry cocktail according to the manufacturer's instructions for 30 min at room temperature and protected from light. For H3K27me3 staining, U2OS cells were seeded into 6-well plates and transfected with the corresponding GFP tagged TET1 expression plasmid. After 16 h, cells were fixed with 4% paraformaldehyde, permeabilized with saponin-based permeabilization reagent and incubated with rabbit polyclonal antibody against H3K27me3 (Sigma-Aldrich, 07-449) for 1 h at 1/100 dilution in the permeabilization reagent. After three cycles of washes/centrifugation with the permeabilization buffer, the pellet was resuspended and stained for 30 min using a 1/1000 dilution of Alexa Fluor 647 secondary antibody (ThermoScientific, A315573). Finally, the samples were washed/centrifuged three times with PBS and the pellets were resuspended with PBS for cytometry analysis.

All the conditions were analyzed on a LSRFortessa X-20 (BD Bioscience) using BD FACSDiva Software v8.0.1. GFP was excited at 488nm laser and detected within 30 nm window centred on 530 nm. For the Click-iT® staining using Alexa Fluor 594, fluorescence was excited with a 575nm laser and detected within a 30 nm window centred on 620 nm. Alexa Fluor-647 was excited at 640 nm laser and detected within a 30 nm window centred on 670 nm. At least 10 5 events were collected for each condition. Quantifications were performed with the FACSDiva software (BD Biosciences).

Statistics

All the experiments were performed in at least three replicates with reproducible results for each replicate. When a single experiment is shown, a typical replicate was chosen. For the curves displaying time courses, means±SD are shown. Boxplots were generated using a web-tool developed by the Tyers and Rappsilber labs (http://boxplot.tyerslab.com/). The box limits correspond to the 25th and 75th percentiles and the bold line indicates the median value. The whiskers extend 1.5 times the interquartile range and outliers are shown by dots. The numbers in parentheses below each condition refer to the number of cells, for the imaging data, or to the number of experiments, for the cytometry data. p values were calculated using either two-sided unpaired Student's t-test assuming unequal variances for data showing normal distributions, or ACCEPTED MANUSCRIPT -CLEAN COPY unpaired two-samples Wilcoxon test. On the boxplots, * refers to p<0.05, ** to p<0.01, *** to p<0.001, **** to p<0.0001 and n.s. to non significant.

DECLARATIONS

Ethics approval and consent to participate

Not applicable

Consent for publication

Not applicable

Competing interests

The authors declare no competing interests.

Authors' contributions


Théo Gnassounou

Pierre Humbert

Laurent Oudre

ADAPTIVE SUBSAMPLING OF MULTIDOMAIN SIGNALS WITH PRODUCT GRAPHS

Keywords: Graph signal processing, graph product, subsampling, graph learning, Laplacian matrix estimation

In this paper, we propose an adaptive subsampling method for multidomain signals based on the constrained learning of a product graph. Given an input multidomain signal, we search for a product graph on which the signal is bandlimited, i.e. have limited spectral occupancy. The subsampling procedure described in this article is composed of two successive steps. First, we use the input data to learn a graph that will be optimized to favor efficient sampling. Then, we derive an algorithm for choosing the best nodes and provide a sampling strategy for multidomain signals. Experiments on synthetic data and two real datasets show the efficiency of the proposed method and its relevance for multidomain data compression and storing.

INTRODUCTION

In standard signal processing, the sampling theory states that a bandlimited signal sampled above its Nyquist rate can be perfectly reconstructed. This important property is at the cornerstone of the sampling of Euclidean signals. However, when signals are defined over a more complex domain, the design of adaptive sampling strategies is still an active topic of interest. In order to deal with signals residing on irregular domains, Graph Signal Processing (GSP) [START_REF] David I Shuman | The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains[END_REF][START_REF] Ortega | Graph signal processing: Overview, challenges, and applications[END_REF] has emerged as a powerful alternative to standard approaches. In this formalism, the graph defines a support, and the signals, now called graph signals, are defined on this support. This allows to capture the structure on which a signal evolves, hence providing more information than considering the signal alone. By generalizing concepts and tools of signal processing to signals recorded over graphs, GSP has proven its success in many tasks such as filtering [START_REF] Sandryhaila | Discrete signal processing on graphs[END_REF], reconstruction [START_REF] Chen | Signal recovery on graphs: Variation minimization[END_REF], and sampling [START_REF] Antonio G Marques | Sampling of graph signals with successive local aggregations[END_REF].

For the later, one idea proposed in the univariate case is to reconstruct a graph signal from its measurement at some nodes by taking advantage of its underlying graph. This approach known as graph sampling set selection (or subset sampling) is now well investigated [START_REF] Chen | Discrete signal processing on graphs: Sampling theory[END_REF][START_REF] Chen | Signal recovery on graphs: Random versus experimentally designed sampling[END_REF][START_REF] Tsitsvero | Signals on graphs: Uncertainty principle and sampling[END_REF]. For instance, (in the noiseless setting) assuming that the graph signal is bandlimited one can show that a random selection of a reasonable number of samples/nodes is sufficient to enjoy a perfect reconstruction with a high probability [START_REF] Shomorony | Sampling large data on graphs[END_REF]. Unfortunately, there exist some major limitations of such methods.

First, to date, most articles have focused on univariate signals. However, recent publications in GSP have argued the need for multidomain graph signal processing that allows to deal with tensor data instead or vector data [START_REF] Sandryhaila | Big data analysis with signal processing on graphs[END_REF][START_REF] Varma | Random sampling for bandlimited signals on product graphs[END_REF]. Indeed, in several application contexts such as sensor networks, data streams are recorded as multivariate time series that evolve on a network, thus defining at least three domains: space, time and dimension. Second, these methods are mostly based on the computation of the first eigenvectors of the Laplacian. This step is computationally prohibitive for large graphs. Third, as in many computational tasks such as spectral clustering, and semi-supervised learning, the availability of the underlying graph is a core assumption. However, in most situations, no natural graph can be derived or defined and this graph must be inferred from the available data.

In this paper, we propose to address these issues in order to provide an efficient adaptive subsampling method for multidomain signals. Our approach extensively relies on tensor algebra to take into account the multidomain property of signals. To circumvent the complexity problem, we exploit the product structure of large graphs by factorizing them as the product graph of smaller graphs. This is a standard technique used to lower the complexity of GSP algorithms [START_REF] Sai | Learning product graphs from multidomain signals[END_REF][START_REF] Asad | Learning product graphs underlying smooth graph signals[END_REF]. The adaptive property of our subsampling procedure is based on a graph learning step that will search for a graph on which the signals can be efficiently sampled. As this graph learning problem is inherently an ill-posed problem, most state-of-the-art approaches use additional constraints such as smoothness (minimal signal variability between adjacent nodes) in order to retrieve an adequate graph [START_REF] Dong | Learning graphs from data: A signal representation perspective[END_REF][START_REF] Kalofolias | How to learn a graph from smooth signals[END_REF]. In the context of adaptive sampling, which is the core of this article, a less-studied constraint naturally emerges: the bandlimitness [START_REF] Sardellitti | Graph topology inference based on sparsifying transform learning[END_REF][START_REF] Humbert | Learning laplacian matrix from graph signals with sparse spectral representation[END_REF][START_REF] Batiste | Learning laplacian matrix from bandlimited graph signals[END_REF]. Indeed, Graph Fourier Transform (GFT) theory insures that if a graph signal if K-bandlimited (i.e. has only K non null spectral coefficients), then a perfect recovery is possible from its values on only K well-chosen nodes [START_REF] Chen | Discrete signal processing on graphs: Sampling theory[END_REF]. In order to take advantage of this property, we will seek for a graph where the signals have a sparse spectral representation, which will in turn ensure that the signals will be easy to sample efficiently.

BACKGROUND

We now introduce the notations of this paper on tensor algebra and graph product (see [START_REF] Tamara | Tensor decompositions and applications[END_REF][START_REF] Hammack | Handbook of product graphs[END_REF] for more complete introductions).

Tensor algebra. For d1, d2, . . . , dP ∈ N+, let denote a tensor of order P by Y ∈ R d 1 ×d 2 ×...×d P . This tensor can be seen as a discretized multidomain signal with each of its entries indexed over P domains. The matrix Y (m) in R dm×(d 1 ...d m-1 •d m+1 ...d P ) represents the tensor Y unfolded along the dimension m . The mode-m matrix product between a matrix X ∈ R j×dm and Y is defined as Y ×m X ⇐⇒ XY (m) . The operator ⊗ represent the Kronecker product. When multiple products are necessary, we use upper version of the notations.

Graph notations. Throughout the paper, let consider a weighted and undirected graph G = (V, E), where V is a set of N nodes and E is a set of edges. The combinatorial Laplacian matrix of this graph is a N by N matrix defined as L = D -W , where D is the degree matrix and W the weight matrix. Since G is an undirected graph, L is a symmetric and positive semi-definite matrix verifying L = XΛX T with Λ the diagonal matrix of non-negative eigenvalues of L and X the matrix of the corresponding eigenvectors as columns. Assuming that G has one connected component, L has λ1 = 0 for first eigenvalue associated with the eigenvectors X:,1 = 1N / √ N with 1N the unitary vector of size N . The space of such Laplacian matrices is given by:

L = {L ∈ S N + | L1N = 0N , Lij = Lji ≥ 0, ∀i = j} , (1)
where S N + is the set of real, symmetric and positive semi-definite matrices of size N × N .

Graph Signal Processing. A graph signal is a function y : V → R N assigning a scalar value to each node of a graph G. This function can be represented as a vector y ∈ R N , where yi is the function value at the i-th node. The eigenvectors of the Laplacian of G provide a Fourier-like basis for graph signals, allowing to decompose any signal into its spectral components. From this formalism, the Graph Fourier Transform (GFT) of y is defined by h = X T y. A K-bandlimited graph signal is thus a signal for which hi = 0 in K entries [START_REF] David I Shuman | The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains[END_REF][START_REF] Ortega | Graph signal processing: Overview, challenges, and applications[END_REF].

Cartesian graph product. Consider two weighted and undirected graphs G1 = (V1, E1) and G2 = (V2, E2) with N1 and N2 nodes. According to the foregoing there are L1 = X2Λ1X T 1 and

L2 = X2Λ2X T 2 .
A graph G is said to be a Cartesian graph when it is the result of the Cartesian product between G1 and G2. The Laplacian matrix of G is then given by:

L = (X1 ⊗ X2)(Λ1 ⊕ Λ2)(X T 1 ⊗ X T 2 ) , (2)
where ⊕ denotes the Kronecker sum i.e. for

A ∈ R N 1 ×N 1 and B ∈ R N 2 ×N 2 , A ⊕ B = I2 ⊗ A + B ⊗ I1 with
In the identity matrix of size Nn.

The notion of GFT and K-bandlimitedness can be extended in this multidomain context. The multidomain GFT of a graph signal Y ∈ R N 1 ×N 2 , is given by:

H = X1YX T 2 ⇐⇒ H = Y ×1 X T 1 ×2 X T 2 . (3)
Furthermore, given a vector K = (K1, K2), a K-bandlimited graph signal is a signal for which H has K1 rows and K2 columns different than zero. In other word, a signal is K-bandlimited when it is simultaneously Ki-bandlimited respectively in each domain.

It is possible to generalize equations ( 2) and (3) to P graphs (Gp) P p=1 with (Np) nodes. Denoting by Lp = XpΛpX T p their respective Laplacian matrices, the Laplacian matrix of their Cartesian products is [START_REF] Hammack | Handbook of product graphs[END_REF]:

L = ( P p=1 Xp) P p=1 Λp( P p=1 X T p ) . (4)
The multidomain GFT for the graph signal

Y ∈ R N 1 ×N 2 ×...×N P is then given by H = Y P × p=1 X T p .
Recall that by definition of the mode product, we also have

Y = H P × p=1
Xp .

SUBSAMPLING OF MULTIDOMAIN SIGNALS

Let consider a multidomain signal Y ∈ R N 1 ×...×N P and a vector K = (K1, . . . , KP ). We want to efficiently subsample Y by using only P p=1 Kp samples. This problem can be decomposed into two main tasks. First, learn an adaptive multidomain graph representation where the signal is K-bandlimited. Second, select the P p=1 Ki nodes (i.e. samples) that allow to achieve the best reconstruction errors.

Learning the graph Fourier basis

In this step, our goal is to learn the P graph Fourier basis {Xp} P p=1 of each subgraph composing the Cartesian graph together with the sparse spectral representation H of the signal on this graph. The associated optimization problem writes as:

minimize {Xp}p,H Y -H P × p=1 Xp 2 F s.t. X T p Xp = I, xp 1 = 1N p / Np, ∀p H (p) 2,0 ≤ Kp, ∀p , (5)
where • F is the Frobenius norm and H (p) 2,0 simply counts the number of nonzero rows in the matrix H (p) . In this problem, the first constraints insures that the matrices {Xp}p are proper graph Fourier basis and the second that the input signal is K-bandlimited in the spectral domain induced by the learned basis. Indeed, the norm imposes zeros on the rows of the unfolding H which make it adapted for the bandlimitedness assumption. To solve [START_REF] Antonio G Marques | Sampling of graph signals with successive local aggregations[END_REF] we use an alternating procedure on X1, . . . , XP and H.

Xp computation. To learn Xp with fixed ({Xp} p = p , H), we reformulate the objective function of ( 5) into:

Y -H P × p=1 Xp 2 F = Y ×p X T p -H P × p=1,p = p X T p 2 F = X T pY ( p) -H(p) 2 F , (6)
with ) . This manipulation leads to the following optimization problem: minimize

H(p) (H P × p=1,p = p X T p ) ( p
X p∈R N p ×N p X T pY ( p) -H(p) 2 F s.t. X T pX p = Ip, xp 1 = 1N p / Np , (7)
where a solution can be derived from [START_REF] Humbert | Learning laplacian matrix from graph signals with sparse spectral representation[END_REF]: let X be any matrix that verifies the constraint of [START_REF] Chen | Signal recovery on graphs: Random versus experimentally designed sampling[END_REF]. Introducing the matrix M = ( XT Y ( p) H(p) 1:,1: ) the submatrix containing everything except the first column and the first row. Let PDQ T be the singular value decomposition of M. The problem accepts the following closed form:

Xp = X 1 0 T N p -1 0N p -1 PQ T . ( 8
)
This can be computed in O(N 3 p ) leading to a complexity of O( p N 3 p ) to compute all {Xp}p. Note than the use of graph products is particularly interesting in this context, since the standard complexity (when the full graph is not supposed to be a Cartesian graph) is O(( P p=1 Np) 3 ). H computation. To compute H with fixed {Xp} p , we use the strategy proposed in [START_REF] Humbert | Subsampling of multivariate time-vertex graph signals[END_REF]. By fixing all dimensions except the p-th, the solution of the subproblem is obtained by sorting the rows of the matrix Y P × p=1 X T p ( p) by their 2-norm and then by selecting the Kp rows with larger norms (row/column-wise hard thresholding) [START_REF] Richard G Baraniuk | Model-based compressive sensing[END_REF]. The complexity of this sorting process is O( P p=1 Np + Np log(Np)) for each p. In the following, we will denote by Kp the set of Kp indexes corresponding to the selected rows. The matrix Xp(:, Kp), is thus Xp where only the column indexed by Kp are kept.

Fig. 1. Evolution of the RMSE (± standard deviation) with the percentage of removed nodes for the three dimensions.

Selection of the best subset of nodes

In this second step, we want to choose the best subset of nodes. This choice is based on the frequency components kept in the previous step. Considering a graph from a graph product, selecting the best nodes of each graph is equivalent to selecting a subset of rows from each associated Xp. To this end, we introduce the selection matrices {Sp} P p=1 , respectively in R Kp×Np , with elements satisfying: Sp(i, j) ∈ {0, 1}, ∀i, j Sp(i, j) = 1, and ∀j, i Sp(i, j) ≤ 1 [START_REF] Antonio G Marques | Sampling of graph signals with successive local aggregations[END_REF]. The sampled graph signal Ys lies in

R K 1 ו••×K P instead of R N 1 ו••×N P and is defined as Ys = Y P × p=1
Sp. Y can be approximated from Ys according to the following equation:

Ŷ = Ys P × p=1 Xp(:, Kp)(SpXp(:, Kp)) -1 . (9)
The only thing left to do is to properly choose the {Sp}p. For each dimension, one related optimization problem is:

minimize

Sp∈R Kp ×Np Y -Ys P × p=1 Xp(:, Kp)(SpXp(:, Kp)) -1 2 F . (10)
This optimization program is computationally expensive because of the particular structure of Sp. For multi-domain graphs, all methods in the literature for the selection of the rows are approximate: these iterative methods select one node per iteration in a greedy procedure. There are different existing works such as D-optimality [START_REF] Shamaiah | Greedy sensor selection: Leveraging submodularity[END_REF] and frame potential [START_REF] Krause | Nearoptimal sensor placements in gaussian processes: Theory, efficient algorithms and empirical studies[END_REF]. Following these approaches, we use a lowcomplexity greedy algorithm to sample signals that reside on the nodes of a product graph [START_REF] Ortiz-Jiménez | Sparse sampling for inverse problems with tensors[END_REF].

Learning a Cartesian graph structure

Although the recovery of the underlying graph is not mandatory for subsampling, it can help to visualize the data in a compact form and to assess the relevance of the selected nodes in an application perspective. In this context, since the Xp have been computed in the optimization problem ( 5), the only missing information is linked to the eigenvalues Λp. In the univariate case, a procedure for learning a Laplacian matrix from the graph Fourier basis and a sparse spectral representation has been introduced by Sardellitti et al. [START_REF] Sardellitti | Graph topology inference based on sparsifying transform learning[END_REF]. We propose here to learn a Cartesian graph with the {Xp}p as eigenvectors, hence extending this strategy to the case of multidomain signals. For the multidomain framework, we cast this problem as:

where H K is the tensor H where the "non-selected rows" of each dimension are removed, and each Cp ∈ R Kp×Kp is a matrix relative to the sparsity of H. Note that, H K , P p=1

H K ×p Cp = vec(Y) T P p=1
Lpvec(Y), where vec(•) denotes the vectorization.

Hence, we look for a Cartesian graph on which vec(Y) is smooth.

For more details on the optimization program and the additional matrices {Cp}p, the readers shall refer to the aforementioned paper. We reformulate the first term of the objective function as:

H K , P p=1 H K ×p Cp = P p=1 tr(H (p) T K CpH (p) K ) , (12)
where

H (p)
K is the tensor H K unfolded on the dimension p. Thus, it turns out that because of the use of a Cartesian graph, the full objective function naturally splits into P separable functions i.e. P separable optimization problems for (Lp, Cp) with fixed {(Lp, Cp)} p = p. Therefore, we can simply solve each independent problem:

minimize L p ,C p ∈R K p ×K p tr(H ( p) T K CpH ( p) K ) + µ Lp 2 F s.t.
Lp ∈ L, tr(Lp) = Np LpXp(:, Kp) = Xp(:, Kp)Cp, Cp 0 .

(

) 13
This last formulation is the same as in Sardellitti et al. [START_REF] Sardellitti | Graph topology inference based on sparsifying transform learning[END_REF] (univariate case) where they proposed to solve it with CVXPY [START_REF] Diamond | CVXPY: A pythonembedded modeling language for convex optimization[END_REF].

EXPERIMENTS

Our method is now tested on synthetic and real data and compared to other non-adaptive sampling strategies. For each simulation, we compute the Root-Mean-Square-Error (RMSE) of the reconstruction as a function of the percentage of removed nodes. The code to reproduce all the subsequent results is available at https://github.com/ Tgnassou/adaptive subsampling multidomain graph.

Synthetic data

We consider a Cartesian graph constructed from three sub-graphs G1, G2, and G3 of different sizes. As the bandlimited property is relative to the clusters of a graph [START_REF] Sardellitti | Graph topology inference based on sparsifying transform learning[END_REF], we generate them using the Stochastic Block Model (SBM), one of the most widely studied generative model that exhibits community structure [START_REF] Paul W Holland | Stochastic blockmodels: First steps[END_REF]. The sizes of the clusters, and the edge densities connecting clusters are parameters of SBM. In this experiment, each graph is composed of 3 clusters of respectively 4, 5 and 6 nodes. Dimensions of the three graphs are thus respectively N1 = 12, N2 = 15 and N3 = 18. Two nodes of the same group have a probability of 0.9 to be connected together and two nodes of different group have a probability of 0.2 to be connected together. Then, we generate noisy compressed bandlimited signals. They are known to provide more flexibility for the low-frequency component, hence reflecting better the data structure exhibited in several application contexts [START_REF] Chen | Signal recovery on graphs: Random versus experimentally designed sampling[END_REF]. Graph signals are generated as Y = H ×1 X1 ×2 X2 ×3 X3 + E where E ∼ N (0, 0.01) and H(k1, k2, k3) ∼ N (1, 0.5) for all k1 ≤ K1 = 3, k2 ≤ K2 = 3 Fig. 2. Learned Graphs of the representatives and the votes, where the larger yellow nodes are the selected nodes.

and k3 ≤ K3 = 3. For the others indexes, H(k1, k2, k3) is equal to the outer product of three vectors w1,w2 and w3, where wi(ki) = 1 for all ki ≤ Ki and wi(ki) = (Ki/ki) 2β for ki > Ki with β = 2 and i = {1, 2, 3}. With this construction, the values on frequencies upper than (K1, K2, K3) drastically decrease to zero. Method. We compare the reconstruction performances to those obtained by a correlation graph, the real graph (that was used for data generation), and a random graph (Erdős-Rényi). The quality of the reconstruction is measured by the RMSE between the learned signal and the noiseless real signal (Y -E). For each method, we reconstruct the signal with a fixed number of nodes and compare the RMSE performances. For sake of clarity, we display the results on Figure 1 dimension per dimension. Note that when the graph is known (which is the case for the three competitor methods), the optimization problem (5) reduces to the learning of H.

Results. Evolution of the RMSE with respect to the percentage of removed nodes is displayed in Figure 1 (average on 10 simulations). The importance of the structure of the graph clearly appears. The random graph and the correlation graph do not allow the suppression of nodes since they do not capture the structure of the data. For these graphs, the RMSE grows linearly with the percentage of removed nodes. On the other hand, our adaptive procedure efficiently learns a graph structure that is appropriate for efficient subsampling. Thanks to the learned graph, it is possible to remove more than 60% and still have a good reconstruction with a RMSE close to zero. In fact, for each dimension our learned graph allows to select only 3 nodes (reflecting the 3 clusters in the graphs and the Ki parameters), and still have a RMSE equals to 0.03. These curves underline another surprising property: the RMSE obtained with our graph is sometimes actually better than the one obtained with the true graph. This may be due to the fact that the bandlimitness constraint used in our method does not necessary learns the original graph but rather a graph structure that is optimised for sampling. Note that during our experiments, we also tried with other graphs such as Erdős-Rényi or Barabási graphs and we obtained similar results.

Real data

Illinois House of Representation. In this part, we measure the effectiveness of our method on a real dataset consisting on the votes of the Illinois House of Representatives during 2015-2017. This dataset is available at voteview.com. There are N1 = 18 seats at the Illinois House of representatives, one for each district. Each seat corresponds to a US Representative belonging either to The Democratic or The Republican parties. In our study, we used the N2 = 50 first available votes. Each vote is represented by a vector of size 18 with values +1 (Yes) and -1 (No).

In this case we have two dimensions, one for the US Representatives, another for the votes. We use our algorithm to compute the Fourier basis, select the nodes and learn the graph structure. Figure (2) shows the two graphs (one per dimension). Two clusters appear in the US representatives graph corresponding to Democrats and Republicans. For the votes graph, three clusters appear corresponding to questions with the same answer distribution (e.g. questions where all Republicans answer Yes and all Democrats answer No). As expected, we see that the selected nodes (in yellow) are coherent because we selected a node in each cluster: two nodes for the first graph and three nodes for the second graph.

Brain Computer Interface. This data set contains ElectroEn-cephaloGram data (EEG) recorded during a Brain Computer Interface (BCI) experiment consisting on 4 different motor imagery: movement of left hand, movement of right hand, movement of both feet, and movement of tongue. Each action is imagined and repeated 10 times (leading to N1 = 4 • 10 = 40 trials). The EEG is recorded with N2 = 22 electrodes and each time signal contains N3 = 188 samples. For a much more complete description of the data set see [START_REF] Brunner | BCI competition 2008-Graz data set A[END_REF]. The signal is therefore a three dimensional tensor.

Figure 3 displays the evolution of the RMSE as a function of the percentage of removed nodes for each dimension for our method and for a random graph (naive sampling). The sampling performances differ according to the dimension of interest (trials, space/electrodes, time samples). For the trial dimension, subsampling is tricky since as soon as a node is removed, the reconstruction is largely degraded, which is probably due to the fact that we work on the raw signals that do not allow to correctly represent the different types of movements. In contrast, it is possible to remove 82% of the electrodes and 75% of the time samples and to obtain a global RMSE of 2.8 (which corresponds to the inflexion point on the plots). The graph of electrodes learned in this configuration is displayed on the bottom right of Figure 3. The nodes which are selected for sampling are in yellow.

CONCLUSION

In this paper, we proposed a method to adaptively sample multidomain signals. By merging results from GSP, tensors, and graph products, we showed that it is possible to learn a graph structure that allows sampling and a useful data modelling and visualization. Using the learned structure we were able to deduce the number of nodes that could be deleted, select the best remaining nodes and accurately reconstruct the graph signal. This method also offers nice perspectives for the processing of datasets of N multidomain signals or size N1 × . . . × NP . Indeed, by constructing a new tensor of dimension N × N1 × . . . × NP (where the first dimension of the tensor serves to store the different signals), our algorithm can directly be used.

  s.t.Lp ∈ L, tr(Lp) = Np ∀p LpXp(:, Kp) = Xp(:, Kp)Cp, Cp 0 ∀p ,

Fig. 3 .

 3 Fig.3. Evolution of the RMSE with the percentage of removed nodes for the dimension of trials, the spatial dimension and the time dimension (top and bottom right). Electrodes learned graph (bottom right). In yellow, the four selected nodes/channels.


Yvon Maday
email: maday@ann.jussieu.fr

Julien Salomon
email: salomon@ceremade.dauphine.fr

CEREMADE † J Salomon

A greedy algorithm for the identification of quantum systems

The control of quantum phenomena is a topic that has carried out many challenging problems. Among others, the Hamiltonian identification, i.e, the inverse problem associated with the unknown features of a quantum system is still an open issue. In this work, we present an algorithm that enables to design a set of selective laser fields that can be used, in a second stage, to identify unknown parameters of quantum systems.

Introduction

The possibility to use coherent light to manipulate molecular systems at the nanoscale has been demonstrated both theoretically [START_REF] Beauchard | Local controllability of a 1D Schrödinger equation[END_REF] and experimentally [START_REF] Weinacht | Controlling the shape of a quantum wavefunction[END_REF]. Different types of methods have proven their relevancy for various settings, ranging from electron to large polyatomic molecules [START_REF] Bris | Towards efficient numerical approaches for quantum control[END_REF][START_REF] Levis | Selective bond dissociation and rearrangement with optimally tailored, strong-field laser pulses[END_REF][START_REF] Rice | Optimal Control of Quatum Dynamics[END_REF][START_REF] Shenvi | Nonlinear kinetic parameter identification through map inversion[END_REF][START_REF] Tadi | Explicit method for parameter identification[END_REF]. At the same time, the ability to generate a large amount of quantum dynamics data in a small time frame can also be used to extract from experiments the values of unknown parameters of quantum systems. The corresponding inverse problem, usually called Hamiltonian identification has recently been subject to significant developments through encouraging experimental results [START_REF] Geremia | Optimal identification of Hamiltonian information by closed-loop laser control of quantum systems[END_REF]. Various formulations in an optimization settings have been studied. Because of the nature of the available data, zero order methods were first tested, see e.g. the technique of map inversion [START_REF] Shenvi | Nonlinear kinetic parameter identification through map inversion[END_REF]. The use of optimal control techniques was then introduced [START_REF] Bris | Hamiltonian Identification for Quantum Systems: Well-posedness and Numerical Approaches[END_REF][START_REF] Geremia | Optimal Hamiltonian identification: The synthesis of quantum optimal control and quantum inversion[END_REF]. Contrary to this last class of methods, we present in this work a methodology that enables to handle situations where the experimental measurements are provided only at a given time. Our approach is based on a precomputation that provides a family of selective laser fields. Roughly speaking, these laser fields are designed iteratively to highlight variations in the parameters that are subject to the identification. In a second stage, these fields and the experimental measurements are used to assemble a nonlinear system satisfied by the to-be-identified parameters. The paper is organized as follows: the optimization framework and the assumptions we use are presented in Sec. 2. In Sec. 3, the structure of our algorithm is given. The procedures used in the two parts of this algorithm are described in Sec. 4 and Sec. 5. The identification step is explained in Sec. 6. Details about practical implementation and some numerical results are given Sec. 7. We conclude with some remarks in Sec. 8. Throughout this paper, Ω is a spacial domain in R d , d = 1, 2, 3, L 2 denotes the space of complex valued square integrable functions over Ω, and ., . the usual Hermitian product associated to L 2 . The following standard convention is used:

a|O|b := a, O(b) , a ∈ L 2 , b ∈ L 2 , O ∈ L(L 2 ; L 2 )
the set of all linear operator from L 2 into L 2 . Finally, we use ℜ(z), ℑ(z) to denote respectively the real and the imaginary part of a complex number z.

The identification problem

We first introduce the model and the framework used in this paper.

Control of the Schrödinger Equation

Consider a quantum system ψ ∈ H 1 , with norm ψ L 2 = 1, evolving according to the Schrödinger equation

i ψ = [H 0 + V + ε(t)µ]ψ ψ(0) = ψ 0 , (1)
where H 0 is the kinetic energy operator, V ∈ L(L 2 ; L 2 ) the potential operator and µ ∈ L(L 2 ; L 2 ) the dipole moment operator coupling the system to a timedependent external laser field ε(t). In this context, ε reads as a control since it can be chosen by the experimenter.

In the settings we consider here, we assume that the internal Hamiltonian H = H 0 + V is known so that the goal is to identify the dipole moment operator µ. The generalization to the identification of V should not give rise to any particular problem and is left to a future contribution.

The basic hypothesis made on µ is that it belongs to (or actually can be conveniently approximated by) a finite dimensional space spanned by some basis set B µ = (µ ℓ ) ℓ=1,...,L .

Experimental measurements and controllability

In order to perform the identification, we assume that given a time T and a laser field ε ∈ L 2 (0, T ), the experimenter can measure, for some fixed state ψ 1 ∈ L 2 , with norm ψ 1 L 2 = 1, the value ϕ(µ, ε) := ψ 1 , ψ(T ) . Note that all what follows still holds when considering several measurements a time T , i.e., in the case where a set of measurement ( ψ ℓ , ψ(T ) ) ℓ=1...,p , with p > 1, is known. Finally, we assume that the system under consideration is wavefunction controllable, i.e., that ε ∈ L 2 (0, T ) → ψ(T ) is surjective.

Formulation of problem

Our identification method is based on a particular formulation of the identification problem that we now briefly introduce. Denote by µ ⋆ the actual dipole moment operator of a given system. The solution µ = µ ⋆ of our problem also solves the minimization problem:

inf µ∈L(L 2 ;L 2 ) sup ε∈L 2 (0,T ) |ϕ(µ, ε) -ϕ(µ ⋆ , ε)| 2 .
(

) 2
This settings highlights the fact that as long as µ = µ ⋆ , a selective laser field should be designed so that the difference between µ and µ ⋆ is discerned through the measurement ϕ(µ, ε).

Structure of the algorithm

Our algorithm consists in designing, through a finite iterative procedure, a set of selective laser fields. We start with the general structure of our algorithm. Details about its steps are given in the next sections.

The selective laser fields computation greedy algorithm

Starting from the basis set B µ = (µ ℓ ) ℓ=1,...,L , the algorithm builds up iteratively a set of L selective laser fields as follows.

Algorithm 1 (Selective laser fields computation greedy algorithm) Let us define ε 1 a laser field that solves the problem:

sup ε∈L 2 (0,T ) |ϕ(µ 1 , ε)| 2 .
Suppose now that at the step k, with 1 < k ≤ L, a laser field ε k-1 is given. The computation of ε k is performed according to the two following sub-steps:

1. Fitting step : Find (α k j ) j=1,...,k-1 that solves the problem:

                 ϕ( k-1 j=1 α k j µ j , ε 1 ) = ϕ(µ k , ε 1 ) . . . ϕ( k-1 j=1 α k j µ j , ε m ) = ϕ(µ k , ε m ) . . . ϕ( k-1 j=1 α k j µ j , ε k-1 ) = ϕ(µ k , ε k-1 ), (3)
in the minimum mean square error sense.

2. Discriminatory step : Find ε k that solves the problem:

ε k = argmax ε∈L 2 (0,T ) |ϕ(µ k , ε) -ϕ( k-1 j=1 α k j µ j , ε)| 2 .
The initialization of the algorithm is somehow arbitrary, the only requirement is that ε 1 has a link with the type of measurement. In our case, we decide to maximize it.

Remark 1 Note that, in opposition to usual approaches (see e.g. [START_REF] Geremia | Optimal Hamiltonian identification: The synthesis of quantum optimal control and quantum inversion[END_REF][START_REF] Bris | Hamiltonian Identification for Quantum Systems: Well-posedness and Numerical Approaches[END_REF]), our method plays the role of a precomputation step since the actual measurements ϕ(µ ⋆ , ε) are not required at this stage.

Intuitive interpretation of the algorithm

In the first sub-step of an iteration of Algorithm 1, one looks for a defect of selectivity of the current laser fields ε 1 , . . . , ε k-1 : in the case the minimum reaches zero, two distinct dipole moment operators give rise to two identical measurements when exited with the laser fields ε 1 , . . . , ε k-1 . On the contrary, the second sub-step aims at computing a laser field that compensates this defect. These two sub-steps corresponds respectively to the minimization part and to the maximization part of the formulation (2).

Remark 2 Even if no hierarchy is assumed in the basis B µ , this algorithm should be viewed as a first step towards future works that handle infinite dimensional systems. In such a framework, the sum k-1 j=1 α k j µ j would read as an asymptotic expansion of the dipole moment operator.

This algorithm belongs to the class of greedy algorithms, since it follows the problem-solving's heuristic of making the locally optimal choice (in the second sub-step) at each stage with the hope of finding the global optimum that solves [START_REF] Bris | Towards efficient numerical approaches for quantum control[END_REF].

Let us first focus on the first sub-step of the algorithm. Consider an integer k such that 1 < k ≤ L and denote by K k the functional (defined on R k-1 ):

K k (α) = k-1 m=1 |ϕ(µ k , ε m ) -ϕ( k-1 j=1 α j µ j , ε m )| 2 .
During this sub-step, one has to find the minimum of the cost functional K k . To do this, a standard global minimization algorithm applied to this minimum mean square error associated problem. Note that, for small values of L, the gradient of the functional K k can be computed thanks to the formula:

∇K k (α).δα = k-1 m=1 2ℜ ψ α ε m (T ) -ψ k ε m (T ), ψ 1 δψ α ε m (T ), ψ 1 ,
where ψ α ε m and ψ k ε m are the solutions of Eq. ( 1) with ε = ε m as laser field, and µ = k-1 j=1 α j µ j and µ = µ k respectively as dipole moment operator. The variation δψ α is computed thanks to:

       iδ ψα ε m = ε k-1 k-1 j=1 α j µ j δψ α ε m +[H 0 + V + ε k-1 (t) k-1 j=1 δα j µ j ]ψ α ε m δψ α ε m (0) = 0.
In this way the computation of the components of ∇K k (α) can be parallelized to make the use of gradient methods feasible.

Discriminatory step

To achieve the second sub-step of Algorithm 1, we adapt an efficient strategy usually used in in quantum control. This strategy has given rise to a large class of algorithms often called "monotonic schemes". For a general presentation of these algorithms, we refer to [START_REF] Salomon | A monotonic method for solving nonlinear optimal control problems[END_REF].

Improvement of the selectivity of a given laser field

Let us present in more details how this strategy applies in our case. Note first that, given a laser field ε ∈ L 2 (0, T ), and two dipole moment operators µ and µ, one has:

|ϕ( µ, ε) -ϕ( µ, ε)| 2 = ψ(T ) -ψ(T )|O ψ1 | ψ(T ) -ψ(T ) ,
where O ψ1 = ψ 1 .ψ T 1 , ψ and ψ are the solutions of Eq. ( 1) with respectively µ = µ and µ = µ as dipole moment operator.

In order to compare the selectivity of ε and ε ′ , we introduce the functional:

J(ε) = ψ(T ) -ψ(T )|O ψ1 | ψ(T ) -ψ(T ) -β T 0 ε 2 (t)dt,
which has to be maximized. For sake of simplicity, we omit the dependence of J with µ and µ in the notations. The additional term β T 0 ε 2 (t)dt, is introduced for two complementary reasons: first, as it penalizes the L 2 -norm of the laser field, it enables to obtain feasible laser fields and secondly, it improves the convergence of Algorithm 2 below. Consider now another laser field ε ′ ∈ L 2 (0, T ), and denote by ψ ′ and ψ ′ the corresponding solutions of Eq. ( 1) with µ = µ and µ = µ respectively. We introduce the two adjoints states defined by:

i ˙ χ = [H 0 + V + ε(t) µ] χ χ(T ) = O ψ1 ψ(T ) -ψ(T ) , (4)
and i ˙ χ = [H 0 + V + ε(t) µ] χ χ(T ) = O ψ1 ψ(T ) -ψ(T ) . (5)
One has:

J(ε ′ ) -J(ε) = δψ ′ (T ) -δψ(T )|O ψ1 |δψ ′ (T ) -δψ(T ) +2ℜ δψ ′ (T ) -δψ(T ), χ(T ) -χ(T ) -β T 0 ε ′2 (t) -ε 2 (t)dt = δψ ′ (T ) -δψ(T )|O ψ1 |δψ ′ (T ) -δψ(T ) + T 0 (ε ′ (t) -ε(t)) 2ℑ χ(t)| µ| ψ ′ (t) -2ℑ χ(t)| µ| ψ ′ (t) -β (ε ′ (t)+ε(t)) dt, (6)
where we denote δψ ′ (T ) = ψ ′ (T )ψ ′ (T ) and δψ(T ) = ψ(T )ψ(T ). Identity [START_REF] Levis | Selective bond dissociation and rearrangement with optimally tailored, strong-field laser pulses[END_REF] gives a criterion to guarantee that ε ′ is more selective than ε. Indeed, suppose that ε ′ satisfies for all t ∈ [0, T ] the condition:

(ε ′ (t) -ε(t)) 2ℑ χ(t)| µ| ψ ′ (t) -2ℑ χ(t)| µ| ψ ′ (t) -β(ε ′ (t) + ε(t)) ≥ 0, (7)
then J(ε ′ ) ≥ J(ε).

Various ways to ensure that (7) holds. For example [START_REF] Tannor | Control of photochemical branching : Novel procedures fornding optimal pulses and global upper bounds[END_REF], one can define ε ′ at each time t as the solution of the equation:

ε ′ (t) -ε(t) = θ β 2ℑ χ(t)| µ| ψ ′ (t) -2ℑ χ(t)| µ| ψ ′ (t) -β (ε ′ (t) + ε(t)) , (8)
where θ is a given strictly positive number. In this case, one has:

J(ε ′ ) -J(ε) = δψ ′ (T ) -δψ(T )|O ψ1 |δψ ′ (T ) -δψ(T ) + β θ T 0 (ε ′ (t) -ε(t)) 2 dt ≥ 0,
which is the desired conclusion. In Sec. 7.1, we present an alternative that can be obtained in a time discretized settings.

Discriminatory sub-algorithm

We derive form the previous considerations the following iterative procedure to define a laser field ε k that maximizes J(ε):

Algorithm 2 (Discriminatory sub-algorithm) Let T ol be a positive number. Consider an initial guess ε k 0 and compute the corresponding solutions of Eq. (1) with µ and µ, say ψ 0 and ψ 0 . Set err = 2.T ol. While err > T ol, do:

1. Use Eqs. [START_REF] Geremia | Optimal identification of Hamiltonian information by closed-loop laser control of quantum systems[END_REF][START_REF] Geremia | Optimal Hamiltonian identification: The synthesis of quantum optimal control and quantum inversion[END_REF] with ε = ε k ℓ , ψ = ψ ℓ and ψ = ψ ℓ , to compute χ ℓ and χ ℓ , respectively.

Compute simultaneously the laser field ε k

ℓ+1 and the states ψ ℓ+1 and ψ ℓ+1 the solutions of coupled system composed of Eq. ( 8) with χ = χ ℓ , χ = χ ℓ and Eq. (1) with µ = µ and µ = µ respectively.

3. ℓ ← ℓ + 1, err = |ε k ℓ+1 -ε k ℓ |.
In [START_REF] Maday | Monotonic time-discretized schemes in quantum control[END_REF], one shows that Eq. ( 8) has a solution and presents some efficient numerical nonlinear solvers to compute it.

Identification procedure

Once the L selective fields ε 1 , ..., ε L have been computed, one can use them experimentally to obtain the corresponding measurements ϕ(µ ⋆ , ε 1 ), ..., ϕ(µ ⋆ , ε L ).

The identification procedure consists then in finding the linear combination (α 1 , ..., α L ) that solves the following nonlinear system:

                 ϕ( L j=1 α j µ j , ε 1 ) = ϕ(µ ⋆ , ε 1 ) . . . ϕ( L j=1 α j µ j , ε k ) = ϕ(µ ⋆ , ε k ) . . . ϕ( L j=1 α j µ j , ε L ) = ϕ(µ ⋆ , ε L ). (9)
in the mean square sense. In this view, the standard global optimization procedure used for the first sub-step of algorithm can be applied to the associated problem.

Note that, in a finite-dimensional settings, the existence of a solution is guaranteed.

Numerical implementation and results

We give here details about the practical implementation of Algorithm 1, and show its efficiency on an example.

Numerical solvers

In order to solve numerically Eq. ( 1), we use the second order Strang operator splitting [START_REF] Strang | On the construction and comparison of difference schemes[END_REF]. Given M > 0, a time step ∆t such that M.∆t = T and an approximation ψ j of ψ(j.∆t) with j < M , this method leads in our case to the following iteration: ψ j+1 = e iH ∆t 2 e iεj µ∆t e iH ∆t 2 ψ j .

In the second sub-step of Algorithm 1, Discriminatory sub-algorithm 2 is adapted to this discrete settings. In this way, we consider the time-discretized version of the cost functional J:

J ∆t (ε) = ψ M -ψ M |O ψ1 | ψ M -ψ M -β∆t M-1 j=0 ε 2 j ,
where ε ∈ R M-1 . Fix now two discrete laser fields ε and ε ′ , one can then repeat the computation done in Sec. 5.1 to obtain:

J ∆t (ε ′ ) -J ∆t (ε) = δψ ′ M -δψ M |O ψ1 |δψ ′ M -δψ M +∆t M-1 j=0 ε ′ j -ε j 2ℑ χ j | µ ∆t (ε ′ j , ε j )| ψ ′ j -2ℑ χ j | µ ∆t (ε ′ j , ε j )| ψ ′ j -β ε ′ j + ε j , (11)
where the vectors χ, χ, ψ ′ and ψ ′ are computed using the iteration [START_REF] Salomon | Convergence of the time-discretized monotonic schemes[END_REF] with µ = µ ∆t (ε ′ j , ε j ) and µ = µ ∆t (ε ′ j , ε j ). These matrices are the approximations of µ and µ respectively defined by:

µ ∆t (ε ′ j , ε j ) = e -iH ∆t 2 e iε ′ j e µ∆t -e iεj e µ∆t i∆t(ε ′ j -ε j ) e iH ∆t 2 µ ∆t (ε ′ j , ε j ) = e -iH ∆t 2 e iε ′ j b µ∆t -e iεj b µ∆t i∆t(ε ′ j -ε j ) e iH ∆t 2 .
For the sake of simplicity, instead of solving the discrete version of Eq. ( 8), we compute ε ′ j using one step of a Newton optimization method applied to its corresponding term in the sum of Eq. ( 11). This strategy, and the one corresponding to Eq. ( 8) are presented in more details in [START_REF] Maday | Monotonic time-discretized schemes in quantum control[END_REF]. Their convergence are proven in [START_REF] Salomon | Convergence of the time-discretized monotonic schemes[END_REF].

Numerical test

Settings

To illustrate the ability of our approach, we consider a simple finite dimensional settings where H 0 , V and µ are 3 × 3 Hermitian matrices with entries in C and ψ(t) ∈ C 3 . The internal Hamiltonian we consider is:

H = 10 -2   1 0 0 0 2 0 0 0 4   .
Since Eq. ( 1) with such an internal Hamiltonian is generically controllable, we choose to define the basis B µ randomly so that the systems handled by our algorithm are almost surely controllable. In order to work in a general framework, we chose µ ⋆ also randomly. In our example, we consider: The states ψ 0 and ψ 1 are

µ ⋆ =   2.
ψ 0 =   1 0 0   , ψ 1 =   0 0 1   .
We choose T = 4000π, which corresponds to 20 periods of the transition associated to the smallest frequency of the system.

Algorithm parameters

The minimum mean square error problems (3)-( 9) are solved by standard pseudo-Newton solvers. In order to make a global search, we repeat the minimization 10 times with random initialization. The parameter β is adapted to make Algorithm 2 converge. In our case β = 10 -2 .

Numerical results

The precomputation is achieved by our algorithm in approximately 80 min CPU. The dipole moment operator is regained with a relative error

Concluding remarks

The Selective laser fields computation greedy algorithm presented in this paper shows a good efficiency in a general settings. However, there is some room for improvement of our strategy. First, the choice of the basis B µ could be improved, e.g. through an iterative procedure. Secondly, the experimental measurements could be used during the computation of the selective fields in order to design an online procedure. Lastly, some work has to be done to design a more specific approach to treat the first sub-step of the algorithm. The identification procedure presented in Sec. 6 would also certainly take advantage of such a study.

µ ⋆ -µ 2 µ ⋆ 2 ≈ 9 .Figure 1 :

 291 Figure 1: Selective laser fields obtained by Algorithm 2.

ACKNOWLEDGMENTS

The problem of identification in this context was raised during discussions with H. Rabitz from Princeton University and G. Turinici from Dauphine University, we thank them for helpful inputs. This work was supported by the french A.N.R, "Programme blanc C-Quid" and PICS CNRS-NSF collaboration between the Department of Chemistry, Princeton University, and University Paris Dauphine.