Dr Arnaud Gagneur
email: arnaud.gagneur@usherbrooke.ca

Geneviève Héry-Arnaud

Séverine Croly-Labourdette

Gisèle Gremmo-Feger

Sophie Vallet

Jacques Sizun

Roland Quentin

Didier Tandé

Infected breast milk associated with late-onset and recurrent group B streptococcal infection in neonatal twins: a genetic analysis

Keywords: group B streptococcus, breast milk, recurrent infection, genetic analysis, late-onset disease

Asymptomatic excretion of group B streptococcus (GBS) in breast milk may be an under-recognised cause of neonatal and recurrent infection. We report the case of late-onset and recurrent infection in newborn twins resulting from ingestion of maternal breast milk infected with GBS. Genetic analysis of isolates is equally presented.

Introduction

Group B streptococcus (GBS) is the most frequent cause of neonatal sepsis and meningitis. Most cases occur in the first week of life and are related to vaginal carriage in the mother (early-onset disease). Conversely, late-onset disease (between 1 week and 3 months of age) is less common and is handtransmitted by nursery personnel or via other nosocomial or community pathways [START_REF] Klein | Bacterial sepsis and meningitis[END_REF]. Late-onset and recurrent disease have also been reported with the ingestion of infected mother's milk, with some cases confirmed using molecular techniques [2-4, 7, 9, 13]. We report in this article the case of lateonset and recurrent infection in newborn twins resulting from ingestion of maternal breast milk infected with GBS. In this case genetic analysis demonstrated that all GBS isolates (maternal breast milk and vaginal isolates; twin CSF and blood isolates) were identical, but additional genetic analysis also revealed that the GBS isolates were of a particularly virulent clone belonging to serotype III, as are 90% of the strains responsible for late-onset disease [START_REF] Baker | Group B streptococcal infections[END_REF].

Case report

Premature twins were delivered by spontaneous vaginal delivery at 31 weeks gestation, 48 hours following membrane rupture. The 26 year-old mother, gravida 1, para 2, received intrapartum antibiotic prophylaxis (single dose of penicillin G) due to positive vaginal GBS culture at 30 weeks gestation. Twin 1, female, weighed 1600g with an Apgar score of 10/10 at 1 and 5 minutes.

No respiratory failure was noted and early enteral feeding was started with raw breast milk at day 1. Total enteral alimentation with breast milk was obtained at day 6. Twin 2, male, weighed 1720g with an Apgar score of 4/9 at 1 and 5 minutes. He was intubated shortly after birth due to respiratory failure and received one dose of surfactant. By the 3 rd hour of life, he was extubated and nasal CPAP was initiated. Enteral alimentation with raw breast milk was introduced at day 4. Antibiotic therapy with cefotaxime and amoxicillin was prescribed for both infants due to incomplete prepartum antibiotic prophylaxis and stopped at day 2 due to negative C-reactive protein as well as negative gastric and blood cultures. Cerebral ultrasound examination was normal for both infants at day 4.

On day 13, Twin 2 developed cardio-respiratory instability and blood culture tested positive for group B streptococcus. Meningitis was suspected due to elevated CSF protein concentration. Until day 16 Twin 1 was asymptomatic with negative C-reactive protein control. On day 16, she developed respiratory distress and subsequent blood and cerebrospinal fluid cultures tested positive for GBS. Antibiotic treatment with amoxicillin at 200mg/day was prescribed for 14 days, in association with an aminoglycoside during the first 48 hours of treatment. Control blood cultures were negative after day 1, 3 and 5 of treatment. Cerebral ultrasound examination controls were normal for both infants. At day 41 of life Twin 1 developed septic syndrome with parotitis and was transferred to the NICU. Blood culture was positive for GBS. Cardiac ultrasound examination was normal. Antibiotic therapy with cefotaxime (200mg/kg/day) and tobramycin (5mg/kg/day) was initiated. At day 7 of infection tomodensitometry examination identified cerebral microabcess and modification of the antibiotherapy ensued, with the administration of ciprofloxacin (30 mg/kg/day) associated with cefotaxime (250 mg/kg/day) for 3 weeks. Then oral amoxicillin was initiated for 3 additional weeks. Mastitis was diagnosed in the infants' mother 24 hours following discovery of GBS infection in Twin 1 (day 42). Milk culture tested positive for GBS and the maternal infection was treated with amoxicillin for 10 days. Breastfeeding was suspended and a 10-day preventive oral amoxicillin treatment given to the non-infected twin (confirmed via negative blood and CSF cultures as well as negative CRP controls). Following this infection, the infants remained on pasteurized breast milk. Follow-up at one year showed no cerebral anomalies upon ultrasound examination in association with normal neurological examinations at 1 year of life.

Analysis revealed all strains as belonging to serotype III. Epidemiological relationships between maternal and neonatal GBS isolates were investigated by pulsed-field gel electrophoresis (PFGE) of DNA restricted with SmaI [11].

Analysis was conducted on maternal vaginal and raw breast milk isolates (2 isolates), a single Twin 2 blood culture isolate (1 isolate), and on Twin 1 CSF and first-and -second blood culture isolates (3 isolates). All six isolates displayed identical PFGE patterns, revealing their genetic relationship (figure 1A).

Characterization of isolate virulence was conducted by multiplex PCR according to primers and method previously described [11,12]. First, amplification of the tRNA gene clusters at the 3' end of rRNA operons produced a unique fragment of 1.2 Kb (figure 1B); second, hylB amplification produced a 0.3 Kb fragment, showing no IS1548 insertion within the gene (figure 1B). This pattern was correlated to the invasive phylogenetic division I defined by Musser et al. in multilocus enzyme electrophoresis (MLEE) analysis [START_REF] Musser | Identification of a highvirulence clone of type III Streptococcus agalactiae (group B Streptococcus) causing invasive neonatal disease[END_REF][START_REF] Quentin | Characterization of Streptococcus agalactiae strains by multilocus enzyme genotype and serotype: identification of multiple virulent clone families that cause invasive neonatal disease[END_REF].

Discussion

Cases reporting neonatal late and recurrent group B streptococcal disease associated with raw maternal milk are rare, and few are the studies in which genetic evidence is proposed for this scenario [START_REF] Atkins | Recurrent group B streptococcal disease in infants: Who should receive rifampin?[END_REF][START_REF] Kotiw | Late-onset and recurrent neonatal Group B streptococcal disease associated with breast-milk transmission[END_REF]13]. In our case, not only was total DNA macrorestriction analysis conducted, showing indistinguishable patterns for the six isolates, but additional genetic analysis also revealed that the GBS isolates belonged to a particularly virulent clone shown to produce more extracellular neuraminidase [START_REF] Musser | Identification of a highvirulence clone of type III Streptococcus agalactiae (group B Streptococcus) causing invasive neonatal disease[END_REF]. These isolates, as 90% of the strains responsible for late onset disease [START_REF] Baker | Group B streptococcal infections[END_REF], belonged to serotype III.

If the physiopathology of early onset GBS disease is well-documented, little is known about late or recurrent GBS infections. This case offers a novel hypothesis explaining how GBS can cause neonatal infection from 7 days to 3 months following delivery.

In most of the rare cases described [START_REF] Kotiw | Late-onset and recurrent neonatal Group B streptococcal disease associated with breast-milk transmission[END_REF][START_REF] Olver | Neonatal group B streptococcal disease associated with infected breast milk[END_REF], there were no signs of maternal mastitis, indicating a silent maternal duct colonization. Moreover, National Committee of Hygiene guidelines do not systematically screen for GBS in the raw maternal milk supply [START_REF] Anonymous | Breastfeeding and the use of human milk[END_REF]. These may be two reasons for the underestimation of maternal milk as a source of GBS infection. In the present case, genetic analysis affords evidence for maternal milk as the source of neonatal GBS infection. A circular process was hypothesized by Kotiw et al [START_REF] Kotiw | Late-onset and recurrent neonatal Group B streptococcal disease associated with breast-milk transmission[END_REF]. GBS initially colonizes the neonate's oropharynx mucosa from perinatal or other sources, infecting maternal ducts during breastfeeding. The organism multiplies in the milk ducts. As the microbial concentration increases in the milk, the infant is re-infected during breastfeeding. Mastitis may or may not be present [START_REF] Kotiw | Late-onset and recurrent neonatal Group B streptococcal disease associated with breast-milk transmission[END_REF]. However Olver et al described cases of GBS infection in preterm infants fed with maternal milk via nasogastric tube alone [START_REF] Olver | Neonatal group B streptococcal disease associated with infected breast milk[END_REF]. Prematurity is a recognized predisposing factor to GBS infection although breast milk transmission was also described in term infants [START_REF] Byrne | Neonatal group B streptococcal infection related to breast milk[END_REF].

In our unit, expressed mother's milk is systematically pasteurized and frozen for conservation in our lactarium, and each specimen of milk is screened for bacteria before administration to preterm infants. However, preterm infants might also receive raw, freshly-expressed breast milk from their mothers present in the unit. When a preterm infant falls clinically ill while the mother is breastfeeding, the mother's milk should be cultured to rule out or to document possible breast milk transmission. Mother's milk feedings should be suspended while providing banked milk pending culture result. If breast milk is positive for GBS, adequate antibiotherapy should also be prescribed for the mother. Byrne et al reported that it is possible to give the mother the opportunity to continue breastfeeding as desired ; she can be encouraged to maintain her milk supply by pumping and discarding milk until appropriate treatment is completed and negative breast milk cultures are obtained [START_REF] Byrne | Neonatal group B streptococcal infection related to breast milk[END_REF].

Asymptomatic excretion of GBS in breast milk may be an under-recognised cause of neonatal and recurrent infection. Recommendations should be established to prevent these infections, notably in the case of multiple births:

1/ treatment of both twins; 2 / recognition of the possibility of GBS breast milk infection in late onset or recurrent infection; 3/ suspension of breastfeeding upon suspicion of GBS breast milk infection in both of the children; 4/ search for GBS colonization in both mother and children.

However, it should also be noted that use of human milk in the intensive care nursery decreases the incidence of nosocomial sepsis [START_REF] El-Mohandes | Use of human milk in the intensive care nursery decreases the incidence of nosocomial sepsis[END_REF] and breastfeeding should still be considered as the most appropriate nutrition for babies and preterm infants [START_REF] Anonymous | Breastfeeding and the use of human milk[END_REF].

Acknowledgements

We wish to thank Tracey Montagnon for reviewing the manuscript.

Conflict of interest statement

All authors, no conflict of interest.


Lucian Petrescu

