Jean-Marie Ramirez

Qiang Bai

Marie Péquignot #

Fabienne Becker

Alboukadel Kassambara

Alexandra Bouin

Vasiliki Kalatzis #

Marilyne Dijon-Grinand

John De Vos
email: john.devos@inserm.fr

Side scatter intensity is highly heterogeneous in undifferentiated pluripotent stem cells and predicts clonogenic self-renewal

Keywords: pluripotent stem cells, side scatter, pluripotency, clonogenic selfrenewal, heterogeneity, stem cell hierarchies Human-induced pluripotent stem cells -hFF: Human foreskin fibroblasts -PSCs: Pluripotent stem cells -hESCs: Human embryonic stem cells -

In culture, human pluripotent stem cells (PSCs) are phenotypically (for instance, SSEA3 expression level) and functionally (capacity to survive after single-cell dissociation) heterogeneous. We report here that the side scatter (SSC) signal measured by flow cytometry, a variable correlated with membrane irregularity and cell granularity, is very high in PSCs, even higher than in blood polymorphonuclear cells, and markedly heterogeneous.

Moreover, SSC intensity rapidly and strongly decreases upon PSC differentiation into any of the three germ layers. PSCs with high SSC (HSSC cells) or low SSC (LSSC cells) values both express pluripotency markers, but HSSC cells are characterized by more frequent simultaneous expression of the membrane pluripotency factors SSEA3, SSEA4, TRA-1-81, TRA-1-60 and CD24 and by higher mitochondrial content. Functionally, HSSC cells are more likely to generate colonies upon single cell passage than LSSC cells. SSC monitoring might provide a simple, but robust and rapid method to estimate pluripotency variations in culture and unveils a new phenotypic and functional heterogeneity in PSCs.

Introduction

Pluripotent stem cells (PSCs), either human embryonic stem cells (hESCs) or human induced pluripotent stem cells (hiPSCs), retain many features of the inner cell mass (ICM), particularly the expression of pluripotency transcription factors (such as POU5F1 and NANOG) [START_REF] Yu | Pluripotent stem cell lines[END_REF]. In addition, expression fluctuations [START_REF] Chambers | Nanog safeguards pluripotency and mediates germline development[END_REF][START_REF] Guo | Resolution of cell fate decisions revealed by single-cell gene expression analysis from zygote to blastocyst[END_REF][START_REF] Halley | Self-organizing circuitry and emergent computation in mouse embryonic stem cells[END_REF], not only at the mRNA but also at the protein level [START_REF] Gu | Aire regulates the expression of differentiation-associated genes and self-renewal of embryonic stem cells[END_REF][START_REF] Gu | Global expression of cell surface proteins in embryonic stem cells[END_REF][START_REF] Gu | Proteomic analyses reveal common promiscuous patterns of cell surface proteins on human embryonic stem cells and sperms[END_REF],

have been reported for pluripotency transcription factors, surface markers [START_REF] Gu | Aire regulates the expression of differentiation-associated genes and self-renewal of embryonic stem cells[END_REF][START_REF] Gu | Global expression of cell surface proteins in embryonic stem cells[END_REF][START_REF] Gu | Proteomic analyses reveal common promiscuous patterns of cell surface proteins on human embryonic stem cells and sperms[END_REF][START_REF] Cui | Spatial distribution and initial changes of SSEA-1 and other cell adhesion-related molecules on mouse embryonic stem cells before and during differentiation[END_REF][START_REF] Enver | Cellular differentiation hierarchies in normal and culture-adapted human embryonic stem cells[END_REF][START_REF] Hong | Cell fate potential of human pluripotent stem cells is encoded by histone modifications[END_REF][START_REF] Ramirez | Brief report: benchmarking human pluripotent stem cell markers during differentiation into the three germ layers unveils a striking heterogeneity: all markers are not equal[END_REF][START_REF] Stewart | Clonal isolation of hESCs reveals heterogeneity within the pluripotent stem cell compartment[END_REF], epigenetic marks [START_REF] Hong | Cell fate potential of human pluripotent stem cells is encoded by histone modifications[END_REF][START_REF] Hayashi | Dynamic equilibrium and heterogeneity of mouse pluripotent stem cells with distinct functional and epigenetic states[END_REF][START_REF] Tanasijevic | Progressive accumulation of epigenetic heterogeneity during human ES cell culture[END_REF] and lineage-specific genes [START_REF] Hough | A continuum of cell states spans pluripotency and lineage commitment in human embryonic stem cells[END_REF]. PSCs might thus permanently shift between cell compartments characterized by different levels of pluripotency [START_REF] Enver | Cellular differentiation hierarchies in normal and culture-adapted human embryonic stem cells[END_REF][START_REF] Stewart | Clonal isolation of hESCs reveals heterogeneity within the pluripotent stem cell compartment[END_REF]. The apparently stochastic fluctuation in the expression of pluripotency genes, cell surface glycoproteins or epigenetic marks might reflect the unstable pluripotency equilibrium of cultured PSCs and a better control of culture conditions could help minimizing these variations [START_REF] Chambers | Nanog safeguards pluripotency and mediates germline development[END_REF]. Alternatively, this heterogeneity may be inherent to cell fate decisions as it is also observed in state transitions during neurogenesis or hematopoiesis [START_REF] Arias | Gene expression heterogeneities in embryonic stem cell populations: origin and function[END_REF] and may even be a common biological feature [START_REF] Eldar | Functional roles for noise in genetic circuits[END_REF].

Flow cytometry is a powerful, quantitative method to study PSC diversity at the single cell level [START_REF] Ho | Flow cytometric analysis of human pluripotent stem cells[END_REF] and has significantly contributed to the characterization of the fluctuations of pluripotency cell surface markers in PSCs [START_REF] Enver | Cellular differentiation hierarchies in normal and culture-adapted human embryonic stem cells[END_REF][START_REF] Ramirez | Brief report: benchmarking human pluripotent stem cell markers during differentiation into the three germ layers unveils a striking heterogeneity: all markers are not equal[END_REF][START_REF] Stewart | Clonal isolation of hESCs reveals heterogeneity within the pluripotent stem cell compartment[END_REF]. In addition to the quantification of fluorescence emission (biomarker detection), other types of information can be derived from the analysis of flow cytometry data. Specifically, as the laser beam passes through the cell, light is deflected and refracted. This scattered light is collected by the forward scatter (FSC) and by the side scatter (SSC) photodiodes. FSC intensity mainly correlates with the cell size, whereas SSC is a measure of the cell refractive index that depends on the cell granularity or internal complexity (for instance, cytoplasmic membrane wrinkling, number and shape of vesicles and mitochondria, development of the endoplasmic reticulum and nucleus structure).

We recently reported that polychromatic flow cytometry can provide information on the expression kinetics of the cell surface pluripotency markers SSEA3, SSEA4, TRA-1-81 and TRA-1-60 [START_REF] Ramirez | Brief report: benchmarking human pluripotent stem cell markers during differentiation into the three germ layers unveils a striking heterogeneity: all markers are not equal[END_REF]. During that work, we noted that human PSCs had very heterogeneous SSC values. Here, we thoroughly investigated SSC intensity variations in different PSC lines and show that PSCs fluctuate between extreme SSC states. Moreover, PSCs with high SSC (HSSC cells) express more often many cell surface pluripotency markers simultaneously than cells with low SSC (LSSC). HSSC and LSSC cells show differences in cell cycle distribution. HSSC cells are more prone to apoptosis and display higher mitochondrial content. However, despite their higher apoptosis rate, HSSC cells have a much higher clonogenic self-renewal capacity than LSSC cells.

Material and Methods

Cell culture

The following hES cell lines were used: HUES1, from Douglas Melton's laboratory [START_REF] Cowan | Derivation of embryonic stem-cell lines from human blastocysts[END_REF] and HD83 and HD291 cells that we derived in our laboratory [START_REF] Ramirez | Brief report: benchmarking human pluripotent stem cell markers during differentiation into the three germ layers unveils a striking heterogeneity: all markers are not equal[END_REF][START_REF] Assou | A gene expression signature shared by human mature oocytes and embryonic stem cells[END_REF]; The hiPSC cell line M4C7 is one of a series of clones obtained by reprogramming the human fibroblast foreskin hFF1 cell line (ATCC) by using lentiviral vectors to express human OCT4/POU5F1, SOX2, NANOG and LIN28 [START_REF] Yu | Induced pluripotent stem cell lines derived from human somatic cells[END_REF]. M4C7 cells grow as typical PSCs and display all the features of PSCs: OCT4 and ABCG2 expression, phosphatase alkaline activity, a typical PSC expression signature by microarray analysis and can differentiate into cells of the three germ layers. All cell lines were routinely cultured in 35 mm wells in KO-DMEM medium supplemented with 20% Knockout serum replacement (KO-SR) (Invitrogen, Carlsbad, CA), 0.1mM non-essential amino acids, 2 mM L-Glutamine, 100 U/ml penicillin, 100 µg/ml streptomycin, 50 µM β-mercaptoethanol and 10 ng/ml bFGF (PSC culture medium) and were passaged mechanically (HD83 and HD291) or semi-mechanically using Collagenase IV (HUES1) or mechanically and then enzymatically using TrypLE Select (Invitrogen) (M4C7). Mouse embryonic fibroblasts (MEFs) or human foreskin fibroblasts (hFF1) were irradiated and used as feeders at 34 000 cells/cm 2 . HFF1 cells and MEFs were maintained in DMEM medium supplemented with 15% heat-inactivated fetal bovine serum (FBS) (Invitrogen). For low-serum condition, PSCs were grown in PSC culture medium with 0.1% instead of 20% KO-SR.

For in vitro differentiation, PSCs were pre-treated with 10 µM Y27632, a p160 Rho-associated kinase (ROCK) inhibitor (Sigma-Aldrich, St. Louis, MO), for 1 hour. Cells were then dissociated with 1X TrypLE Select (Invitrogen) at 37°C for 10 min, seeded at 5 000 cells/cm 2 in 96-well plates pre-coated with 1mg/ml Matrigel (endoderm and ectoderm differentiation) or with feeder cells (mesoderm differentiation). PSCs were again incubated with 10 µM Y27632 for 24 hours and then PSC differentiation into endodermal (HUES1 cells), mesodermal (HD291 cells) or ectodermal cells (M4C7 cells) was induced as previously described [START_REF] Ramirez | Brief report: benchmarking human pluripotent stem cell markers during differentiation into the three germ layers unveils a striking heterogeneity: all markers are not equal[END_REF]. Briefly, endodermal differentiation was induced by addition of activin A to low-serum medium [START_REF] D'amour | Efficient differentiation of human embryonic stem cells to definitive endoderm[END_REF]; ectodermal differentiation was induced by adding Noggin and SB431542, two specific inhibitors of SMAD signaling [START_REF] Chambers | Highly efficient neural conversion of human ES and iPS cells by dual inhibition of SMAD signaling[END_REF]; and mesoderm was obtained by co-culture with OP9 cells and VEGF [START_REF] Klimchenko | A common bipotent progenitor generates the erythroid and megakaryocyte lineages in embryonic stem cell-derived primitive hematopoiesis[END_REF]. Differentiation was accompanied by the sharp down-regulation of the pluripotency specific transcription factor OCT4 and by the induction of lineage-specific markers, such as CD45 (mesoderm), FOXA2 and CXCR4 (endoderm), Nestin and CD56 (ectoderm). In addition, hematopoietic differentiation (mesoderm model) was documented by the appearance of white blood cells, as indicated by May-Grünwald-Giemsa staining [START_REF] Ramirez | Brief report: benchmarking human pluripotent stem cell markers during differentiation into the three germ layers unveils a striking heterogeneity: all markers are not equal[END_REF].

