Thomas Pollet

D Rémy

Tadonléké

Jean-F Humbert

Comparison of primer sets for the study of Planctomycetes communities in lentic freshwater ecosystems

In search of a primer set that could be used to study Planctomycetes dynamics in lakes and especially via fingerprinting methods, e.g. denaturing gradient gel electrophoresis (DGGE), three existing specific primer sets, developed for marine and soil systems, have been tested on water samples from four freshwater ecosystems. The first primer set (PLA46F/PLA886R) allowed PCR amplification of Planctomycetes sequences in only one of the four ecosystems, whereas the second primer set (PLA40F/P518R) amplified Planctomycetes sequences in all the studied ecosystems but with a low specificity, since sequences belonging to Verrucomicrobiales and Chlamydiales clades were also amplified. Finally, the third primer set (PLA352F/PLA920R) allowed amplification of Planctomycetes sequences in the four ecosystems with a very high specificity. It amplified all known Planctomycetes genera and yielded the highest Operational Taxonomic Unit (OTU) richness and diversity estimates. In silico analyses supported these results. Further experiments comparing PLA352F/PLA920R to PLA46F/P1390R (a primer set generating a longer PCR fragment, also used to study Planctomycetes) yielded very similar results. Our findings suggest that the primer set PLA352F/PLA920R provides good estimates of Planctomycetes richness and diversity compared with other, and can thus be used to study Planctomycetes dynamics in lentic freshwater ecosystems.

Introduction

The order of Planctomycetales belongs to a separate phylum in the domain bacteria on the basis of their 16S rRNA gene sequences [START_REF] Schlesner | Assignment of the genera Planctomyces and Pirella to a new family Planctomycetaceae fam. nov. and description of the order Planctomycetales ord[END_REF][START_REF] Woese | Bacterial evolution[END_REF]. It is now known that these bacteria are ubiquitous and display several specificities in terms of cell structure (e.g. cell compartmentalization), genetics and physiology [START_REF] Fuerst | The Planctomycetes: emerging models for microbial ecology, evolution and cell biology[END_REF][START_REF] Lindsay | Cell compartmentalisation in planctomycetes: novel types of structural organisation for the bacterial cell[END_REF]. These bacteria have been isolated from soils, marine and freshwaters, hot springs, wastewater treatment plants and invertebrate animals (e.g. [START_REF] Schlesner | Pirella marina sp. nov., a budding, peptidoglycan-less bacterium from brackish water[END_REF][START_REF] Schlesner | The development of media suitable for microorganisms morphologically resembling Planctomyces spp., Pirellula spp., and other Planctomycetales from various aquatic habitats using dilute media[END_REF][START_REF] Giovannoni | Isosphaera pallida gen. and comb. nov., a gliding, budding eubacterium from hot springs[END_REF][START_REF] Bornemann | Molecular microbial diversity in soils from Eastern Amazonia: evidence for unusual microorganisms and population shifts associated with deforestation[END_REF][START_REF] Fuerst | Isolation and molecular identification of Planctomycete bacteria from postlarvae of the giant tiger prawn, Penaeus monodon[END_REF][START_REF] Axelrood | Molecular characterization of bacterial diversity from British Columbia forest soils subjected to disturbance[END_REF][START_REF] Pimentel-Elardo | Isolation of Planctomycetes from Aplysina sponges[END_REF][START_REF] Gade | Identification of Planctomycetes with Order-, Genus-, and Strain-specific 16S rRNA-targeted probes[END_REF][START_REF] Köhler | Novel lineages of Planctomycetes densely colonize the alkaline gut of soil-feeding termites (Cubitermes spp)[END_REF]. Since the discovery of the implication of Planctomycetes in the anaerobic ammonia oxidation (anammox) process [START_REF] Strous | Missing lithotroph identified as new Planctomycete[END_REF][START_REF] Kuypers | Anaerobic ammonium oxidation by bacteria in the Black Sea[END_REF], the number of studies dealing with these bacteria has considerably increased over years. However these studies generally concern marine ecosystems, soils, wastewater treatment plants and springs [START_REF] Strous | Missing lithotroph identified as new Planctomycete[END_REF][START_REF] Schmid | 16S-23S rDNA intergenic spacer and 23S rDNA of anaerobic ammonium-oxidizing bacteria: implications for phylogeny and in situ detection[END_REF][START_REF] Chouari | Molecular evidence for novel Planctomycete diversity in a municipal wastewater treatment plant[END_REF]Elshahed et al., 2007;[START_REF] Rich | Anaerobic ammonium oxidation (anammox) in Chesapeake Bay sediments[END_REF][START_REF] Shu | Different Planctomycetes diversity patterns in latitudinal surface seawater of the open sea and in sediment[END_REF]Woebken et al., 2008) and only a few studies have been carried out on these bacteria in lentic freshwater ecosystems. In the latter, Planctomycetes have been found to be involved in anammox [START_REF] Schubert | Anaerobic ammonium oxidation in a tropical freshwater system (Lake Tanganyika)[END_REF] and in change in the quality of dissolved organic matter [START_REF] Tadonléké | Strong coupling between natural Planctomycetes and changes in the quality of dissolved organic matter in freshwater samples[END_REF], two key processes in the functioning of aquatic systems. However while the first cited study focused on particular Plancmotycetes, the second was based on fluorescence in situ hybridization, with no information on Planctomycetes composition. In order to increase our knowledge on the structural and functional diversity of the Planctomycetes communities in freshwater ecosystems, we need to have a set of primers allowing specific amplification of this group and an evaluation of its diversity. We wanted that this primer set generates a small size fragment, which can be used in a fingerprinting approach such as denaturing gradient gel electrophoresis (DGGE), although the latter is not restricted to short fragments. In this goal, we tested three existing primer sets developed in marine and soil systems PLA40F/P518R [START_REF] Derakshani | Novel bacterial lineages at the (sub)division level as detected by signature nucleotide-targeted recovery of 16S rRNA genes from bulk soil and rice roots of flooded rice microcosms[END_REF],PLA46/PLA886R [START_REF] Pyneart | Characterization of an autotrophic nitrogen-removing biofilm from a highly loaded lab-scale rotating biological contactor[END_REF] and PLA352F/PLA920R [START_REF] Mühling | Improved group-specific PCR primers for denaturing gradient gel electrophoresis analysis of the genetic diversity of complex microbial communities[END_REF], and described as Planctomycetes-specific (Table 1). Our study was performed on water samples from four contrasting freshwater ecosystems (two French peri-alpine lakes, Annecy and Bourget, and two reservoirs in Burkina Faso, Bagre and Pouytenga) and was complemented by in silico analyses, which compared the theoretical performance of these primers in the context of the current Ribosomal Database Project (RDP) database. Finally, we compared the results provided by the best primer set with those obtained with a primer set (PLA46F/P1390R) generating a longer fragment, and used by other authors to study Planctomycetes [START_REF] Chouari | Molecular evidence for novel Planctomycete diversity in a municipal wastewater treatment plant[END_REF].

Results and discussion

Efficacy and ability of each primer set to amplify Planctomycetes sequences

The in silico analyses (theoretical approach) using the RDP database showed that the forward and the reverse of the primer set PLA352F/PLA920R matched with a higher number of Planctomycetes sequences (60% and 52% respectively), and also exhibited lower percentages of match outside this bacterial group than the forward and the reverse from the two others primer sets (Table 2). It is, however, important to note that in spite of a good specificity of primers obtained by theoretical analysis, experimental analyses are necessary and can reveal a high degree of non-specificity [START_REF] Morales | Empirical testing of 16S rRNA gene PCR primer pairs reveals variance in target specificity and efficacy not suggested by in silico analysis[END_REF]. Very contrasting results were found in the PCR efficiency and specificity using the first three primer sets tested.

The primer set PLA46F/PLA886R allowed amplification of Planctomycetes sequences (100% of sequences obtained belonged to Planctomycetales order) but only in Lake Bourget. No amplification was obtained with this primer set for samples from the three other ecosystems (Table 3). This primer set has previously been used to estimate the Planctomycetes diversity in marine environments and in digesters of wastewater treatment plants [START_REF] Pyneart | Characterization of an autotrophic nitrogen-removing biofilm from a highly loaded lab-scale rotating biological contactor[END_REF][START_REF] Shu | Different Planctomycetes diversity patterns in latitudinal surface seawater of the open sea and in sediment[END_REF]. Previous studies indicated that a low specificity of the primers and a poor efficacy of the PCR reaction can yield poor amplification [START_REF] Becker | PCR bias ecological analysis: a case study for quantitative Taq nuclease assays in analyses of microbial communities[END_REF], as found with this primer set in our study. Moreover, recent studies [START_REF] Bru | Quantification of the detrimental effect of a single primer-template mismatch by real-time PCR using the 16S rRNA gene as an example[END_REF][START_REF] Wu | Quantitative effects of position and type of single mismatch on single base primer extension[END_REF] have shown that single mismatch in the regions where primers bind onto the template strands can also greatly affect the PCR efficacy. This set of primers might thus be able to detect and amplify only a restricted number of Planctomycetes sequences among those present in freshwater ecosystems. In a recent review, Amann and Fuchs (2008) reported that the coverage of the probe PLA46 was 44% (559 of 1271) of the Planctomycetes 16S rRNA sequences in the SILVA database.

Our in silico analysis in the context of a larger database (RDP) indicates a lower coverage (30%) for the PLA46F as a primer for this bacterial group (Table 2).

Table 3. Analysis of the 16S rDNA diversity (Chao1 estimator, total number of sequences, number of Planctomycetes sequences and OTUs and percentage of Planctomycetes sequences) for each primer set and lake.

Primers

Lakes Total number of sequences

Number of

Planctomycetes