Essaïd Bilal
email: bilal@emse.fr

Luisa Elena Iatan

The impact of a uranium mining site on the stream sediments (Crucea mine, Romania)

Keywords: impact of mine wastes, stream sediments, uranium, thorium

XRF methods were used to evaluate the impact of uranium mine dumps on the stream sediments from Crucea region (Romania). In order to estimate the natural and anthropogenic inputs of radioactive and heavy metals in the sediments, normalization to Al was applied. The pollution degree of the bottom sediments show that U, Th and Pb reach medium and punctual high values, while the rest of the elements appears in concentrations close to the background or lower. The measurements carried out in the surroundings of a local uranium mine show that the impact of Crucea mine on water quality downstream of mining area is insignificant.

I. Introduction

Uranium and thorium are omnipresent in our human environment. Various anthropogenic activities involving the processing or use of uranium rich materials in may modify the natural abundance of uranium in stream sediments. These activities include various mining operations and the industrial processing for the manufacturing of nuclear fuel and of other products. Between 1962 and2004, the National Company of Uranium (former Romanian Autonomous Administration for Rare Metals), mined over 1,200,000 tones of pitchblende ore in the East Carpathians. Traditional mining methods rely on the relocation of large volumes of mineralized and waste rock. These intensive mining operations have created over 30 solid radioactive mining wastes in Crucea-Botusana area (Bistrita Mountains). In this paper we present data from four main tributaries of Bistrita River (i.e. Crucea, Troaca Gavanului, Gavan and Cracul Rau creeks) and from two sampling points locating in the upper course of this main stream (one point located downstream and the other located upstream of the confluence with Crucea creek). The aim of this paper is to investigate the impact of mining activities on the stream sediments in the area.

II. The study area

I II I. .1 1. . G
Ge eo og gr ra ap ph hy y The study concerns the uranium mineralization located within Crucea ore deposit, in the East Carpathians, Romania. The uranium ore deposit found in the eastern part of the Bistrita Mountains (40 km southeast of the town of Vatra Dornei) in the headwaters of Crucea, Lesu and Livezi valleys. At present, this is the largest uranium mine in the country. In the past, the mining area covered 18 km 2 , but was gradually overtaken by logging activities. The area has a steep mountainous relief cut by narrow valleys. The surface streams are tributaries of Bistrita River. I II I. .2 2. . G Ge eo ol lo og gi ic ca al l b ba ac ck kg gr ro ou un nd d The host rocks belong to the Crystalline-Mesozoic zone of the northern part of the East Carpathians. The uranium mineralization is confined to a body of brecciated retrometamorphic micaschists, containing carbonates and clay minerals. The richest uranium mineralization is associated with four NW-SE faults with E or W strike. Most mineralized areas occur along fractures, thus denoting a significant structural control over the emplacement of the ore-bodies. Uranium ore bodies form darkgrey lenses or veinlets with thicknesses ranging from decimetres to meters. Primary mineralization consists of pitchblende U 3 O 8 included in bitumen. The intimate association between pitchblende and bitumen is characteristic to this uranium ore deposit. The sulphide content is lower than 10% and consists of pyrite, marcasite, chalcopyrite, and minor galena, sphalerite, arsenopyrite and tetrahedrite in a carbonate gangue. The mineralization is considered to be of hydrothermal origin (Petrescu, 2004). The geochemical background of uranium in the surrounding area is lower than 3-mg.g -1 , whereas in the mining area, it reaches 5.39-mg.g -1 in the waste rocks (rocks with subeconomic mineralization) and 7.26-mg.g -1 in soils (Petrescu and Bilal, 2007). The exploration and mining works include thirtytwo adits, situated between 780 and 1040 m above sea level. Radioactive waste resulted from mining are disposed next to the mining facilities (Petrescu and Bilal, 2007). The waste rock was disposed in piles of variable sizes that are spread over an area of 364,000 m 2 . The average geochemical content of the waste rocks is similar to that of the ore. The average particle size of the waste material is large (>5 cm). Older dumps (18) have been already naturally reclaimed by forest vegetation, which played an important role in stabilizing the waste dump cover and in slowing down the uranium migration processes.

III. Material and methods

The sediment samples were collected at seventeen sampling points (e.g. 1, 3, 6, 11, 16, 19, 23, 25, 27, 29, 32, 33, 35, 41, 43 and 45) (Figure 1) from the bottom of the stream waters. To ensure the representativity of samples, bulk sediments were obtained by mixture of several cores collected. The pool of sediments were collected by manual coring (0-5-cm) from the moist shore and stored in glass jars. They were dried in an oven at 50°C for 48 hr, ground using an agate grinder to pass a 100-µm sieve and analyzed for the following: I II II I. .1 1. . T To ot ta al l e el le em me en nt ta al l c co on nc ce en nt tr ra at ti io on n The total metal concentrations of bottom sediments were determined using a SRS3400 Bruker AXS wavelength dispersive X-Ray fluorescence spectrometer (WD-XRF). Thus, 1.3gram of calcinated powder was added to 1.3-g of oxides (lithium nitrate) and 6.5-g of flux (lithium tetraborate 20%, lithium metaborate 80 %) and then mixed in a crucible that is placed in an automated fusion-fluxing device. The fusion is performed at 1100°C. The cooling is forced to avoid crystallization, and then fused beads with a vitreous, homogeneous structure and flat surface were obtained. The detection limits, in mg.g-1, were: 4 for U, 5 for Th, 5 for Pb and 10 for V. I II II I. .2 2. . S St ta at ti is st ti ic ca al l d da at ta a a an na al ly ys si is s The reliability of the chemical analyses was calculated from the four duplicate samples with the analysis of variance. The F-ratio, which is a measure of the ratio of total data variance to error variance (sample inhomogeneity, analytical error) and which has to be > 4 to be significantly high at the 99% level of confidence, was >1000 for Ca, Cr, Cu, Mn, Pb, V, Na and Zn, >100 for U, Th, Sr, Al, Co, Fe, K, P and Mg and 54 and 32 for Si and Ni respectively.

Hence, any pattern or absence of pattern of any of these elements in the data set is not due to analytical or sampling error. Single correlation analysis was performed to assess relationship among heavy elements concentration. Linear correlation coefficients (r) were calculated between the different samples parameters. For each correlation, statistical significance (p-value) was determined. Traditionally, results that yield p-value < 0.05 are considered borderline statistically significant, i.e. the probability of concluding that two variable are different when they are actually the same is ≤5%. The STATISTICA 7.1 for Windows (StatSoft ® ) program was used for the statistical analyses mentioned above.

IV. Results and discussion

Table 1 shows levels of heavy metals, aluminium oxides and total organic matter (lost on ignition =LOI) found in bottom sediments sampled at the sixteen stations. In all the areas, similar variations in the concentrations of uranium and thorium in sediment samples were observed. The bottom sediment samples collected from mining area show a wide range of total uranium content from 17.20 to 5023.96-mg.g -1 with a mean value of 246.81-mg.g -1 (dry wt.) (Figure 2). Total thorium concentration in sediments collected from Crucea site ranged from 18.70 to 6643.92-mg.g -1 with an average of 391.17-mg.g -1 (dry wt.) (Figure 2). This indicates that the adsorption of the radioactive elements by the sediments is high and variable, influenced by the ore dump -sample relationship. There is a decrease in the concentrations of U and Th in the stream sediments downstream from the Crucea uranium ore (Figure 2). The fate and transport of uranium is governed by its oxidation state, which is either hexavalent (U(VI)) or tetravalent (U(IV)). When mine wastes containing low-grade uranium mineralization are exposed to surface waters containing carbonate ions, the pitchblende is altered, uranium (IV) oxidizing to uranyl carbonates, without producing any megascopically observable changes. This can be illustrated by the following equation: U3O8 + 2H2O + 9 CO3 2-→ 3[UO 2 (CO3)3] 4-+ 4OH + 2e - A previous study performed in this area (Petrescu and Bilal, 2006) showed that in the pH range of waters from Crucea area, more than 85% of the uranium was present as the carbonate complexes UO 2 (CO 3 ) 2 2-, UO 2 CO 3 (OH) 3-, UO 2 CO 3 and UO2(CO3)3 4-and 17% as the hydroxide complexes UO2(OH)2. The anionic carbonate complexes are extremely important because they increase uranium mobility by limiting the extent of uranium adsorption in oxidized waters and by increasing the solubility of uranium minerals [START_REF] Langmuir | Aqueous environmental chemistry[END_REF]. Strong correlation of uranium with thorium in the sediments (r=1; n=17; p<0.0001) suggests that unspecific co-precipitation dominates over selective adsorption according to valence and ion sizing. The lack of correlation between uranium and vanadium in sediments samples (r=0.16; n=17; p<0.0001) is related to ore type. In the Crucea uranium mineralization the vanadium concentration range between 60.7 -335.9-mg.g -1 and we can state that the surrounding rocks are the source of this element. This affirmation confirms the absence of the uranium minerals with vanadium (such as carnotite). Through the weathering of the rocks, vanadium migrates as vanadic acid and in supergene conditions this metal easily concentrates in clay minerals and oxy-hydroxides. In order to estimate the natural and anthropogenic inputs of actinides and heavy metals into sediments, their total concentrations were correlated with Al concentration. Since Al originated from the natural background it was chosen as the normalizing factor. These data are presented in Table 2. It can be seen that only vanadium shown a linear relationship with the concentration of Al, fact indicates his natural background origin. For the rest of selected elements their origin was anthropogenic. This indicates the impact of mining waste disposal on the environment.

The total concentrations shown in Table 1 have been used for the determination of the degree of pollution with radioactive and heavy metals. At this time, no regulations exist for levels of heavy metals in sediments; however several sets of sediment quality criteria have been reported. Thus, the pollution degree of the sediments was classified using the geoaccumulation index (Igeo), which is applied for separate element characteristic in the examined sediments and determines the contribution of each element to overall sediment contamination. Igeo was calculated as follows: et al., 2004), were: Cn = total element concentration (mg.g -1 ); Bn = background concentration for each element (mg.g -1 ). Therefore, based on the concentration levels of these metals, the stream sediments are categorized into seven classes, with the "0" class corresponding to the values close to the geochemical background. The background concentrations were taken from Petrescu and Bilal (2007) (Table 3). The pollution degree of the bottom sediments is shown in Table 4. The Igeo of U, Th and Pb presents medium and punctual high values that represent sediments belonging to classes of strong to extreme pollution (Igeo > 6). A proportion of 47.05% of Th and U concentration values are categorized into the highest pollution class (Igeo class 6) (Figures 3 a andb), while the rest of the elements presents concentration close to the background values or lowers (Figure 3d). These high concentrations of U, Th and Pb are distributed preferentially in the proximities of 5, 8, 9 and 1bis mine dumps. Table 4 also shows that Igeo of U and Th from the sediments taken from downstream of the mineralized zone (P27 and P29) are lower than those taken from the mineralized areas. This indicates that uranium and thorium are not very mobile in the waters and that these elements are precipitated close to the mining dumps.