Quantification of stem cell markers, mitochondrial content, apoptosis rate and cell cycle distribution by flow cytometry

For flow cytometry analyses, cells were dissociated with TrypLE at 37°C for 10 minutes. Cell surface pluripotency markers were quantified using a cocktail of five antibodies against CD24-PB (EXBIO), SSEA3-PE (BD), SSEA4-PERCP (R&D), TRA-1-60-FITC, (BD) and TRA-1-81-APC (BD) and the Sytox Blue Dead Cell Stain (Invitrogen).

Mitochondrial content was determined using MitoTRacker Red CMXRos probes (Invitrogen). PSC culture medium was removed from the 35 mm dishes and cells were incubated with 100 nM MitoTRacker probes in prewarmed (37°C) PSC medium without serum at 37°C for 30 minutes. Cells were then washed with PBS, dissociated with TrypLE, suspended in 500 µl PBS and incubated with 6 µl anti-CD44-APC antibody (BD) to exclude feeder cells.

To evaluate the cell cycle distribution, dissociated cells were fixed in cooled 70% EtOH at 4°C for 15 min, permeabilized using 0.1% Triton X-100 (v/v) and stained with 10 µg/ml propidium iodide plus 250 µg/ml RNase A for 30 min.

To quantify the apoptosis rate, TrypLE-dissociated cells were incubated with the anti-CD44 antibody for 30 min and then stained with Annexin-V-FITC (BD Biosciences) in Annexin-V binding buffer (BD) for 20 min, as recommended by the manufacturer. Sytox Blue stain was added to exclude dead cells.

Flow cytometry data acquisition was carried using a CyAn instrument (Beckman Coulter, Fullerton, CA) and data analyzed using the KALUZA software (Beckman Coulter).

FACS-sorting of HSSC and LSCC cells

PSCs were pretreated with 10 µM Y27632 for 1 h, stained with Sytox Blue (Invitrogen) and incubated with an anti-CD44 antibody (Pharmingen) to label feeder cells. LSSC and HSSC cells were then sorted using an Astrios instrument (debris and aggregated cells were eliminated by specific gating), plated in 96-well plates pre-coated with irradiated hFF1 cells and cultured in the presence of 10 µM Y27632 in PSC culture medium for 24 hours.

Teratoma formation

Animals were handled in strict accordance with the ARVO Statement for the Use of Animals in Ophthalmic and Vision Research and the EU directives (Council Directive 86/6009/EEC). NOD/SCID/IL2RG -/-mice (Charles River, France) were maintained in standard conditions (12 h 90 lux light and 12 h dark) and fed ad libitum with a standard rodent diet. 300 000 sorted HSSC and LSSC cells were centrifuged and pellets resuspended in 30% Matrigel diluted in KO-SR in a final volume of 200 µl per injection. Before subcutaneous injection of the cell mixture, mice were lightly anesthetized with 35 mg/kg ketamine (Merial, France) and 14 mg/kg xylasine (Bayer Healthcare, Germany), shaved on the left flank and disinfected with 70% ethanol. Teratoma appearance and growth were followed weekly.

Approximately eight weeks after the graft, teratomas were dissected, fixed in 3.7% formaldehyde, dehydrated and embedded in paraffin; 4-μm sections were processed for Hematoxylin-Eosin-Safran staining.

Transmission Electron Microscopy (TEM)

FACS-sorted HSSC and LSSC cells were immersed in a solution of 2.5% glutaraldehyde in 0.1M Sorensen's buffer (pH 7.4) at 4°C overnight. They were then rinsed in Sorensen's buffer and post-fixed in 0.5% osmic acid in the dark, at room temperature for 2 h. After two rinses in Sorensen's buffer, cells were dehydrated in a graded series of ethanol solutions (30-100%) and embedded in EmBed 812 using an Automated Microwave Tissue Processor for Electron Microscopy (Leica EM AMW). Thin sections (70 nm; Leica-Reichert Ultracut E) were collected at different levels of each block, counterstained with uranyl acetate and observed using a Hitachi 7100 transmission electron microscope at the Centre de Ressources en Imagerie Cellulaire de Montpellier (France)

Statistical analysis

Statistical comparisons were made using the Student's t-test and p-values <0.05 were considered to be significant. The SSC distribution of the different cell subpopulations was analyzed with the R software (R.2.13, http://www.rproject.org/) and Bioconductor, version 2.8 [START_REF] Gentleman | Bioconductor: open software development for computational biology and bioinformatics[END_REF]. The flowCore package [START_REF] Hahne | flowCore: a Bioconductor package for high throughput flow cytometry[END_REF] was used to import the flow cytometry data in R and the SSC distribution of each cell subpopulation is shown as dot plots.

Results

SSC intensity values in PSCs are very heterogeneous

To thoroughly evaluate the variations in SSC intensity in human PSCs, we assessed SSC intensity distribution in hES (Figure 1 andFigures 2B-D) and hiPS cell lines (Figure 2A) by flow cytometry. In this study, to quantify the PSCs, we excluded dead cells, aggregated cells and debris (Figure 1A). All PSC lines showed a wide SSC distribution, including a sub-population of cells with very high SSC intensity. Comparison with the SSC values of blood cells (assessed using the same settings) (Figure 1C) indicated that SSC intensity in PSCs was 6.5-fold higher than in blood mononuclear cells and even higher than in granulocytes, which are known to have very high SSC values (Figure 1D). In addition, the SSC values of PSCs were markedly heterogeneous (from very low to very high SSC intensity), suggesting that each PSC population includes cells with variable morphological features (internal and/or membrane complexity). To understand why PSCs display such a huge SSC variability, we investigated several features (morphology, expression of pluripotency markers, clonogenic self-renewal, apoptosis and cell cycle distribution) in PSCs with either high SSC (HSSC) or low SSC values (LSSC).

The size of the HSSC population is highly correlated with the pluripotency state

First, we tested whether loss of pluripotency and cell differentiation could be associated with a decrease in SSC intensity by measuring the size of the HSSC population in the four PSC lines following induction of differentiation [START_REF] Ramirez | Brief report: benchmarking human pluripotent stem cell markers during differentiation into the three germ layers unveils a striking heterogeneity: all markers are not equal[END_REF], or culture in low-serum conditions (0.1% KO-SR). The size of the HSSC populations was drastically decreased after only five days of differentiation into ectodermal (M4C7 hiPSC; 6-fold reduction), endodermal (HUES1 ES cells; 7-fold decrease) or mesodermal cells (HD291 ES cells; 7.8-fold reduction) (Figures 2A, B andC), and also when HD83 cells were induced to differentiate in low-serum medium (Figure 2D). As the decrease of the HSSC populations was a very early event during PSC differentiation, we investigated whether HSCC reduction was temporally associated with changes in the expression levels of the cell surface pluripotency markers SSEA3, SSEA4, TRA-1-60, TRA-1-81 and CD24 (Figures 2E, F andG). The decrease of the HSSC cell subset was the most remarkable change between day 0 and day 5 of differentiation in the ectoderm and endoderm differentiation models, whereas reduction in SSEA3 expression, one of the most sensitive pluripotency markers [START_REF] Enver | Cellular differentiation hierarchies in normal and culture-adapted human embryonic stem cells[END_REF][START_REF] Ramirez | Brief report: benchmarking human pluripotent stem cell markers during differentiation into the three germ layers unveils a striking heterogeneity: all markers are not equal[END_REF], came first only during mesoderm differentiation.

Hence, in two of the three differentiation pathways tested, the decrease of the HSSC population significantly exceeded in magnitude loss of SSEA3 expression. Conversely, SSEA4, TRA-1-81 and CD24 expression were not changed or only slightly decreased in all the differentiation models at day 5 of differentiation. These results indicate that SSC could represent a new parameter for monitoring loss of pluripotency during PSC differentiation.

HSSC cells express simultaneously different pluripotent markers

We previously showed that expression of five common cell surface pluripotency markers (SSEA3, SSEA4, TRA-1-81, TRA-1-60 and CD24) is highly variable within a PSC population and that only a portion of these cells expresses simultaneously all five markers [START_REF] Ramirez | Brief report: benchmarking human pluripotent stem cell markers during differentiation into the three germ layers unveils a striking heterogeneity: all markers are not equal[END_REF][START_REF] Assou | A meta-analysis of human embryonic stem cells transcriptome integrated into a web-based expression atlas[END_REF]. We thus asked whether HSSC and LSSC cells had distinct expression profiles of pluripotent markers.

HSSC cells from the hiPS cell line M4C7 contained a significantly higher percentage of SSEA3-positive cells than LSSC cells. In addition, 85-87% of HSSC cells expressed simultaneously SSEA3, SSEA4, TRA-1-81, TRA-1-60 and CD24, whereas only 51-54% of LSSC cells did (Figure 3A). Similar results were obtained with the hES cell lines HUES1 and HD291 (Figures 3B and3C), indicating that this finding was not limited to a specific PSC line and suggesting that high SSC may be correlated with pluripotency.

The clonogenic self-renewal capacity of PSCs is correlated with high SSC values.

To examine more closely the link between SSC intensity and pluripotency, we compared the developmental potential of HSSC and LSSC cells. To quantitate their clonogenic self-renewal capacity, PSCs were FACS-sorted into HSSC and LSSC cells (Figure 4A) and cultured at low density in 96-well plates pre-coated with irradiated hFFs for four days. Microscopy analysis revealed that HSSC cells had higher self-renewal capacity and generated larger colonies than LSSC cells (49 +/-6 and 9 +/-5 colonies per 500 cells, respectively) (Figures 4B andC). Moreover, the number of SSEA3-positive cells per well was higher in colonies generated from HSSC than from LSSC (Figure 4D). Despite their lower clonogenic self-renewal capacity, LSSC cells could form bona fide PSC colonies that contained both HSSC and LSSC cells, as indicated by flow cytometry analysis of the colonies after 5 days of culture (Figure 4E). This finding demonstrates that PSCs fluctuate between these extreme SSC states. Finally, to confirm the pluripotency status of the LSSC population, we grafted immuno-compromised mice with FACS-sorted LSSC cells or with a mixture of LSSC and HSSC cells [START_REF] Lensch | Teratoma formation assays with human embryonic stem cells: a rationale for one type of human-animal chimera[END_REF][START_REF] Prokhorova | Teratoma formation by human embryonic stem cells is site dependent and enhanced by the presence of Matrigel[END_REF]. Both cell preparations produced teratomas of similar size that included patches of endoderm, mesoderm and ectoderm (Figures 4G andH), indicating that LSSC cells are nevertheless pluripotent. Taken together, these results demonstrate that HSSC cells have a higher clonogenic self-renewal capacity, but that each SSC subset contains pluripotent cells that can generate both HSSC and LSSC cells.