Chao1 estimates (±SD) Sequences OTUs 1. After amplification, fragments were cloned in the pGEM-T system II vector (Promega) and then transformed into JM109-competent cell (Promega) according to the manufacturer's instructions. Ninetysix positive clones (white colonies) obtained from each clone library from each lake and each primer set were randomly selected, verified by PCR using the commercial primers SP6/T7 and finally sequenced (GATC Biotech). Sequences were then edited, aligned with Genedoc [START_REF] Nicholas | Genedoc: a tool for editing and annotating multiple sequence alignments[END_REF] and finally checked for chimera using Bellerophon [START_REF] Huber | Bellerophon: a program to detect chimeric sequences in multiple sequence alignments[END_REF] and the Ribosomal Database Project (RDP, [START_REF] Cole | The Ribosomal Database Project (RDP-II): sequences and tools for high-throughput rRNA analysis[END_REF]. Operational Taxonomic Units (OTUs) were defined on the basis of a ≥ 98% sequence identity. The Chaos-1 and abundance-based coverage estimators of species richness were calculated using the software 'EstimateS' (http://viceroy.eeb.uconn.edu/estimates The second primer set, PLA40F/P518R, allowed to amplify Planctomycetes sequences in the four lakes but the specificity of the amplification varied considerably (Table 3). The proportion of sequences belonging to Planctomycetes ranged from a value as low as 13% in Lake Annecy to 96% in the reservoir Pouytenga, intermediate values being found in Lake Bourget (51%) and in the reservoir Bagre (67%). This primer set amplified members of other bacterial orders, which are phylogenetically close to Planctomycetes, namely Verrucomicrobiales and Parachlamydiales (13% and 9% of the total of sequences respectively) in Lake Annecy, and Chlamydiales (6% of all sequences) in the reservoir Bagre. This lack of Planctomycetes specificity has previously been reported for the primer PLA40F in soils and rice roots [START_REF] Derakshani | Novel bacterial lineages at the (sub)division level as detected by signature nucleotide-targeted recovery of 16S rRNA genes from bulk soil and rice roots of flooded rice microcosms[END_REF]. These authors pointed out that PLA40F revealed no mismatch for almost all Planctomycetes and only a very few mismatches for members of the Verrucomicrobia and some Chlamydia. Our results suggest that the primer set PLA40F/P518R might not be suitable for study of freshwater planktonic Planctomycetes, especially in the framework of fingerprinting approach on these communities.

The third primer set tested (PLA352F/PLA920R) provided the best results in term of specificity and efficiency (Table 3). Indeed, a positive PCR amplification was obtained in all sampling sites and almost all the sequences (> 97%) belonged to the Planctomycetes phylum. The remaining sequences (< 3%) belonged to a group of 'unclassified bacteria'.

Operational Taxonomic Unit (OTU) richness found in the Planctomycetes communities using the three sets of primers

The number of Planctomycetes Operational Taxonomic Units (OTUs) identified in the samples from the four ecosystems with the different primer sets did not differ significantly from the values of the Chao1 estimator (Sign test, P = 0.25, Table 3), showing that we obtained a sufficient number of sequences to evaluate the Planctomycetes OTU richness in our samples. This finding was confirmed by the rarefaction curves which were asymptotic in most cases (not shown). Finally it appeared for the same primer set that there were large variations in the richness index between the samples from the four ecosystems (Table 3). Compared with the two others, the primer set PLA352F/PLA920R allowed to obtain higher numbers of OTUs and higher values of the Chao1 index in almost all cases (Table 3). suggests that the primer set PLA352F/PLA920R was able to amplify a larger set of Planctomycetes sequences than the two others, and that the latter underestimated the richness of the Planctomycetes community. These results were supported by the in silico analyses, which revealed only a few sequences detected by PLA352F/PLA920R outside the group of Planctomycetes in the context of RDP (Table 2). As pointed out in previous works [START_REF] Lueders | Evaluation of PCR amplification bias by terminal restriction fragment length polymorphism analysis of small-subunit rRNA and mcrA genes by using defined template mixtures of methanogenic pure cultures and soil DNA extracts[END_REF], our findings confirmed that it is necessary to perform a careful evaluation of the PCR primers used to assess the diversity in microbial communities. Hence, the primer set PLA352F/PLA920R seemed to provide the best compromise considering the ability to obtain a positive PCR amplification, the specificity of the amplification and the ability to detect a large number of OTUs within the Planctomycetes communities.

Planctomycetes diversity obtained using PLA352F/PLA920R versus PLA46F/P1390R

Based on the above results, further experiments were performed to compare results obtained with PLA352F/PLA920R with those found with the primer set PLA46F/P1390R [START_REF] Chouari | Molecular evidence for novel Planctomycete diversity in a municipal wastewater treatment plant[END_REF]. PLA352F/PLA920R generates a short enough PCR fragments (568 bp) to be used in DGGE [START_REF] Mühling | Improved group-specific PCR primers for denaturing gradient gel electrophoresis analysis of the genetic diversity of complex microbial communities[END_REF], one of the most used fingerprinting approach. PLA46F/P1390R allows amplification of a 1344 bp fragment, and has been used to study Planctomycetes communities in waste water treatment plants [START_REF] Chouari | Molecular evidence for novel Planctomycete diversity in a municipal wastewater treatment plant[END_REF]. The comparison was performed on a sample from Lake Bourget (Bourget 2). Results showed that the two primer sets presented a high specificity to Planctomycetes, as 99% and 98% of sequences obtained with PLA352F/PLA920R and PLA46F/P1390R, respectively, belonged to Planctomycetales. The number of OTUs (5 and 6) and the composition of the Planctomycetes community amplified by these two primer sets were very similar (Table 3). In addition, these two primer sets yielded similar values for the Chao1 (Table 3). This indicated that, despite the smaller size of PCR fragments generated by PLA352F/PLA920R, the latter did not underestimate Planctomycetes diversity, in our samples, compared with PLA46F/P1390R.

Phylogenetic analysis of the sequences obtained with the primers PLA352F/PLA920R

Based on the above results, a phylogenetic analysis was performed with data obtained with PLA352F/PLA920R. To date, seven cultured genera of Planctomycetes (Planctomyces, Gemmata, Isosphaera, Pirellula, Schlesneria, Singulisphaera and Zavarzinella) have been described [START_REF] Staley | The order Planctomycetales and the genera Planctomyces, Pirellula, Gemmata and Isosphaera[END_REF][START_REF] Kulichevskaya | Schlesneria paludicola gen. nov., sp. nov, the first acidophilic member of the order Planctomycetales from sphagnum-dominated boreal wetlands[END_REF][START_REF] Tadonléké | Strong coupling between natural Planctomycetes and changes in the quality of dissolved organic matter in freshwater samples[END_REF]2009) but several additional Candidatus genera (Brocadia, Kuenenia, Scalindua, Annamoxoglobus and Jettenia) playing a key role in the anammox process [START_REF] Strous | Missing lithotroph identified as new Planctomycete[END_REF][START_REF] Kuypers | Anaerobic ammonium oxidation by bacteria in the Black Sea[END_REF][START_REF] Kartal | Candidatus 'Anammoxoglobus propionicus' a new propionate oxidizing species of anaerobic ammonium oxidizing bacteria[END_REF][START_REF] Schmid | Anaerobic ammonium-oxidizing bacteria in marine environments: widespread occurrence but low diversity[END_REF][START_REF] Quan | Diversity of ammoniumoxidizing bacteria in a granular sludge anaerobic ammonium-oxidizing (anammox) reactor[END_REF] have recently been proposed. Among the 437 Planctomycetes sequences obtained with the primer set PLA352F/PLA920R in our four sampling sites, 39 OTUs displaying a ≥ 98% sequence identity were identified. Thirty-eight of these OTUs were distributed in four clusters containing all Planctomycetes genera [cluster I (Pirellula/Blastopirellula), cluster II (Planctomyces/Schlesneria), cluster III (Isosphaera/Singulisphaera) and cluster IV (Gemmata/Zavarzinella)]. This indicated that the primers PLA352F/PLA920R were able to amplify members of all the Planctomycetes genera currently described (Fig. 1). Among these 38 OTUs, 20 (∼53%) belonged to cluster I, which comprises Pirellula/Blastopirellula genera. A high diversity of the Pirellula genus and its dominance within Planctomycetes communities have been reported by previous studies in other freshwater environments, namely wastewater treatment plants and river biofilms [START_REF] Chouari | Molecular evidence for novel Planctomycete diversity in a municipal wastewater treatment plant[END_REF][START_REF] Brümmer | Diversity and seasonal changes of uncultured Planctomycetales in river biofilms[END_REF]. Our results, together with those from these previous studies, suggest that members of this genus are probably well adapted to freshwater ecosystems and might play a particular role in these environments.

One OTU (OTU 28), amplified in the reservoir Bagre, showed a high degree of similarity with sequences belonging to anammox Planctomycetes (Fig. 1). This suggests that the primer set PLA352F/PLA920R is potentially able to detect and amplify freshwater Planctomycetes involved in anaerobic ammonium oxidizing processes. This may be of interest even though obtaining of a better representation of anammox organisms when they are present might require primer specific to these organisms [START_REF] Schmid | Biomarkers for in situ detection of anaerobic ammonium-oxidizing (Anammox) bacteria (Minireview)[END_REF]. The existence of anoxic microzones in detrital aggregates [START_REF] Shanks | Reducing microzones and sulfide production in marine snow[END_REF]) might be a factor allowing the presence of anammox bacteria in oxygenated water columns, as found in this study. All the sequences sharing ≥ 98% sequence identity with sequences from GenBank were closely related to environmental sequences and not to cultivated organisms. In addition, 46% of the OTUs displayed sequence similarity < 98% with sequences from GenBank® database (Table S1), which confirmed that the Planctomycetes from freshwater ecosystems are still poorly known. As previously found for dominant phyla in freshwater communities [START_REF] Humbert | Comparison of the structure and composition of bacterial communities from temperate and tropical freshwater ecosystems[END_REF], the fact that best blast hits were all obtained with sequences from freshwater habitats (lakes, rivers and reservoirs) seems to support the conclusions of Brümmer and colleagues (2004) suggesting that Planctomycetes species from freshwaters are different from those from marine environments.

In conclusion, we propose that at the current state of knowledge, the primer set PLA352F/PLA920R could be used to study pelagic Planctomycetes in freshwater ecosystems specifically via fingerprinting methods.

Table 1 .

 1 Planctomycetes-specific 16S rRNA gene PCR primer sets used in this study and conditions of each PCR reaction.

	Primers

Target group Sequenc es (5′-3′) PCR conditions Referenc es No. of cycle s Denaturation Annealing Elongation Temperatu re (°C) Tim e (min ) Temperatu re (°C) Tim e (min ) Temperatu re (°C)

										Tim
										e
										(min
										)
	1. For DNA extraction and PCR amplification, 250 ml of water sample from each lake was first filtered
		through a 2-µm-pore-size polycarbonate membrane filter (Nuclepore) to eliminate larger eukaryotes.
		The bacterioplankton remaining in the filtrate was then collected, through a gentle filtration, on a 0.2 -
		µm-pore-size polycarbonate membrane filters (Nucleopore). Nucleic acid extraction was carried out as
		described in Dorigo and colleagues (2006). PCR amplifications were performed with the PTC-100 TM
		Thermal Cycler (MJ Research) using three sets of primers (PLA46F/PLA886R, Pyneart et al., 2003;
		PLA40F/P518R, Derakshani et al., 2001; and PLA352F/PLA920R, Mühling et al., 2008). To compare
		the two sets of primers PLA352F/PLA920R and PLA46F/P1390R, PCR amplifications have been
		performed only on water samples from the Lake Bourget (referred to as 'Bourget 2'). For each primer
		set, a range of annealing temperatures was tested to determine empirically the temperature that resulted
		in the most specific PCR product.
			CGG
			CRT
	PLA40F P518R Eubacteria Planctomycetal es	GGA TTA GGC ATG ATT ACC GCG	30	94	0,5 60	1	72	2	Boon et al. (2000); Deraksha ni et al. (2001)
			GCT
			GCT GG
			GGA
			TTA
	PLA46F	Planctomycetal es	GGC ATG
			CAA GTC GCC TTG	30	95	1	56	1	72	2	Pyneart et al. (2003)
	PLA886	Planctomycetal	CGA
	R	es	CCA
			TAC
			TCC C
			GGC
			TGC
	PLA352	Planctomycetal	AGT
	F	es	CGA
			GRA TCT TGT GTG	30	96	1	58	1	74	1	Mühling et al. (2008)
	PLA920	Planctomycetal	AGC
	R	es	CCC
			CGT
			CAA
	PLA46F	Planctomycetal es	GGA TTA GGC	30	94	1	57	1	72	2	Chouari et al. (2003)