2 log [ /(1.5 )] geo n n I C B = (Pinto

V. Conclusions

The knowledge concerning uranium and thorium concentration in stream sediments is essential for predicting natural actinide migration and for drafting strategies concerned with the rehabilitation of contaminated sites. The pollution degree of the bottom sediments shows that U, Th and Pb are characterized by medium and punctual high values, while the rest of the elements are in concentration close to the background values or lower. Although neither U nor Th have an appreciable "exchangeable" fraction, the isolation of specific U-and Th-rich sediment fractions helped to identify connections between bioavailability and genesis of sediments, which control the natural cycling of U and Th. However, if the geochemical conditions it happens to be changes, the surface sediments could act as a future source of actinides and heavy metals to the above water body. The measurements carried out in the surroundings of a local uranium mine show that the impact of Crucea mine on sediment quality downstream of mining area is insignificant.

FiguresF

  Figures

Acknowledgements

This work was supported by MIRA Project 050094501 of the Rhône-Alpes Region (France). The inorganic chemical components were determined at Department GENERIC of the Ecole Nationale Supérieure des Mines, Saint Etienne (France). The authors thank Dr. Jacques Moutte and Frédéric Gallice for technical assistance. Thanks are due to Dr. G. Ilinca (University of Bucharest) for the reading of the manuscript and for the helpful suggestions. Also, we would like to thank anonymous reviewers for their constructive inputs.

. The impact of the environment of the uranium mine dumps from Crucea and Botusana ore deposit (Bistritei Mountains, Romania): PhD thesis, Department of Mineralogy, University of Bucharest, 196 p. (in Romanian with English abstract). Petrescu L., Bilal E., 2006. Speciation of natural uranium and thorium in surface waters around a uranium mine (Bistrita Mts., Romania). Geochimica & Cosmochimica Acta, vol. 70, issue 18, supplement 1, pp. A486 Petrescu L., Bilal E., 2007. Environmental impact assessment of a uranium mine, East Carpathians, Romania: Metal distribution and partitioning of U and Th. Carpathian Journal of Earth and Environmental Sciences, vol. 2/1, pp. 39-50. Pinto M.M.S.C., Silva M.M.V.G., Neiva A.M.R., 2004. Pollution of water and stream sediments associated with the Vale De Abrutiga Uranium Mine, Central Portugal. Mine, Water and the Environment, v. 23, pp. 66-75.


Cécile Straub

Philippe Quillardet

Julia Vergalli

Nicole Tandeau De Marsac

Jean-François Humbert
email: humbert@biologie.ens.fr

A day in the life of Microcystis aeruginosa strain PCC 7806 as revealed by a transcriptomic analysis

Keywords: Cyanobacteria, Microcystis aeruginosa, transcriptome, microarray, day/night cycle, metabolism

The cyanobacterium, Microcystis aeruginosa, is able to proliferate in a wide range of freshwater ecosystems and to produce many secondary metabolites that are a threat to human and animal health. The dynamic of this production and more globally the metabolism of this species is still poorly known. A DNA microarray based on the genome of M. aeruginosa PCC 7806 was constructed and used to study the dynamics of gene expression in this cyanobacterium during the light/dark cycle, because light is a critical factor for this species, like for other photosynthetic microorganisms. This first application of transcriptomics to a Microcystis species has revealed that more than 25% of the genes displayed significant changes in their transcript abundance during the light/dark cycle and in particular during the dark/light transition. The metabolism of M. aeruginosa is compartmentalized between the light period, during which carbon uptake, photosynthesis and the reductive pentose phosphate pathway lead to the synthesis of glycogen, and the dark period, during which glycogen degradation, the oxidative pentose phosphate pathway, the TCA branched pathway and ammonium uptake promote amino acid biosynthesis. We also show that the biosynthesis of secondary metabolites, such as microcystins, aeruginosin and cyanopeptolin, occur essentially during the light period, suggesting that these metabolites may interact with the diurnal part of the central metabolism.

Introduction

Microcystis aeruginosa is one of the most common bloom-forming cyanobacteria in freshwater ecosystems, and is widely distributed on the five continents. Members of this genus are known to be able to synthesize numerous secondary metabolites, including microcystins, which are hepatotoxins frequently involved in animal and human poisoning [START_REF] Briand | Health hazards for terrestrial vertebrates from toxic cyanobacteria in surface water ecosystems[END_REF][START_REF] Kuiper-Goodman | Human health aspects[END_REF]. The genomes of two M. aeruginosa strains (PCC 7806 and NIES-843) have recently been sequenced [START_REF] Frangeul | Highly plastic genome of Microcystis aeruginosa PCC7806, a ubiquitous toxic freshwater cyanobacterium[END_REF][START_REF] Kaneko | Complete genomic structure of the bloom-forming toxic cyanobacterium Microcystis aeruginosa NIES-843[END_REF], revealing their high plasticity, due in part to the presence of numerous, long, repeated sequences and transposase genes.

The annual life cycle of M. aeruginosa is characterized by two distinct phases. The first is a benthic phase, during which it lives on sediment during the Winter and the early Spring, and the second a planktonic phase, during which it lives in the water column during Summer and Autumn [START_REF] Reynolds | Seasonal variations in the vertical distribution and buoyancy of Microcystis aeruginosa Kiitz. Emend. Elenkin in Rostherne Mere[END_REF]. During this second phase, M. aeruginosa cells are organized in colonies, which carry out daily vertical migrations in the water column [START_REF] Rabouille | Simulation of carbon reserve dynamics in Microcystis and its influence on vertical migration with Yoyo model[END_REF]. These vertical migrations are due to hydrodynamic processes, and also to the ability of these microorganisms to regulate their buoyancy [START_REF] Reynolds | Water-blooms[END_REF][START_REF] Ibelings | Microcystis changes its buoyancy in response to the average irradiance in the surface mixed layers[END_REF]. One of the most spectacular effects of this process is the accumulation of M. aeruginosa colonies in the top few centimeters of the water column in the early morning, and their subsequent progressive disappearance from the surface over the course of the day. The light exposure of M. aeruginosa cells can be very limited due to these vertical migrations and to low light levels in the water column when high cyanobacterial biomass reduces water transparency.

This limited exposure to light can have very important consequences for the metabolism of this photosynthetic microorganism, and may have resulted in adaptative mechanisms to ensure the optimal use of the light available. For example, two recent publications have suggested that M. aeruginosa is better able to tolerate lack of light than green algae [START_REF] Tan | Seasonal variation of Microcystis in Lake Taihu and its relationships with environmental factors[END_REF][START_REF] Wu | Physiological comparison between colonial and unicellular forms of Microcystis aeruginosa Kütz[END_REF]. However, little firm information is available about the metabolism of this species during the light/dark cycle. Microarray technology, which provides an overall view of gene expression over time, is one of the best ways to obtain such data. Several recent publications on three cyanobacterial species, Synechocystis sp. strain PCC 6803 [START_REF] Hihara | DNA microarray analysis of cyanobacterial gene expression during acclimation to high light[END_REF][START_REF] Gill | Genome-wide dynamic transcriptional profiling of the light-to-dark transition in Synechocystis sp. strain PCC 6803[END_REF], Cyanothece sp. strain 51142 [START_REF] Stöckel | Global transcriptomic analysis of Cyanothece 51142 reveals robust diurnal oscillation of central metabolic processes[END_REF] and Prochlorococcus marinus strain MED4 [START_REF] Zinser | Choreography of the transcriptome, photophysiology, and cell cycle of a minimal photoautotroph, Prochlorococcus[END_REF] have demonstrated the effectiveness of this approach.

In the present work we set out to construct a microarray based on the genome of M. aeruginosa strain PCC 7806 in order to study the dynamics of gene expression in this cyanobacterium during the light/dark cycle.

Methods

Strain and culture conditions

Microcystis aeruginosa strain PCC 7806 (hereafter designated M. aeruginosa), provided by the Pasteur Culture Collection of Cyanobacteria (Institut Pasteur of Paris), was grown in BG11 medium [START_REF] Rippka | Genetics assignments, strain stories and properties of pure cultures of cyanobacteria[END_REF] containing 2 mM of NaNO 3 and 10 mM of NaHCO 3 . Cultures were placed in an incubator (INFORS AG, Switzerland) continuously sparged with 1% (v/v) CO 2 in air, at 22°C, shaken (90 rpm), and exposed to 12h-12h light-dark cycles. The photosynthetic photon flux density during the light period (50 µmol photons.m -2 .s -1 ; cool white OSRAM L 18/640) was measured using a LICOR LI-185B quantum/radiometer/photometer equipped with a LICOR LI-193SB spherical sensor.

Experimental protocol

Two independent experiments, each with three biological replicates, were conducted. in 1L Erlenmeyer flasks. To train the circadian clock before sampling, the cultures were incubated for 10 days under the same experimental conditions. Cell growth was followed by monitoring the OD 750 using a Uvikon 933 VIS-US spectrophotometer (Kontron Instruments).

Samples were taken in the middle of the exponential-growth phase (OD 750 0.6 to 0.8)

of the cultures over a period of 24h (all times are expressed in terms of a 24-h clock): at T0

(00:00, the middle of dark period), T1 (05:30, i.e. 30 min before the transition from dark to light), T2 (06:30, i.e. 30 min after the transition from dark to light), T3 (07:30, i.e. 90 min after the transition from dark to light), T4 (12:00, the middle of the light period), T5 (18:30, i.e. 30 min after the transition from light to dark) and finally T6 (24:00, the middle of the dark period) (Supporting Information S1). T6 of this cycle corresponds to T0 of the next 24-h cycle. For each time point, 50 ml of the cultures of each replicate were collected, and centrifuged (15000 g, 2 min, 4°C). The cell pellets were then frozen in liquid nitrogen, and stored at -80°C until used for RNA extraction. In the first experiment, samples were collected only at times T0, T2, T4, T5 and T6, whereas in the second experiment all the time points were sampled.

RNA isolation

Total RNA was isolated and purified using the kit RiboPure TM -Bacteria (Ambion) according to the manufacturer's instructions. Cell pellets were resuspended in 750 µl RNAwiz, and 500 µl Zirconia Beads were added. The bacterial cells were lysed for 200 s (10 x 20 s with 1 min intervals on ice) in a FastPrep FP120 (Qbiogene) at speed 5 m.s -1 . DNase I treatment was not performed. RNA concentrations were quantified using a NanoDrop ND-1000 Spectrophotometer (NanoDrop Technologies), and the quality of the RNA extracts was checked using the Agilent 2100 Bioanalyzer, RNA 6000 Nano reagents and RNA Nano Chips (Agilent Technologies) according to the Manufacturer's instructions.

cDNA synthesis and Cy-dye labeling Five µg of the purified total RNA was used for the random priming of cDNA synthesis and fluorescence labeling (Cy3 or Cy5 dye, Amersham, GE Healthcare) using the SuperScript TM Indirect cDNA Labeling System (Invitrogen) according to the manufacturer's instructions.

Indirect labeling relies on the incorporation of aminoallyl-modified nucleotides into cDNA derived from the target RNA. The aminoallyl side chains were then used to label the cDNA with either Cy3 or Cy5. The efficiency of the labeling procedure was analyzed by a NanoDrop ND-1000 spectrophotometer (NanoDrop technologies). The labeled cDNA was used immediately for hybridization to a microarray.

Construction of Microcystis aeruginosa microarrays

The transcriptomic study was performed using an Agilent 4X44K two-color microarray (Agilent technologies). Each slide contained 4 arrays each with 44,000 spots. Specific oligonucleotides (60-mers) were designed for 5085 protein-coding genes corresponding to 96% of the total predicted protein-coding genes in the annotated genome of the M. aeruginosa [START_REF] Frangeul | Highly plastic genome of Microcystis aeruginosa PCC7806, a ubiquitous toxic freshwater cyanobacterium[END_REF]. Oligonucleotide probes for each gene were designed by the eArray protocol (Agilent technologies). Five, four, three, two and one oligonucleotide probes were defined for 68%, 9%, 11%, 9% and 3% of genes, respectively. Each oligonucleotide was randomly printed in duplicate across one array by means of the SurePrint technology (Agilent Technologies).

Hybridization and scanning of DNA microarrays

Hybridization was performed according to the protocol for the two-color microarray-based gene expression analysis recommended by Agilent Technologies. Briefly, Cy3-labeled cDNAs (150 pmol) from the control culture (T0) was mixed with Cy5-labeled cDNAs (150 pmol) from other time points (T1 to T6), and they were hybridized to one array. The color-swapped cDNA samples (i.e. the same samples with the reverse combination of dye labels) were hybridized to the other array on the same glass slide. The resulting mixture was hybridized to an array at 65°C for 17h. Following hybridization, each slide was washed as recommended by Agilent. After the slide had been air dried, it was scanned using a GenePix 4200AL scanner (Axon Instruments, Union City, Calif.) at 5µm resolution. The signal intensity of each spot and its local background were determined using GenePix Pro software (Version 6, Axon Instruments). The net signal intensity was calculated by subtracting the median signal intensity of all the pixels within the local background area from the median signal intensity of all the pixels within the spot area. All the spot areas recognized by the automatic alignment function of the GenePix Pro were confirmed visually. The flagged spots

were not used for subsequent data analysis.

Microarray data analyses

Any bias in signal intensity between the two fluorescent dye channels in a microarray was normalized by locally weighted linear regression analysis or Lowess normalization [START_REF] Quackenbush | Microarray data normalization and transformation[END_REF] using MIDAS (freely available from http://www.tm4.org/midas.html). For all normalizations, the smoothing parameter was set to 0.33. Low-intensity spots (Cy3 or Cy5 intensity <2500) were filtered and then, pairs of flip-dye data were subjected to a flip-dye consistency check in which data outside of a two standard-deviations cut-off were removed. The relative expression level of a gene at a given time point was calculated as log2 (Cy5/Cy3), where Cy5 and Cy3 were normalized signal intensities from experimental and control cDNAs.

Several probes (usually, five) were designed for each protein-coding gene, and as two independent experiments were performed in triplicate (biological replicates) and swapping two colors (CY-3/CY-5, CY-5/CY-3, technical replicates), it follows that the transcript abundance obtained for each gene at each sampling time was based on 12 to 60 measurements.

Statistical analyses of the data were performed using the Significance Analysis of Microarrays (SAM) algorithms [START_REF] Tusher | Significance analysis of microarrays applied to the ionizing radiation response[END_REF] in MeV (http://www.tm4.org/mev.html), on the median of log2-transformed signal ratios (fold change) of the replicate spots. This analysis was done gene by gene (time point TX/T0 tested against time point T0/T0; One class analysis), but also by taking into account all the genes in the same analysis (Multi class analysis). SAM analysis uses repeated permutations of the data to find out whether the expression of any genes is significantly related to the response. The data for each gene are permuted, and a t-test is computed for both the original and the permuted data for each gene. One valuable feature of SAM is that it gives estimates of the False Discovery Rate (FDR), which is the proportion of genes likely to have been picked out by chance as being significant. Stringent parameters were applied for this criteria because we chose a delta value so that estimated FDR was equal to zero. Finally, a gene was considered as presenting a significant variation in its transcript abundance when it was found significant by SAM analysis (One class or Multi class) for at least one sampling time.

The time point T0 consisted in examining one sample against itself (T0/T0), so that the log2-transformed signal ratio for each gene is expected to be near the value 0. As a control, we have checked using SAM that at time T0, no gene displayed a mean variation of their expression level that was significantly different from 0.

Results

General overview

During the normalization step, a strong dye bias (CY-3/CY-5 data very different from CY-5/CY-3 data) was found for the genes displaying very high expression values (Table S1).

These values were eliminated from the analysis by applying a flip dye consistency check after the normalization step, a function provided by the MIDAS software. As a result, 70 genes were eliminated from the analysis. For 258 genes, data were only available at certain time sampling points, whereas a complete time series data was obtained for 4757 genes. Very good overall reproducibility of the results was found for these genes, as illustrated for example for hypD in Supporting Information S1.

To obtain an overview of the global pattern of gene expression during the 24-h light/dark cycle, we performed a principal component analysis (PCA) on the transcript abundance values of the 4757 genes (displaying a complete time series data) obtained for each replicate at each sampling time. The projection of these points onto the plane defined by the first two axes accounting for 45.8% of the total variance (Fig. 1), displayed a spatial ordination of the samples during the light/dark cycle. For each sampling time, all the replicates of the two independent experiments were grouped together and clearly separated from those for the other sampling times. Moreover, the spatial ordination of these points described a clockwise circle beginning at 00:00 and finishing at 24:00.

Temporal dynamics of the transcriptome

From the statistical analyses, it appeared that the transcript levels of 1344 genes varied significantly at least at one sampling time during the light/dark cycle, which is equivalent to 25% of the M. aeruginosa protein-coding genes identified in the genome. The complete data set is provided in Table S2. We found also that among these 1344 genes, 40% displayed a zenith (peak-up) or a nadir (peak down) in their transcript abundance at 07:30 (i.e. 1h30 after the transition from dark to light) knowing that most of them were characterized by a maximum transcript abundance (Fig. 2A). For most of these genes, changes could already be observed at 06:30, 30 min after the lights had been switched on. At 12:00, the number of genes differentially expressed steadily decreased, and then remained stable untill 18:30 (i.e.

30 min after the transition from light to dark). In the middle of the night (24:00), around 100 genes still displayed significant variations in their expression relative to T0. Overall, the maximum or minimum values of the expression of all the genes usually occurred at 07:30, with a second peak at 18:30 (Fig. 2B). This indicates that light/dark (or dark/light) transition triggers major changes in the expression of the M. aeruginosa genome.

Transcriptome analysis and clusters of orthologous genes (COG) categories

The 1344 genes displaying significant variations during the light/dark cycle were further sorted using COG (cluster of orthologous genes) categories [START_REF] Tatusov | The COG database: a tool for genome-scale analysis of protein functions and evolution[END_REF]. It was possible to assign a COG category to 853 (63%) of these 1344 genes. On the basis of their pattern of expression during the light/dark cycle, the COG-assigned genes could be divided into 4 groups (Fig. 3).

The first group (29% of the COG-assigned genes) contained genes that peaked at both 07:30 and 18:30 (Fig. 3A). These genes belonged to five COG categories (C: Energy production and conversion; E: Amino acid transport and metabolism; F: Nucleotide transport and metabolism; G: Carbohydrate transport and metabolism and K: Transcription), and among them, the highest transcript levels were obtained for genes involved in photosynthesis and energy metabolism.

The second group (36% of the COG-assigned genes) contained genes that peaked only at 07:30 (Fig. 3B). Genes in this group belonged to nine COG categories (H: Coenzyme metabolism; I: Lipid metabolism; J: Translation, ribosomal structure and biogenesis; M: Cell envelope biogenesis, outer membrane; N: Cell motility and secretion; P: Inorganic ion transport and metabolism; T: Signal transduction mechanisms; U: Intracellular trafficking, secretion, and vesicular transport and V: Defense mechanisms;) and among them, the highest transcript levels were found for genes involved in lipid metabolism and translation processes.

The other two groups included genes representing only a small number of COG categories. The third group (10% of the COG-assigned genes) contained genes displaying high transcript levels throughout the light period (Fig. 3C). These genes belonged to three COG categories (L: DNA replication, recombination and repair; O: Posttranslational modification, protein turnover, chaperones and Q: Secondary metabolites biosynthesis, transport and catabolism). The fourth group contained genes belonging to only one COG category (D: Cell division and chromosome partitioning), and displaying high transcript levels in the middle of the light period (Fig. 3D).

Finally, 25% of the COG-assigned genes belonged to categories R and S, which are poorly characterized.

Carbon transport and energy metabolism

Most of the genes with expression peaks at both 07:30 and 18:30 are involved in carbon transport and energy metabolism processes (COG categories C and G). A more detailed analysis shows that transcript abundance of the ccm (carbon dioxide concentrating mechanism) and rbc (ribulose bisphosphate carboxylase) genes peaked soon after the transition from dark to light and then rapidly decreased after the transition from light to dark (Fig. 4A). Similar patterns of expression were found for most of the genes encoding the ndh (NADH dehydrogenase subunits) and for the atp genes (ATP synthase subunits) (Figs. 4C and4D). On the other hand, the opposite pattern of expression was found for genes involved in the respiratory terminal oxidase system (cta genes), such as the genes encoding three cytochrome c oxidase subunits and the protoheme IX farnesyltransferase (Fig. 4B).

Only four pet genes (cytochrome b 6 f complex) displayed significant variations in their transcript abundance. For the petCA gene cluster, the maximum transcript level occurred at 07:30 (Supporting Information S1), while the petBD gene cluster peaked at 12:30 and 18:30 (Supporting Information S1). Like petBD, but unlike petAC, cydA and cydB (cytochrome bd-I quinol oxidase) displayed their lowest transcript levels in the morning, and their highest ones in the afternoon (Fig. 4B). Finally, the hoxEFUYH (bidirectional hydrogenase) and the hypCD genes (maturation of the enzyme) displayed similar patterns of expression, with transcript levels peaking at the beginning of the dark period (Table S2).

Eleven photosystem I genes (psa) and 20 photosystem II genes (psb) have been identified in M. aeruginosa. These include psbD, psbD2, psbQ, psbW, psb27 and psb29, which displayed significant changes in their transcript patterns with maximum values at 06:30

or 07:30 (Fig. 5). On the other hand, psaK2 transcript abundance was very low at 05:30 and 06:30, but peaked at 12:00. The psbU gene displayed a single transcript abundance peak, which occurred at the beginning of the dark period (18:30).

Finally, when considering the TCA cycle, four genes (isocitrate dehydrogenase(icd), citrate synthase (gltA), succinate dehydrogenase flavoprotein subunit (sdhA) and succinate dehydrogenase iron-sulfur protein subunit (sdhB)) among the seven genes described in the M. aeruginosa genome for this pathway, displayed significant changes in their transcript abundance, with maximum transcript abundance values during the night and minimum values just at the beginning of the light period (Fig. 6).

Nitrogen and phosphate transport and metabolism

The transportation and metabolism of nitrate and phosphate are very important for freshwater cyanobacteria as these nutrients are generally growth limiting in continental aquatic ecosystems. While the transcript levels of narB (nitrate reductase) did not significantly change during the light-dark cycle, the expression of the nirA (nitrite reductase) and nrt (nitrate/nitrite transport) genes both peaked soon after the light was switched on (Fig. 7A). In contrast, the expression of the three ammonium permease genes occurred later (amt1 and amt2: max value at 12:00; amt3: max value at the beginning of the dark period) (Fig. 7B). The same pattern of expression as for amt3 was found for the genes involved in the transport and synthesis of glutamate (gltS and gltBD, respectively) or of glutamine (glnH and glnA, respectively), as well as for the gltA (citrate synthase) and icd (isocitrate dehydrogenase), two genes required for the synthesis of 2-oxoglutarate, and for the glnB gene that codes for the signal transduction PII protein, which is known to coordinate nitrogen and carbon metabolism (Fig. 7B and Table S2).

It is worth noticing that gltBD (NADH-dependent glutamate synthase subunits) and glsF (ferredoxin-dependent glutamate synthase) display opposite patterns of expression, with a high level of glsF transcripts observed from 06:30 to 12:00 followed by a decrease during the dark period, which is when the gltBD transcripts reached their maximum value. Like the nirA and nrt genes, the transcript levels of argJ (a bifunctional enzyme involved in the first step of the arginine biosynthetic pathway from glutamate) display a sharp increase followed by a decrease during the light-dark cycle, peaking at 07:30 (Fig. 7B and Table S2).

Three ABC-type phosphate transporter genes (pst) displayed only slight variations in their transcript abundance, which was characterized by a slight increase between 00:00 and 07:30, followed by a decrease. In contrast, the two phosphate regulator genes phoB (a twocomponent system response regulator of a phosphate sensing system) and phoU (a phosphate uptake regulator) displayed contrasting patterns of expression, with a marked decrease in phoB transcript levels at 12:00 and 18:30, whereas those of phoU increased during the night (18:30 and 24:00) (Table S2).

Cell division and chromosome partitioning

The transcript levels of most COG category D genes (D: Cell division and chromosome partitioning) were high in the middle of the light period. However, their pattern of expression was essentially due to variations in the transcript levels of genes either with unknown functions or only indirectly related to cell division. The transcript abundance of only one gene that may be involved in chromosome partitioning (MIC_1479) peaked in the middle of the light period (Table S2). The transcript abundance of ftsH, which encodes a protein essential for cell division and classified as COG category O, peaked at 5:30, just before the dark/light transition, whereas ftsK/spoIIIE, which encodes a putative cell division protein, and minD which encodes a septum site-determining protein peaked at 6:30, just after the dark/light transition. Similarly, genes encoding rod shape-determining proteins (mreB; mreC, belonging to the COG category M) displayed a transcript abundance peak at 07:30 (Table S2).

Secondary metabolites COG analysis of the genes belonging to category Q (Q: Secondary metabolites biosynthesis, transport and catabolism) has shown that most of the genes involved in the biosynthesis of secondary metabolites displayed high transcript levels throughout the light period. More detailed analysis showed that all the genes involved in the biosynthesis of microcystins (mcy genes), with the exception of mcyB and mcyD, were characterized by an increase in their transcript abundance shortly after the light was switched on (at 06:30 and 07:30, respectively), then remained stable or declined slightly until night began (18:30), and finally displayed a sharp decrease when the light was turned off (Fig. 8A). The mcyB and mcyD genes displayed the same overall transcription patterns, but the changes were not statistically significant due to wide variations in the values obtained for the different replicates. The transcription patterns of the genes involved in the biosynthesis of both aeruginosin and cyanopeptolin were generally similar to those of the microcystin genes, with only a slight difference for the mcnA and mcnB genes, which displayed a single transcript abundance peak just after the transition from dark to light (06:30), followed by a decrease (Fig 8B and8C).

Similarly, the transcript levels of the genes required for the biosynthesis of the patellamidelike microcyclamide, also rose just after the transition from dark to light (06:30 and 07:30), and then fell during the light period (Fig. 8D). Finally, the two PKS I gene clusters (PKSI modular/PKSIII and PKSI iterative) did not show any significant variation in their transcript levels during the light/dark cycle (Table S2).

Regulation of the light/dark cycle (Circadian clock and sigma factors)

The day/night variations in the abundance of transcripts of numerous genes suggest that they are influenced by both light and an endogenous circadian rhythm. This prompted us to look at the dynamics of the genes involved in the circadian clock (kai genes). Whereas the transcript levels of the kaiA gene showed no significant modifications during the light/dark cycle, the transcription patterns of the kaiB and kaiC genes and of the two-component sensor histidine kinase sasA gene (KaiC-interacting protein), were very similar (Fig 9A). In contrast, the rpaA gene, which together with sasA forms a two-component regulatory system involved in circadian clock regulation, did not display any significant changes in transcript abundance (Table S2). Finally, kaiB1, a paralog of the kaiB gene, displayed an inverse pattern of expression to kaiB and kaiC (Fig. 9A), with a low value at the onset of the light period and a high value at the onset of the dark period.

Sigma factors could also significantly influence the light/dark expression of numerous cyanobacterial genes. Among the sig genes identified in the genome of M. aeruginosa PCC 7806, four sigma factors, sigB, sigD, sigE and sigH, displayed significant changes in their transcript levels during the light/dark cycle (Fig. 9B). Very similar patterns of transcript abundance were observed for sigB and sigE, with minimum values at 07:30 and maximum values at 18:30. The transcript profiles for sigD and sigH differed both from each other, and from those of sigB and sigE (Fig. 9B).

DISCUSSION

As far as we are aware, this study is the first to provide data about the overall expression of the M. aeruginosa genome during a light/dark cycle. More than 25% of the genes displayed significant variations in their transcript abundance during the 24h cycle, most of which occurred following the dark/light transition. Similar results were reported by Stöckel et al.

[13], who found that nearly 30% of the genes of Cyanothece sp. strain 51142 exhibited cycling behavior, and that most of them were involved in core metabolic processes (photosynthesis, respiration, energy metabolism…). On the other hand, for Prochlorococcus marinus strain MED4, which has a much smaller genome (1700 kb) than either Cyanothece or Microcystis (5000 kb), Zinser et al. [START_REF] Zinser | Choreography of the transcriptome, photophysiology, and cell cycle of a minimal photoautotroph, Prochlorococcus[END_REF] estimated that 80% of the genes displayed cyclic expression during the diel cycle. The core genome of Prochlorococcus has been estimated to consist of >1200 genes [START_REF] Kettler | Patterns and implications of gene gain and loss in the evolution of Prochlorococcus[END_REF], which corresponds to more than 70% of the MED4 genome, which means that most of the genes expressing cyclic variations in MED4 belong to the core genome. Similarly, only 176 of the 1344 genes displaying significant variations in the study reported here were specific to Microcystis (i.e. were not found in any of the other twelve cyanobacterial genomes in the study of Frangeul et al. [START_REF] Frangeul | Highly plastic genome of Microcystis aeruginosa PCC7806, a ubiquitous toxic freshwater cyanobacterium[END_REF]). Collectively, these findings indicate that in cyanobacteria most of the genes displaying cyclic variations during a dark/light cycle belong to the core genome, and so are involved in basic metabolic processes.

This also suggests that in cyanobacterial genera with large genomes, such as Microcystis or

Cyanothece, numerous genes are likely to be involved in light-independent adaptation processes. However, as found in Cyanothece [START_REF] Stöckel | Global transcriptomic analysis of Cyanothece 51142 reveals robust diurnal oscillation of central metabolic processes[END_REF][START_REF] Toepel | Transcriptionnal analysis of the unicellular, diazotrophic cyanobacterium Cyanothece sp. ATCC 51142 grown under short day / night cycles[END_REF], our findings showed that light/dark transitions play an important role in the metabolism of Microcystis.

Cyanobacteria are known to be able to perform oxygenic photosynthesis and respiration simultaneously [START_REF] Vermaas | Photosynthesis and respiration in cyanobacteria[END_REF], but a contrast between night and day was obvious in the transcript abundance of the genes involved in photosynthesis and of those involved in respiration. Indeed, the ccm and rbc genes (oxygenic photosynthesis) displayed maximum transcript abundance after the switch to light (06:30 and 07:30), whereas the cta genes (respiration), displayed maximum transcript abundances when the dark period began (18:30).

Similar results were found for rbcL in Prochlorococcus [START_REF] Zinser | Choreography of the transcriptome, photophysiology, and cell cycle of a minimal photoautotroph, Prochlorococcus[END_REF] and Cyanothece [START_REF] Toepel | Transcriptionnal analysis of the unicellular, diazotrophic cyanobacterium Cyanothece sp. ATCC 51142 grown under short day / night cycles[END_REF] during a light-dark cycle. In the same way, most of the genes coding for components of oxidative phosphorylation, such as NADH dehydrogenase (ndh genes) and ATP synthase (atp genes), were up-regulated at the beginning of the light period, as reported in Synechocystis sp. strain PCC 6803 [START_REF] Gill | Genome-wide dynamic transcriptional profiling of the light-to-dark transition in Synechocystis sp. strain PCC 6803[END_REF][START_REF] Kucho | Global analysis of circadian expression in the cyanobacterium Synechocystis sp. PCC 6803[END_REF]. In cyanobacteria, ATP and NADPH are used in many metabolic processes including CO 2 fixation. In particular, NADPH appears to be used preferentially for carbon fixation processes, rather than for respiration [START_REF] Vermaas | Photosynthesis and respiration in cyanobacteria[END_REF][START_REF] Cooley | Succinate:quinol oxidoreductases in the cyanobacterium Synechocystis sp. PCC 6803: presence and function in metabolism and electron transport[END_REF]. On the other hand in cyanobacterial thylakoids, the main respiratory electron transport activity into the plastoquinone pool involves the succinate dehydrogenase genes (sdh) [START_REF] Cooley | Succinate:quinol oxidoreductases in the cyanobacterium Synechocystis sp. PCC 6803: presence and function in metabolism and electron transport[END_REF]. These genes, and others from the TCA branched pathway, displayed the same pattern of transcription as those involved in respiration (data not shown). Finally, the responses of bidirectional hydrogenase genes (HoxFUYH and hypCD) support the hypothesis that these enzymes may be involved in fermentation or respiration, rather than in photosynthesis (see the review of Tamagnini et al. [START_REF] Tamagnini | Cyanobacterial hydrogenases: diversity, regulation and applications[END_REF]).

With regard to nutrient uptake, if, as expected, carbon metabolism occurs at the beginning of the light period, this makes it even more surprising to find ammonium assimilation occurring at the end of the light period and during the first few hours of the night.

Ammonium assimilation has been clearly shown to occur in cyanobacteria (e.g. [START_REF] De Marsac | Control of nitrogen and carbon metabolism in cyanobacteria[END_REF][START_REF] Muro-Pastor | Ammonium assimilation in cyanobacteria[END_REF]). In Spirulina platensis, Boussiba [START_REF] Boussiba | Ammonia uptake in the alkalophilic cyanobacterium Spirulina platensis[END_REF] found that ammonium uptake proceeded at the same rate in light and in darkness. In contrast, our result suggests that the assimilation of carbon and of ammonium occur at different times, like CO 2 and N 2 fixation in Cyanothece sp. ATCC 51142 [START_REF] Stöckel | Global transcriptomic analysis of Cyanothece 51142 reveals robust diurnal oscillation of central metabolic processes[END_REF] and in Trichodesmium sp. strain IMS-101 [START_REF] Finzi-Hart | Fixation and fate of C and N in the cyanobacterium Trichodesmium using nanometer-scale secondary ion mass spectrometry[END_REF]. The preferential assimilation of ammonium during the night could be explained by the ecological conditions that drive the cell biology of Microcystis. This species is known to perform vertical migrations in the water column, and to be located lower down in the water column during the first few hours of the dark period [START_REF] Rabouille | Simulation of carbon reserve dynamics in Microcystis and its influence on vertical migration with Yoyo model[END_REF]. However, while an increase of phosphorus concentrations with depth has been well documented for many lakes, nothing suggests that ammonium concentrations display similar changes.

Altogether, these data suggest that light is not the sole trigger for the transcription of genes involved in photosynthesis and respiration. Indeed their transcription may also be under the control of an endogenous circadian clock, which makes it possible to anticipate their regulation in order to ensure an optimum efficiency of the processes when the changes in environmental conditions arise. The circadian clock is an endogenous biological timing mechanism found in almost all organisms [START_REF] Pittendrigh | Temporal organization: reflections of a Darwinian clock-watcher[END_REF], however Cyanobacteria are the simplest organisms known to exhibit circadian rhythms. The kaiABC genes have been identified as essential genetic components for circadian oscillation in most cyanobacteria [START_REF] Brunner | Transcriptional and post-transcriptional regulation of the circadian clock of cyanobacteria and Neurospora[END_REF]. We did indeed find oscillations in the transcript abundance of the kaiB and kaiC genes, whereas there was no such significant variation for the kaiA gene, a regulator of kaiBC transcription [START_REF] Ishiura | Expression of a gene cluster kaiABC as a circadian feedback process in cyanobacteria[END_REF],

during the light/dark cycle. This is consistent with previous findings showing that the transcript levels of kaiBC displayed high amplitude oscillations, whereas the oscillations of the kaiA mRNA were characterized by a lower amplitude [START_REF] Ishiura | Expression of a gene cluster kaiABC as a circadian feedback process in cyanobacteria[END_REF]. In addition to the kaiABC locus, M. aeruginosa also possesses two paralogs of kaiB, and one of kaiC located at other loci. One of the kaiB paralogs exhibited strong periodicity, with an inverse pattern to that of the kaiABC locus, but it is not known to what extent this kaiB paralog is involved in the circadian clock process. These data do not allow us to pinpoint which specific genes of M. aeruginosa are light induced or circadian controled. A recent paper [START_REF] Elvitigala | Effect of continuous light on diurnal rhythms in Cyanothece sp. ATCC 51142[END_REF] on Cyanothece suggests on the basis of new analytical methods that most of the genes displaying peak expression in the middle of the light period are light responsive, whereas those that are upregulated at the onset of the dark period are circadian controled. Whether this holds true for M. aeruginosa or not remains to be experimentally demonstrated.

With regard to the sigma factors, we found interesting conflicting results in the transcript abundances of the sigD and sigE genes. As reported by Imamura and Asayama [START_REF] Imamura | Sigma factors for cyanobacterial transcription[END_REF],

in cyanobacteria sigD contributes to the light-induced transcription of the psbA2 and psbA3

genes. This is consistent with the finding that psbA2 transcript abundance peaks at 05:30.

More generally, the transcription pattern of the sigD gene, which peaks at 05:30, seems to suggest that the SigD function could be crucial for all genes displaying a light-induced transcription in cyanobacteria [START_REF] Imamura | Sigma factors for cyanobacterial transcription[END_REF]. With regard to the sigB and sigE genes, the fact that the transcripts peak at 18:30 is consistent with the putative function of these sigma factors in the regulation of nitrogen metabolism [START_REF] Imamura | Sigma factors for cyanobacterial transcription[END_REF]. Indeed, as we show in the present study, many genes involved in nitrogen metabolism (amt3, glnA, gltS…) display similar transcript abundance patterns to the sigB and sigE genes. However, our findings on the transcription of mcy genes and of some other genes involved in the production of secondary metabolites, do not confirm the hypothesis that SigD and SigE are involved in regulating these NRPS genes.

From a more general standpoint, this first transcriptomic study of a Microcystis species has shown that the metabolism of this cyanobacteria is clearly compartmented between the light period, during which carbon uptake, photosynthesis, and the reductive pentose phosphate pathway permit the synthesis of glycogen, and the dark period during which glycogen degradation, the oxidative pentose phosphate pathway, the TCA cycle and ammonium uptake permit amino acid biosynthesis. This organization of the central metabolic pathways is very similar to that described in Cyanothece sp. strain ATCC 51142 [START_REF] Stöckel | Global transcriptomic analysis of Cyanothece 51142 reveals robust diurnal oscillation of central metabolic processes[END_REF], a unicellular marine diazotrophic species, suggesting that many cyanobacterial species may share the same central metabolic organization. More specifically, Microcystis is able to produce secondary metabolites, particularly microcystins. The fact that all the genes involved in the biosynthesis of microcystin, and also that of aeruginosin and to some extent that of cyanopeptolin, displayed the same expression pattern with clear upregulation during the light period suggests that these metabolites are essentially produced during the light period and that they may interact with the diurnal part of the central metabolism.

Another question to address is the representativeness of all the results presented in this study compared to the conditions encountered by the Microcystis cells in natural environments. Even if we have chosen to apply a moderate light irradiance to the cultures, to reproduce at the best the natural irradiance conditions, two main differences persist. The first one is the absence of progressive increase and decrease in the light intensity during the transition from dark to light and light to dark, respectively. Such experimental conditions could provoke a stress and explain that for several genes (atpA and atpG, for example) a decrease in transcript abundance was temporarily observed at 06:30. Another difference with natural conditions concerns the continuous exposition of the cells to light for 12 hours, whereas in natural aquatic ecosystems, Microcystis cells are only located at the top of the water column during the first hours of the day, before to be mixed in this water column and to have a discontinuous access to the light during the rest of the day. Interestingly, in spite of these limitations, our study clearly shows that a large part of the carbon and energy metabolism is performed in the first hours of the light period, suggesting the existence of an adaptative process controlling the daily cycle of M. aeruginosa in natural ecosystems.

Legend of Figures

Figure 1 .

 1 Figure 1. Principal component analysis performed on the transcript abundance values

Figure 2 .

 2 Figure 2. Temporal distribution of genes displaying significant changes in their

Figure 3 .

 3 Figure 3. Different patterns in peaking genes classified with regard to their COG

Figure 4 .

 4 Figure 4. Significant variations in the transcript abundance of genes involved in carbon

Figure 5 .

 5 Figure 5. Significant variations in the transcript abundance of genes involved in

Figure 6 .

 6 Figure 6. Significant variations in the transcript abundance of genes involved in the

Figure 7 .

 7 Figure 7. Significant variations in the transcript abundance of genes involved in

Figure 8 .

 8 Figure 8. Significant variations in the transcript abundance of genes involved in the

Acknowledgments

We would like to thank the Génopole of Institut Pasteur for their technical assistance. Monika Ghosh is acknowledged for improving the English version of the article.

Legend of Supplemental materials

Table S1. List of genes with missing data. The second column, gene or gene info, shows either the gene name or brief information about the gene, and the fourth column gives the number of missing data (i.e. the value is one when data are missing data at one sampling time; this value is two when data are missing at two sampling times…). Except for the 70 genes, which have been eliminated by the flip dye consistency check, all the other genes were taken into account in the analyses.

Table S2. Complete data set of our transcriptomic results.

The second column, gene or gene info, shows either the gene name or brief information about the gene. The fold change (Tx/T0) is given for each gene at each sampling time. The last column (STAT) gives the result of the statistical analysis (sig: significant variation of the gene transcript abundance; ns: non significant variation of the gene transcript abundance).

Supporting Information S1

1. Design of our experiment. The figure shows the various sampling times during the light (yellow)/dark (gray) cycle. T1:0:00; T2: 5:30; T3: 6:30; T4: 7:30; T5: 12:00; T6: 18:00; T7: 24:00.

Variations in the transcript


Roland Badeau

Nancy Bertin

Emmanuel Vincent

STABILITY ANALYSIS OF MULTIPLICATIVE UPDATE ALGORITHMS FOR NON-NEGATIVE MATRIX FACTORIZATION

Keywords: Optimization methods, non-negative matrix factorization, multiplicative update algorithms, stability, Lyapunov methods

Multiplicative update algorithms have encountered a great success to solve optimization problems with non-negativity constraints, such as the famous non-negative matrix factorization (NMF) and its many variants. However, despite several years of research on the topic, the of their convergence properties is still to be improved. In this paper, we show that Lyapunov's stability theory provides a very enlightening viewpoint on the problem. We prove the stability of supervised NMF and study the more difficult case of unsupervised NMF. Numerical simulations illustrate those theoretical results, and the convergence speed of NMF multiplicative updates is analyzed.

INTRODUCTION

Non-negative matrix factorization (NMF) is a popular technique allowing the decomposition of two-dimensional non-negative data as a linear combination of meaningful elements in a dictionary [START_REF] Lee | Learning the parts of objects by nonnegative matrix factorization[END_REF]. Given an F × T data matrix V having non-negative entries, NMF consists in computing a rank-K truncated approximation V of matrix V (with K < min(F, T )) as a product V = W H, where both the F ×K matrix W and the K×T matrix H have non-negative entries. The columns of matrix W form the elements of the dictionary, and the rows of H contain the coefficients of the decomposition. NMF can be considered either as a supervised, or as an unsupervised learning tool. In the case of supervised learning [START_REF] Vincent | Musical source separation using time-frequency source priors[END_REF], the dictionary W is estimated from training data in a preprocessing stage, and matrix H only has to be computed given the data in matrix V . In the case of unsupervised learning [START_REF] Lee | Learning the parts of objects by nonnegative matrix factorization[END_REF], both matrices W and H have to be computed given V . Several algorithms have been proposed in order to compute an NMF. The most popular are the multiplicative update algorithms initially proposed by Lee and Seung [START_REF] Lee | Algorithms for non-negative matrix factorization[END_REF]. These algorithms can be applied both to supervised and to unsupervised NMF.

A curious point is that to the best of our knowledge, the convergence properties of multiplicative update algorithms for unsupervised NMF have not been clearly identified. Indeed, Lee and Seung proved that the objective function decreases at each iteration [START_REF] Lee | Algorithms for non-negative matrix factorization[END_REF]. However, this proves neither that the limit value of the objective function is a local minimum, nor that the successive iterates converge to a limit point. In other respects, some numerical examples have been presented in [START_REF] Gonzalez | Accelerating the Lee-Seung algorithm for nonnegative matrix factorization[END_REF], where the Karush-Kuhn-Tucker conditions are not fulfilled after a high (but finite) number of iterations, but this does not contradict possible asymptotic convergence to a local minimum. Finally, since multiplicative updates involve ratios, numerical problems can be encountered if the denominator becomes arbitrarily small. In order to circumvent this problem, it is proposed in [START_REF] Lin | On the convergence of multiplicative update algorithms for nonnegative matrix factorization[END_REF] to add a small positive quantity to the denominator, and it is proved that any accumulation point of the sequence of the iterates computed in this way is a stationary point. However there is no guarantee that such a stationary point is a local minimum, nor that the algorithm converges to this accumulation point.

In this paper, we intend to analyze the convergence properties of NMF multiplicative update algorithms. We apply Lyapunov's first and second methods [START_REF] Lasalle | The stability and control of discrete processes[END_REF] to find some criteria which guarantee the exponential or asymptotic stability of the local minima of the objective function. This approach is applied to prove the stability of supervised NMF multiplicative updates, and we finally show how Lyapunov's first method provides some interesting insights into the convergence properties of unsupervised NMF multiplicative updates. The theoretical results presented in the paper are confirmed by numerical simulations involving both supervised and unsupervised NMF, and the convergence speed of NMF updates is investigated.

MULTIPLICATIVE UPDATE ALGORITHMS AND NMF

Given a matrix V ∈ R F ×T + and an integer K < min(F, T ), NMF consists in computing a reduced-rank approximation of V as a product V = W H, where W ∈ R F ×K + and H ∈ R K×T + . This probcan be formalized as the minimization of an objective function

D(V |W H) = F f =1 T t=1 d v f t K k=1 w f k h kt , (1)
where d a scalar divergence (i.e. a function such that ∀x, y ∈ R+, d(x|y) ≥ 0, and d(x|y) = 0 if and only if y = x). βdivergences [START_REF] Eguchi | Robustifying maximum likelihood estimation[END_REF] are defined for all β ∈ R\{0, 1} as

d β (x|y) = 1 β(β -1) x β + (β -1)y β -βxy β-1 . (2)
The Euclidean distance corresponds to β = 2, and Kullback-Leibler Itakura-Saito divergences are obtained when β → 1 and β → 0, respectively. The generalization of Lee and Seung's multiplicative updates to the β-divergence takes the following form [START_REF] Kompass | A generalized divergence measure for nonnegative matrix factorization[END_REF]:

W ← W ⊗ (V ⊗ (W H) β-2 )H T (W H) β-1 H T η (3)
H ← H ⊗ W T (V ⊗ (W H) β-2 ) W T (W H) β-1 η (4)
where η = 1, the symbol ⊗ and the fraction bar denote entrywise matrix product and division respectively, and the exponentiations must also be understood entrywise. In the case of unsupervised NMF, updates (3) and ( 4) are computed alternately, whereas in the case of supervised NMF, the update (4) only is computed at each iteration, W being kept unchanged. We focus in this paper on the generalization of this approach to an exponent step size η > 0 (possibly different from 1) and we analyze the convergence properties of the multiplicative updates (3) and (4) in section 4. As will be shown in section 5, η actually permits to control the convergence rate, and in particular to outperform the standard case η = 1.

In [START_REF] Kompass | A generalized divergence measure for nonnegative matrix factorization[END_REF], it was proved that if η = 1 and β ∈ [1, 2], then the objective function is non-increasing at each iteration of ( 3) and ( 4). The following proposition proves that ( 3) and ( 4) actually satisfy the same decrease property for all η ∈]0, 1] (i.e. 0 < η ≤ 1).

Proposition 1. Consider the objective function

D(V |W H) de- fined in equation (1), involving the β-divergence (2), with β ∈ [1, 2]. If η ∈]0, 1]
, if all entries in the numerator and denominator in recursion (3) are non-zero and if (W , H) is not a fixed point of (3), then (3) makes the objective function strictly decrease. Similarly, if η ∈]0, 1], if all entries in the numerator and denominator in recursion (4) are non-zero and if (W , H) is not a fixed point of (4), then (4) makes the objective function strictly decrease.