HSSC and LSSC cells have distinct apoptosis rates and cell cycle profiles

To further characterize the biological features of HSSC and LSSC PSCs, we evaluated their apoptosis rate and cell cycle distribution. Indeed, during FACS sorting of the SSC populations, we noticed that HSSC cells were less viable than LSSC cells (data not shown). To quantify apoptosis and cell death, M4C7 cells were stained with Annexin-V and Sytox Blue, respectively, and then the percentages of early and late apoptotic cells, live cells and dead cells were quantified in the HSSC and LSSC populations (Figure 5A). Early and late apoptosis as well as cell death were significantly higher (2.4-, 4-and 6.7fold, respectively) in HSSC cells than in LSSC cells (Figure 5B).

Next, we investigated cell cycle progression in HSSC and LSSC cells. The percentage of cells in the G2/M phase was 2.1-fold higher, whereas the percentage of cells in G0/G1 and S phase was reduced by 1.5-and 2-fold, respectively, in HSSC in comparison to LSSC cells (Figures 5C andD).

HSSC and LSSC cells show distinct ultrastructural characteristics and mitochondrial content

As differences in SSC intensity are related to irregularities in the cell membrane, cytoplasm or nucleus, we compared the ultrastructure of HSSC and LSSC PSCs by transmission electron microscopy (TEM). The most notable differences between the two cell populations were that HSSC cells were characterized by many rounded mitochondria that were distributed around the nucleus, whereas LSSC cells had many lipid vesicles in the cytoplasm (Figures 6A andB). The higher mitochondrial content of HSSC cells was confirmed also by quantification of the mitochondrial-specific dye MitoTracker by flow cytometry (see Material and Methods). Indeed, MitoTracker accumulation was significantly higher in HSSC cells than in LSSC cells (Figure 6C).

Plural pluripotent stem cell compartments

The SSC parameter can thus delineate different PSC cell populations that do not overlap with the differential expression of the cell surface pluripotency marker SSEA3 (Figure 3). Indeed, both HSSC and LSSC populations contained SSEA3-positive and -negative cells. As differential expression of SSEA3, SSEA4, TRA-1-81, TRA-1-60 and CD24 can delineate several PSC subsets that co-exist in variable percentages in the different PSC lines 8 , we assessed the association of these different markers and the SSC parameter (Figures 7A, B andC). Like for SSEA3 alone, there was no overlap of the SSEA3, SSEA4, TRA-1-81, TRA-1-60 and CD24 compartments with the HSSC and LSSC populations, suggesting that SSC might represent another level of PSC heterogeneity. Moreover, many sub-populations expressed few pluripotency markers.

Discussion

The molecular definition of stem cells is an old goal that has been only partially met. Single adult stem cells cannot be isolated and they can only be enriched based on specific factors the expression of which has been reported to be heterogeneous [START_REF] Gu | Aire regulates the expression of differentiation-associated genes and self-renewal of embryonic stem cells[END_REF][START_REF] Gu | Global expression of cell surface proteins in embryonic stem cells[END_REF][START_REF] Gu | Proteomic analyses reveal common promiscuous patterns of cell surface proteins on human embryonic stem cells and sperms[END_REF][START_REF] Arias | Gene expression heterogeneities in embryonic stem cell populations: origin and function[END_REF]. This is also true for PSCs. For instance, SSEA3 is a cell surface marker strongly associated with pluripotency, but a PSC colony will contain both SSEA3-positive and -negative cells. Although SSEA3negative cells are less prone to self-renewal and appear primed for differentiation, they can still generate SSEA3-positive cells in vitro and are pluripotent [START_REF] Enver | Cellular differentiation hierarchies in normal and culture-adapted human embryonic stem cells[END_REF][START_REF] Stewart | Clonal isolation of hESCs reveals heterogeneity within the pluripotent stem cell compartment[END_REF]. Our work reveals that also SSC, a basic parameter that is routinely monitored when cells are analyzed by flow cytometry, is very heterogeneous in human PSCs and many cells have very high or very low SSC values. SSC heterogeneity was comparable in different human PSC lines (three ES and one hiPS cell line), although hiPSCs are often thought to be more heterogeneous than ES cells and thus less stable [START_REF] Narsinh | Single cell transcriptional profiling reveals heterogeneity of human induced pluripotent stem cells[END_REF]. As a matter of fact, expression of cell surface pluripotency markers in the hiPS cell line M4C7 appeared to be less variable than in the three ES cell lines analyzed.

The wide SSC distribution could not be explained by major morphological or phenotypically differences, as cycling and non-cycling cells, apoptotic and non-apoptotic cells, SSEA3-positive and -negative cells were detected in both HSSC and LSSC populations, although some variations were observed.

Neither it was explained by culture biases, such as the use of the Rock inhibitor Y27632. Indeed, cell lines that were cultured using mechanical passaging (HD83 cells) or semi-enzymatic passaging (HUES1cells) and that thus were not exposed to this inhibitor were also highly heterogeneous (see Figure 2D and2B). The only consistent difference was the higher number of round-shaped mitochondria near the nucleus and the 2.1-fold increase in mitochondrial content in HSSC in comparison to LSSC cells, suggesting that HSSC cells have a higher mitochondrial energy production than LSSC cells. This morphological feature could be related to random segregation of mitochondria at cell division, as previously suggested for the cell-to-cell variability in mitochondrial mass observed in a the Hela cell line [START_REF] Das Neves | Connecting variability in global transcription rate to mitochondrial variability[END_REF].

We then show that the size of the HSSC population diminishes dramatically and rapidly (day 5) following directed PSC differentiation or serum deprivation. SSC could thus be used as a simple, rapid and inexpensive tool to follow loss of pluripotency during cell culture or differentiation of PSCs, although it cannot substitute the analysis of pluripotency marker expression or functional assays, such as teratoma formation.

Although LSSC cells express lower amounts of pluripotency markers (suggesting that differentiation takes place essentially in this compartment), both HSSC and LSSC populations can form new PSC colonies, which contain HSSC and LSSC cells, and teratomas. In the model proposed by Enver et al, the SSEA3-positive ES compartment gives rise also to SSEA3-negative cells that can revert to the SSEA3-positive compartment or retain their multilineage differentiation potential and then differentiate [START_REF] Enver | Cellular differentiation hierarchies in normal and culture-adapted human embryonic stem cells[END_REF][START_REF] Stewart | Clonal isolation of hESCs reveals heterogeneity within the pluripotent stem cell compartment[END_REF]. Here we show that both HSSC and LSSC compartments contain SSEA3-positive andnegative cells (Figure 7C) and that LSSC cells can generate HSSC cells, albeit at a lower rate than HSSC cells. PSCs might thus fluctuate between compartments, as suggested by our findings (see Figure 4) and similar experiments based on sorting for SSEA3 [START_REF] Enver | Cellular differentiation hierarchies in normal and culture-adapted human embryonic stem cells[END_REF]. Overall, there might be a permanent shuffling between compartments, with a preferential exit toward apoptosis for cells with the highest SSC and toward differentiation for cells with the lowest SSC values (Figure 7D).

Given the well-known PSC heterogeneity for some specific markers, it is surprising that SSC heterogeneity was never previously reported. As PSC heterogeneity may be a fundamental feature of undifferentiated cells [START_REF] Chang | Transcriptome-wide noise controls lineage choice in mammalian progenitor cells[END_REF] and especially of pluripotent stem cells [START_REF] Arias | Gene expression heterogeneities in embryonic stem cell populations: origin and function[END_REF], their high and variable SSC values may have important biological implications, as suggested by the higher mitochondrial content and clonogenic self-renewal capacity of the HSSC population.

Figure 1 .

 1 Figure 1. Analysis of the side scatter parameter in pluripotent stem cells

Figure 4 :

 4 Figure 4: Clonogenic self-renewal efficiency of HSSC and LSSC PSCs.

Figure 5 :

 5 Figure 5: Cell cycle distribution and apoptosis rate in HSSC and LSSC

Figure 6 :

 6 Figure 6: Ultrastructural analysis and mitochondrial content of HSSC

(

  A and B) Analysis of the ultrastructure of both HSSC and LSSC cells by transmission electron microscopy reveals specific characteristics for each subset: high level of rounded mitochondria in HSSC cells (asterisks) and many lipid vesicles in the cytoplasm of LSSC cells (arrows). Cells were FACSsorted and aggregates, debris, dead cells and feeder cells were excluded as detailed in Figure 4. (C) Flow cytometry using MitoTracker probes indicates that HSSC cells have a significantly higher mitochondrial content than LSCC cells. Data are the mean ± SEM of three independent experiments.

Figure 7 :

 7 Figure 7: Plural pluripotent stem cell compartments.

Acknowledgement