Table 1 .

 1 Planctomycetes-specific 16S rRNA gene PCR primer sets used in this study and conditions of each PCR reaction.

				PCR conditions
	Primers Target group	Sequenc es (5′-3′)	No. of cycle s	Denaturation Tim Temperatu e re (°C) (min )	Annealing Tim Temperatu e re (°C) (min )	Elongation Tim Temperatu e re (°C) (min )	Referenc es
		ATG
		CAA
		GTC
		GAC
		GGG
	P1390R Eubacteria	CGG TGT
		GTA
		CAA

Table 2 .

 2 Planctomycetes-specific 16S rRNA gene PCR primers used in this study and their specificity revealed by the in silico analysis.

					Identical	Percentage matches	Percentage
	Primer Target group	Sequence (5′→3′)	Escherichia coli position	matches within the target	within the target group compared with all Planctomycetes	matches within the target group compared with
					group	sequences a	total hits
	 a.
			GGA TTA
	PLA46F Planctomycetales	GGC ATG	46-63	2838	30	77
			CAA GTC
			GCC TTG
	PLA886R Planctomycetales	CGA CCA TAC TCC	886-904	4574	48	95
			C
			CGG CRT
	PLA40F Planctomycetales	GGA TTA	40-57	2446	26	86
			GGC ATG
			ATT ACC
	P518R	Eubacteria	GCG GCT	518-534	1784	18	0.24
			GCT GG
			GGC TGC
	PLA352F Planctomycetales	AGT CGA	350-367	5684	60	91
			GRA TCT
			TGT GTG
	PLA920R Planctomycetales	AGC CCC	920-937	4973	52
			CGT CAA

A total of 9488 Planctomycetes sequences were available in the RDP database.  The theoretical performance of each PCR primer was tested (in silico analyses) using the 'Probe Match' function within the Ribosomal Database Project (RDP)

[START_REF] Cole | The Ribosomal Database Project (RDP-II): sequences and tools for high-throughput rRNA analysis[END_REF]

.

Acknowledgements