Proposition 1 is proved in [START_REF] Badeau | Stability analysis of multiplicative update algorithms and application to non-negative matrix factorization[END_REF]. Note that this property does not guarantee that the limit value of the criterion is a local minimum, nor that the successive values of W and H converge to a limit point.

STABILITY DEFINITIONS

Let us recall a few definitions in Lyapunov's stability theory of discrete dynamical systems [START_REF] Lasalle | The stability and control of discrete processes[END_REF], which aim to characterize the convergence properties of the general recursion x (p+1) = φ(x (p) ) (where p denotes the iteration index, x (p) is the iterate at iteration p, and function φ is called a mapping), in a neighborhood of a fixed point x (such that φ(x) = x). Notation . denotes any vector norm.

Definition 1 (Lyapunov stability). A fixed point

x ∈ R n + of the recursion x (p+1) = φ(x (p) ), where mapping φ : R n + → R n + is continuous in a neighborhood of x, is said to be Lyapunov stable if ∀ε > 0, ∃δ > 0 such that ∀x (0) ∈ R n + , x (0) -x < δ ⇒ x (p) -x < ε ∀p ∈ N.
property means that initializing the recursion close enough to x guarantees that the subsequent iterates remain in a bounded domain around x. it does not guarantee local convergence. A fixed point which is not Lyapunov stable is called unstable.