Supported by grants from the Association Française contre les Myopathies (AFM), the Région Languedoc-Roussillon (Chercheur d'Avenir 09-13198 01) and the Agence Nationale de la Recherche (ANR-07-BLAN-0076-01, ANR-08-BIOT-012-02). We thank E. Andermarcher for critical reading of this manuscript. We thank Montpellier RIO Imaging (MRI) for access to the flow cytometry department, the Centre de Ressources en Imagerie Cellulaire (CRIC) for TEM and the Plateau Animalerie Inserm U583, Institut des Neurosciences de Montpellier (INM) for mice husbandry. The authors are grateful to Chantal Cazevieille and Cécile Sanchez for their technical assistance and interpreting TEM images. We thank the "Réseau d'Histologie Expérimentale de Montpellier (RHEM)" histology facility for processing our teratoma samples and Pascal Roger for the interpretation of the Hematoxylin-Eosin-Safran stained sections.

Author Disclosure Statement

All authors disclose any potential conflicts of interest. SSC: side scatter; FSC: forward scatter; PMN: granulocytes; Mono: mononuclear cells; Lympho: lymphocytes; Leu: leukocytes; MNC: mononuclear cells.


Olivier Meste

Frédéric Brau

Alice Guyon

A New Method for the Robust Estimation of the Motile Cilia Beating Frequency

A New Method for the Robust Estimation of the Motile Cilia Beating Frequency O. MESTE 1 , F. BRAU 2 , A. GUYON 2 1 Lab I3S, CNRS University of Nice-Sophia Antipolis 2 IPMC, CNRS University of Nice-Sophia Antipolis

Abstract-Motile cilia play a major role in fluid circulation in many biological systems such as the oviduct, the airway and the cerebral ventricles and defects in cilia beating lead to pathologies called ciliopathies. The estimation of the cilia beating frequency (CBF) is of great interest to understand how the CBF modulates the liquid fluxes and how CBF is controled by the intraand/or extra-cellular medium composition of the ciliated cell in a physiological process. Such estimation is almost impossible to be computed on a single cilium because of overlapping patterns. In addition, motion artifacts and camera defaults may hinder or bias the computation of the frequency variations during long lasting experiments. In this paper is provided a methodological development matching the cilia properties. It is first based on a grid removal and a shift compensation. Second, an harmonic model of the observed cilia through image recordings is proposed within a Maximum Likelihood Estimator framework. It is proven that in contrast to classical approaches the best estimation of the frequency is obtained by averaging the squared Fourier transform of individual pixels followed by a particular summation over the frequency, namely the Compressed Spectrum. Illustrative examples show the performances of this robust method compared to classical ones.

I. INTRODUCTION

As mentioned in [START_REF] Salathe | Regulation of mammalian ciliary beating[END_REF], motile cilia play a crucial role in clearing mucus and debris from the airways under normal conditions, as can be seen in patients with abnormal airway ciliary beating, namely patients with primary ciliary dyskinesia. Because of its presence in third ventricle of the brain cavities, motile cilia also play a role in circulating cerebro-spinal fluid (CSF), where abnormal ciliary beating has recently been linked to hydrocephalus and other developmental cerebral abnormalities. In addition, they are speculated to be a part of the central glucose sensing system [START_REF] Garcia | Hypothalamic ependymal-glial cells express the glucose transporter GLUT2, a protein involved in glucose sensing[END_REF]. In both case, the circulation and the clearing are performed by the coordination of the beating cilia producing a wave. The coordination involving several cells is still speculative although it is well accepted that cilia on a single cell beat together. Beating properties are usually addressed by the inverse of the duty cycle, namely the Ciliary Beating Frequency (CBF). This frequency is strongly affected by the molecular components of the cells and the submerging fluid, thus contributing to local homeostatic processes. By this way CBF changes may allow the indirect assessment of metabolic variation. It is clear that the more accurate the estimation of this change will be, the more subtle variation of CBF in response to changes of the medium will be measured.

Unfortunately, cilia are densely packed and close together on all naturally occurring ciliated surfaces. Because this density is observed in the three dimensional case, any imaging system will fail in tracking a single cilium. Instead, the recorded image will be a combination of overlapping cilia, eventually beating at slightly different frequencies. Indeed, because of the long duration experiments and the time constant of the CBF changes according to the molecular components variations, the recording conditions are not optimal and makes difficult the CBF estimation. Furthermore, because of the fluid flux generated by the cilia and the perfusion system, the sample being recorded through the image acquisition system is shifted from frame to frame. Lastly, permanent patterns such as the artifacts due to the acquisition system could be recorded during the experiment and hinder a correct image registration.

The listed limitations can be overcome by using adequate processings usually applied in sequence. Anyhow, this could explain with only few papers cope with the difficult estimation of this frequency while it is of great interest in the Biology field. This problem cannot be efficiently approached using image processing tools such as patterns detection and tracking because of the overlapping cilia. However, attempts based on optical flow [START_REF] Kim | An Automated Measurement of Ciliary Beating Frequency using a Combined Optical Flow and Peak Detection[END_REF] and correlation methods [START_REF] Yi | Directional disorder of ciliary metachronal waves using two-dimensional correlation map[END_REF] failed with our specific application and need high quality records or manual interventions. There is a need for a robust CBF estimator that uses all the prior available to describe advantageously the temporal and spatial properties of the cilia. Many approaches are indeed signal processing methods such as the popular Fourier analysis [START_REF] Yi | Correlation between ciliary beat frequency and metachronal wave disorder using image analysis method[END_REF] applied on temporal series. Along the same line, the underlying model is not fully exploited according to the cilia properties when considering classical Fourier analysis.

The present report proposes a methodological approach that will cope with all these limitations and lack of adequate models by first removing grid artifacts, as described in section III. In section IV, the compensation of the geometrical sample shift is explained and illustrated. In section V, a full derivation of our approach is provided, including theoretical development and comparison to Fourier based analysis and evaluated on illustrative examples. This work will be ended by a conclusion.

II. IMAGE ACQUISITION

A. Brain slices

All the protocols were carried out in accordance with French standard ethical guidelines for laboratory animals (Agreement N 75-178, 05/16/2000) and were approved by the IPMC care committee. Attention was paid to use only the number of animals requested to generate reproducible results. 12-80 day old male wild type C57BL6 mice were bred in the local animal facilities and maintained on a 12h dark/light cycle (7am/7pm) with food and water ad libitum. Mice were decapitated and brains were rapidly removed into ice-cold gassed medium (95% O2/5% CO2) containing (in mM): 125 NaCl, 2.5 KCl, 1 MgCl2, 0.4 CaCl2, 1.25 NaH2PO4, 26 NaHCO3 and 25 glucose. Coronal slices of the median part of the hypothalamus (250m), were cut with a HM650V vibratome (Microm, Walldorf, Germany) and placed in a holding chamber at 34C for 30 min. Slices were then transferred into another incubating chamber at room temperature (22-25C) continuously perfused at a rate of 1-2ml.min-1 with Phosphate Bicarbonate Buffer Saline (PBBS), pH 7.4 when bubbled with 95% O2/5% CO2, containing additional CaCl2 (final concentration: 2mM) and used for recordings. Solutions were applied to the slices through the perfusion system.

B. Acquisition Systems

Two types of acquisition system are typically available for CBF studies. In this context, the confocal microscope technology is mostly used for linescan analysis as described in section V. Xt image sequences were acquired using a TCS SP5 Laser Scanning Confocal Microscope (Leica Microsystems, Nanterre, France) [START_REF] Doyle | Use of confocal linescan to document ciliary beat frequency[END_REF]. The brain slice was placed on the stage of the DMI6000 inverted microscope and first, a bright field image was done (XY scan) with the 488 nm laser line using the transmission detector (image a), through a HCX APO 40X/0.80 Water immersion objective. Then a line was drawn perpendicular to the axis of the cilia and a scanning was done in Xt mode along this line with a sampling frequency of 1400 Hz . Synthetic images (512x1024 pixels corresponding to a duration of 0.73 sec) of the traces of the beating of the cilia were recorded each 10 sec during the experiment. ! Fig. 1. A line is drawn parallel to the surface of the ciliate cells where the cilia are moving. The light intensity is determined along this line and plotted against time to build a synthetic image.

The second system acquires bright field image sequences on an upright Axioskop microscope (CarlZeiss, Le Pecq, France) using an ApoChromat 63X/1 water immersion objective. 200 images were captured every minute at a frame rate of 200 frames/ second with a GigE Vision color CMOS camera (PicSight, Leutron Vision, Switzerland). To use the camera at its maximum speed, the image size was reduced to acquire a 200x200 pixels sub area of the image. The acquisition and automatic recording of all the successive 1sec sequences was done using StreamPix software (Norpix, Montreal, Canada). This second acquisition system was prefered to be used in routine analysis in place of the confocal linescan [START_REF] Lechtrek | High-speed digital imaging of ependymal cilia in the murine brain[END_REF] strategy more sensitive to the focus drifts and the movements of the sample. However, the image recorded through this second system will be also processed offline with a linescan analysis similar to the one used with confocal microscope, for comparison.

III. ARTIFACTS REMOVAL

Grid pattern in the CCD sensor is one artifact that should be removed prior to any image processing. The image frame from the second acquisition system displayed in fig. 3 and the averaged image AV from a sequence of 200 images displayed in fig. 4 clearly exhibit such frame in addition to the cilia. Note that while the cilia don't appear in the average image because of the moving pattern, the grid is still visible because of its permanent location. We consider now the acquisition of a background sequence of 200 images (same illumination, without brain slice) and its computed average AVG displayed in fig. 5 10µm Fig. 3. A single frame containing ciliary beating and interferences It is assumed now that the grid AVG is added to each frame, and consequently to the average image AV, up to a constant multiplicative factor. This factor remains to be estimated in 10µm Fig. 4. Averaged frames AV containing ciliary beating and interferences order to efficiently subtract the grid. A single line of the frames is considered to compute this estimation. Let's call av(n) and avg(n) the corresponding lines of AV and AVG, whose means have been subtracted and where n stands for the column index. Lines that are free of cilia samples will produce better estimation. Considering the illustrative examples of fig. 4 and fig. 5, a line parallel to the border of the sample will be selected for the following computation.

The statistical properties of the line av(n) do not allow an efficient use of a least square method for the estimation. Indeed, with respect to the periodical spatial properties of the grid, it can be proven that the autocorrelation function of avg(n) exhibits a large spike for lag m = 0, i.e. column shift equal 0. Thus, for a guessed value α of the multiplicative coefficient that has to be estimated, the correlation function is computed such as:

r(m, α) = N -m-1 n=0 [av(n) -α.avg(n)] avg(n + m) (1)
where N stands for the length of the two lines. the criteria:

c(α) = |r(0, α) -r(1, α)| + |r(0, α) -r(-1, α)| |r(1, α)| + |r(-1, α)| (2)
with respect to the parameter α. This criteria tends to be high as long as the grid remains in the line of interest, and small when the grid is perfectly removed. This criteria is given in fig. 7 where a minimum is visible for the value α = 0.55.

The corresponding function r(m, 0.55) is displayed in fig. 8, where compared to fig. 6 the spike has been perfectly substracted, meaning that the introduced criteria is adequate and the coefficient accurately estimated.

Once the coefficient has been estimated, the weighted average image AVG containing the grid can be subtracted to the average image that contains the beating cilia. The figure fig. 9 shows that compared to the contaminated one (see fig. 4) the grid has been perfectly removed. This corrected average image will be used in the following as the reference for the registration process that will allow the image realignment. the optimal two-dimensional ( li, ĉo) alignment or registration is estimated by using the minimization:

( li, ĉo) = arg min li,col i j

(IM k (i, j) -RIM (i + li, j + col)) 2

(3) An exhaustive search is computed within a range of possible values using the full images. This criteria, similar to a crosscorrelation function, produces for each image IM k a pair of values that will allow the alignment of the frames sequence. In [START_REF] Cabasson | Time Delay Estimation: A New Insight Into the Woody's Method[END_REF], in addition to the optimal estimation of the reference image an iterative scheme is proposed for realignment in the one dimensional case that could be extended to this application. However, an exhaustive search guarantees to reach the optimal values and can be speeded up by reducing the search domain according to the alignment data estimated from the previous image. Indeed, because of the low strength exerted by the cilia and the slow perfusion flow, the spatial shift observed frame by frame is weak allowing an efficient relative realignment. An example of the registration performed of a sequence of 6000 frames is given in figure 10,11 where it can be noticed that the shift is weak and continuous along the acquisition. In fig. 12, it is clear that the spatial shift is almost diagonal in this experiment.

Note that optimality of the estimation is not addressed in this study. However, because it is based on mean square minimization, asymptotically optimal when the noise is gaussian, the estimation will be robust with respect to interferences. To this respect, the image alignment is indeed influenced by the static part of the sample and not by the cilia. It means that the registration will not track the cilia motion because it plays a very little role in the global image variation considered frame by frame. Furthermore, the complex nature of the cilia observation, including overlapping and multiplicity, does not allow a direct tracking frame by frame of a single cilium. In spite of these limitations, a global modeling of the sub image containing only cilia will be proposed as shown in the following V. FREQUENCY ESTIMATION Two classic measurements of the Cilia Beating Frequency (CBF) consist in a visual inspection of concatenated lines, named confocal linescan and the computation of the Fast Fourier Transform (FFT) of averaged pixels. These two methods will be briefly introduced, followed by the development of a new method named MLE.

A. Confocal Linescan

From the first system described in section II a synthetic image (fig. 2) is created where the periodic behavior of the cilia pattern should be visible. The peaks of the sinusoids were manually determined for each file (arrows in fig. 2) and the frequency was calculated from the average of 5-6 in between Fig. 12. The estimated pairs of spatial shifts (pixels). They are mostly superimposed because of the small variations estimated as integers peak measurements using an homemade ImageJ (Wayne Rasband, National Institute of Mental Health, Bethesda, Maryland, USA) macro. In addition to demanding acquisition conditions, this computation needs manual interventions leading to a large variance in the CBF estimation.

B. Average of pixel intensity

As the name implies, first the pixel intensities in the entire sub-image was averaged. This averaged value gave the average pixel intensity of the sub-image. This process was done for the entire set of images. The Fast Fourier Transform (FFT) of the average value was then obtained. We plotted the frequency vs. magnitude, the frequency with the maximum magnitude giving the CBF. This method was rather simple and not much computationally costly.

C. Average of the FFT of individual pixels

Here first the squared modulus of FFT of the individual pixels of the sub image were obtained. Once this process was completed, they were stored. This process repeats for all the pixels in all the sub images in the given set of images. Then, the average of the processed FFT is obtained and plotted. Thus a plot of frequency vs. magnitude is obtained. The frequency with the maximum magnitude is the CBF. This second method was much more computationally costly but quite effective. Note that in contrast to the average of pixels, the modulus of the FFT is invariant with respect to delay in the time domain. It is expected that due to its own motion, one cilia affects pixels in a bounded area with a small delay function of the velocity and the inter-pixel distance. This fact is in favor of any spectral approach based on modulus averaging.

Although the behavior model definition of single cilia may be well addressed [START_REF] O'callaghan | Analysis of ependymal ciliary beat pattern and beat frequency using high speed imaging: comparison with the photomultiplier and photodiode methods[END_REF], the estimation of the CBF based on this model is a difficult task due to the observation itself. This difficulty is explained not only by the presence of artifacts such as acquisition noise and spatial shift of the active ciliated cells but also because of the overlapping multiple cilia. This phenomenon makes tedious the frame by frame tracking of geometrical features and consequently hinders the accurate estimation of the CBF. In contrast to this spatial approach, robustness and accuracy are improved by considering temporal and spatial coherence. In other words, we consider that within a region of interest, all pixels share a common temporal feature, maybe shifted in time, that is the frame by frame intensity variation. Eventually, pixels could also capture adjacent cilia beating at the same frequency because of the metachronal waves. Assuming that the variation is periodical and that each pixel records an abrupt intensity variation, we propose a harmonic model that combines the fundamental frequency (f ) and the first and second harmonics. Considering that the model should be invariant with respect to time delays, each pixel intensity variation p i (t) from the total I within the subimage, is modeled by:

p i (t) = 3 j=1 [θ i,2j-1 cos(2πf jt)+θ i,2j sin(2πf jt)]+n i (t) (4)
Where n i (t) is a white gaussian error composed by the modeling error and observation noise. In this model, j and θ stand for the harmonic number and weight of the decomposition, respectively. Its is worth noticing that pairs of values (θ i,2j-1 , θ i,2j ) convey information of both the amplitude of the harmonic and its temporal shift although the analysis is focused on the frequency (f ). Among the available estimators provided by the estimation theory, the Maximum Likelihood Estimator (MLE) [START_REF] Kay | Fundamentals of StatisticalSignal Processing[END_REF] is selected for its properties with respect to the error. It is assumed here that the error is gaussian, stationary and uncorrelated spatially and temporally. Considering the I pixels from the sub-image, maximizing the likelihood function is equivalent to minimize the following criteria:

J 1 = I i=1 p i -M f θ i 2 (5)
Where p i and θ i stand for the vector formed by the temporal values of pixel i extracted from the successive frames and the parameters vector, respectively. Note that temporal delays introduced by the sequential recording of the motile cilia by adjacent pixels are inherent to the model. In this model, the parameter f is the fundamental frequency and only the first and second harmonics are taken into account. Then, M f matrix contains cosines and sines functions defined over the time t k (with t 1 t k t K ) corresponding to the acquisition of the frames such that:

M f = [M 1.f M 2.f M 3.f ] (6)
with

M j.f =    cos(2πf jt 1 ) sin(2πf jt 1 ) . . . . . . cos(2πf jt K ) sin(2πf jt K )   
Indeed, the criteria J 1 has to be minimized with respect to the set of parameters f and θ i 's. Developing [START_REF] Yi | Correlation between ciliary beat frequency and metachronal wave disorder using image analysis method[END_REF], the solution is thus given by the minimization: θ1 , . . . , θI , f = arg min θ1,...,θ I ,f

I i=1 p T i p i -2p T i M f θ i + θ T i M T f M f θ i (7)
The derivation of the criteria with respect to the vectors θ i provides the estimated θi 's:

θi = [ θi,1 , θi,2 , θi,3 , θi,4 , θi,5 , θi,6 ] T = (M T f M f ) -1 M T f p i (8)
It can be proven that after substituting these solutions in [START_REF] Doyle | Use of confocal linescan to document ciliary beat frequency[END_REF], the minimization is finally a maximization with respect to f :

f = arg max f I i=1 p T i M f (M T f M f ) -1 M T f p i (9)
Whatever the value of f , the matrix (M T f M f ) -1 can be approximated by kI, where I stand for the identity matrix and k a constant. Note that for specific values of f this approximation is indeed the equality. Thus ( 8) is simplified by: θi = kM T f p i [START_REF] Kay | Fundamentals of StatisticalSignal Processing[END_REF] and ( 9) can be replaced by:

f = arg max f I i=1 θT i θi (11)
In fact, a sum ( θ2 i,2j-1 + θ2 i,2j ) with j = 1, 2, 3 is equivalent to the squared modulus of the Fourier Transform (FT) at the frequency j.f , namely |F T (p i )| 2 (j.f ). Thus the maximization in ( 9) is equivalent to:

f = arg max f 3 j=1 I i=1 |F T (p i )| 2 (j.f ) (12)
The sum of the squared modulus of the fundamental frequency and its harmonics is also called the Compressed Spectrum CS(f ). Then, ( 12) is replaced by:

f = arg max f I i=1 CS i (f ) (13)
This last expression allows an easy interpretation of this maximization, using the mentioned approximation. It is finally equivalent to the problem addressed by the average of squared modulus of FFT where the harmonic modeling is accounted in the estimation. In summary, this last solution uses all the priors related to the observation of the motile cilia. Note that the optimization ( 12) should be performed over the single parameter f . Sophisticated optimization algorithms are available but with the risk of converging toward a local maxima. We will prefer in the sequel an exhaustive search over a given range of the f values, allowing the addition of a priori knowledge and fixing the frequency resolution.

D. Methods evaluation

Albeit the linescan technique has been designed to work on image captured by confocal microscope it has been tested on the second acquisition system because it will be preferred for routine acquisition.

In fig. 13 is given the linescan image from this system where the extraction of the periodic pattern seems to be feasible (see white arrows), however confounding. Nevertheless, such method fails in providing accurate CBF estimates with such image. Instead, the following methods will be preferred for their good performances with such image acquisition modality. Synthetic image built with the same method than for the linescan, applied to bright field images. The light intensity is determined along a line and plotted against time In order to focus the analysis on a small set of cilia's, it is defined a 3x17 pixels sub image in the vicinity of the cilia base (black rectangles in fig. 14). In order to analyze the time-varying CBF, sequences of 200 frames are processed every 10 seconds. So, for a total amount of 6000 frames, 30 estimated CBF are provided by the mentioned methods for two sub images (in fig. 14, left and right rectangles). In fig. 17 and fig. 18, the results from the averaged FFT, FFT of averaged pixels, MLE are given for the two sub images, respectively. It can be noticed that results are almost similar whatever the sub images and the method. However, the MLE methods slightly outperforms the others because the frequency variations tends to be smoother, as what is expected.

VI. CONCLUSION

Motivated by the great interest in the Biology field for CBF measurement, an estimator has been proposed according to the priors available in the time and spatial domain. Robustness has been addressed by introducing a harmonic model of the temporal variations of pixels intensities. It is worth noting that amplitudes and frequencies have been assumed constant in the sequences of 200 frames. A more general extended model could be introduced to overcome this limitation but with the counterpart of a lower robustness and a higher computational cost.

It has been shown that the resulting method has a clear interpretation in term of time-invariant estimation. In parallel with an image artifacts removal, the method outputs exhibit lower variability as expected.

Fig. 2 .

 2 Fig. 2.Synthetic Xt image obtained by scanning along a line perpendicular to the axis of the cilia. The CBF is determined by the average of the period calculated on 5 peaks of the sinusoids obtained in the plot of light intensity against time.

10µmFig. 5 .

 5 Fig. 5. Average of frames (from a blank field) AVG containing only interferences

Fig. 6 .

 6 Fig. 6. The correlation function r(m, α) for the value α = 0

Fig. 7 .Fig. 8 .Fig. 9 .

 789 Fig. 7. The criteria c(α) with a clear minimum

Fig. 10 .Fig. 11 .

 1011 Fig. 10.The estimated line shift (pixels) for each frame (k = 1 . . . 6000)

  100ms

Fig. 13 .

 13 Fig. 13.

  10µm

Fig. 14 .

 14 Fig. 14. Left and right sub images bounded by black rectangles defining regions of interest

Fig. 15 .Fig. 16 .

 1516 Fig. 15. Average of squared FFT (solid line), squared FFT of averaged pixels (dotted line)

Fig. 17 .Fig. 18 .

 1718 Fig. 17. CBF within the right sub image: averaged FFT (dotted line), FFT of averaged pixels (dashed line), MLE (solid line)

ACKNOWLEDGEMENT

Authors would like to thank N. Ramkumar, G. Conductier and J.L. Nahon for their contributions to this work. This work has been partly funded by grant PEPS 2010 and supported by INSB.


Daisuke Ishii

Guillaume Melquiond
email: guillaume.melquiond@inria.fr

Shin Nakajima

Inductive Verification of Hybrid Automata with Strongest Postcondition Calculus

Safety verification of hybrid systems is a key technique in developing embedded systems that have a strong coupling with the physical environment. We propose an automated logical analytic method for verifying a class of hybrid automata. The problems are more general than those solved by the existing model checkers: our method can verify models with symbolic parameters and nonlinear equations as well. First, we encode the execution trace of a hybrid automaton as an imperative program. Its safety property is then translated into proof obligations by strongest postcondition calculus. Finally, these logic formulas are discharged by state-of-the-art arithmetic solvers (e.g., Mathematica). Our proposed algorithm efficiently performs inductive reasoning by unrolling the execution for some steps and generating loop invariants from verification failures. Our experimental results along with examples taken from the literature show that the proposed approach is feasible.

Introduction

Hybrid systems, transition systems with continuous dynamics, are a good model for embedded systems that have a strong coupling with the physical environment. Achieving the desired reliability levels of such systems has brought a challenging and important problem in formal methods research.

To date, verification of hybrid systems has been extensively studied with two prominent approaches: model checking and logical analysis. The model-checking approach has been successfully applied to practical examples with tools such as HyTech [START_REF] Henzinger | HyTech: A model checker for hybrid systems[END_REF], PHAVer [START_REF] Frehse | PHAVer: Algorithmic verification of hybrid systems past HyTech[END_REF], and HybridSAL [START_REF] Tiwari | HybridSAL relational abstracter[END_REF]. The approach is said algorithmic: tools numerically over-approximate a certain class of hybrid automata (HA) to have piecewise-linear systems, and apply model-checking methods [START_REF] Alur | Discrete abstractions of hybrid systems[END_REF]. The second approach is based on logical analysis [START_REF] Platzer | Logical Analysis of Hybrid Systems[END_REF]. While the theory of logical analysis has been studied extensively, there are few practical tools. A notable and successful exception is KeYmaera [START_REF] Platzer | KeYmaera: A hybrid theorem prover for hybrid systems[END_REF]. The logical analytic approach can be applied to the class of hybrid programs which generalize the automata handled by model checking. Indeed, this class includes systems with symbolic parameters and nonlinear dynamics. There is, however, a major drawback: the larger the class of systems is, the less automatic its verification becomes. Engineers thus have to apply some proof strategies during the interactive verification process, which requires understanding the target model.

In this paper, we propose a partly automated tool for the logical analysis of HA that makes heavy use of state-of-the-art arithmetic solvers. Our goal is to prove safety properties. First, our method encodes executions of HA into straight-line imperative programs. This formalism allows us to construct a lassoshaped structure based on induction: after exhibiting at most m steps of continuous evolution and discrete transition, any execution of the system forms a loop with a length of at most n steps between some specific regions of the state space. Then, the imperative program is transformed into a conjunction of verification conditions as a result of strongest postcondition (SP) calculus. The resulting logic formula involves real-arithmetic predicates and ordinary differential equations (ODEs). The generated conditions can be discharged using solvers such as Mathematica for some nonlinear HA.

The contribution of this work is as follows. The use of an imperative language and SP calculus gives a straightforward justification of the soundness of our method for generating the finite-length verification conditions from HA. The algorithm we propose realizes an automated verification process, although some user interactions are needed to determine efficiently (a) correct numbers m and n of steps to unroll the execution and (b) the loop invariant that represents the initial region of the loop. Computer algebra techniques, however, are employed to automate most of the work in generating loop invariants.

This paper is organized as follows. Section 2 introduces the class of hybrid automata. Section 3 describes a simple imperative language for simulating HA and the corresponding SP calculus. In Section 4, we present the concept of induction and loop unrolling, and describe an algorithm for automated verification. Section 5 describes an implementation using Mathematica. Section 6 reports how our implementation behaves on several examples and provides a comparison of the results with existing tools. Section 7 describes some related studies.

Hybrid Automata

In this paper, we model hybrid systems as hybrid automata (HA) [START_REF] Henzinger | Verification of Digital and Hybrid Systems[END_REF]. Definition 1. A hybrid automaton is a tuple HA = L, V , Init, G, R, F, I that consists of the following components:

-A finite set L = {l 1 , . . . , l p } of locations. -A finite set V = {x 1 , . . . , x q } of real-valued variables. R V is the set of all of
the valuations of the system. -An initial condition Init in L×R V that specifies the initial states. In the first rule, dφ dt = F l is an abbreviation of dφ( t) dt = F l (φ( t)). We say that an execution is length-k canonical when of the form σ 0

-A family G = {G l,l ′ } l∈L,l ′ ∈L of guard conditions G l,l ′ in R V . -A family R = {R l,l ′ } l∈L,l ′ ∈L of reset functions R l,l ′ : R V → R V . -A family F = {F l } l∈L of vector fields F l : R V → R V . -A family I = {I l } l∈L of location invariants I l in R V . t > 0 φ(0) = ν ∀ t ∈ [0, t] dφ dt = F l ∧ I l [φ( t)] l, ν t -→ l, φ(t) G l 1 ,l 2 [ν1] ν2 = R l 1 ,l 2 (ν1) I l 2 [ν2] l1, ν1 0 -→ l2, ν2
t1 -→ σ 1 0 -→ σ 2 t3 -→ • • • 0 -→ σ 2k
that alternates continuous and discrete phases.

In this paper, we assume that multiple discrete transitions do not occur in an instant. We also assume that no discrete transition occurs initially. Thus, any execution can be expressed as a canonical execution. Verification of noncanonical executions can be considered as future work. Infinite-length canonical executions are supported; yet in presence of Zeno points (infinite number of transitions in finite time), HA executions are handled only up to the first point.

Example 1. Water-level monitor (WLM) [START_REF] Alur | The algorithmic analysis of hybrid systems[END_REF][START_REF] Manna | Deductive verification of hybrid systems using STeP[END_REF]. A controlled water tank is modeled as a four-location constant-rate HA, as illustrated in Figure 2. It supplies water at a constant rate rate out , whereby, in location off (and sw-on), the water level y decreases as dφy dt = rate out . In location on (and sw-off), the system pumps water to refill the tank, which results in the water level changing as dφy dt = rate in . A sensor observes y and switches between the locations on and off when the level reaches the thresholds low or high. However, it takes delay seconds for switching, hence the locations sw-on and sw-off. In this paper, we constrain the values for the constant parameters as follows:

min ≤ low ∧ high ≤ max ∧ low < high ∧ delay > 0 ∧ max ≥ high + rate in • delay ∧ min ≤ low + rate out • delay. (1)
Because the discrete transition edges in the automaton form a single cycle, the trace of locations that were reached is the same for all of the executions.

Definition 2. A safety property (or an inductive invariance) is expressed by a formula P , where P is a predicate on L×R V . HA |= P denotes that HA satisfies P , that is, predicate P holds initially and is preserved by every discrete transition and continuous evolution.

Example 2. In the following sections, we will prove that the level stays between a lower and an upper limit, which is expressed by the following safety property: 3 Modeling HA Executions with Programs

(min ≤ y ≤ max ). on ẏ = rate in y ≤ high y=low sw-off ẋ = 1, ẏ = rate in 0 ≤ x ≤ delay off ẏ = rate out y ≥ low sw-on ẋ = 1, ẏ = rate out 0 ≤ x ≤ delay y = high, x ′ = 0 x = delay y = low , x ′ = 0 x = delay
In this section, we introduce the theoretical foundation of our method. It analyzes finite and infinite executions of a HA by reusing traditional tools in program verification. We first introduce a simple imperative language in which the statements simply sketch the executions of the HA (Section 3.1). Then, we provide a notion of strongest postcondition for each program statement given a precondition, and we prove that this calculus derives the safety of the HA (Section 3.2).

Imperative Language

Given a HA, we define an untyped imperative language Imp HA . This language is basic, since it does not even provide loops. For the purpose of this work, it has only sequences and the commands evolve and trans. The command evolve expresses a continuous evolution of the HA for a given duration, while trans expresses a discrete transition.

Definition 3. The language Imp HA is given by the following syntax:

s ::= skip | s; s | evolve t | trans Definition 4. A program state (denoted S or S i
) is a map from variable names to program values. A special variable x s is associated to the "current" state (∈ L × R V ) of the HA execution. For the sake of readability, pseudo-variables are introduced to access part of the HA state as follows:

x s = x l , • = •, x v =
•, (x 1 , . . . , x i , . . . , x q ) . We assume this equivalence is always maintained automatically when a new value is assigned to a pseudo-variable.

Figure 3 describes the operational semantics of the language. [[e]]

S denotes the term obtained by replacing each free variable of an expression e by its associated value in the program state S. S{x → v} denotes the program state obtained by adding to S that variable x is associated to the value v. The rules for skip and sequence are the usual ones. The rules for evolve and trans are derived from the operational semantics of a HA execution. Note that we allow the statement evolve 0, so that the theorems presented in this paper have a simple way to check the safety property for the initial state or after a discrete transition. Lemma 1. For any execution σ 0

t1 -→ σ 1 0 -→ • • • t k -→ σ 2k-1
of the HA, assuming that σ = σ 0 holds for the initial program state, there is an execution of the following Imp HA program that does not block (that is, it reduces to skip) and such that the final program state satisfies

x s = σ 2k-1 . evolve t 1 ; trans; • • • ; evolve t k
Note that this program might also have either blocking executions or executions that end on a different HA state; the former are made irrelevant by our SP-based approach, while the latter are expected due to the non-deterministic nature of HA. For the programs above, the execution is canonical only for the first k -1 continuous steps; the last duration t k can be arbitrarily short. It can also be arbitrarily large, if the HA stays infinitely long in that continuous evolution.

Since we can now express any partial execution of a HA as a program, we can state the safety property of the HA as a property that every non-blocking program must satisfy in its final state. Lemma 2. If, for all non-blocking programs of Imp HA of the above form starting from an initial program state σ ∈ Init, property P holds in the final program state, then P is a safety property for the HA (up to the first Zeno point, if any).

Strongest Postconditions

In this section, we instantiate the principles of program verification [START_REF] Hoare | An axiomatic basis for computer programming[END_REF][START_REF] Dijkstra | Guarded commands, nondeterminacy and formal derivation of programs[END_REF] with Imp HA . We are not interested in manual verification, so we will skip over the definition of Hoare triples and directly go to the topic of verification conditions (VCs). Moreover, since we are not dealing with reachability but only safety, we do not have to prove that programs are non-blocking, we can just assume they are. Therefore, weakest preconditions (WPs) and strongest postconditions (SPs) are dual from each other for our purpose. Should we have to perform backward reachability analysis, WP computation would be better suited. This is not the case though, so we choose SP, so as to follow the direction of time.

Lemma 3 (Soundness of SP). For any program s in Imp HA , if the initial state satisfies a given property P , the final state satisfies SP (P, s) (assuming s terminates) with SP inductively defined as follows. 3 SP (P, skip) := P SP (P, s 1 ; s 2 ) := SP (SP (P, s 1 ), s 2 )

SP (P, evolve t) := ∃φ P [x v ← φ(0)] ∧ φ(t) = x v ∧ (∀ t ∈ [0, t] dφ dt = F x l ∧ I x l [φ( t)]) SP (P, trans) := ∃ l ′ , x ′ v P [x s ← l ′ , x ′ v ] ∧ G l ′ ,x l [x ′ v ] ∧ x v = R l ′ ,x l (x ′ v ) ∧ I x l [x v ] Proof.
SP ((x l = on ∧ y = low ), evolve t) ⇒ y ≤ max .
Let us assume that we are in program state such the left-hand side holds, and we prove that y ≤ max holds. From the definition of SP , we know that there exists a function φ such that

(x l = on ∧ φ y (0) = low ) ∧ φ(t) = (x, y) ∧ (∀ t ∈ [0, t] dφ dt = F x l ∧ I x l [φ( t)])
As a consequence, we have y = low +rate in •t (by solving the ODE) and y ≤ high (by unfolding the location invariant I x l ). The latter property, in conjunction with Constraint (1) of Example 1, proves the goal y ≤ max by linear arithmetic.

Remark 1. As we will later pass the verification conditions to automated tools, it is important to eliminate as many quantifiers as possible beforehand. For instance, SP (P, trans) has the form

∃l ′ Q[l ′ ]. This is equivalent to the disjunction Q[l 1 ] ∨ . . . ∨ Q[l p ]
with l 1 , . . . , l p all the locations. In the case of SP (P, evolve t), Example 3 shows how one can get rid of ∃φ if the ODE admits a closed form.

4 Inductive Verification Method

Induction Strategy

We now present an algorithm derived from Lemma 2 that performs safety verification of a HA. The statement of Lemma 2 is unpractical, as it requires verifying infinitely-many programs. This section describes how we can build weaker yet more practical variants of it, by only considering a bounded number of programs. The approach is as follows. Let us assume that there is a predicate P + such that P + ⇒ P and from an initial state, any execution of HA reaches a state satisfying P + after alternating at most m continuous evolutions and m discrete transitions, -from any state satisfying P + , any execution of HA reaches a state satisfying P + after alternating at most n continuous evolutions and n discrete transitions.

Verifying the safety property is therefore simple:

-For the initial m-step execution, we check that every intermediate state is safe and that the execution finally reaches the region represented by predicate P + (base case). -For the n-step execution from the region P + , we check that every intermediate state is safe and that the execution finally reaches the region P + .

The success of our approach depends on whether we can exhibit some lengths m and n and some predicate P + for a given HA. We first show the simplest case (m = 0 and n = 1): the base case is the verification of the initial states, and the induction is performed on a continuous phase followed by a discrete phase.

Theorem 1 (Simplest case). Given a predicate P + such that P + ⇒ P holds in any state, the following inference rule is correct: VC 0 : Init ⇒ P + VC 1 : ∀t ≥ 0 SP (P + , evolve t) ⇒ P VC -1 : ∀t ≥ 0 SP (P + , evolve t; trans) ⇒ P + HA |= P Proof. VC 0 checks that the initial states satisfy the property P + . VC -1 inductively verifies that all of the possible two consecutive continuous and discrete phases σ i ti+1 --→ σ i+1 0 -→ σ i+2 evolve for the arbitrary duration t i+1 from a state σ i that satisfies P + to a state σ i+2 that again satisfies P + . VC 1 ensures that the safety property was not broken during the continuous phase.

⊓ ⊔

We now extend the above induction to a more generic case.

Theorem 2 (Unrolled case).

SP 1 ≡ SP (Init ∧ ¬P + , evolve t 1 ) VC 1 : ∀t 1 ≥ 0 SP 1 ⇒ P SP 2 ≡ SP (SP (SP 1 , trans) ∧ ¬P + , evolve t 2 ) VC 2 : ∀t 1 , t 2 ≥ 0 SP 2 ⇒ P . . . SP m ≡ SP (SP (SP m-1 , trans) ∧ ¬P + , evolve t m ) VC m : ∀t 1 . . . t m ≥ 0 SP m ⇒ P SP 0 ≡ SP (SP m , trans) VC 0 : ∀t 1 . . . t m ≥ 0 SP 0 ⇒ P + SP m+1 ≡ SP (P + , evolve t 1 ) VC m+1 : ∀t 1 ≥ 0 SP m+1 ⇒ P . . . SP m+n ≡ SP (SP (SP m+n-1 , trans) ∧ ¬P + , evolve t n ) VC m+n : ∀t 1 . . . t n ≥ 0 SP m+n ⇒ P SP -1 ≡ SP (SP m+n , trans) VC -1 : ∀t 1 . . . t n ≥ 0 SP -1 ⇒ P + HA |= P
Proof. This theorem is an extension of Theorem 1. It verifies that a state satisfying P + can be reached in at most m steps initially (from VC 1 to VC 0 ), and then inductively that P + can always be reached again in at most n steps (from VC m+1 to VC -1 ).

⊓ ⊔ Remark 2. Only VC 0 and VC -1 check that P + holds after an execution; all the other VCs check the safety property P only. Moreover, except for VC m+1 , all these other conditions compute the SP by assuming that P + does not hold. Indeed, there might be less than n transitions before reaching again a state satisfying P + (or m transitions initially).

Verification Algorithm

Given a HA, a safety property P , and the maximal numbers m max and n max of steps to unroll, the algorithm in Figure 4 tries to check that all the hypotheses of Theorem 2 hold, and thus that HA |= P holds too. 4 The algorithm performs the inductive verification with every m ≤ m max and n ≤ n max (line 1). It iteratively computes the base case (line 4) and then the induction step (line 7). Procedure Validate returns true if the given logic formula holds, false if it cannot conclude. The verification succeeds if all the verification conditions are successfully validated (line 10).

When the verification fails during the induction step, we strengthen the loop invariant (line 8) so that the failing condition holds, and we perform the verification anew. Possibly, procedure Learn strengthened the invariant so much that we detect it is now useless (line 3). In this case, we leave from the inner recursion and try the verification with another m and n, or otherwise return false.

Note that the algorithm does not specify how to enumerate m and n. Typically, we enumerate from m = 0 and n = 1, but for certain models, we can guess the values, e.g., from the size of a lasso-shaped automaton. In the algorithm, the verification of the base case (lines 4-6, named BaseCase) and the induction step (lines 7-11, named Induction) are independent, thus we can also reverse the order of the two verification processes.

Loop Invariant Generation

In the following, we present the loop invariant generation method implemented in procedure Learn. Let us assume that the verification of a condition VC i ≡ ∀t 1 . . . t i ≥ 0 SP (P + , s) ⇒ P has failed in the induction step. Then, Learn(VC i ) generates a lemma Q from the failed verification. Specifically, Learn generates a formula Q such that VC i becomes valid after we update the loop invariant as P + := P + ∧ Q. Basically, Learn searches for Q such that SP (Q, s) ⇒ VC i holds by applying algebraic transformations to VC i . Note that all the occurrences of variable x s (the current state) in Q refer to the time P + holds, while the ones in VC i refer to the state at the end of the execution of s. end while 13: end for 14: return false Learn computes Q by using a quantifier elimination (QE) method, such as the Resolve procedure of Mathematica:

Q := QE(∀x s ∀t 1 . . . t i (SP ((P + ∧ x 0 = x s ), s) ⇒ P ))[x 0 ← x s ].
To simplify the loop invariant, the other local variables in VC i , i.e., φ, t, l ′ , x ′ v introduced in the SP calculus in Lemma 3 and t i introduced in Theorem 2, should also be removed. Unfortunately, QE with mixed quantifiers and function quantifiers is a hard problem in general. See Remark 1 and the next section for details on how we perform this simplification.

The formula computed for Q is often a large disjunctive formula that is unusable as a loop invariant. For instance, some sub-formulas of Q describe states that are never accepted by the HA. Such sub-formulas are not only useless but make the verification process expensive. So we strengthen Q according to the following strategies:

-Lemma separation. We split Q at the (top-most) disjunction operators and employ one (or several) of the resulting sub-formulas. -Location disabling. When we remove a sub-formula of Q that is related to some location l, we insert the constraint x l = l. The resulting loop invariant might be effective when combined with loop unrolling.

Implementation

We have implemented the method presented in the previous sections using Mathematica 8.0.45 , which can perform the computations in a fully symbolic manner.

Note that the loop invariant generation by Learn (line 8) is not automatic but guided by the user so as to apply the strategies described in Section 4.3. Validate is implemented in three different ways by using the built-in procedures of Mathematica, FullSimplify, Reduce, and FindInstance. We also rely on Mathematica's DSolve to find closed form of ODEs whenever possible.

In the implementation of BaseCase and Induction, we optimize the computation in two ways. First, we do not validate each VC i separately but try to reuse the common assumptions. When validating VC i , the algorithm computes SP i which axiomatizes the state after executing the corresponding program s i , and then validates SP i ⇒ P/P + . If we perform the validation of VCs in ascending order, we can compute SP i from SP i-1 efficiently. Second, we perform locationwise validation of VCs to avoid the inefficiency that occurs when the execution of program s spans multiple locations. So we replicate the SP and instantiate each copy with a different location (cf. Remark 1). Throughout the computation, we manage the set of the copies instead of the original SP. Although it causes Validate to be called more often, the computation is more efficient in general.

Example 4. We verify the safety property of Example 2 for the HA in Example 1 with this implementation. Following the main algorithm, we first compute with m = 0 and n = 1. We run BaseCase to check that Init entails P + ≡ P , and it returns true. Next, we simulate a continuous and discrete change by running Induction. It computes the SP separately for each of the locations, on, sw-off, off , and sw-on, and validates VCs. For VC 1 , the validation for locations on and off succeeds but the validation for sw-off and sw-on fails. Procedure Learn generates the following lemmas for these two locations.

Q sw-off ≡ min + x • rate in ≤ y + delay • rate in ≤ max + x • rate in ∨ x = delay ∨ y + delay • rate in < low + x • rate in , Q sw-on ≡ min + x • rate out ≤ y + delay • rate out ≤ max + x • rate out ∨ x = delay ∨ high + x • rate out < y + delay • rate out .
Here, we can use either of the two presented strategies for improving the loop invariant. For instance, location disabling appends

Q 1 := x l = sw-off ∧ x l = sw-on
to P + . The VCs are then successfully validated with m = 0 and n = 2. The lemma-separation strategy makes use of the additional lemmas generated by Learn. Here, we divide each lemma into three parts at the top-most disjunction operator. Then, the first part of each lemma (denoted Q sw-off,1 and Q sw-on,1 ) makes the verification successful. More precisely, if we append

Q 2 := (x l = sw-off ⇒ Q sw-off,1 ) ∧ (x l = sw-on ⇒ Q sw-on,1 )
to P + , the validation succeeds with m = 0 and n = 1.

Experiments

To confirm the feasibility of our method and to compare it with existing tools, we applied it to several examples taken from the literature. We also verified the examples using the existing tools, HyTech, PHAVer, and KeYmaera, for comparison. The encoded models for the implementation is available at http:// www.ueda.info.waseda.ac.jp/ ~ishii/pub/mathybrid/. Table 1 reports the results of verifying the examples using our implementations. The columns are: the number of locations; the number of variables; the way loops are unrolled (i.e., m/n); how many times P + had to be improved by the main algorithm; the computational time taken by the BaseCase and Induction procedures implemented in Mathematica; the time taken by HyTech (version 1.04f, indicated by "H") or PHAVer (version 0.38, indicated by "P"); and the time taken by KeYmaera (version 3.0). The notation "-" means that the verification failed. The experiments were run on a 3.4GHz Intel Xeon processor with 4GB of RAM. Note that the computational time for our method only measures the process after we found the loop invariants, since their generation requires some human interaction.

Considered Examples

WLM. Example 1 could be verified with our proposed method in a reasonable time, as explained in Example 4. In [START_REF] Manna | Deductive verification of hybrid systems using STeP[END_REF], the same instance was handled by using a mathematical solver manually, whereas our Mathematica implementation verified the instance by simply following the algorithm. The model-checking (MC) tools could not handle this instance because of the nonlinear terms caused by the parameterized flow rate. KeYmaera verified this example but the model had to be given a loop invariant beforehand [START_REF] Platzer | Guide for KeYmaera hybrid systems verification tool[END_REF].

Leaking gas burner (LGB) [START_REF] Alur | The algorithmic analysis of hybrid systems[END_REF]. Our implementation verified this rectangular HA consisting of two locations L = {leaking, non-leaking} as follows: Induction failed in the verification of the first continuous evolution in the two locations. The lemma generated for leaking was successful. For non-leaking though, we had to resort to our location-disabling strategy. Then, the verification succeeded with m = 4 and n = 2. This model was verified efficiently by the MC tools. KeYmaera could not verify the model, even with the loop invariant.

Temperature control [START_REF] Alur | The algorithmic analysis of hybrid systems[END_REF]. Our implementation verified this problem after some preliminary transformations. First, we verified that location shutdown of the HA is never reached. In order to get a loop invariant, we strengthened the safety property by appending the negation of the guard condition of the transition edge to shutdown. The failure of Induction led to a lemma of the form Q 1 ∨ Q 2 ∨ Q 3 , but setting each sub-lemma as a loop invariant did not make the verification successful. After some trials, we found that the lemma Q 1 ∧ (Q 2 ∨ Q 3 ) was a necessary loop invariant. Finally, the verification succeeded for m = 1 and n = 1. This model was also verified efficiently by the MC tools. KeYmaera could not verify the model, even with the loop invariant we had found.

Bouncing ball. This simple nonlinear HA describes a ball with a constant acceleration. As exemplified in [START_REF] Platzer | Logical Analysis of Hybrid Systems[END_REF], we verified that the height of the ball never exceeds the initial energy level of the ball, assuming that the reflection coefficient is smaller than 1. We first attempted the verification under a simple constraint that specified only the sign of each parameter and generated a lemma equivalent to the energy consumption constraint in [START_REF] Platzer | Logical Analysis of Hybrid Systems[END_REF]. We succeeded by setting this lemma as the initial condition and the loop invariant. KeYmaera verified the model given the energy consumption constraint as the initial condition.

European train control system (ETCS) [START_REF] Platzer | Logical Analysis of Hybrid Systems[END_REF][START_REF] Hasuo | Exercises in nonstandard static analysis of hybrid systems[END_REF][START_REF] Abrial | Formalizing hybrid systems with Event-B[END_REF]. The simple model borrowed from [START_REF] Platzer | Logical Analysis of Hybrid Systems[END_REF] is about a train at a position z that should not exceed a limit m. The original model does not have guard conditions so we set them manually based on the analysis in [START_REF] Abrial | Formalizing hybrid systems with Event-B[END_REF]. We attempted to verify the safety property z < m by running the algorithm with m = 0 and n = 1. Verification succeeded after we obtained a loop invariant from the failure in the validation of VC 1 . This model was also verified in [START_REF] Platzer | Logical Analysis of Hybrid Systems[END_REF][START_REF] Hasuo | Exercises in nonstandard static analysis of hybrid systems[END_REF] by using several strategies for the model transformation and loop invariant generation. MC tools could not verify the model because of the nonlinear constraints. KeYmaera verified the model by setting a specific parameter constraint as described in [START_REF] Platzer | Logical Analysis of Hybrid Systems[END_REF].

Highway [START_REF] Jha | Reachability for linear hybrid automata using iterative relaxation abstraction[END_REF]. This model concerns an autonomous highway with n vehicles. We solved instances for n = 9 and n = 19, which were also computed by the specific method in [START_REF] Jha | Reachability for linear hybrid automata using iterative relaxation abstraction[END_REF]. PHAVer verified the instance of n = 9 but the computation for n = 19 failed after consuming the available memory. KeYmaera could not verify this example.

Discussions

The MC tools verified three examples quite efficiently. However, our method was better for the other examples. First, it can handle uncertain parameters. Example 1 involves such parameters, as described in Equation (1). In the bouncing ball example, the initial height, velocity, and reflection coefficient are parameterized. Although HyTech and PHAVer verify the same problems with constant values given to the parameters, they cannot verify the instances that involve uncertain parameters. Second, our method scales better: for the highway example, PHAVer can handle only the instances up to n = 15 [START_REF] Jha | Reachability for linear hybrid automata using iterative relaxation abstraction[END_REF].

Although KeYmaera handles various hybrid programs automatically, it did not succeed on most hybrid programs that were translated from hybrid au-tomata. Users often need to annotate models with a loop invariant that might be difficult to extract from the original problem [START_REF] Platzer | Guide for KeYmaera hybrid systems verification tool[END_REF]. Otherwise, users need to interact with the underlying theorem prover to investigate the correct derivation tree with various deduction rules. Our approach is limited in verification strategies, i.e., induction and loop unrolling, but the results show that the approach is effective for various examples in practice.

Although our method requires that the executions are lasso shaped (from the point of view of the loop invariant P + ), many examples in the literature can be handled. It, however, requires other verification strategies for the case of compositional and distributed hybrid automata.

Related Work

Various tools for the logical analysis of hybrid systems have been proposed. These methods translate hybrid systems into an underlying verification framework, such as STeP [START_REF] Manna | Deductive verification of hybrid systems using STeP[END_REF], PVS [START_REF] Ábrahám-Mumm | Verification of hybrid systems: Formalization and proof rules in PVS[END_REF], SAL [START_REF] Ghosh | Automated symbolic reachability analysis: with application to delta-notch signaling automata[END_REF], Fluctuat [START_REF] Bouissou | HybridFluctuat: A static analyzer of numerical programs within a continuous environment[END_REF], and Event-B [START_REF] Abrial | Formalizing hybrid systems with Event-B[END_REF][START_REF] Su | Complementary methodologies for developing hybrid systems with Event-B[END_REF]. However, neither the translation nor the verification is fully automated, because some invariants must be added manually, and the theorem provers require some interactions.

Another tool, KeYmaera [START_REF] Platzer | KeYmaera: A hybrid theorem prover for hybrid systems[END_REF][START_REF] Platzer | Logical Analysis of Hybrid Systems[END_REF], developed by Platzer et al., has been successful in recent years. This tool supports hybrid programs that are annotated using differential (algebraic) dynamic logic. A dedicated theorem prover verifies the programs by using a set of proof strategies [START_REF] Platzer | Logical Analysis of Hybrid Systems[END_REF]. With its imperative language, which is more expressive than HA, and its corresponding logic (which depends on 141 inference rules [START_REF] Platzer | KeYmaera: A hybrid theorem prover for hybrid systems[END_REF]), KeYmaera is able to perform various logical analysis through a variety of strategies, including induction, and can serve as a basis for a complete verification framework [START_REF] Platzer | Logical Analysis of Hybrid Systems[END_REF]. In contrast, our framework consists of a light imperative language that is sufficiently expressive to encode HA executions and a logical framework that is introduced to pursue automated verification with the induction strategy.

Recently, a logical analysis tool based on the framework of Hoare logic and relying on infinitesimal variables was proposed [START_REF] Hasuo | Exercises in nonstandard static analysis of hybrid systems[END_REF]. Although its verification scheme comes with several strategies and an invariant generation technique, its practical uses are still unclear.

There are techniques for hybrid systems that generate polynomial invariants by analyzing the executions of a HA via Gröbner basis manipulations [START_REF] Sankaranarayanan | Constructing invariants for hybrid systems[END_REF][START_REF] Rodriguez-Carbonell | Generating polynomial invariants for hybrid systems[END_REF]. These methods could be integrated in the Learn procedure of our framework.

The proposed method also relates to BMC methods. BMC of infinite executions based on induction has been proposed (e.g., [START_REF] De Moura | Bounded model checking and induction: From refutation to verification[END_REF]), but this approach is applied to discrete systems with continuous states. Most of the BMC tools for hybrid systems handle only finite executions. This is not the case for Hybrid SAL Relational Abstracter [START_REF] Tiwari | HybridSAL relational abstracter[END_REF]. This tool is a translator from hybrid systems to discrete systems with a specific abstraction method. Our method directly handles HA without the abstraction. This paper presents a tool for logical analysis of safety properties of HA, which is able to deal with a large class of linear and nonlinear HA, in contrast with the model-checking approach found in major existing tools.

Rather than introducing various derivation rules to automatically verify HA, we are using a simple process inspired from deductive program verification: strongest postcondition calculus. It allows us to compute logical formulas that, once proved, guarantee the safety of the HA. Our experiments show that our method succeeds in a reasonable time on some example HA from literature, including some that were not solvable with existing tools. The verification process amounts to finding loop invariants, as is the case for program verification. This search for sufficient invariants is guided by the responses from the decision procedures assisted by Mathematica.

A limitation of our approach is that the invariant generation process still requires some human interaction. Efficient automated search of invariant generations is the next challenge for us to tackle. Another direction for further research would be to explore the relation between our approach and some methods from model checking, e.g., verification of an over-approximated model [START_REF] Alur | Discrete abstractions of hybrid systems[END_REF].

Fig. 1 .-→ σ 1 t2-

 11 Fig. 1. Operational semantics of HA.

Fig. 2 .

 2 Fig. 2. Water-level monitor.

Fig. 3 .

 3 Fig. 3. Operational semantics of Imp HA .

3 P

 3 [x ← e] denotes the substitution of all the occurrences of variable x in P with e.

Fig. 4 .

 4 Fig. 4. Algorithm for inductive verification.

  Let us assume that there are S and S ′ such that [[P ]] S holds and S, s * S ′ , skip. We just have to prove that [[SP (P, s)]] S ′ holds. The proof is performed inductively on the structure of the statement s by checking that every case of SP is implied by the operational semantics of Imp HA . This is a consequence of the operational semantics of HA given on Figure1.

⊓ ⊔

Example 3. Let us prove that, if the HA of Figure

2

is in a state satisfying x l = on ∧ y = low , then any continuous evolution of duration t leads to a state satisfying y ≤ max . By Lemmas 2 and 3, it is sufficient to prove that the following implication holds in any program state:

  To fix this discrepancy, Input: HA; P ; m max ∈ N ≥0 ; n max ∈ N >0 Output: true: HA |= P ; false: cannot decide P within m max + n max steps 1: for m ∈ {0, . . . , m max }; n ∈ {1, • • • , n max } do

	2:	P + := P
	3:	while P + ≡ false do
	4:	if ¬∀i ∈ {0, . . . , m} Validate(VC i ) then
	5:	break
	6:	end if
	7:	if ∃j ∈ {m + 1, . . . , m + n, -1} ¬Validate(VC j ) then
	8:	P + := P + ∧ Learn(VC j )
	9:	else
	10:	return true
	11:	end if
	12:

Table 1 .

 1 Experimental results.

	example	locs vars unroll lemmas Mathematica MC tool KeYmaera
	WLM (Ex. 1) 4	2	0/1	2	0.85s	-	1.8s
	LGB	2	3	4/2	3	2.22s	0.004s (H)	-
	temp. control 4	3	1/1	4	2.82s	0.012s (H)	-
	bouncing ball 1	2	0/1	1	0.49s	-	0.9s
	ETCS	2	3	0/1	1	4.48s	-	3.1s
	highway 9	10 9	0/2	1	0.22s	0.22s (P)	-
	highway 19	20 19 0/2	1	3.64s	-	-

http://www.wolfram.com/mathematica/

Acknowledgments The authors are indebted to the anonymous referees for their helpful comments. This work was partially funded by JSPS (KAKENHI 23-3810).