This work is part of the fulfilment of the requirements for a PhD (Doctorat d'Université) degree by T.P., who was supported by a French ministry of research doctoral grant to R.D.T. The study was partly supported by funding from the 'Cluster Environnement' of the Rhône-Alpes Region (France) to R.D.T. We are grateful to all the colleagues involved in the sampling campaign in Burkina Faso (P. Cecchi, M. Bouvy and B. Le Berre) and the IRD (Research Unit 167) for funding this campaign. We would like also to thank Ursula Dorigo who has been involved in the sampling and in the DNA extraction for Lakes Bourget (Bourget 1) and Annecy.


Jean-Louis Barrat
email: jean-louis.barrat@ujf-grenoble.fr

David Rodney
email: david.rodney@grenoble-inp.fr

David Rodney Simap

Portable implementation of a quantum thermal bath for molecular dynamics simulations

Keywords:

 proposed that the quantum statistics of vibrations in condensed systems at low temperature could be simulated by running molecular dynamics simulations in the presence of a colored noise with an appropriate power spectral density. In the present contribution, we show how this method can be implemented in a flexible manner and at a low computational cost by synthesizing the corresponding noise 'on the fly'. The proposed algorithm is tested for a simple harmonic chain as well as for a more realistic model of aluminium crystal. The energy and Debye-Waller factor are shown to be in good agreement with those obtained from harmonic approximations based on the phonon spectrum of the systems. The limitations of the method associated with anharmonic effects are also briefly discussed. Some perspectives for disordered materials and heat transfer are considered.

Introduction

Standard constant-energy molecular dynamics (MD) simulations explore the phase space of a system by following the classical equations of motion, therefore producing a classical micro canonical distribution function. Very early it was realized that these equations could be modified in several different manners to produce a canonical phase space distribution at fixed temperature. Among the possible ways to achieve this is the implementation of a "Langevin thermostat" [START_REF] Grest | Molecular-dynamics simulation for polymers in the presence of a heat bath[END_REF] that couples each particle independently to a fictitious thermal bath by introducing in the equation of motion a random force and a friction term related by the classical fluctuation-dissipation theorem. This is illustrated by the following equation for a single degree of freedom X of mass M submitted to an external force F (X):

M dV dt = -M γV + F (X) + 2M γ θ(t) (1)
where V = dX/dt is the velocity, and θ is a white noise verifying:

θ(t)θ(t ′ ) = k B T δ(t -t ′ ). (2)
While this type of thermostat (and various other more sophisticated versions developed by Hoover, Nosé and others, see for example Ref. [START_REF] Frenkel | Understanding Molecular Simulation[END_REF]) provides a convenient control of the temperature, it is intrinsically classical, and in particular the equipartition theorem is strictly verified. Any quadratic degree of freedom will on average correspond to an energy k B T /2. While this is a minor problem in liquid systems at room temperature, in which quantum effects are effectively rather small, the situation is quite different in solids.

The Debye temperature, that signals the onset of quantum effects, is often of the order of a few hundred Kelvins, meaning that quantum effects are indeed important. The most prominent manifestation of this fact is the deviation from the Dulong and Petit law of heat capacity, which is manifest in most solids already at room temperature. Obviously quantum dynamics in solids is well handled by the usual harmonic theory of solids, possibly including anharmonic corrections. Using this theory however requires a full calculation of the normal mode spectrum, which may be cumbersome for large or aperiodic systems. Calculations of transport properties such as heat conductivity also represents a real challenge, as well as static properties in strongly anharmonic systems. While static aspects can be dealt with using path integral methods, dynamical effects can not. Moreover, path integrals represent at low temperatures a rather costly alternative to much simpler classical simulations due to the large number of replicas to include in the path integrals in order to converge to the quantum limit [START_REF] Müser | Simulation of material properties below the debye temperature: a path-integral molecular dynamics case study of quartz[END_REF].

It is therefore highly desirable to develop methods that would allow one to describe correctly quantum effects at a relatively low computational cost, even if the representation is approximate. Recently two very interesting steps in this direction were taken by Dammak et al [START_REF] Dammak | Quantum thermal bath for molecular dynamics simulation[END_REF] and by Ceriotti et al [START_REF] Ceriotti | Nuclear quantum effects in solids using a colored-noise thermostat[END_REF]. The fundamental idea introduced in both papers is that the correct distribution of positions and impulsions for a quantum harmonic oscillator can be reproduced by using an equation similar to Eq. 1 but with a colored noise replacing the white noise θ(t). In Ref. [START_REF] Ceriotti | Nuclear quantum effects in solids using a colored-noise thermostat[END_REF], the friction term was moreover substituted by a Non-Markovian term, thus permitting a greater flexibility in the types of thermostats used. This non-Markovian term, and the corresponding noise, were replaced by a coupling to a finite set of classical degrees of freedom forming the "quantum bath", and the parameters characterizing the evolution of the quantum bath were fitted to obtain the correct energy distribution for a quantum harmonic oscillator. In Ref. [START_REF] Dammak | Quantum thermal bath for molecular dynamics simulation[END_REF], a slightly different approach was used. The friction term γ was maintained Markovian, and the noise was replaced by a colored noise with power spectrum Θ(ω) = hω 1 2 + (exp(hω/k B T ) -1) -1 . In this manner, the correct energy distribution for the harmonic oscillator was directly obtained, and a very satisfactory agreement with low-temperature experimental data was demonstrated. The only drawback of the approach proposed in Ref. [START_REF] Dammak | Quantum thermal bath for molecular dynamics simulation[END_REF] is practical and concerns the colored noise generation. The authors used a method which is classical in the physics community [START_REF] Billah | Numerical method for colored-noise generation and its application to a bistable system[END_REF][START_REF] Lü | Numerical simulation of generalized langevin equation with arbitrary correlated noise[END_REF][START_REF] Rodriguez | Efficient generation of correlated random numbers using chebyshev-optimal magnitude-only iir filters[END_REF], but which requires to perform an inverse Fourier transform of the power spectrum of the noise over the entire time duration of the trajectory. The complete noise signal must therefore be generated and stored in the computer memory before starting the simulation. The corresponding memory requirement, for large numbers of particles or for long trajectories, is important, making the implementation of the method rather cumbersome. Our contribution in the present paper is to show how this problem can be circumvented by using a synthesis of the noise employed in the signal processing community [START_REF] Oppenheim | Discrete-time signal processing: second edition[END_REF][START_REF] Rodriguez | Efficient generation of correlated random numbers using chebyshev-optimal magnitude-only iir filters[END_REF], based on an appropriate linear filter applied to a white noise, thus making the quantum thermal bath method much easier to implement in practice.

The manuscript is organized as follows. Section 2 recalls the basis of the method and discusses some physical issues related to convergence. Section 3 describes the synthesis of the noise using a linear filter. Applications to two selected examples are considered in Section 4.

The quantum thermal bath method

The principle of the quantum thermal bath is readily illustrated on the simple case of an harmonic oscillator. Using a unit mass and frequency ω 0 , the equations of motion are

Ẋ = V (3) V = -ω 2 0 X -γV + 2γ θ(t) (4)
where X and V are the position and velocity of the oscillator, γ the friction and θ(t) the external noise. The Fourier components are easily obtained as

X(ω) = 2γ θ(ω) ω 2 0 -ω 2 + iωγ (5)
and

V (ω) = 2γ iωθ(ω) ω 2 0 -ω 2 + iωγ . ( 6
)
And finally the average energy of the oscillator is given by:

E = dω 2π 1 2 ω 2 0 |X(ω)| 2 + 1 2 |V (ω)| 2 = dω 2π γ ω 2 + ω 2 0 (ω 2 -ω 2 0 ) 2 + ω 2 γ 2 Θ(ω) (7)
where Θ(ω) is the power spectral density (PSD) of the noise. In the following, we will use capital letters and a˜to denote spectral densities associated with a given random noise. For a constant Θ(ω) (white noise), integration of Eq. 7 shows that the energy is constant and independent of γ. If moreover Θ(ω) = k B T as required by the classical fluctuation-dissipation theorem, the energy is just equal to the thermal energy k B T .

For a colored noise with a finite frequency support [-Ω max , Ω max ], Eq. 7 yields, provided Ω max > ω 0 +γ, E ≃ Θ(ω 0 ) with an accuracy that depends on γ (typically the actual value of the energy will average Θ(ω) in the frequency range [ω 0 -γ/2, ω 0 + γ/2]). Hence, by choosing

Θ(ω) = h|ω| 1 2 + 1 exp(h|ω|/k B T ) -1 (8)
and a value of γ small compared to ω 0 , one ensures that the energy of the oscillator is given by the Bose-Einstein formula including the zero point energy. The reasoning is readily extended to a set of linearly coupled harmonic oscillators (since the friction and covariance of the noise are diagonal matrices). In this case, each eigenmode will have an energy equal to that of an harmonic oscillator with the corresponding frequency.

A slight word of caution is in order here. The spectral density given in Eq. 8 does not have a finite support, and moreover it diverges proportionally to ω for high frequencies. The integrand in Eq. 7 thus decays as 1/ω for large ω and the integral diverges at high frequency for any finite value of γ. Clearly this divergence is associated with the fact that the quantum thermal bath includes fluctuations at arbitrary high frequencies, which are picked up by the broadening of the oscillator response function due to the damping γ. These high frequency fluctuations, however, are not expected to be present in any real system: in a crystal, the actual bath of oscillators that a given oscillator couples with is physically limited to frequencies of the order of the Debye frequency. No higher frequency is present in the system. Hence, it appears reasonable to modify Eq. 8 by introducing an upper cutoff frequency Ω max of the order of a few times the highest frequencies observed in the system. We will see in the following that such a cutoff frequency appears most naturally in a numerical implementation of the colored noise. The divergent contribution to the energy in Eq. 7 then behaves as γ ln(Ω max ), and can be made arbitrarily small by choosing a small enough value of γ, once Ω max has been specified. In general, and for a given cutoff frequency, the independence of the results on the damping parameter γ should be checked to ensure that the relevant spectrum of frequencies is properly thermostated.

On the other hand, it is also clear from Eq. 7 that, for a fixed cutoff Ω max , the result will be a decreasing function of the damping γ. In practice, γ should be kept small enough (compared to ω 0 ) that the integrand actually behaves like a δ function, while the equilibration time remains reasonable. We will see in the following that values of γ/ω 0 ∼ 10 -2 provide a good compromise between these various requirements.

Noise synthesis

The numerical implementation of the quantum bath as described above reduces to generating a colored noise with a prescribed PSD, Θ(ω). We use here a classical method of signal processing based on filtering a white noise signal [START_REF] Oppenheim | Discrete-time signal processing: second edition[END_REF]. For a continuous noise, we introduce the filter

H(ω) = Θ(ω) (9)
and its inverse Fourier transform H(t) ( Θ is given in Eq. 8). The noise θ(t) is then obtained by convoluting H(t) with a random white noise r(t) with PSD R(ω) = 1. We have:

θ(t) = ∞ -∞ H(s)r(t -s)ds (10)
The spectral density of the resulting noise is then | H(ω)| 2 R(ω) = Θ(ω), which is the target PSD. This scheme can be implemented in a discrete algorithm suitable for computer calculation. The filter H(ω) is first discretized in 2N f values with steps δω over an interval [-Ω max , Ω max [:

Hk = H(kδω), k = -N f ...N f -1. ( 11
)
The cutoff frequency Ω max can be chosen on physical grounds as discussed above, and will typically be taken to be a few times the maximum eigenmode frequency of the system. A discrete Fourier transform yields H(t) over the interval [-π/δω, π/δω[ with a timestep h = π/Ω max . H(ω) being even and real, so is H(t):

H n = 1 2N f N f -1 k=-N f Hk cos π N f kn . ( 12
)
The convolution is then performed discretely in the range [-π/δω, π/δω[, yielding a discrete colored noise with a time step h:

θ n = θ(nh) = N f -1 m=-N f H m r n-m , (13)
where r m is a white noise drawn from a gaussian or a uniform distribution with a variance √ h in order to yield a unity PSD. In practice, {H m } m=-N f ...N f -1 is first computed using Eqs. 11 and 12 and stored. The noise generator is then initialized by drawing 2N f initial values {r m } m=0...2N f -1 and time is initialized to N f h. Then, at each subsequent time step, with time nh and n ≥ N f , r n-N f is discarded and a new random number r n+N f is drawn and stored. The correlated noise is then obtained by computing the convolution in Eq. 13. The present algorithm has a CPU computational cost close to a classical Langevin thermostat since the only difference is the convolution in Eq. 13. The main difference is the memory storage because the present algorithm requires to store 4N f values per degree of freedom, corresponding to the filter and white noise. This requirement can however be easily managed and is much less than in the method used by Dammak et al [START_REF] Dammak | Quantum thermal bath for molecular dynamics simulation[END_REF], where an inverse Fourier transform is performed to create a signal over the complete duration of the simulation with a discretization that corresponds to the time step used by the MD integrator, i.e. a required storage (number of time steps) × (number of degrees of freedom).

In principle, there is no specific reason to set the time step of the MD integrator (δt) equal to the time step of the noise generator h corresponds to the inverse of the cutoff frequency Ω max ). Taking again the harmonic oscillator as an example, a reasonable time step for integration of the equation of motion is δt = 0.01ω -1 0 . The cutoff frequency for the noise, on the other hand, should be physically limited to a few times ω 0 , let us say for concreteness Ω max = 2ω 0 . Hence the noise time step can be taken to be M times larger than the MD time step, with M = 50 for this specific example. The simplest way to manage the two time steps is to generate the noise according to the procedure outlined above, with a time step h = M δt, and to keep the noise value constant for M integration time steps δt. This amounts to generating a noise with an actual power spectrum

Ξ(Ω) = H(ω) 2 C(ω) 2 , ( 14
)
where

C(ω) = sin(ωh/2) ωh/2 (15)
is the Fourier transform of the square function of amplitude unity and duration h. Obviously the function Ξ(ω) is not the desired PSD Θ(ω), but this problem is easily cured by replacing H(ω) in Eq. 11 by a corrected spectrum

H1 (ω) = H(ω)C(ω) -1 (16)
so that the proper spectral density is generated over the interval [0, Ω max ].

As an illustration, Fig. 1 compares on a specific example the PSD of the noise generated by the above procedure with or without the final correction. Before correction, the intensity of the noise is below Θ and decreases rapidly at high frequencies. The reason is the suppression of the high frequencies due to the noise plateaus between successive updates. On the other hand, after correction, the PSD of the noise matches very well Θ up to the cutoff frequency Ω max . The results of the quantum thermal bath are illustrated below for two simple examples.

Examples

Linear chain of oscillators

We consider first the very simple example of a linear chain of oscillators coupled by Hookean springs between nearest neighbors, with periodic boundary conditions. The frequency of an individual spring is ω 0 . In Fig. 2 the energy obtained from a MD run using a thermostat with Ω max = 2ω 0 and γ = ω 0 /50 is compared to the exact result for this system.

In Fig. 3, we show the influence of the friction parameter γ. As expected, a too large γ results in a broadening of the oscillator response function which does not pick up the correct range of frequency spectrum. As the spectrum has been cut off to a rather low frequency here, the result is an artificial decrease of the energy. However, a value of γ of the order of ω 0 /50 provides both a good approximation to the spectrum and reasonable equilibration time.

Aluminium crystals

We model an aluminum perfect crystal using the many-body embedded atom method potential developed by Ercolessi and Adams [START_REF] Ercolessi | Interatomic potentials from first-principles calculations: the force-matching method[END_REF]. The simulation cell contains a perfect periodic crystal made of 6 × 6 × 6 repeating cells with crystallographic orientation X = [100], Y = [010] and Z = [001]. The cell comprises N=864 atoms. We compute the evolution of the crystal energy and Debye-Waller factor as a function of temperature using the quantum thermal bath thermostat. In both cases, the results of the MD simulations are compared to classical mechanics calculations performed using a classical Langevin thermostat, to path-integral molecular dynamics (PIMD) calculations (see Ref [START_REF] Müser | Simulation of material properties below the debye temperature: a path-integral molecular dynamics case study of quartz[END_REF] for details) and to quantum and classical harmonic approximations

Energy

Fig. 2 (Color online) Reduced energy (in units of hω0) as a function of reduced temperature hω0/kB, for a chain of 50 oscillators. The dots are the result from a numerical calculation using the implementation of the quantum bath described in the text, with Ωmax = 2ω0, γ = 0.02ω0, and δω = 0.02ω0, i.e. N f = 50. The full line is the exact energy. based on the phonon frequency distribution of the crystal, {ω i } i=0;3N -1 . The latter is computed after diagonalization of the Hessian matrix of the crystal.

In the harmonic approximation, the quantum average energy is given by:

E harmo = 3N -4 i=0 Θ(ω i ), (17)
where the 3 translational modes have been omitted and Θ is the harmonic oscillator energy given in Eq. 8. The classical harmonic approximation is obtained by replacing Θ by k B T . The comparison is shown in Fig. 4. The quantum thermal bath allows to recover very well the quantum harmonic approximation at both frictions γ = 0.1 THz and 1 THz and the figure shows the convergence to the classical regime and the gradual slight departure from the harmonic approximation at high temperatures. The highest frequency in the crystal being ω max = 59.7 THz, we conclude again that a ratio γ/ω max in the range 1/50 to 1/500 is sufficient to ensure the convergence of the energy. In this example, N f = 100, meaning that the memory requirement is about 100 times larger than would be the case for a classical molecular dynamics simulation, independently of the duration of the simulation. The present algorithm is therefore much less demanding than the previously proposed implementation [START_REF] Dammak | Quantum thermal bath for molecular dynamics simulation[END_REF], in which memory storage was proportional to the duration of the simulation. The second example is the Debye-Waller factor computed here as the atomic mean-square displacement

u 2 = 1 N N i=1 (r i (t) -r i (0)) 2 .
In the harmonic approximation, the Debye-Waller factor is expressed as:

u 2 = 1 N -1 3N -4 i=0 h 2M ω i coth hω i 2k B T = 1 N -1 3N -4 i=0 Θ(ω i ) M ω 2 i ( 18
)
The result is shown in Fig. 5, where again a good agreement is seen between the QTB MD simulations and the harmonic approximation at low temperatures. It should be noted however that, in contrast with the PIMD results, the QTB simulations tend to slightly overestimate the harmonic approximation, a feature that will be discussed in section 5 below.

Some words of caution

In this section, we point out some differences between this method and the use of a classical Langevin thermostat. With a classical thermostat, the equations of motion can be transformed into a Fokker-Planck equation for the probability distribution of the positions and momenta. This Fokker-Planck equation admits as a stationary solution the classical Gibbs-Boltzmann distribution, independently of the manner in which each specific coordinate is coupled to the thermostat. In particular, with a set of weakly coupled (eg. trough a small anharmonic term in the energy) harmonic oscillators, the coupling of a finite number of oscillators to the thermostat is sufficient to ensure the thermalization of the whole system at the same temperature T .

Energy

Fig. 6 (Color online) Reduced energies (measured in units of hω1) of two coupled quantum oscillators with frequencies ω1 and ω2, as a function of the reduced temperature (measured in units of hω1/kB). The oscillator with the smaller frequency (ω1) is directly coupled to a thermal bath with the same parameters as used in figure 2. The second oscillator has a frequency ω2 = √ 2ω1 and is not coupled to the thermal bath. The coupling between the two oscillators is ensured by a linear (squares) or quartic (dots) term Hc in the energy (see text). In both cases the strength of the coupling is ǫ = 0.05 in reduced units. The points correspond to simulation data, the lines are the exact energies for decoupled oscillators.

In the case of colored noise, the transformation to a Fokker Planck equation is not possible, and there is to our knowledge no general expression available for the stationary phase space distribution. The quantum thermal bath presented above is coupled to all degrees of freedom, and ensures that each independent harmonic mode acquires the correct quantum energy, independently of the coupling between modes. One may however question the influence of such a coupling, and in particular the manner in which a "quantum" mode would transfer energy to other modes that are not directly coupled to the quantum bath. We have performed an exploratory numerical study of this situation by considering two weakly coupled oscillators, X 1 and X 2 , with frequencies ω 1 and ω 2 and identical masses, and coupled trough a weak coupling Hamiltonian H c (X 1 , X 2 ). The first oscillator is coupled to the quantum thermal bath, while the second one evolves according to the classical equations of motion, M Ẍ2 = -M ω 2 2 X 2 -∂ X2 H c . We have considered the case of a linear coupling H c = ǫX 1 X 2 , and the case of a nonlinear coupling

H c = ǫX 2 1 X 2 2 /2.
In both cases, we find (see figure 6) that both oscillators equilibrate, within numerical accuracy, at the desired energy. The equilibration time in the case of a nonlinear coupling is however extremely long, and the error on the energy of the second oscillator remains significant even after 10 7 periods of oscillation.

In the more practical case of a weakly anharmonic system with a large number of modes, such as a crystal at low temperatures, a related question arises. Here each mode is coupled to the quantum bath with a coupling strength γ. However, the equation of motion is also affected by the coupling to other modes through the anharmonic terms of the energy. Classically, such anharmonic couplings give rise to a finite lifetime of the mode, characterized by an effective damping Γ anh . In order for the quantum thermal bath to operate properly, it is desirable that the coupling strength γ is stronger than this anharmonic coupling Γ anh . On the other hand, as discussed previously, a large γ also leads to inaccurate results for a single mode. One therefore has to find a reasonable compromise, in which the coupling is strong enough so that the mode can be thermalized within its lifetime Γ -1

anh , but weak enough to keep the thermalisation accurate. In fact, the slight deviation observed for the Debye-Waller factor in Fig. 5 can be ascribed to the failure in finding a completely satisfactory compromise. By computing the lifetime of the short wavelength phonons using a classical simulation at a temperature that corresponds to the quantum zero point energy, we find indeed that the order of magnitude of the lifetime is of a few picoseconds for the present interatomic potential. This implies that the friction γ should be at least in the Terahertz range, which indeed yields the best results for the Debye-Waller factor. It is found that smaller values of γ yield a repartition of energy among the normal modes of the crystal that tends to a classical distribution, independent of the mode frequency.

Conclusions and perspectives

In this manuscript, we have presented an "on the fly" implementation of a quantum thermal bath for molecular dynamics simulations, suitable for long calculations in large systems. Two simple test cases were considered as validations. The present algorithm is easily portable with limited memory requirements, and therefore opens the possibility of studying complex systems at a relatively modest computational cost.

Direct applications of interest in this context include, for static properties, all cases in which the use of the harmonic approximation requires a numerically demanding matrix diagonalisation. This is the case in particular of amorphous or nanocrystalline systems at low temperatures. Using the quantum thermal bath allows one to access directly thermodynamic quantities from a simple molecular dynamics calculation.

A more involved potential application, which may require further developments, is the investigation of thermal transport properties in nanostructures and nanostructured materials, for which the use of classical approximations based on phonon spectra and lifetimes is not practical. The study presented in section 5 suggests that modes that are only indirectly coupled to the thermal bath also equilibrate correctly at their quantum energy. Therefore it should be possible, as is done in classical heat transfer simulations, to equilibrate a temperature profile by connecting a sample to two quantum thermal baths operating at different temperatures and in different regions of space.

Finally, the fundamental properties of the quantum thermal bath, and in particular of the phase space distribution created by the colored noise in a system which is not strictly harmonic, should also be investigated in the future.

1 (

 1 Fig.

Fig. 3 (Fig. 4 (

 34 Fig.3(Color online) Reduced energy (in units of hω0) for reduced temperatures hω0/kB = 0.2 and hω0/kB = 1.0 for a single oscillator, as a function of the friction coefficient γ/ω0. The point at γ = 0 is the exact energy.

2 )Fig. 5 (

 25 Fig. 5 (Color online) Debye-Waller factor, or average mean square displacement, in an Aluminum crystal as a function of temperature. The same parameters are used here as in Fig. 4.

Acknowledgements We thank Hichem Dammak, Yann Chalopin, and Jean-Jacques Greffet for useful correspondence. We are particularly grateful to Olivier Michel (Gipsa-Lab, Grenoble) who suggested the method used for the synthesis of the noise. The work was initiated during a stay at the Kavli Institute for Theoretical Physics, Santa Barbara. The authors are members of Institut Universitaire de France.


Nicolas Normand

Robin Strand

Pierre Evenou

Aurore Arlicot

Pierre Arlicot

Aurore}

Path-Based Distance with Varying Weights and Neighborhood Sequences

This paper presents a path-based distance where local displacement costs vary both according to the displacement vector and with the travelled distance. The corresponding distance transform algorithm is similar in its form to classical propagation-based algorithms, but the more variable distance increments are either stored in look-up-tables or computed on-the-fly. These distances and distance transform extend neighborhood-sequence distances, chamfer distances and generalized distances based on Minkowski sums. We introduce algorithms to compute, in Z 2

, a translated version of a neighborhood sequence distance map with a limited number of neighbors, both for periodic and aperiodic sequences. A method to recover the centered distance map from the translated one is also introduced. Overall, the distance transform can be computed with minimal delay, without the need to wait for the whole input image before beginning to provide the result image.

Introduction

In [START_REF] Rosenfeld | Sequential operations in digital picture processing[END_REF] discrete distances were introduced along with sequential algorithms to compute the distance transform (DT) of a binary image, where each point is mapped to its distance to the background. These discrete distances are built from adjacency and connected paths (path-based distances): the distance between two points is equal to the cost of the shortest path that joins them. For distance d 4 ("d" in [START_REF] Rosenfeld | Sequential operations in digital picture processing[END_REF]), defined in the square grid Z 2 , each point has four neighbors located at its top, left, bottom and right edges. Similarly, for distance d 8 ("d * " in [START_REF] Rosenfeld | Sequential operations in digital picture processing[END_REF]), each point has four extra diagonally located neighbors. In both cases, d 4 and d 8 , the cost of a path is defined as the number of displacements. These simple distances have been extended in different ways, by changing the neighborhood depending on the travelled distance ( [START_REF] Rosenfeld | Distances functions on digital pictures[END_REF][START_REF] Borgefors | Distance transformations in arbitrary dimensions[END_REF]), by weighting displacements [START_REF] Montanari | A method for obtaining skeletons using a quasi-Euclidean distance[END_REF][START_REF] Borgefors | Distance transformations in arbitrary dimensions[END_REF], or even by a mixed approach of weighted neighborhood sequence paths [START_REF] Strand | Weighted distances based on neighbourhood sequences[END_REF].

Section 2 presents definitions of distances, disks and some properties of nondecreasing integer sequences that will be used later. Section 3 introduces a new generalization of path-based distances where displacement costs vary both on the displacement vector and on the travelled distance. An application is presented in section 4 for the efficient computation of neighborhood-sequence DT in 2D.

Preliminaries

Lambek-Moser inverse of a integer sequence [START_REF] Lambek | Inverse and complementary sequences of natural numbers[END_REF]. Let the function f define a non-decreasing sequences of integers (f (1), f(2), . . . ) For the sake of simplicity, we call f a sequence. The inverse sequence of f , denoted by f † , is a non-decreasing sequence of integers defined by:

f (m) < n ⇔ f † (n) < m . (1)
An interesting property of a sequence f and its inverse f † is that, by adding the rank of each term to these two sequences, we obtain two complementary sequences f (m) + m and f † (n) + n [START_REF] Lambek | Inverse and complementary sequences of natural numbers[END_REF]. This property extends the results given by Ostrowski et al. [START_REF] Ostrowski | Solutions to problem 3173[END_REF] about Beatty sequences [START_REF] Beatty | Problem 3173[END_REF]. From [START_REF] Lambek | Inverse and complementary sequences of natural numbers[END_REF], we deduce that the inverse of the sequence

f (m) = τ m with a scalar τ , is f † (n) = n τ -1 so f (m) + m = (1 + τ )m and f † (n) + n = (1 + 1 τ
)n -1 are two complementary sequences. If τ is irrational, these sequences are Beatty sequences and, for any positive n, (1

+ 1 τ )n -1 is equal to (1 + 1 τ )n as given in [1]. Proposition 1. f † (f (m) + 1) + 1 is the rank of the smallest term greater than m where f increases. Proof. f † (f (m) + 1) + 1 = m ⇔ f † (f (m) + 1) < m f † (f (m) + 1) ≥ m -1 ⇔ f (m ) ≥ f (m) + 1 f (m -1) < f(m) + 1 ⇔ f (m ) > f(m) and f (m -1) ≤ f (m) .
If we extend f with f (0) = 0, and define g by g(0) = 0, g(n + 1) = f † (f (g(n)) + 1) + 1, then f (g(n)) takes, in increasing order, all the values of f , each one appearing once. Definition 1 (Discrete distance). A function d :

Z n ×Z n → N is a translation- invariant distance if the following conditions holds ∀x, y, z ∈ Z n , ∀λ ∈ Z: 1. translation invariance d(x + z, y + z) = d(x, y) , 2. positive definiteness d(x, y) ≥ 0 and d(x, y) = 0 ⇔ x = y , 3. symmetry d(x, y) = d(y, x) ,
In the following sections, we will drop definiteness and symmetry to define "asymmetric pseudo-distances". Definition 2 (Disk). The disk D(p, r) of center p and radius r and the symmetrical disk Ď(p, r) are the sets:

D(p, r) = {q : d(p, q) ≤ r} , Ď(p, r) = {q : d(q, p) ≤ r} . (2)
Table 1. Example of a non-decreasing sequence f and its Lambek-Moser inverse. f is the cumulative sequence of the periodic sequence (1, 2, 0, 3), f † its inverse. f † (f (n) + 1)+1 locates the rank of the next f increase. For instance, f (6) = 9, f † (f (6)+1)+1 = 8 is the rank of appearance of the first value greater than 9, which is 12 in this case. By definition, any disk of negative radius is empty and the disk of radius 0 only contains its center (D(p, 0) = p ).

Definition 3 (Distance transform). The distance transform DT X of the binary image X is a function that maps each point p to its distance from the closest background point:

DT X : Z n → N DT X (p) = min d(q, p) : q ∈ Z n \ X . (3)
Alternatively, since all points at a distance less than DT X (p) to p belong to X ( Ď(p, DT X (p) -1) ⊂ X) and at least one point at a distance to p equal to

DT X (p) is not in X ( Ď(p, DT X (p)) ⊂ X) then: DT X (p) ≥ r ⇔ Ď(p, r -1) ⊂ X . ( 4
)
The DT is usually defined as the distance to the background which is equivalent to the distance from the background by symmetry. The equivalence is lost with asymmetric distances, and this definition better reflects the fact that DT algorithms always propagate paths from the background points.

In this paper, we consider path-based distances, i.e. distance functions that associate to each couple of points (p, q), the minimal cost of a path from p to q. For a simple distance, a path is a sequence of points where the difference between two successive points is a displacement vector taken in a fixed neighborhood N , and the cost (or length) of a path is the number of its displacements. The cost of the path (p 0 , . . . , p n , p n + v) derives from the cost of the path (p 0 , . . . , p n ):

L(p 0 , . . . , p n ) = r ⇒ ∀v ∈ N , L(p 0 , . . . , p n , p n + v) = r + 1 . (5)
Rosenfeld and Pfaltz specifically forbid paths where a point appears more than once [START_REF] Rosenfeld | Sequential operations in digital picture processing[END_REF]. This restriction has no effect on the distance because a path where a point appears more than once can not be minimal. In a similar manner, they exclude the null vector from the neighborhood, forbidding a point to appear several times consecutively. As before, it has no effect on the distance. Notice that, in terms of distance, forbidding a path is equivalent to giving it an infinite cost, so that it can not be minimal. ( 5) can be rewritten as:

L(p 0 , . . . , p n ) = r ⇒ ∀v, L(p 0 , . . . , p n , p n + v) = r + c v ,
where

c v = 1 if v ∈ N ∞ otherwise .
For a NS-distance characterized by the sequence B:

L(p 0 , . . . , p n ) = r ⇒ ∀v, L(p 0 , . . . , p n , p n + v) = r + c B v (r) , (6)
where the displacement cost c B v (r) is 1 for a displacement vector in the neighborhood B(r + 1) and infinite otherwise:

c B v (r) = 1 if v ∈ N B(r+1) ∞ otherwise (7)
For a weighted distance with mask

M = {(v k ; w k ) ∈ Z n × N * } 1≤k≤m
, the distance increment only depends on the displacement vector, but not on the distance already travelled:

L(p 0 , . . . , p n ) = r ⇒ ∀v, L(p 0 , . . . , p n , p n + v) = r + c v , (8)
c v = w if (v; w) ∈ M ∞ otherwise
Briefly, the displacement cost for a vector v and the travelled distance r, is 1 or ∞, independently of r for simple distances, is equal to 1 or ∞ whether v belongs or not to N B(r) for a NS-distance, is in N * ∪ {∞} according to the chamfer mask and independently of r for a weighted distance.

In the following, we propose to use a displacement cost, denoted by c v (r), with values in N * ∪ {∞}, that depends both on the displacement vector v and on the travelled distance r. According to the previous remarks, the cost associated to the null displacement will always be unitary:

∀r ∈ N, c 0 (r) = 1 . (9)
3 Path-based Distance with Varying Weights Definition 4 (Path). We call path from p to q, any finite sequence of points P = (p = p 0 , p 1 , . . . , p n = q) with at least one point, and denote by P(p, q), the set of these paths.

Notice that this definition of a path is not related to any adjacency relation.

The sequence P = (p) is allowed as a path from p to itself. It is distinct from P = (p, p), the path from p to itself with a null displacement.

Definition 5 (Partial and total costs of a path). Let N be a set of vectors containing the null vector 0 and the positive displacement costs c v (with c 0 (r) = 1 and c v ∈N (r) = ∞). The total cost of the path P = (p 0 , p 1 , . . . , p n ) is:

L(P ) = L n (P ) , (10)
where L i (P ) is the partial cost of the path truncated to its i + 1 first points (i.e., to its i first displacements):

L 0 (P ) = L(p 0 ) = 0 , ( 11
) L i+1 (P ) = L(p 0 , . . . , p i+1 ) = L i (P ) + c pip i+1 (L i (P )) . ( 12
)
Definition 6. We use the notation C vk (r) = r + c vk (r). c vk (r) is the relative cost of the displacement v k when the distance travelled so forth is r. C vk (r) represents the partial cost of the path after this displacement (the absolute cost of this displacement):

L i+1 (P ) = L i (P ) + c pip i+1 (L i (P )) = C pip i+1 (L i (P )) . ( 13
)
Definition 7. The pseudo-distance induced by v k , c vk is defined by: d(p, q) = 0 ⇔ p = q d(p, q) = min P ∈P(p,q) L(P ) .

Definition 8. We call minimal relative (resp. absolute) cost of displacement, denoted by ĉ (resp. Ĉ), the quantity ĉv (r) = min c v (s) + sr, ∀s ≥ r (resp.

Ĉv (r) = min C v (s), ∀s ≥ r ).
Proposition 2 (Preservation of cost order by concatenation). Appending the same displacement to existing paths preserves the relation order of their costs. Let P = (p 1 , • • • , p n P ) and Q = (q 1 , • • • , q n Q ) be two paths with costs L(P ) and L(Q), v a vector and P = (p 1 , • • • , p n P , p n P +v), Q = (q 1 , • • • , q n Q , q n Q +v) the extended paths with costs L(P ) and L(Q ) measured with minimal displacement costs. Then:

L(P ) ≤ L(Q) ⇒ L(P ) ≤ L(Q ) .
(14) Proof. From (13), L(P ) = Ĉv (L(P )) and L(Q ) = Ĉv (L(Q)). By definition of Ĉv , s ≤ r ⇒ Ĉv (s) ≤ Ĉv (r), which gives (14). Proposition 3. Let N = v k be a set of vectors and, c v (r), the displacement costs for these vectors. There exists a path P from p to q of cost L(P ) = r measured with costs c v (r) if and only if there exists a path P from p to q of cost L (P ) = r measured with the minimal displacement costs ĉv (r). Proof. Consider the cost of P after i displacements, L i (P ) = L i (p 0 , p 1 , . . . , p i ), we note m 0 = 1, m 0<i≤n = 1 + L i (P ) -L i-1 (P )ĉp i-1 pi (L i-1 (P )) = 1 + c p i-1 pi (L i-1 (P ))ĉp i-1 pi (L i-1 (P )) and M i = i j=0 m i the cumulated sum of m i . Clearly, if L(P ) is finite then each m i is finite and positive because ĉv (r) is less than or equal to c v (r) by construction. Let P be the (finite) path obtained by m i occurrences of each point p i :

P = (p 0 , p 1 . . . p 1 m1 , . . . , p i . . . p i mi , . . . , p n . . . p n mn ) .
We take as an induction hypothesis that the partial cost of P after m i occurrences of p i , L

Mi-1 (P ), is equal to L i (P ). It holds for i = 0 because L M0-1 (P ) = L m0-1 (P ) = L 0 (P ) = 0 = L 0 (P ). If the hypothesis holds for i -1, then the partial cost of P after the first occurrence of p i is L

Mi-1 (P ) = L i-1 (P ) + ĉp i-1 pi (L i-1 (P )), and after m i -1 repeats of p i , equals:

L Mi-1+mi-1 (P ) = L Mi-1 (P ) = L i-1 (P )+ĉ p i-1 pi (L i-1 (P ))+m i -1 = L i-1 (P )+c p i-1 pi (L i-1 (P )) = L i (P )
and the hypothesis is true at rank i. Therefore, for every path of finite cost r measured with L, there exists a path with the same cost measured with L . This is shown in Fig. 2a.

Conversely, let P be a path with finite cost measured by L . We build a path P where each point of P appears m i times consecutively with m i such that m

i -1 + c pip i+1 (L i (P ) + m i -1) = ĉpip i+1 (L i (P ))
. By definition of ĉ, ∀r, ∃s : ĉv (r) = c v (s)+s-r, so m i exists. Let M 0 = 0 and M 0<i≤n =

i-1 j=0 m j , be the cumulated sum of the previous terms of m i . The induction hypothesis is that the partial cost of P , measured with L, at the first occurence of p i , L M i (P ), is equal to L i (P ). It holds for i = 0 with a null partial cost L M 0 (P ) = L 0 (P ) = 0 = L 0 (P ). If the hypothesis holds at rank i, the partial cost of P , after m i -

1 repetitions of p i , if L M i +m i -1 (P ) = L M i (P ) + m i -1 = L i (P ) + m i -1,
and at the first occurence of p i+1 , equals

L i (P )+m i -1+c pip i+1 (L i (P )+m i -1) = L i (P )+ĉ pip i+1 (L i (P )) = L i+1 (
P ) and the hypothesis also holds at rank i + 1. An example of such a path is shown on Fig. 2b.

Corollary 1. Displacement costs c v and ĉv induce the same pseudo-distance.

According to [START_REF] Rosenfeld | Distances functions on digital pictures[END_REF], any path from p to q of cost less than r can be extended with null displacements to reach cost r:

L(p 0 , . . . , p n = q) = s < r ⇒ L(p 0 , . . . , p n = q, . . . , q 1+r-s ) = r (15)
Proposition 4. There exists a path of cost r from p to q if and only if d(p, q) ≤ r.

Proof. If a path of cost r from p to q exists then by definition of the distance, d(p, q) = r if P cost is minimal, d(p, q) < r otherwise. Conversely, if d(p, q) = s then there exists a path of cost s from p to q that, according to (15), can be extended to cost r ≥ s. measured with displacement costs cv . P = (p0, p1, p1, p1, p2, p2), solid lines, is built in such a way that its cost L (P ) measured with minimal displacement costs ĉv , is equal to L(P ) = 8. (b) Given P = (p0, p1, p2), shown with solid lines, has a total cost L (P ) = 7 measured with displacement costs ĉv . P = (p0, p0, p1, p1, p1, p1, p2), dashed lines, is built in such a way that L(P ) = L (P ) = 7.

Corollary 2. For any value of r greater than or equal to d(p, q), there exists a path from p to q which cost is exactly r. The closed disk centered in p with radius r is the set of points for which a path from p of cost equal to r exists: q ∈ D(p, r) ⇔ ∃P ∈ P(p, q), L(P ) = r .

(16)

An iterative construction rule of disks is deduced from (16):

∀r > 0, D(p, r) = v∈N q : ∃P ∈ P(p, qv) and C v (L(P )) = r

= v∈N s : Cv (s)=r D(p + v, s) (17)
4 Minimal Delay Distance Transform

In [START_REF] Wang | An algorithm for a generalized distance transformation based on Minkowski operations[END_REF], Wang and Bertrand, proposed a single scan asymmetric generalized DT based on a neighborhood for which there exists a scanning order such that when a point p in the image is scanned, all neighbors of p have already been scanned (forward scan condition). Then, they extended this result to a sequence where two neighborhoods with forward scan condition are alternated (i.e., B = (1, 2)) [START_REF] Wang | Some sequential algorithms for a generalized distance transformation based on Minkowski operations[END_REF]. In the following we propose a method to compute an asymmetric generalized DT based on any number of neighborhoods having forward scan condition used in an arbitrary order defined by a sequence B, either periodic or not. For our purpose, we will use translated versions of regular NS-distances neighborhoods, in order to meet the forward scan condition for each of them.

The resulting translated distance map can easily be transformed back into a regular, symmetrical, NS-distance map.

Proposition 5. The DT of an image X with the distance induced by the neighborhood N and the displacement costs C v is such that:

DT X (p) = 0 if p ∈ X min Ĉv (DT X (p -v)), v ∈ N * otherwise ( 18
)
where Ĉv represents the minimal absolute displacement costs corresponding to C v (definition 8).

Proof. Case p ∈ X directly results from definitions 3 and 7. Suppose now that p ∈ X so any path from q ∈ X to p has at least one displacement. Prop. 3 states that distances induced by v k , C vk and v k , Ĉvk are equal so we consider the latter cost increments for which prop. 2 holds. According to prop. 2, if P = (q = p 0 , . . . , p n = pv) is a minimal path from q to pv then P = (q = p 0 , . . . , p n , p + v) has a minimal cost -among paths from q to p with second last point pv -equal to Ĉv (L(P )). So Ĉv (DT X (pv)) is the shortest distance from a point q ∈ X to p via p-v. Since all paths which last displacement v does not belong to N have an infinite cost and can not be minimal, (18) holds.

Generalized Distance Transform

When all vectors in N * are directed forward relatively to the scan order, (18) propagates paths from background pixels in a single scan. As a consequence, a generalized DT using any number of neighborhoods N 1 . . . N n , selected by a sequence B, B(i) ∈ [1, n], derives directly from [START_REF] Ostrowski | Solutions to problem 3173[END_REF]18) and minimal costs given by: Ĉv (r) = min s : s > r and v ∈ N B(s) .

let χ v (r) denote the characteristic function of the set N B(r) (i.e., χ v (r) = 1 if v ∈ N B(r) ; 0 otherwise) and χ Σ v (r) its cumulative sum (χ Σ v (r) = s≤r χ v (r)).

Then according to prop. 1:

Ĉv (r) = [χ Σ v ] † (χ Σ v (r) + 1) + 1 . (20)
Algorithm 1 produces a generalized DT using any sequence of neighborhoods (N represents their union) in forward scan condition, using displacement costs given by (20). A similar algorithm was already presented for the decomposition of convex structuring polygons [START_REF] Normand | Convex structuring element decomposition for single scan binary mathematical morphology[END_REF]

Fig. 1 .

 1 Fig.1. Total cost of a path P = (p0, p1, p2). Costs of displacements p0p1, p1p2 and p2p3 depend on the partial costs L0(P ) = 0, L1(P ) = cp 0 p 1 (0) + 0 and L2(P ) = cp 1 p 2 (L1(P )) + L1(P ). The total cost of P is cp 2 p 3 (L2(P )) + L2(P ).

Fig. 2 .

 2 Fig. 2. (a) Given P = (p0, p1, p2), shown with dashed lines, has a total cost L(P ) = 8measured with displacement costs cv . P = (p0, p1, p1, p1, p2, p2), solid lines, is built in such a way that its cost L (P ) measured with minimal displacement costs ĉv , is equal to L(P ) = 8. (b) Given P = (p0, p1, p2), shown with solid lines, has a total cost L (P ) = 7 measured with displacement costs ĉv . P = (p0, p0, p1, p1, p1, p1, p2), dashed lines, is built in such a way that L(P ) = L (P ) = 7.

  .

4. 2

 2 Translated NS-distance transform The sequence of disks for a NS-distance induced by a sequence B is produced by iterative Minkowski sums of neighborhoods: D(p, 0) = p , D(p, r) = D(p, r -1) ⊕ N B(r) .

For each neighborhood N j , we apply a translation vector t j such that the translated neighborhood N j = N j ⊕ t j is in forward scan condition. In a translation preserved scan order, t j translates the first visited point in N j to the origin. Assuming a nD standard raster scan order:

The translated neighborhoods N 1 and N 2 obtained with t 1 = (0, 1) and t 2 = (1, 1) are depicted in Fig. 3a and Fig. 3b. Characteristic functions for 3c-e) are respectively 1 B , 2 B and the constant value 1 resulting in the following minimal displacement costs:

Periodic sequence. When B is a periodic sequence, minimal relative costs ĉv are also periodic sequences. Take the periodic sequence of the octagonal distance B = (1, 2), then 1 B (r

Rate-based sequence. Suppose now that the sequence of neighborhoods is defined as a Beatty sequence (as in [START_REF] Hajdu | Approximating the Euclidean distance using non-periodic neighbourhood sequences[END_REF]):

This allows to compute Ĉ1

v and Ĉ2 v on the fly. For the octagonal distance,

An exemple result of algorithm 1 for the translated octagonal distance (with displacement costs obtained either from sequence B = (1, 2) either from τ = 3 2 ) is shown in Fig. 4b.

Symmetric DT from asymmetric DT

Let t(r), r ∈ N * be a sequence of translation vectors such that the translated disks D (p, r) = D(p + t(r), r) and Ď (p, r) = Ď(pt(r), r) are increasing according to the set inclusion. For a sequence of disks produced by translated neighborhoods defined in (21), the translation vectors are:

In particular, for the 2D case:

DT

X has equivalence with values of DT X :

Consequently:

Knowing DT X (p) and DT X (p+t), we can deduce the values of DT X (p-t(r-1)) for all values of r between DT X (p + t) and DT X (p) for which t(r) = t(r -1) + t, i.e., t = t B(r) . Algorithm 2 recovers the values r of the centered DT by selecting all r in the interval [DT X (p+t j ), DT X (p)] such that B(r) = j. Iterating through values r with B(r) = j is achieved using prop. 1. Values of DT X become available before the whole image is computed. For instance, in a standard raster scan, as soon as line y is processed, all lines of DT X above yr max are fully recovered (where r max denotes the maximal value of DT in that line).

Conclusion

In this paper, a path-based pseudo-distance scheme where displacement costs vary both with the displacement vector and with the travelled distance was presented. This scheme is generic enough to describe neighborhood-sequence distances, weighted distances as well as generalized distances produced by Minkowski sums. It was shown that a set of displacement costs can be provided in a minimal form, where each displacement vector is assigned a non-decreasing sequence of costs, without altering the distance function. These non-decreasing sequences are directly applied in the distance transform algorithm to keep track of the costs of minimal paths from the background. An application to a translated neighborhood-sequence distance transform in a single scan was presented along with a method to recover the proper, centered, distance transform. Combined methods provide partial result with a minimal delay, before the input image is fully processed. Their efficiency can benefit all applications where neighborhoodsequence distances are involved, particularly when pipelined processing architectures are involved, or when the size of objects in the source image is limited.

The pseudo-distance presented here is strongly linked to the properties of non-decreasing integer sequences studied by Lambek and Moser. An implementation in C language is publicly available at http://www.irccyn.ec-nantes. fr/~normand/LUTBasedNSDistanceTransform.


Moussa Chnani
email: moussa.chnani@utbm.fr

Hattab Maker

Denis Candusso

Marie-Cécile Péra

Daniel Hissel

ELECTRICAL ANALOGY MODELLING OF PEFC SYSTEM FED BY A COMPRESSOR

The PEFC generator for automotive application is expected to have a low cost, a low weight and a low size, to compete with more and more efficient combustion engine. To reach this aim, the complete system has to be taken into account, not only the stack itself, but also the fluid circuit ancillaries. A model is developed with the aim in view to be integrated in the simulation of an electrical vehicle power train. As many components have to be modelled, a macroscopic approach has been chosen.

A general scheme of the system is proposed, which structure is representative of an embedded generator, i.e. few sensors, few actuators. Gaseous hydrogen is stored in a tank. The anode output is closed by a solenoid valve. It is opened when the fuel cell voltage reaches a minimum threshold, allowing flushes of the channels. Air is provided by a compressor which flow is controlled by the motor speed.

The modelling of the fuel cell electrical response is developed, based on semi-empirical approach. The decrease of the output voltage which can be attributed to the anodic dead end operation is taken into account. The fluidic behaviour of the gas circuits has been dealt through an electrical analogy to facilitate the implementation in a usual electrical engineering software (Matlab/Simulink ® ). Each component of volume V and fluidic resistance R f is represented by a RC cell. In a formal approach, the flow is related to the current and the pressure is related to the voltage. The model is validated with experimentations carried out on a low power test bench (100W). Then, it is simulated on a load cycle, compatible with vehicle application dynamics.

Introduction

Polymer Electrolyte Fuel Cell is one of the promising technologies under development for embedded power supply. The PEFC generator for automotive application is expected to have a low cost, a low weight and a low size, to compete with more and more efficient combustion engine. To reach this aim, the complete system has to be taken into account, not only the stack itself, but also the fluid circuit ancillaries. The behaviour of the last can have a great impact on the global performance. Moreover, the generator architecture should allow to minimise the sensor number and to optimise the reactant consumption.

Simulation is a binding milestone to fulfil these requirements. As many components have to be modelled, a macroscopic approach has been chosen. The model is developed with the aim in view to be integrated in the simulation of an electrical vehicle power train.

In the second section, a general scheme of the system is proposed, which structure could be implemented on board. The modelling of the electrical response is developed, based on a semi-empirical approach. The hydraulic behaviour of the gas circuits has been dealt through an electrical analogy to facilitate the implementation in usual electrical engineering software (Matlab/Simulink ® has been used). In the third section, the model is validated with experimentations carried out on a low power test bench, under constant current and also on a load cycle, compatible with vehicle application.

Component modelling

System structure

An embedded fuel cell system has to fulfil several requirements as compactness, low cost, reliability and durability. Then, at least, the system structure should low the number of actuators and sensors, as much as possible, and guarantee anyway an acceptable performance level of the generator. Figure 1 is an example of such an architecture. Gaseous hydrogen is stored in a tank. A mechanical valve limits the pressure up to 1500mbars absolute. The anode output is closed by a solenoid valve. It is opened when the fuel cell voltage reaches a minimum threshold, allowing flushes of the channels. Air is provided by a compressor. The air is saturated by a humidifier. The laboratory test bench [START_REF] Harel | Experimental investigation on PEFC : a test Bench Design[END_REF] allows us to carry out measurements useful to characterise such a structure with reservations concerning the compressor, which will be developed in the section 3.

Electrical response modelling

The electrical response of the fuel cell is modelled by the following equation [START_REF] Maker | Modélisation d'une pile à combustible de type PEM[END_REF] :

rI I cT P bT P aT U U H air pol - + + + = ln ln ln 2 0 (1)
With: P air : the inlet air pressure P H2 : the inlet hydrogen pressure I : the fuel cell current U pol : the fuel cell voltage r : the fuel cell apparent resistance U 0 , a, b, c are parameters to be determined This expression is used to determine the voltage at the initial time and at the end of the flush. It is defined from a minimal current (because of the logarithmic curve term) below which the voltage is fixed at the value of U0. When the stack is operated in an anodic dead end mode, the stack voltage decreases with time, even if all the operating parameters and the current are kept constant. This phenomenon can be attributed to the accumulation of water and nitrogen, which have migrated from the cathode. This voltage drop is modelled by an exponent function of time and current u(t, I) (2). [3] 2), the voltage of the FC can be written in the following expression :

) ). ( exp( ) ( ) ( ) , ( t I f I e I d I t u ∧ ∧ ∧ - ⋅ + = δ (2)
) (t u U U pol FC δ ⋅ =
(3)

Test bench modelling by electrical analogy

The fluidic behaviour of the gas lines was modelled using an electric analogy, for the integration in Matlab/Simulink®. Each component (pipe, valve, compressor...) of volume V and fluidic resistance R f is represented by a RC cell (figure 2).

P inlet Volume V i 1 i C C U 1 U 2 i 2 R P outlet P inlet -P outlet Flow inlet Flow outlet P inlet Volume V i 1 i C C U 1 U 2 i 2 R P outlet P inlet -P outlet Flow inlet Flow outlet
Capacity ∫ ⋅ = dt i C 1 U ∫ ⋅ = ⋅ = dt F V RT n V RT P

Table1

equations used in the electrical analogy Such an electric representation allows applying methods of electricity and automation. The evolutions of the outlet air flow and air pressures have the typical look of a RC circuit response [START_REF] Moussa | Modelling of a PEFC Generator[END_REF].

Model identification and validation of FC generator

Structure of test bench

The parameters of the model are identified using experimental measurements on a low power test bench (<1kW) developed in L2ES laboratory [START_REF] François | Design of a test bench for a low power PEMFC[END_REF].

The figure 5 shows the simplified structure of the laboratory test bench : The figures 6 and 7 show the two models adopted for the gas lines (air and hydrogen). The identification of each element of the air and hydrogen lines is based on analytical calculations in static and dynamic modes. In statics, resistances are calculated, as in electricity.

In dynamics, by identifying on tests the time-constant RC (and/or also by knowing the volume of each element), the capacities C can be calculated. Two types of current profile were used to identify the model parameters and to validate the results of simulation: a profile of constant current (I=40A) and a standard profile of transport (J227).

Constant current (I=40A)

A load current of 40A was applied to a 3 cell stack, with a current density of 0.4A/cm2. The hydrogen line operates in dead end mode and using the valve placed at the FC outlet, periodic flushes are carried out.

Time [s]

Fuel cell voltage [V]

Experimental voltage Simulated voltage Figure 8 experimental and simulation fuel cell voltages

The voltage change with time is well reproduced by the exponential function (2). In the simulation, flushes have not been started on voltage threshold. The reason why is that in the first phase of the experiment, the stack voltage curve is not repeatable from one flush to an other as long as the membrane hydration is not stabilised, so a time delay between simulation and experiment cannot be avoided. Then the simulation flush starts have been synchronised on the experimental times. Figure 12 zoom -experimental and simulation hydrogen pressures at stack inlet

The inlet hydrogen pressure is constant in the normal operating mode as the flow is constant; it drops when the solenoid valve is opened, and then recovers after the valve shutting. This trend can be seen well on the simulated curve; however there is an error on the pressure level reached during the flush. It can be easily explained as fluidic resistances are taken constant in the model despite the possible strong non-linearity of such a phenomenon. As during the flush, the flow is much higher than during the constant operating point, the linear model resistance cannot take into account the pressure drop anymore. It could have been possible to introduce a non linear resistance but the final accuracy gain is not high enough to justify the complexity raise of the model. The air flow circuit is constant and allows the tuning of fluidic resistances. The air pressure upstream is maintained constant with 1350mbars. Using the flowregulator, the air flow was controlled with 10.1Nl/min.

Transportation current profile (J227)

A standard speed-time profile for personal vehicle (J227) has been used to define a current profile for our 3 cell stack. It has been calculated on the basis of a software developed by the Laboratory on Transportation and Environment (INRETS LTE).

In order to run the stack in the best conditions, the air flow is kept constant at a value corresponding to the maximum current value.

Figure16

experimental and simulation hydrogen pressures at stack inlet (J227 current profile)

Motor-compressor modelling

The air supply circuit based on a motor-compressor group is the most power hungry. It represents up to 80% of the overall ancillary consumption. Consequently, in order to perform both the energy optimization and the energy control of an embedded fuel cell system, a simple model of the motor-compressor has been developed. The choice of each compressor depends on its outlet pressure and the range of power of the fuel cell. A pallet compressor has been chosen (simple construction and lower cost). A permanent magnet synchronous motor is used for the drive of the pallet compressor.

The permanent magnet synchronous motor

The speed regulation of the motor drive is done by varying the I q with I d =0 [START_REF] Mestan | Contribution à l'optimisation énergétique d'un système pile à combustible embarqué[END_REF].

The following table give the equations of the permanent magnet synchronous motor :

Driving torque

Φ ⋅ ⋅ ⋅ = q m I p C 2 3 Load torque ω c r P C = Mechanical equation • • • ⋅ + ⋅ = - θ θ f J C C r m
Motor speed

• = = θ θ ω dt d Electrical power ω ⋅ ⋅ = q phi elec I K P

Tableau 1 mechanical equations of permanent magnet synchronous motor

The inputs of this model are the current I q and the compressor load torque. The outputs of this model are the consumed electrical power, the mechanical power and the motorcompressor group rotational speed.

Compressor modeling

In the work presented in [START_REF] Mestan | Contribution à l'optimisation énergétique d'un système pile à combustible embarqué[END_REF], an accurate model of compressor has been developed, for the study of optimized control laws. It has been validated on an experimental bench. This modelling was used as reference for the identification of our simple model. The compressor being element of the chain, it is modelled in the same manner as the other components of the system (on the principle of a cell RC) by a simple first order transfer function.

Results of simulation of the fuel cell generator

A-Constant current (I=40A)

The simulation of the fuel cell generator (motor-compressor /stack/ hydrogen supply) is based on the same tests as in the case of the test bench structure (without compressor). Figures 1718show the results simulation with constant current solicitation : The strategy of air supply applied in the case of the test bench (constant air flow) is not viable for a real application (the fuel cell does not even provide the sufficient power to feed its own auxiliaries). In this simulation, the motor speed is controlled in order to have a variable flow and to keep a constant air stoichiometry. The motor-compressor consumes 16.5% of the power of the fuel cell, by neglecting the consumption of the other ancillaries.

C-Electrical efficiency of the generator (stack and motor-compressor)

The hydrogen flow is often measured in standard litres per minute (slpm). If this is the case the electrical efficiency should be derived as follows :

.

2 0 sup , ) ( )) ( ( g O H H n P I U f plied hydrogen compressor motor elec stack stack el ∆ × - × = - η (4) . 2 0 273 , ) ( )) ( ( 60 1000 g O H H T R p V P I n U f hydrogen slpm compressor motor elec stack cell cell el ∆ × × × - × × = - η (5)
Where p=101325Pa, R=8.3145(J/mol K), T 273 =273K, The average global efficiency of the generator is approximately 45%.

For transport current profile (J227), the instantaneous global efficiency is periodic and for a null current the global efficiency cannot be calculated, so figure 23 shows one period of generator efficiency with fuel cell current : The average global efficiency of the generator is approximately 60%. This is an upper limit as the consumption of the other ancillaries (mainly the humidifier) has not been taken into account. Furthermore, the average current is rather low so that hydrogen flushes have not been performed.

Conclusion

PEFC generator is a promising technology for the future of electrical vehicles. Simulation is a binding milestone to evaluate the potential interest for propulsion or auxiliary power unit. The whole system has to be taken into account, that includes the stack and the ancillaries needed to manage the fluid feeding of the stack. The proposed model is based on a homogeneous approach of the system components, from the fluidic point of view. As it is based on an electrical analogy, it can be easily implemented in electric engineering softwares, particularly if the whole power train (converter, motorisation and electrical buffers) has to be considered. The model parameters have been tuned around a nominal operating point and linearity has been assumed around this point. It is clear that it reduces its validity area, however simplicity has been favoured and it is very efficient in the first step of the system sizing. Further improvements will be investigated, concerning mainly two aspects: first a better model of polarisation resistance, and the modelling of the humidifier consumption.

Figure 1 diagram

 1 Figure 1 diagram of the setup, consisting of air and hydrogen supplies

∧I:

  Average value of the fuel cell current. d, e et f are parameters to be determined. According to the relation (1) and (

Figure 2 principle

 2 Figure 2 principle of the analogy used between pipe and electric RC cell

Figures 3 and 4

 4 Figure 3 experimental air flows at stack inlet and outlet

  Figure 5 structure of test bench

  Figure 6 model adopted for the air line

  Figure 9 experimental and simulation hydrogen flows at stack inlet

  Figure13fuel cell current (J227 current profile)

  Figure15experimental and simulation hydrogen flows at stack inlet (J227 current profile)

  Figure 17experimental and simulation fuel cell voltages

  Figure20instantaneous electrical power of the fuel cell (J227 current profile)

  Figure 21instantaneous global efficiency of generator (I=40A)

  Figure 23instantaneous global efficiency of the generator (J227 current profile)

5

References