Definition 2 (Asymptotic stability). A fixed point

x ∈ R n + of the recursion x (p+1) = φ(x (p) )
, where mapping φ : R n + → R n + is continuous in a neighborhood of x, is said to be asymptotically stable if it is Lyapunov stable and there exists δ > 0 such that ∀x

(0) ∈ R n + , x (0) -x < δ ⇒ x (p) -→ p→+∞ x.
property means that initializing the recursion close enough to x guarantees the convergence to x. A fixed point which is Lyapunov stable but not asymptotically stable is marginally stable.

Definition 3 (Exponential stability). A fixed point x ∈ R n

+ of the recursion x (p+1) = φ(x (p) ), where mapping φ : R n + → R n + is continuous in a neighborhood of x, is said to be exponentially stable if there exists δ, α > 0 and ρ ∈]0, 1[ such that ∀x (0) ∈ R n + ,

x (0) -x < δ ⇒ x (p) -x ≤ α x (0) -x ρ p ∀p ∈ N. (5)
In this case, the minimum value of ρ such that equation ( 5) stands is called the rate of convergence of the recursion.

property ensures a linear speed of convergence; it also implies asymptotic stability. A fixed point which is asymptotically stable, but not exponentially stable, is generally characterized by a sublinear speed of convergence (depending on the initialization). Note that all the stability properties defined above are local, which means that those properties hold in a neighborhood of the fixed point x.

STABILITY ANALYSIS OF NMF ALGORITHMS

In this section, we show how Lyapunov's stability theory applies to the particular problem of NMF introduced in section 2.

Supervised NMF

In supervised NMF, matrix W is kept unchanged and matrix H only is updated by the multiplicative update (4). Below, we remap the entries of W and H into vectors w and h of dimensions F K and KT , respectively. The following propositions are proved in [START_REF] Badeau | Stability analysis of multiplicative update algorithms and application to non-negative matrix factorization[END_REF][START_REF]Supporting document for the paper "stability analysis of multiplicative update algorithms and application to non-negative matrix factorization[END_REF].

Proposition 2. Let w ∈ R F K + . Let h ∈ R KT + be a local minimum of the objective function h → D(w, h) defined in equation (1). Let [∇ 2
h D(w, h)]+ be the matrix extracted from the Hessian matrix of D w.r.t. h, by selecting the rows and columns whose index i is such that ∂D ∂h i (w, h) = 0. Finally, let h → φ h (w, h) denote the mapping defined by equation (4), and ∇φ h T (w, h) denote its Jacobian matrix 1 w.r.t. h. Then ∀β ∈ R, ∃η ⋆ ∈]0, 2] such that

• If η ∈]0, η ⋆ [ and if matrix [∇ 2
h D(w, h)]+ is positive definite, then h is an asymptotically stable fixed point of mapping φ h .

• If η = 0, then h is a marginally stable fixed point.

• If η / ∈ [0, η ⋆ ],
then h is an unstable fixed point.

• h is exponentially stable if and only if η ∈]0, η ⋆ [, matrix [∇ 2 h D(w, h)]+ is positive definite, and ∀i such that hi = 0, ∂D ∂h i (w, h) > 0. In this case, the rate of convergence is equal to the spectral radius ⋆ = ρ ∇φ h T (w, h) < 1.

If moreover β ∈ [1, 2], then η ⋆ = 2.
Proposition 2 proves the exponential or the asymptotic stability of the local minima of (1) with respect to H under mild conditions. The following proposition proves that conversely, exponentially stable fixed points of recursion (4) are local minima of (1).

Proposition 3. Let w ∈ R F K

+ . For all β ∈ R, if η > 0 and if h is an exponentially stable fixed point of mapping φ h , then h is a local minimum of function h → D(w, h).

Unsupervised NMF

Analyzing the stability of the algorithm which alternates multiplicative updates (3) and (4) happens to be particularly difficult. Indeed, it is well known that unsupervised NMF admits several invariances: the product W H is unchanged by replacing matrices W and H by 1 ∀1 ≤ i, j ≤ KT , the (i, j) th coefficient of ∇φ h T (w, h) is

∂φ h j ∂h i .
the non-negative matrices W ′ = W D and H ′ = D -1 H, where D is any diagonal matrix with positive entries. Consequently, the local minima of the objective function are never isolated (any local minimum is reached on a continuum of matrices W ′ and H ′ whose product is equal to W H). For this reason, the minima of the objective function can never be exponentially stable. Below, we remap the entries of w and h into a vector x, of dimension F K + KT . Proposition 4. Let x be a local minimum of the objective function D defined in equation [START_REF] Lee | Learning the parts of objects by nonnegative matrix factorization[END_REF]. Let φ x denote the mapping defined by the composition of (3) and (4), and ∇φ xT (x) denote its Jacobian matrix w.r.t. x. Then ∀β ∈ R, ∃η ⋆ ∈]0, 2] such that

• If η ∈]0, 2[, ρ ∇φ xT (x) ≥ 1. If moreover η ∈]0, η ⋆ [, ρ ∇φ xT (x) = 1.
• If η = 0, then x is a marginally stable fixed point.

• If η / ∈ [0, 2]
, then x is an unstable fixed point.

• ∀η ∈ R, x is not an exponentially stable fixed point.

If moreover β ∈ [1, 2], then η ⋆ = 2.
This proposition is proved in [START_REF] Badeau | Stability analysis of multiplicative update algorithms and application to non-negative matrix factorization[END_REF][START_REF]Supporting document for the paper "stability analysis of multiplicative update algorithms and application to non-negative matrix factorization[END_REF]. If η ∈]0, 2[, Proposition 4 does not permit to conclude about the possible stability of the local minimum. This is because 1 is always an eigenvalue of the Jacobian matrix ∇φ xT (x), which is generally multiple (we suppose that its multiplicity accounts for the invariances of the factorization).

SIMULATION RESULTS

In this section we propose some numerical simulations which illustrate the theoretical results presented in section 4.

Supervised NMF

First, we study the stability of the multiplicative update (4) applied to matrix H, while keeping matrix W unchanged.

Example of sub-linear convergence speed

In this first experiment, the dimensions are F = 3, T = 3 and K = 2. The multiplicative update ( 4) is applied to the Kullback-Leibler divergence (β = 1) with a step size η = 1 (which corresponds to the standard multiplicative update). Let

V =   1 2 3 2 3 4 3 4 5   and W =   1 1 2 1 3 1   . ( 6
)
It can be noticed that V is singular, and that it can be exactly factorized as the product V = W H, where

H = 1 1 1 0 1 2 . ( 7
)
Thus we know that the lowest value of D(V |W H) w.r.t. H is 0, and that this global minimum is reached for H in [START_REF] Eguchi | Robustifying maximum likelihood estimation[END_REF]. This example was chosen so that matrix [∇ 2 h D(w, h)]+ is positive definite, but ∃i/hi = 0 and ∂D ∂h i (w, h) = 0. Thus according to Proposition 2, this global minimum is an asymptotically stable, but not exponentially stable fixed point. Therefore the speed of convergence of the multiplicative update (4) may be sub-linear. Fig. 1 shows the results obtained by initializing (4) with a matrix H having all coefficients equal to 2. As can be noticed in Fig. 1-(a), the objective function D converges to 0 (its global decrease was proven in Proposition 1). Besides, Fig. 1-(b) represents the sequence

1 H (p) -H F - 1 H (p+1) -H F
(where . F denotes the Frobenius norm, H (p) is the matrix computed at iteration p and H is the matrix defined in equation ( 7)) as a solid blue line. This sequence converges to a finite negative value (represented by the dashed red line), which shows that H (p) -H F = O(1/p). As predicted by the theoretical analysis, the convergence speed happens to be sub-linear (at least for the proposed initialization).

Example of linear convergence speed

In this second experiment, all variables are defined as in section 5.1.1, except that the top left coefficient of V is replaced by 0.9. Consequently, this matrix is no longer singular, thus the global minimum of the objective function w.r.t. H cannot be zero. Instead, a local (possibly global) minimum w.r.t. H can be computed by means of multiplicative update (4), initialized as in section 5.1.1. Numerically, we observed that [∇ 2 h D(w, h)]+ is positive definite and ∀i/hi = 0, ∂D ∂h i (w, h) > 0, thus Proposition 2 now proves the exponential stability of this local minimum, with a convergence rate equal to the spectral radius ⋆ . Fig. 2-(a) shows that the objective function D monotonically decreasing. Besides, Fig. 2-(b) represents the sequence H (p+1) -H F H (p) -H F as a solid blue line, and the value ρ ⋆ as a dashed red line. It can be noticed that this sequence converges to ρ ⋆ , which shows that H (p) -H F = O (ρ ⋆p ). As predicted by the theoretical analysis, the convergence speed is linear, with a convergence rate equal to ρ ⋆ .

Optimal step size

In this third experiment, all variables are defined as in section 5.1.2, and we are looking for an optimal step size η. Since β = 1, Proposition 2 proves that the local minimum is exponentially stable if and only if 0 < η < 2. In Fig. 2-(c), the solid red line presents the spectral radius ρ ⋆ as a function of η, for all η ∈] -0.1, 2.1[. This simulation result confirms that ρ ⋆ < 1 if and only if 0 < η < 2, and it shows that there is an optimal value of parameter η, for which the rate of convergence is optimal. In particular, the standard step size η = 1 is not optimal. Additional experiments (not displayed in the Figure ) showed that the objective function D diverges for a value of eta outside the range [0, 2]. Step size η Fig. 2. Example of linear convergence speed

Unsupervised NMF

We now study the case of unsupervised NMF, which alternates multiplicative updates (3) and ( 4) for W and H. In this fourth experiment, all variables are defined as in section 5.1.2, and the algorithm is initialized in the same way. Fig. 3-(a) shows that the objective function D is monotonically decreasing (to a non-zero value). As in section 5.1.3, the solid red line in Fig. 3-(b) represents the spectral radius ρ ∇φ xT (x) as function of the step size η, for all η ∈] -0.1, 2.1[. We note that ρ ∇φ xT (x) > 1 if η / ∈ [0, 2], and ρ ∇φ xT (x) = 1 in the range η =]0, 2[, which confirms that the local minimum is not exponentially stable. Finally, the solid blue line represents the maximum among the magnitudes of the eigenvalues of matrix ∇φ xT (x) which are different from 1 2 This suggests an optimal value η ≈ 1.875, which is again different from the standard step size η = 1. Indeed it can be verified that the lowest value of the objective function D (after 100 iterations) is reached when the algorithm is run with this optimal value of η. Step size η Fig. 3. Unsupervised NMF

CONCLUSIONS

In this paper, we analyzed the convergence properties of NMF multiplicative update algorithms based on β-divergences, where we introduced an exponent step size η. In the case of supervised NMF, we have presented some criteria which guarantee the exponential or asymptotic stability of the multiplicative updates for any η ∈]0, η ⋆ [, where ∀β ∈ R, η ⋆ ∈]0, 2], and if β ∈ [START_REF] Lee | Learning the parts of objects by nonnegative matrix factorization[END_REF][START_REF] Vincent | Musical source separation using time-frequency source priors[END_REF], η ⋆ = 2. We then studied the more complex case of unsupervised NMF. In particular, we proved the unstability of the multiplicative updates if η / ∈ [0, 2]. Finally, the theoretical results presented in the paper were confirmed by numerical simulations involving both supervised and unsupervised NMF. Those simulations showed that the convergence rate depends on the value of η, and that there exists an optimal value of η provides the fastest convergence rate. An algorithmic outlook of this work would be the design of multiplicative update algorithms with an optimal or an adaptive exponent step size. The proposed stability analysis could also be extended to the hybrid case of constrained unsupervised NMF (as suggested in [START_REF]Supporting document for the paper "stability analysis of multiplicative update algorithms and application to non-negative matrix factorization[END_REF]), to non-negative tensor factorization, or to general multiplicative update algorithms, designed for any optimization problem with non-negativity constraints [START_REF] Badeau | Stability analysis of multiplicative update algorithms and application to non-negative matrix factorization[END_REF][START_REF]Supporting document for the paper "stability analysis of multiplicative update algorithms and application to non-negative matrix factorization[END_REF]. Lyapunov's stability theory could also be used to study the stability of other alternating minimization procedures for NMF, such as that presented in [START_REF] Finesso | Approximate nonnegative matrix factorization via alternating minimization[END_REF].

Fig. 1 .

 1 Fig. 1. Example of sub-linear convergence speed

This maximum value discards the eigenvalues equal to 1, which are due to the invariances of the factorization.

The research leading to this paper was supported by the French GIP ANR under contract ANR-06-JCJC-0027-01, Décompositions en Éléments Sonores et Applications Musicales -DESAM, and by the Quaero Programme, funded by OSEO, French State agency for innovation.