Ian M Carr
email: i.m.carr@leeds.ac.uk

Nick Camm

Graham R Taylor

Ruth Charlton

Sian Ellard

Eamonn G Sheridan

Alexander F Markham

David T Bonthron

GeneScreen: a program for high-throughput mutation detection in DNA sequence electropherograms

Keywords: mutation detection, software, DNA sequencing Word count: 3729 2

come   L'archive ouverte pluridisciplinaire

INTRODUCTION

DNA sequence analysis is a central technology in many branches of biology, perhaps most critically where the relationship between genotype and phenotype is under scrutiny. The massive expansion in the catalogue of human disease-associated genetic variation over the past three decades has been parallelled by the development of diverse methods for detection of sequence variants, particularly approaches that are easily scalable for screening large numbers of DNA samples. Landmark techniques for detecting unknown sequence variants (typically in PCR-generated genomic templates) include chemical heteroduplex cleavage [START_REF] Dahl | Detection and localization of base changes in RNA using a chemical cleavage method[END_REF] , denaturing gradient gel electrophoresis 2 3 , single-stranded conformation polymorphism [START_REF] Boothroyd | Single base mutation in the hormone binding domain of the thyroid hormone receptor beta gene in generalised thyroid hormone resistance demonstrated by single stranded conformation polymorphism analysis[END_REF] and high-resolution melting curve analysis [START_REF] Ririe | Product differentiation by analysis of DNA melting curves during the polymerase chain reaction[END_REF] . The driver for the development of most such methods has been the labour-intensiveness of complete DNA sequencing (earlier by manual and more recently by semi-automated versions of the Sanger approach). Although some of these mutation screening methods can offer high sensitivity for detection of unknown variants, this may demand careful optimization of experimental conditions to suit individual amplicons. Consequently, as modern versions of the Sanger sequencing method have improved in speed and quality 6 7 , sequencing has persisted as the commonest method of choice for mutation detection. Its near-universal applicability under standard conditions makes it particularly attractive for laboratories analysing a wide range of rare single-gene disorders (for whom the optimization of other methods for analysis of diverse genomic targets can present major challenges). Despite the current rapid development of clonal sequencing technologies, high-quality Sanger sequencing of medium-sized (100-1000 bp) PCR-amplified templates is likely to remain a benchmark technology for mutation detection and verification for the foreseeable future.

A major strength of Sanger sequencing is that it provides information on each nucleotide sequenced over a read length of up to 1 kb. Currently, its main limitation in many laboratories is that the large quantity of sequence data generated by modern instruments is both tedious and time-consuming to search for variant positions.

Analysis of electropherograms is also highly prone to human error. For large sequencing projects automated analysis of sequence electropherograms has been an absolute requirement. Techniques arising out of the Human Genome Project tended to focus on the detection of sequence variants (including sequencing errors) in sequences of cloned DNA, with the high quality data then used to construct sequence contigs 8 9 . This high quality data allows detection of polymorphisms through the presence of variant nucleotides at the same position in multiple traces [START_REF] Kwok | Comparative analysis of human DNA variations by fluorescence-based sequencing of PCR products[END_REF] . Since individual sequences originate from clones, this approach does not require computer programs to deal with heterozygous trace positions (which in this setting would indicate sequence artefact, register a low quality score and so be omitted from subsequent analysis).

However, as the focus widened to include SNP detection and resequencing projects, software was developed that could genotype heterozygous positions present in PCR products. 10 11 Such applications identify heterozygous positions by the presence of two overlapping peaks of similar height or by a direct comparison to a reference trace file [START_REF] Crowe | SeqDoC: rapid SNP and mutation detection by direct comparison of DNA sequence chromatograms[END_REF] . The presence of a SNP can be verified by the identification within the study population of samples homozygous for either allele.

While these approaches work well for the detection of SNPs, additional problems complicate the discovery of rare pathological variants. For example, for de novo disease-causing mutations, and even for typical recessive diseases that show extreme allelic heterogeneity, it is unlikely to be possible to observe a second example of a heterozygous mutation among a cohort of samples; homozygosity for the rare allele will also not be observable. Consequently, there will be no supporting evidence to validate the variant within the study cohort. Also, while many heterozygous positions show two overlapping electropherogram peaks of similar heights, this is not always true; heterozygous positions where the allelic peaks have different migration rates and/or peak heights are quite common. When screening for SNPs to describe genetic variability within a population, some level of false negative results is acceptable. In contrast, when searching for pathological changes, it is more important to minimize the false negative rate, even at the cost of a false positive rate requiring assessment by a human operator. Consequently, while it is possible to create a highly automated process for SNP detection, most diagnostic and research laboratories searching for rare pathogenic variants have not felt able to rely entirely on computer programs to automate their detection. We have therefore focussed our attention on maximizing the efficiency with which a human operator, with the aid of the computer, can detect mutations within a large number of electropherograms and assess their likely significance. The application we have developed, GeneScreen, automates the process of detecting possible sequence variants and then, critically, presents the data in a manner that allows a user very quickly to categorize and annotate variants as true or false positives. The program then allows the easy export of the annotated variants for reporting and cataloguing. GeneScreen's focus on a highly efficient user interface, and its integrated capability for annotating rare variants in standard nomenclature on genomic, cDNA or protein reference sequences, are points of difference from previous applications, such as inSNP [START_REF] Manaster | InSNP: a tool for automated detection and visualization of SNPs and InDels[END_REF] , varDetect [START_REF] Ngamphiw | VarDetect: a nucleotide sequence variation exploratory tool[END_REF] and novoSNP [START_REF] Weckx | novoSNP, a novel computational tool for sequence variation discovery[END_REF] , all of which offer various capabilities for mutation detection and annotation, but have not been widely adopted by genetic diagnostic laboratories.

MATERIALS AND METHODS

Programming

Programming was done using Microsoft Visual Studio 2005. GeneScreen runs on Microsoft Windows operating systems, and requires the .NET framework 2.0. Freely available downloads and detailed user guides are at http://dna.leeds.ac.uk/genescreen.

Sequencing data

Sequence electropherograms were generated on Applied Biosystems (Foster City, CA, USA) capillary sequencing instruments in the NHS Regional DNA Laboratories in Leeds and Exeter and the CRUK Genomic Services Mutation Detection Facility (Leeds), either as part of internal quality control of diagnostic sequencing, or for external research projects. (For ABCC8, the patients and sequences were a subset of those previously described [START_REF] Ellard | Semi-automated unidirectional sequence analysis for mutation detection in a clinical diagnostic setting[END_REF] .) Only sequence variants in exons or within 20 bp of a splice site were annotated. According to the priorities of each sequencing project, only selected exons were sequenced in some genes, while in others, larger exons were analysed as a number of overlapping sequences. For example, only protein-coding exons were analysed; the alternatively spliced exon 4 from BRCA1 was not sequenced; BRCA1 exon 10 and BRCA2 exons 10, 11 and 27 were amplified as several shorter products. The TP53 sequences were derived from archived tumour samples, while the ABCC8, BRCA1, BRCA2 and GCK DNA were from patients suspected of having germline mutations, based on the presentation of disease either in the subject under analysis or in a close family member.

For a random subset of the ABCC8, GCK and TP53 sequence files, the ABI Sequencing Analysis v5.2 program was used to create phd.1 files (containing Phredlike quality scores). These files were then processed using the program "QVregion.exe" (http://dna.leeds.ac.uk/qvregion/) in order to obtain the quality scores for the ten worst nucleotides, over the region screened for sequence variants (as defined above-exons + 20 bp of flanking intron).

Data analysis workflow

In outline, GeneScreen detects sequence variants by following a multistep process for each ABI file. Initially the fluorescence data from each ABI file is extracted and used to deduce the nucleotide sequence. This is then aligned to the reference sequence enabling any sequence variants to be identified. Finally, any sequence variants identified are annotated at the genomic, cDNA and protein sequence level where appropriate. It is important to note that Genescreen base-calls each trace independently of the reference or other sequences in the dataset and does not attempt to correct sequencing artefacts. The data analysis workflow is described in greater detail in the Supplementary Data.

RESULTS

Method of use

A full illustrated description of GeneScreen usage is at http://dna.leeds.ac.uk/genescreen/guide/. A brief outline of certain aspects follows.

The initial tabbed user interface of GeneScreen (Figure 1) allows access to each of several sequential steps required for scrutiny of a large batch of test sequences.

Reference file creation

To identify the positions of sequence variants, GeneScreen requires a reference file.

This can be either a plain text sequence, or a special gene descriptor reference file created by GeneScreen itself, using as input cDNA and genomic sequences plus an optional tab-delimited list of the names and sequences of the sequencing primers. The gene descriptor file specifies the position of each exon, relative to the start codon and the beginning of the genomic sequence. This information is then used by GeneScreen to annotate sequence variants according to Human Genome Variation Society nomenclature (http://www.hgvs.org/index.html).

Identification and display of sequence variants

Once a reference file is available, a directory containing a large batch of electropherograms can be loaded and analyzed (taking ~15 s per 96 traces aligned to a 250-bp reference sequence on a typical 2.4 GHz desktop processor). GeneScreen can be configured to reject sequences containing more than a pre-specified number of sequence variants, and the minimum peak height ratio for a heterozygous call is also adjustable. GeneScreen then displays the collated information on all detected sequence variants in the form of a data grid (Figure 2

Exporting and archiving sequence variants

Selected sequence variants can be saved as a tab-delimited plain text (with or without linked images of the variant), a LOVD-compatible data import file, or as a web page that includes images of the exported sequence variants (Figure 3). (Note that GeneScreen's LOVD-compatible file is designed to match standard LOVD2 installations.) If a different layout is required, the form can easily be edited using a spreadsheet application. To extend the utility of the web page output format, new variants can be incrementally added to an existing page, such that all the variants are ordered by gene name, position in the genomic reference sequence and patient ID (irrespective of the order in which they were originally selected for export).

GeneScreen also has a web page editing function that allows the variants in the web page to be edited, updated or deleted. As a result of these features, the web page serves the dual functions of a data presentation file and a database of selected sequence variants.

Design decisions: sequencing artefacts and miscalls

GeneScreen is designed for rapid interactive analysis of large sequence batches by a human operator, not for full automation. It leaves the operator to make decisions related to common sequencing problems, such as "dye blobs", low signal, high background and anomalous peak mobility, and does not introduce computational corrections for such anomalies. This decision stems from the desire to avoid computational artefacts, including the overlooking of genuine sequence variants. Some illustrations of these concerns are shown in Figure 4. Firstly, 4A and 4B show forward and reverse sequences around a heterozygous substitution. In 4A, the different mobilities of the allelic peaks have resulted in a miscalled insertion. A computational correction for this mobility shift could be applied, but this might adversely affect other traces, such as that in Figure 4C; here, software "correction" of the overlapping C and A peaks would result in an incorrect heterozygous call. (The reverse sequence in Figure 4D reveals the peak overlap to be another artefact of anomalous migration.)

Resolving doubtful cases by reference to the reverse strand sequence, as in both of these examples, is not always possible, for various reasons. Firstly, it may not be obvious which of two discordant forward and reverse sequences has been miscalled.

Secondly, forward and reverse sequences may not both be available, either for reasons of economy or technical difficulty (as when sequencing through heterozygous indels, or long simple tandem repeats). Finally, forward and reverse sequences may occasionally be wrongly paired through operator error, as a result of which true heterozygous variants might be discarded, because of appearing to be absent from the sequence in the other orientation.

Correcting for variation in peak height can also result in conflicting results. Figure 4E shows a heterozygous C/T call at a position which in the reverse sequence (Figure 4F) is clearly homozygous C. This error results from the weakness of the genuine C peak in Figure 4E, compared both to flanking positions and to the background. If a scaling correction (based on a reference trace) were used to accept the small blue C peak as genuine, while discarding the red T signal as background, a risk would be introduced of miscalling similarly appearing positions where two small superimposed peaks do in fact reflect true heterozygosity.

Similar to Figure 4E, Figure 4G shows a miscall due to a sequencing artefact that lies between two peaks. The artefactual peak is of comparable height to the genuine C in Finally, using peak height ratios as a criterion of heterozygosity can be problematic.

Figure 4I

-J shows a position that is clearly heterozygous in the reverse sequence, but for which the allelic peak heights differ greatly in the forward sequence. A base-caller configured to use peak height ratios to ignore the spurious peak in 4E might also be liable to miscall traces such as 4I.

In most situations, overlooking a true sequence variant has more serious implications than identification of false positives. GeneScreen has therefore been designed with default settings and a user interface that minimize the impact of dealing with a large number of false positive calls. Rapid real-time visual inspection of all suspect positions within a large number of traces (including comparison of paired forward and reverse reads) allows correct decisions to be made with minimal user intervention.

Analysis of test sequences

To test GeneScreen's accuracy in detecting sequence variations, we used it to analyse several thousand sequences from ABCC8, BRCA1, BRCA2, GCK and TP53 (Table 1).

Each set of test data was analysed using both Mutation Surveyor (SoftGenetics, State College, PA) and GeneScreen, by operators familiar with each piece of software.

Each program detected an identical set of 177 rare sequence variants (Table 1). We also found GeneScreen to be at least as fast as Mutation Surveyor. (For example, the 30 amplicons covering the BRCA1 open reading frame of 96 subjects were analysed by one operator in ~5 hours, compared to ~9 hours using Mutation Surveyor.) GeneScreen also correctly genotyped all the common variants (SNPs) present in the test sequence files. It should be noted, though, that the identification of common SNPs is rather easy for a mutation-screening program to perform, because they show up as sequence variants at the same position in multiple sequence traces. They are therefore excluded from Table 1. Apart from common SNPs, none of the sequence variants in ABCC8, BRCA1, BRCA2 or GCK was homozygous. Most of the TP53 variants likewise appeared heterozygous, but in some cases LOH in tumour DNA resulted in variants appearing to be homozygous.

Table 1

Description of test sequences used to validate GeneScreen and mutations identified (excluding common SNPs). *See main text for discrepancy between number of exons and number of PCR products. †Amplicons were only sequenced in the forward direction. "Indels" refers collectively to deletions, insertions and insertion-deletions.

Sequence quality and miscall rate

To assess whether the false positive variants identified by GeneScreen resulted in most part from poor quality sequence, or from poor performance of its base-caller, we We also noted a strongly non-random distribution of the false positive results, according both to the identity of the individual PCR product (i.e. its sequence) and to the quality of the sequencing template (depending in turn on genomic DNA quality).

For this reason, we have not attempted to derive statistical measures of the sensitivity and specificity of GeneScreen; experimental variation, or choice of different genomic targets, would render such figures invalid. For example, our TP53 DNA samples had been extracted from archived tumour samples, and consistently produced higher false positive rates than samples extracted from blood. Similarly, sequences that involved reading through a long mononucleotide run were more susceptible to false positive results.

Platform-dependence is a perennial issue that has undoubtedly affected the adoption by less technically minded users of some highly sophisticated software packages (such as the Staden package, originally developed solely for a Unix/Linux environment). Limited user acceptance can also result from lack of familiarity with the command line environment through which some programs such as PolyPhred 11 are controlled. However, development of a truly cross-platform application inevitably involves greater effort than focussing on one operating system. To keep development effort manageable, GeneScreen utilizes the .NET framework and is hence not crossplatform. This choice, though, does ensure that it is accessible to most non-specialist computer users.

Compared to some other programs for sequence variant detection, GeneScreen does not go to great lengths to maximize the specificity with which variant positions are identified; this is in order to avoid loss of sensitivity for detection of rare variants. Our design consequently places considerable onus on the operator to distinguish true from false positives, and to that end, the grid presentation of GeneScreen maximizes the efficiency with which the user can inspect possible variant positions. The higher the sequence quality, the fewer artefactual variants are identified on the grid, and the quicker the task of analysis; consequently efficient mutation identification requires high quality sequence.

The initial interface window, showing GeneScreen's various tabs for loading reference files and test sequences, and performing analysis for the presence of mutations.

Figure 2

The GeneScreen results grid. Each column of the grid represents a position where one or more test sequences has a variant (red cell). The sequence around the putative variant is displayed by simply hovering the mouse over the cell. The grey cells represent a region that could not be aligned to the reference sequence.

Figure 4

Examples of migration artefacts and other sequencing anomalies that can interfere with mutation calling. See text for details.

Figure 5

Quality values of sample sequences, classified according to their number of identified sequence variants.

Licence for Publication

The Corresponding Author has the right to grant on behalf of all authors and does grant on behalf of all authors, an exclusive licence (or non exclusive for government employees) on a worldwide basis to the BMJ Publishing Group Ltd to permit this article (if accepted) to be published in Journal of Medical Genetics and any other BMJPGL products and sublicences such use and exploit all subsidiary rights, as set out in our licence (http://group.bmj.com/products/journals/instructionsfor-authors/licence-forms).

3 .

 3 ), within which each column represents a position at which a sequence variant exists in one or more trace files, while each row corresponds to one electropherogram. Provided appropriately annotated reference files were used, positions within exons are coloured blue, and introns white. Locations where a sequence deviates from the reference sequence appear as red cells, with the nucleotide substitution in black text. For each cell it is possible to view its underlying sequence data as an image of the electropherogram around the mutation, or a comparison of the mutant electropherogram to a control trace with normal sequence. The grid and displayed image can be dynamically linked, so that simply moving the mouse cursor across the grid displays the local sequence trace corresponding to the underlying cell, allowing very rapid visual inspection (Figure2). Left-clicking on a red cell invokes a floating menu that allows the user to annotate and save the sequence variant, to genotype SNPs, or to analyse and annotate complex mutations. Even for heterozygous insertions or deletions, GeneScreen can usually perform a correct subtractive analysis of the superimposed sequences downstream of the mutation point, correctly calling and annotating the mutation, as shown for the examples below in Figure All the analysis and annotation options are shown in detail in the online program guide.

Figure 4E .

 4E Figure 4E. It is therefore apparent that any filtering aimed at disregarding the artefactual peak in Figure 4G might adversely affect the calling of traces similar to Figure 4E.

Figure 3

 3 Figure 3Example of a web-page output format generated by GeneScreen. The blue and green sequence rows in the image panels show how the program deconvolutes heterozygous deletions and deduces the two allelic sequences. Note that GeneScreen flips the reverse strand sequence to match the orientation of the reference; consequently the superimposed staggered deletion alleles appear at the left of the image on reverse strand reads. The program includes a facility for editing and updating these web-pages to incorporate newer results.

Figure 1 Figure 2 Figure 3 Figure 4 Figure 5

 12345 Figure 1

  subset of our traces using the commercial ABI Sequencing Analysis v5.2 program, which assigns per-residue quality values defined in similar fashion to Phred qualities[START_REF] Ewing | Base-calling of automated sequencer traces using phred. II. Error probabilities[END_REF] . (We reasoned that if the false positives were attributable to a weak base-caller, quality values of reads containing false positives would be similar to those of files containing no variants at all.) Additionally, because sequence artefacts resulting in false positives might reflect either globally poor signal-to-noise or discrete problems such as 'dye blobs' or compressions, we used an approach of categorizing sequence quality according to the average quality values of the ten worst nucleotides in a sequence. This was found empirically to allow detection of both generally poor quality and short regions of poor sequence in otherwise good reads, while not being overly sensitive to the presence of true heterozygous variants (which also lower the quality score). Figure5shows that traces in which GeneScreen identified a false positive variant tend to be those with the lowest quality values as measured in this way.

	Gene ABCC8 31 Subjects Exons Amplicons* Traces Substitutions Indels 39 39 1209 † 20 2 analysed a random Of the 2143 trace files used for this analysis, 152 (6.8%) generated at least one false
	positive variant call, while 123 (5.4%) contained at least one true variant (including
	BRCA1 96 SNPs). Predictably, the latter group tended to contain only one variant per read (133 24 30 5760 15 3
	BRCA2 96 variants in 123 sequences), in contrast to reads with false positives (328 variants in 27 32 6144 5 6
	GCK 152 sequences). Overall, 29% of all variants identified in this analysis were true 160 10 10 1600 † 49 5
	variants. Within the 385 628 nt of sequence coverage, this corresponded to a false
	TP53 positive rate of 1 per 1175 nt (compared to a true positive rate of 1 per 2899 nt). 95 11 10 1900 57 15
	Total	16613	146	31

DISCUSSION

Identification of unknown functionally important rare sequence variants is an important activity in many research and diagnostic laboratories. The accuracy and volume of sequence data available at a given cost have steadily increased over the last decade, exposing the less standardized process of data analysis as increasingly tedious and error-prone. It is not surprising, therefore, that a number of computer programs have been developed with the aim of screening electropherograms for the presence of mutations. The fact that most diagnostic DNA laboratories continue to inspect patient sequences visually presumably reflects perceived limitations of available software.

Aside from the question of sensitivity, these perceived limitations may include aspects such as complexity of the software, platform-dependence, or its targetting primarily at other uses such as de novo sequence assembly or SNP detection. Unlike large genome centres (that may be focussed on the latter areas) diagnostic and small genetics research laboratories have generally not had the expertise to undertake dedicated software development. Ease of use is one reason for the wide popularity in diagnostic laboratories of the commercial Mutation Surveyor package (SoftGenetics LLC, PA, USA).

Programs aimed mostly at identifying SNPs can improve their specificity using population criteria, such as a requirement for each allele of a putative variant to occur a certain number of times in a set of sequence traces. This approach to reducing the false positive background is not applicable to rare disease-causing variants.

Developed primarily with the latter goal, GeneScreen therefore base-calls each sequence independently and then flags variants after comparison to the known reference sequence. No attempt is made to detect and remove false positives by reference to external data, since in some circumstances this would be likely to result in the overlooking of true positives that appear similar to sequencing artefacts (Figure 4). The authors strongly believe that the best way to keep the number of false positives manageable is to produce high quality sequence data.

To compensate for the variability in peak height and mobility that are intrinsic to dyeterminator sequencing, some mutation detection programs employ direct subtractive analysis of the difference between test and control sequence traces. As well as Mutation Surveyor, this approach is used by the tracediff component of the Staden package [START_REF] Bonfield | Automated detection of point mutations using fluorescent sequence trace subtraction[END_REF] and SeqDoC 12 . This method can perform excellently when sequence batches derive from templates of highly reproducible quality and quantity. It can be rather susceptible, though, to differences in quality and signal intensity between test and control samples [START_REF] Crowe | SeqDoC: rapid SNP and mutation detection by direct comparison of DNA sequence chromatograms[END_REF] . Since such sample variability is common and often unavoidable when dealing with patient material, we decided against such an approach in the design of GeneScreen. A brief summary of the features of a number of mutation detection programs is shown in Table 2 Subjectively refers to ease of setup and "learning curve" for those who are not computer specialists.


Gabriella Gaglio
email: ggaglio@unime.it

Simon Allen

Lee Bowden

Mark Bryant

Eric R Morgan

C Gortázar

Parasites of European hedgehogs (Erinaceus europaeus) in Britain: epidemiological study and coprological test evaluation

Keywords: Hedgehog, Helminth, Ectoparasite, Epidemiology, Coprological diagnosis

Seventy-four European hedgehogs (Erinaceus europaeus) that had died in wildlife rehabilitation centres were dissected and their parasite burdens documented. Overall parasite prevalence was 91%, and a total of six helminth species were isolated: five nematodes (Crenosoma striatum, Eucoleus aerophilus, Capillaria erinacei, Capillaria ovoreticulata and Capillaria spp.), one trematode (Brachylaemus erinacei) and one acanthocephalan (Oliganthorhynchus erinacei). The tick Ixodes hexagonus and flea Archeopsylla erinacei were also collected. The effect of parasite infection on body condition was assessed by correlation of burdens with the residuals of weight-skeletal length regression. Tick presence was positively related to body condition; for other parasites, no significant relationship was found. Faecal egg or larval count was closely correlated with adult parasite burden for C. striatum and Capillaria/Eucoleus spp., but not for other species. Coprological analysis should therefore be useful for in vivo studies of nematode parasite infection in hedgehogs. The epidemiology of parasites in hedgehogs and their possible role in recent population declines are discussed.

Introduction

Hedgehogs (Erinaceus europaeus) are one of the most distinctive and familiar groups of mammals and are widely distributed throughout the Old World (Reeve 1994). In Great Britain, the hedgehog is found almost everywhere, but tends to be scarce or absent from wet areas and extensive pine forests [START_REF] Morris | Hedgehogs[END_REF]. In recent years, the number of hedgehogs in Great Britain appears to have markedly declined (The Mammal Society 1999). The cause of this decline is not known. [START_REF] Reeve N ; ) Hedgehogs | Mortality factors affecting wild hedgehogs: a study of records from wildlife rescue centres[END_REF] reported that just over 40% of hedgehog mortalities were thought to be caused by road traffic accidents, garden and pet injuries and predation, with 59% of deaths resulting from 'natural cases' such as disease. Hedgehogs are the most common species of mammal admitted to wildlife hospitals in the UK [START_REF] Bullen | Hedgehog rehabilitation. British Hedgehog Preservation[END_REF][START_REF] Molony | Factors affecting the success of wildlife rehabilitation in Britain[END_REF]. As nonterritorial nocturnal species with a flexible diet consisting of a wide range of invertebrates, as well as small vertebrates and carrion (Reeve 1994), rehabilitated hedgehogs are sometimes released at novel locations under the assumption that they are readily capable of adapting to their new environs. For this reason, [START_REF] Molony | The effect of translocation and temporary captivity on wildlife rehabilitation success: an experimental study using European hedgehogs (Erinaceus europaeus)[END_REF] utilised hedgehogs as a focal species for translocation studies. Their results suggest that direct translocation had a damaging effect on survival, but a short period in captivity prior to release greatly improved the chance of survival. These beneficial effects of a period of retention in captivity were also evident for those individuals that had been rehabilitated following injury or illness.

Given the recent and possibly ongoing population declines, rehabilitation of sick or injured hedgehogs could have conservation as well as welfare value. Anecdotal evidence from rehabilitation centres suggests that parasite infection, especially lungworm, is an important cause of illness in hedgehogs. Parasites are also known to have significant effects on the population dynamics of their hosts in other systems [START_REF] Irvine | Parasites and the dynamics of wild mammal populations[END_REF], and increased levels of infection could contribute to population declines. Unfortunately, clinical and ecological studies on parasites of hedgehogs are held back by a lack of validated tests for infection in the living animal.

The present study aims to record the parasite species present in British hedgehogs, along with typical levels of infection in animals admitted to rehabilitation centres, and to assess impacts on individual fitness through decreased body condition. Coprological test results are compared with burdens of adult parasites in order to assess their value as indicators of infection status.

Materials and methods

Study design and sampling method

Between November 2008 and July 2009, 74 hedgehogs (42 males and 32 females) were examined for parasites. The animals had died of various causes between December 2006 and September 2008 in two wildlife rehabilitation centres in Great Britain, one in south Wales and the other in the east of England. Hedgehogs that had spent more than 1 week in the rehabilitation centre before death were excluded in order to reduce bias arising from weight loss during prolonged captivity. Carcasses were frozen for storage and transport. After defrosting overnight and before dissection, weight and overall body length were noted, as well as the distance between the carpus and the end of the distal phalanx of the second digit of the left forelimb (=manus length).

Adult parasite counts

The stomach and intestines of hedgehogs were removed, and parasites were extracted using the total worm count technique [START_REF] Euzeby | Diagnostic expérimental des helminthoses animales Travaux pratiques d'helminthologie vétérinaire[END_REF]. Briefly, the intestines were separated from the mesentery and expressed before being incised along their length, and the mucosa was rinsed in water with firm digital pressure. The stomach was similarly, but separately, processed. The heart and lungs were removed, the aorta clamped using forceps and the lungs flushed using a technique modified from foxes [START_REF] Morgan | Angiostrongylus vasorum and Eucoleus aerophilus in foxes (Vulpes vulpes) in Great Britain[END_REF]): Tap water was forced into each lobe of the lung using a 10-ml syringe with a 24-gauge needle to disrupt the pulmonary tissue and dislodge any nematodes. The fluid was collected as it flowed out through the trachea. Blood clots and mucus that blocked the trachea were removed using a pipette and deposited with the rest of the lung flush to be examined for nematodes. Extracts from guts and lungs were inspected under the dissection microscope and parasites counted, removed and placed in lactophenol on a standard microscope slide for identification. The integument and the bag in which hedgehogs were frozen were inspected for ectoparasites, which were removed for counting and identification. Parasites were identified to species level according to [START_REF] Yamaguti | Systema helminthum. The digenetic trematodes of animals[END_REF][START_REF] Yamaguti | Systema helminthum. Acanthocephala[END_REF], [START_REF] Skrjabin | Tricocephalidae and Capillaridae of animals and man and the diseases caused by them[END_REF], [START_REF] Hillyard | Synopses of the British fauna. Ticks of North-West Europe[END_REF] and [START_REF] Whitaker | Handbook for the identification of British insects. Fleas (Siphonaptera)[END_REF].

Faecal egg/larva counts Two grammes of faeces was taken from the rectum during dissections. Faecal egg and larval density were estimated using a modified McMaster technique (MAFF 1986), with a sensitivity of 50 eggs/larvae per gramme, using a salt and sucrose flotation solution (specific gravity, 1.280). Where <2 g of faeces was available, dilution factors were modified accordingly.

Statistical analysis

Descriptive statistics, including the prevalence (P), the mean intensity (I m ) and the mean abundance (A m ), were calculated for each parasite species or genus recovered, as proposed by [START_REF] Bush | Parasitology meets ecology on its own terms[END_REF]. Abundance was compared between sexes and seasons (summer/winter) using Student's t test. The relationships between different parasite species, between parasite burden and body condition (as estimated by the residual of mass/manus length regression) and between adult parasite burden and faecal egg or larval count were compared using Pearson correlation, with Bonferroni correction for multiple testing. In each case, overdispersion in parasite burden was dealt with by log(x+ 1) transformation. After checking for normal distribution of residuals, linear regression was used to predict the number of adult parasites present from the density of eggs or larvae in the faeces and body mass from manus length. The sensitivity of coprological tests was estimated as the number of infected animals identified coprologically as a proportion of those found to carry parasites on postmortem examination.

Results

Adult parasites

Helminth parasites were present in 67 out of 74 hedgehogs examined (91%). A total of four species of nematode (Crenosoma striatum, Eucoleus aerophilus syn. Capillaria aerophila, Capillaria erinacei, Capillaria ovoreticulata, plus Capillaria sp. females), one trematode (Brachylaemus erinacei) and one acanthocephalan (Oliganthorhynchus erinacei syn. Echinorhynchus erinacei) were detected. The ectoparasites Ixodes hexagonus and Archeopsylla erinacei were also found. A summary of levels of infection is presented in Table 1. The most common parasite was C. striatum in the lungs, and the most abundant was Capillaria spp. in the stomach and small intestine. Gut flukes (Brachylaemus), Eucoleus in the lungs and ticks (I. hexagonus) were also common, with some high individual burdens, whilst acanthocephalans and fleas were less common.

Correlations between parasite species

Using serial bivariate correlation, positive relationships were found between burdens of C. striatum in the lung and Capillaria spp. in the stomach (r=0.497, p<0.001) and intestine (r =0.646, p <0.001), Capillaria spp. in the stomach and intestine (r=0.809, p<0.001) and B. erinacei and O. erinacei (r=0.473, p<0.001).

Coprological test results

Parasite eggs or larvae were found in 51 of 74 (69%) hedgehogs. Capillaria sp. eggs (including E. aerophilus, which are morphologically similar) were found in 46 (62%), Crenosoma larvae in 38 (51%) and trematode eggs in 10 (13%). Test sensitivity, as defined above, was 78% for Crenosoma, 89% for Capillaria and 24% for Brachylaemus, with no false positives detected for any species. The relationship between adult worm burden and faecal egg/ larval count is shown for Crenosoma and Capillaria in Fig. 1. The density of trematode and acanthocephalan eggs correlated poorly with corresponding adult burdens.

Variation with sex, season and body condition

There were no significant differences in the abundance of any parasite species between sexes or seasons. The regression of body weight and manus length used to estimate body condition is shown in Fig. 2. Tick abundance was positively related to body condition (r=0.417, p<0.001).

Discussion

In Europe, epidemiological studies on the parasite fauna of hedgehogs have been conducted in Germany [START_REF] Laux | Extent and intensity of endoparasitic infection in the hedgehog[END_REF][START_REF] Wieland | Endoparasiten beim Igel[END_REF], Norway [START_REF] Keymer | Zoonoses and other findings in hedgehogs (Erinaceus europaeus): a survey of mortality and review of the literature[END_REF] and Italy [START_REF] Poglayen | Helminthes found in hedgehogs (Erinaceus europaeus) in three areas of Italy[END_REF]. Studies on British hedgehog parasites focused particularly on parasitic lung disease [START_REF] Majeed | Lesions associated with a Capillaria infestation in the European hedgehog (Erinaceus europaeus)[END_REF][START_REF] Majeed | Occurrence of the lungworms Capillaria and Crenosoma spp. in British hedgehogs (Erinaceus europaeus)[END_REF][START_REF] Robinson | Veterinary care of the hedgehog[END_REF][START_REF] Cousquer | Analysis of tracheal sputum for diagnosing and monitoring verminous pneumonia in hedgehogs (Erinaceus europaeus)[END_REF]Molenaar 2008, unpublished data), and few epidemiological studies of hedgehog parasites have been undertaken in Great Britain. [START_REF] Boag | The prevalence of helminth parasites from the hedgehog Erinaceus europaeus in Great Britain[END_REF] found prevalence of 77% for C. striatum, 85% for C. erinacei, 8% for Capillaria spp. and 8% for Vampirolepis erinacei. [START_REF] Bunnell | The incidence of disease and injury in displaced wild hedgehogs (Erinaceus europaeus)[END_REF] found I. hexagonus (14%), C. striatum (11%) and sarcoptic mange (6%). Reeve (1994) and [START_REF] Bexton | Hedgehogs[END_REF] wrote about clinical signs during lung and intestinal Capillaria, trematode, acanthocephalan and ectoparasite infections. Severe infections of intestinal parasites can cause lethargy, weight loss and diarrhoea, and lung parasites can cause lethargy and dyspnoea [START_REF] Bexton | Hedgehogs[END_REF], whilst ectoparasites, especially ticks, can induce anaemia [START_REF] Pfäffle | Tickinduced blood loss leads to regenerative anaemia in the European hedgehog (Erinaceus europaeus)[END_REF]. [START_REF] Sainsbury | Health and welfare of rehabilitated juvenile hedgehogs (Erinaceus europaeus) before and after release into the wild[END_REF] considered parasitism to be involved in mortality of rehabilitated hedgehogs after release. Our results confirm that hedgehogs are commonly parasitised and can carry high burdens particularly of lungworms (especially C. striatum), gastrointestinal Capillaria spp., gut flukes (B. erinacei) and ticks (I. hexagonus).

Flea burdens are likely to be underestimated in this study since fleas leaving the host after death and before storage of the carcase were not counted.

The high levels of parasitism found in hedgehogs, together with apparent increases in the abundance of some parasites in other wildlife species but relevant to hedgehogs, such as ticks [START_REF] Scharlemannn | Trends in ixodid tick abundance and distribution in Great Britain[END_REF]) and slugtransmitted lungworms [START_REF] Morgan | Angiostrongylus vasorum and Eucoleus aerophilus in foxes (Vulpes vulpes) in Great Britain[END_REF][START_REF] Morgan | Canine pulmonary angiostrongylosis: the influence of climate on parasite distribution[END_REF], make increased parasitism a plausible factor in observed population declines. However, the effect of parasites on hedgehogs at individual and population levels is difficult to ascertain. The results reported here are from animals admitted to rehabilitation centres and might not be representative of the general population, but could nevertheless be more sensitive indicators of negative effects of parasitism on individual fitness. The present study, however, failed to find negative correlations between body condition, as estimated by residuals of skeletal length/mass regression, and the abundance of any parasite species or total parasite burden. At face value, this might suggest that parasites at these levels of infection do not have negative effects on hedgehog body condition or fitness. However, such crosssectional correlations are notoriously insensitive indicators of effects on host fitness, which could be subtle, and restricted to certain times of year or population subclasses [START_REF] Irvine | Parasites and the dynamics of wild mammal populations[END_REF]. Moreover, antagonistic relationships between parasitism and resource acquisition could obscure observed effects. For example, many parasite species are acquired from food, so well-fed individuals might accumulate higher parasite burdens. This is compounded by food choice, e.g. between gastropod molluscs (intermediate hosts of Crenosoma and Brachylaemus) and insects/earthworms, which itself is likely to be influenced by age [START_REF] Dickman | Age-related dietary changes in the European hedgehog, Erinaceus europaeus[END_REF] as well as by nutritional state and habitat use [START_REF] Jones | Diet of hedgehogs (Erinaceus europaeus) in the upper Waitaki Basin, New Zealand: implications for conservation[END_REF]. Gastropods form a small percentage of hedgehog diets, and marginalised animals with poor prey choice might be more likely to ingest these, and to acquire infection, as well as being more likely to present at rehabilitation centres. Similarly, the positive correlation observed in the present study between body condition and tick burden could reflect greater range size or habitat quality in fitter animals. Habitat type and vegetation cover are well known to influence tick abundance [START_REF] Medlock | Investigation of ecological and environmental determinants for the presence of questing Ixodes ricinus (Acari: Ixodidae) on Gower, South Wales[END_REF], and ticks have been used as markers of habitat use in hedgehogs [START_REF] Thamm | Ectoparasite infestation of hedgehogs (Erinaceus europaeus) are associated with small-scale landscape structures in an urban-suburban environment[END_REF]. It is also possible that feeding whilst in rehabilitation enabled parasitised animals to gain weight, obscuring negative relationships between burden and body condition. Positive correlations between burdens of different parasite species, as found here amongst helminths, could indicate co-transmission through diet or habitat use or synergistic interaction in the host, for example via immunosuppression.

The use of mass-length regression as an indicator of body condition in hedgehogs could be a useful tool to detect effects of parasites on fitness, especially in living animals. However, this approach is not without problems. Although residuals were normally distributed, smaller hedgehogs tended to be in poorer condition (Fig. 2 better body condition but also more heavily infected with parasites. If this is the case, even strong negative effects of parasitism on body condition could be obscured. This would further be compounded by the tendency to underestimate parasite burdens in older age classes of free-living animals [START_REF] Pacala | The relation between the number of parasites per host and host age-population dynamic causes and maximum likelihood estimation[END_REF]. Unfortunately, few helpful, validated measures of age are available in hedgehogs. Dentition [START_REF] Delecourt | Age determination test in the fox and hedgehog performed by counting cementum rings[END_REF] and skeletal maturity [START_REF] Morris | Epiphyseal fusion in forefoot as a means of age determination in hedgehog (Erinaceus europaeus)[END_REF]) can be useful early in life, but not to reliably distinguish the ages of mature hedgehogs, whilst eye lens weight [START_REF] Augusteyn | On the relationship between rabbit age and lens dry weight: improved determination of the age of rabbits in the wild[END_REF]Gacic et al. 2007;[START_REF] Janova | Is the eye lens method of age estimation reliable in voles?[END_REF]) has yet to be validated for this species. Better methods for age determination in hedgehogs are needed if the effects of parasites on individual fitness and population viability are to be determined. Observation of parasite eggs or larvae in the faeces was found to be a useful indication of infection and, for C. striatum and combined Capillaria and Eucoleus species, a good reflection of adult parasite burden. This suggests that there is little density dependence in egg output, e.g. mediated by acquired immunity. Moreover, the strong correlation between egg or larval density and adult burden validates the use of coproscopy as a measure of parasitism in living animals. This opens the way for more detailed studies of parasite dynamics, especially longitudinal studies in free-living animals, as well as the effects of these parasites and their removal on the clinical performance of animals in rehabilitation. Such studies are needed if the value of antiparasitic treatment in captivity is to be evaluated and will form an important part of ongoing investigations into the possible role of parasites and other diseases in population decline.

Fig. 1

 1 Fig. 1 Relationship between adult worm burden and density of corresponding immature parasite stages in the faeces for a C. striatum larvae: regression equation log (adults+1) = 0.231+0.578 × log(larvae+ 1). Pearson correlation coefficient 0.895, analysis of variance (ANOVA), F 1, 72 = 289.18, p<0.001. b Capillaria spp. eggs, including E. aerophilus: regression equation log(adults+1) = 0.293+0.560 × log (eggs+1). Pearson correlation coefficient 0.897, ANOVA, F 1, 72 = 295.61, p<0.001. Egg and larval density are treated as predictors of adult burden

Table 1

 1 Parasites recovered from 74 British hedgehogs by site of infection Capillaria spp. consisted of C. erinacei and C. ovoreticulata as well as females identified to genus only P percent prevalence, I m mean intensity, A m mean abundance, SD standard deviation, R range

		Lungs		Stomach	Intestine		Intestine and	Skin
							gut serosa
		Crenosoma	Eucoleus	Capillaria spp. Capillaria spp. Brachylaemus	Oliganthorhynchus	Ixodes	Archeopsylla
		striatum	aerophilus			erinacei	erinacei	hexagonus	erinacei
	P	71	32	66	61	55	18	59	8
	I m (±SD)	71 (±153)	9 (±15)	110 (±135)	92 (±126)	85 (±179)	3.5 (±2.8)	24 (±57)	11 (±13)
	A m (±SD) 50 (±132)	3 (±9)	73 (±122)	56 (±108)	47 (±140)	0.6 (±1.7)	14 (±45)	0.9 (±1.7)
	R	1-961	1-56	3-591	2-605	1-991	1-11	1-283	2-33

Pearson correlation coefficient 0.895, analysis of variance (ANOVA), F 1, 72 = 289.18, p<0.001. b Capillaria spp. eggs, including E. aerophilus: regression equation log(adults+1) = 0.293+0.560 × log (eggs+1). Pearson correlation coefficient 0.897, ANOVA, F 1, 72 = 295.61, p<0.001. Egg and larval density are treated as predictors of adult burden

  ), and there is likely to be confounding of this statistic with age. Since parasite accumulation is age dependent, it is possible that older individuals are generally not only in generally

		1400 1400
		1200 1200
		1000
	g / g s / ass ma m	800 800 600 600
		400 400
		200 200
		0
		16 16	18 18	20 20	22 22	24 24	26 26	28 28
					manus length / mm
	Fig. 2 Regression of body mass (g) against manus length (mm).
	Residuals were used as an index of body condition. Pearson
	correlation coefficient = 0.704; regression equation: Mass = -662+
	43.9 × manus length. ANOVA, F 1, 71 = 69.65, p<0.001

Eur J Wildl Res (2010) 56:839-844


Nadra Harchaoui

Slimane Hellal

Jean-Georges G Gasser
email: gasser@univ-metz.fr

Benoit ; Grosdidier

and atomic transport properties of liquid alkaline metals by molecular dynamics simulations

Keywords: atomic structure, liquid metals, atomic transport, molecular dynamic simulations, diffusion model pseudopotentials first-principle model potential, molecular dynamics, liquid metal, structure factor, selfdiffusion

Self consistent Shaw Optimized Model Potential

Application to the determination of structural and atomic transport properties of liquid alkaline metals by molecular dynamics simulations.

I. Introduction

First-principle pseudopotentials remain up to today an important tool of investigation and have attracted several theoretical works [1][START_REF] Harrison | Pseudopotentials in the Theory of Metals[END_REF][START_REF] Bachelet | Strategies for Computer Chemistry[END_REF]. Indeed, to understand the physical properties of simple metals from a fundamental point of view, one needs a quantum mechanic treatment of the metallic bond. This can be achieved by combining the pseudopotential formalism and electronic structure calculations with second-order perturbation recipes. This way leads to the "pair theory" for metals [START_REF] Egelstaff | An Introduction to the Liquid State[END_REF] 2 in which the interactions between ions are described in terms of an effective pairwise potential ( ) r V eff . An interesting but stringent test of ability for the "pair theory" to accurately predict most of liquid state properties lies in the ionic structure calculation [START_REF] Egelstaff | An Introduction to the Liquid State[END_REF][START_REF] Hansen | Theory of Simple Liquids[END_REF][START_REF] Waseda | The Structure of Non-crystalline Material: Liquid and Amorphous Solids[END_REF][START_REF] Hafner | From Hamiltonians to Phase Diagrams[END_REF]. This one is given by the paircorrelation function ( ) r g in real space or by the static structure factor ( ) q S in the reciprocal one. It is closely conditioned by the shape of the ionic pair potential ( ) r V eff . Furthermore, within the "pair theory" of fluids, it was shown that there is a one-to-one relationship between g(r) and ( ) r V eff [START_REF] Egelstaff | An Introduction to the Liquid State[END_REF][START_REF] Gray | Theory of Molecular Fluids[END_REF].

However, a reliable test needs computer simulation method to bypass approximate theories of liquid state [START_REF] Egelstaff | An Introduction to the Liquid State[END_REF][START_REF] Hansen | Theory of Simple Liquids[END_REF]. In the present work, we have firstly determined ( ) r g and ( )

q S
for some alkali metals (lithium at 463 K, sodium at 378 K, and potassium at 338K). These two important functions are computed by using the molecular dynamic (MD) simulation technique [START_REF] Allen | Computer Simulation of Liquids[END_REF] with the effective ionic pair potentials as main ingredients. We can already notice that ( ) r V eff and ( ) q S play a central role in the physic of simple liquids. For instance, the microscopic theories of collective dynamics involve the use of precise form of ( )

q S
and of ( ) r V eff [START_REF] Monaco | [END_REF]11]. At the second stage, we have extended the pseudopotentiel-based calculations in order to investigate the atomic transport properties which are quite sensitive to the interactions. This task can be achieved through the calculated Velocity-Autocorrelation Function (VAF), the corresponding spectral density and the self-diffusion constant.

Their temperature-dependence behavior in liquid state is discussed.

The ionic potentials are built within the framework of the first-principle pseudopotential formalism.

We use the Optimized Model Potential (OMP) of Shaw [1] to describe the electron-ion interaction.

Many-body effects in valence electrons gas are properly included through a full-nonlocal screening [12] of the bare electron-ion model with a suitable dielectric screening function [13,14] which has a good theoretical background. That is not without some difficulties due mainly to the nonlocal and energy-depending nature of the OMP-model [15][16][17]. We point out these difficulties in the pseudopotential context in section 2 where we will bring out the calculation of the OMP parameters as a central problem. For clarity needs, we briefly remind in section 3 some quantitative aspects of the ionic structure and of the self-diffusion for pure liquid metals. In the same section, we will also present the computing conditions with the standard molecular dynamics. In the last section 4, the essential features of the calculated pair potentials ( )

r V eff
and the behavior of the structural quantities ( ) r g and ( )

q S
for the studied alkaline metals will be discussed in connection with the energydependence of OMP parameters. Next to this, we will report and interpret our theoretical results upon other quantities that provide detailed and systematic microscopic understanding of the self-diffusion phenomena in liquid state. Axiomatic principles of the pseudopotential method within the plane-wave-basis formalism are wellknown and are described in a number of monographs and textbooks [START_REF] Harrison | Pseudopotentials in the Theory of Metals[END_REF][START_REF] Bachelet | Strategies for Computer Chemistry[END_REF][START_REF] Heine | Solid State Physics[END_REF][START_REF] Cohen | Electronic Structure and Optical Properties of Semiconductors[END_REF]. The main goal of this method is to remove the core states in the electronic structure calculations so that the strong fullelectron potential in the one-electron Schrödinger-type equation is replaced by an effective much weaker potential, namely the "pseudopotential". There are many ways to construct pseudopotentials (see brief review in [START_REF] Bachelet | Strategies for Computer Chemistry[END_REF][START_REF] Cohen | Electronic Structure and Optical Properties of Semiconductors[END_REF]). During the last three decades they have been generated from all-electron calculations for a free atom so that their main scattering properties are captured [START_REF] Hamann | [END_REF][21][22][23][24]. Such scattering properties can be indeed rephrased in terms of logarithmic derivatives arguments [START_REF] Hamann | [END_REF]21].

However, these atomic calculations involve various underlying problems that can spoil the quality of the pseudopotential. Most of them are mainly associated with non-linear aspects of density-functional exchange-correlation-energy [25]. Much more, "transferability" of the bare ion pseudopotential [21] can be also altered as regard to some other considerations such as the Kleinmann and Bylander factorization [22,26] or the "chemical hardness" criterion [27]. In addition, Hafner [START_REF] Hafner | From Hamiltonians to Phase Diagrams[END_REF] already stressed the fact that «the use of norm-conserving pseudopotentials [START_REF] Hamann | [END_REF]21] in perturbation calculations yields some serious problems». If one now turns to alkali metals, Na is a prototype element in the field of application of the pseudopotential theory. However, unlike Na which is an ideal case, light element Li, having strong pseudopotential, requires great attention in the atomic calculations [23]. Previously, Shenoy and Halder [28] have emphasized the singularity of Li ion. The latter, without p-core states, must be described by a nonlocal pseudopotential. To elude the difficulties which appear when we dealt with ab initio pseudopotential, we performed our electronic structure calculations using OMPmodel potential of Shaw [1,[15][16][17]. This model is built in the spirit of the "Quantum Defect Method" [29] so that it matches to the observed atomic energy levels of the isolated ion [START_REF] Moore | Atomic Energy Levels[END_REF]. In this approach, difficulties about exchange and correlation with the core electrons are somewhat minimized [START_REF] Heine | Solid State Physics[END_REF]. The real-space representation of Shaw's bare-ion OMP-model potential reads: are related to the first ones through the following "optimization" relationship:

l l l l l P ˆ (r) w r z w ˆ0 0 ion V 0 ion ∑ = = - - = , ( 2
( )       - - = r z E A r) Θ(R (r) w V ion l l l , with m , m , m m P ˆl l l l l ∑ = - = = (2.
V z (E) A (E) R = l l (2.3)

B. Model Potential for a metal

We consider now a homogeneous metallic medium at observed density, so that, the valence charge V z , the atomic volume 0 Ω , the atomic radius a R such as

3 a R 4π Ω 3 0 =
, the electronic number density " n ", the Fermi momentum

( ) 1/3 F n 3π k 2 =
, and the electron-sphere radius

3 1 V a s z R r - = are
then well-defined. So, to describe valence-electron-ion interaction in a metal environment using OMPmodel is not an easy task mainly because it depends on the energy on an absolute scale. As a preliminary step, the OMP-model parameters must be evaluated at the valence energy k E r in the metal shifted by an amount

k ∆E r namely at k k k ∆E E E r r r - = *
. Because it is found that model parameters have rough linear energy-dependence, Animalu and Heine [START_REF] Animalu | [END_REF] gave a faithful procedure for the pure metal case. Indeed, let us consider with a high degree of accuracy, that:

( ) ( ) ( ) *       ∂ ∂ + = * F E F k F * k E A E - E E A E A l r l r l (2.4)
In (2.4)

xc xc s 3 2 v 2 F µ ε r 0.6Z 2m k 5 2 BEE MIE E + - +         + - - = (2.5)                 - + = ∆ 2 a M s v xc F r R 4 3 3 2r Z µ E (2.6)
where xc ε is the exchange-correlation energy per electron for a homogeneous electron gas. Within LDA approximation, the exchange-correlation potential can be obtained as: , several interpolation formulas for the correlation energy c ε are available in the literature. Since our present calculations otherwise performed will be compared with previous work of Cowley [32], it may a good thing to recall that this author used the one given by Nozières and Pines [33]

( ) * = F v M E A Z R l . Hence to estimate ( ) * F 0 E A
, equations (2.5) and

(2.6) together are solved iteratively. We think that this scheme is inconsistent since it doesn't lead to a unique determination of * F E . On the other hand, Cowley assuming that

( ) * = F 0 v M E A Z R performed similar iterative scheme to get ( ) * F E A l
with clearly a unique value of * at the shifted Fermi energy, are directly evaluated in the selfconsistent manner within the first order pseudopotential-perturbation theory. This method has only be used for electronic transport properties [37] through the electron ion interaction, never for atomic structure and atomic transport through the ion ion effective potential. It turns out that from this selfconsistent calculation of

( ) * F E A l
, for which we use the abbreviation OMP-Self, the binding energy per electron BEE can be conversely estimated accounting of either the Animalu and Heine procedure or the Ballentine and Gupta ones. Before, equation (2.4) must be solved for:

( ) ( ) E A 0 A E A E 0 F F       ∂ ∂ - = * * l l l (2.
( ) * F E A l
for Li, Na and K are gathered in Table I along with previous theoretical results of Cowley. This author used the three like Animalu and Heine procedures denoted above OMP-AH, OMP-BG1, and OMP-BG2. We note that OMP-Self, OMP-AH, and in some measure OMP-BG1, respective calculations fairly converge. On the other hand, we agree with Cowley that « the inclusion of the particular form of inhomogeneity correction » given by Ballentine and Gupta (OMP-BG2) leads to « a gross shift in the absolute energy of the conduction band ». Further, we will confirm this statement. In the last column of table I, we also report the crude estimate of BEE issued from our accurate OMP-Self calculation of

( ) * F E A l
together with the use of equations (2.5),

(2.6) and (2.7). We wish to point out that our predicted BEE for Li is in a close agreement with experimental data. On the other hand, the result for Na is rather less good and the K one's turns only qualitative.

Since OMP-model is a nonlocal operator that also depends on the energy, the implementation of electronic structure calculation is not straightforward. First, we have dealt with some other cumbersome concepts that are closely linked to the energy dependence of the OMP-model. It is not convenient here to dwell too long upon these. Their signification and their importance are widely stressed in the literature. Shaw [16] introduced effective valence * v z , and effective masses k m r and

( ) k E m r
. However, the depletion hole [15] d n that is also a manifestation of the energy-dependence of the OMP-model, accounts for the difference between the true and pseudo-charge densities in the metals. It is instructive to note that the expression for the depletion hole given by Shaw and Harrison [15], equation (2.8) in this reference, is closely related to the Friedel "sum rule" and so to the scattering properties. Another crucial problem is connected with the many-body effects among electron-valence gas. These ones are incorporated in sitting up the OMP-screened-model potential. To do accurate calculations taking account of the nonlocal nature and energy-dependence of the OMPmodel, one used full-consistent screening [12,17] with the "electron-test-charge dielectric function" ε(q) and the local-field correction function G(q) [13,14]. To include the exchange and correlation, the electric dielectric function ε(q) is related to the Hartree dielectric function (q) ε H following the relationship:

( )

) ε G)(1 (1 ) (m 1 1 ε H 2 E F k - - + = (2.8) where               - + - + + = η 2 η 2 Log 4η η 4 1 η πk 2 1 ε 2 2 F H , and F k q η =
. In (2. 8), ( )

F k E m
specifies the Shaw effective mass evaluated at the Fermi momentum. Among the most well-physically-based forms for G(q) at the present time (see, e. g., ref [14]) we have chosen to use in our calculations the one of Ichimaru and Utsumi [13]. . In the second-order perturbation theory, one has need of the matrix element of the bare model potential that is to say:

k q q 0 ion f v k w q k N r r r r r + = + (2.9)
The so-called unscreened-form factor (2.9) contains two matrix elements q v and k q f r r which correspond to respectively the local and to the nonlocal parts of the bare model potential. If the complex expression for k q f r r can be found in the Shaw paper [1], q v can be deduced from the simple form:

2 0 q q Ω z 4π v V - =
(2.10)

The corresponding matrix element of the screened model potential or screened form factor is also useful. It takes the following expression [START_REF][END_REF]:

( ) ( ) q g f q ε w k w q k N k q 0q ion + + = + r r r r r (2.11)
In (2.11) the term 0q w is a sum of the local potential due to the valence charge q v (2.10) and that due to the renormalized depletion hole dq v . This last is corrected by exchange and correlation through [START_REF][END_REF] ( )

q G so that, ( ) dq q v G 1 v w 0q - + = (2.12)
where:

2 0 E dq dq q Ω m n 4π v = (2.13)
If we know the spatial distribution of the depletion hole ( ) R is a weighted mean of (E) R l model radii [START_REF] Appapillai | [END_REF]. With this choice, the normalized pseudo-wave-function and the normalized eigenfunction have same shape and same amplitude in the region outside the sphere, like that must be. If we define here the modified Coulomb factor as

r
[ ] G 1 q 4π v 2 * c - =
, the nonlocal screening contribution ( ) q g in (2.11) may be then written: The ( )

( ) ( ) k d ) E (E m m ) k , q f( ε 2π v g(q) 3 Ω 3 F 0 q k 0 k q k E k E * c
q g
function accounts for the exchange-correlation effect and energy-dependence of the OMPmodel potential. This definition of ( ) q g [START_REF][END_REF] is otherwise different in comparison with Shaw's one [16]. The first order in expansion for the energy of electronic states in perturbation theory is given in terms of effective masses as follows [16]:

( ) k E k 2 o k m 2m k E r r r = (2.15)

C. Effective interatomic pair potential

From the first-principle model potential presented above, one can determine the effective pair potential (r) V eff between an ion and another one in liquid metal. The second-order perturbation theory together the pseudopotential formalism leads to the familiar expression [START_REF] Hafner | From Hamiltonians to Phase Diagrams[END_REF][START_REF][END_REF]41]:

[ ] dq q sin(qr) (q) F π 2 1 r z (r) V 0 N 2 * V eff         - = ∫ ∞ (2.16)
This pair potential is connected to the model pseudopotential via the so-called normalized energywave-number characteristic (q) F N [START_REF] Hafner | From Hamiltonians to Phase Diagrams[END_REF][START_REF][END_REF]41]. The bare model potential for ions is suitably screened, taking account properly of its nonlocal and energy-dependent nature (full-nonlocal screening). The core shift is self-consistently evaluated according to the procedure that Hallers et al. [37] used for electronic transport calculation. To our knowledge it is the first time that this method is used for effective potentials. Hence, (q) F N has the very complex expression [START_REF] Waseda | The Structure of Non-crystalline Material: Liquid and Amorphous Solids[END_REF]13,[START_REF][END_REF]:

      - - + + + - -         - = 2 2 0q 2 * v 0 N dq 0q 2 2 v G) G(1 h g ε g w 2 w ε ε 1 G 1 1 4ππ q Ω (q) F
(2.17)

Except the term ( ) q h , all quantities in this equation are defined in the precedent subsection. As for ( ) q g defined by (2.14), the ( )

q h
function arises from the nonlocal part of the bare-ion OMP-model.

By including exchange-correlation effect and effective masses, it can be written [START_REF][END_REF]:

( ) ( ) k d ) E (E m m ) k , q ( f 2π v h(q) 3 0 2 F k k q k 0 k q k E k E 3 * c r r r r r r r r r ∫ ≤ + + - = (2.18)
So, within the "pair theory" of liquid metals, the effective pair potential (2.16) can be employed later on to describe the atomic structure and the self-diffusion phenomena. In order to check the accuracy of the pair potentials (2.16) that are built within pseudopotential formalism previously described, we used molecular dynamics simulations to solve the classical equations of motion for a system of N atoms interacting via these potentials. So, provided atomic positions ( ) t r α r and velocities ( ) t v α r at time t for each particle « α », we can examine the atomic structure and transport properties for liquid metal. First we check the pair-correlation functions ( ) r g computed as the time average [START_REF] Allen | Computer Simulation of Liquids[END_REF]:

( ) ( ) ( ) ∑ = + = N 1 α 2 α ∆r r 4π ∆r r r, ∆N ρ 1 r g (3.1)
Above, α ∆N stands for the number of the atoms in the distance between r and ∆r r + from the atom « α », and henceforth ρ will be the average number density at given temperature T. The corresponding static structure factor ( ) q S is obtained through the Fourier transformation [START_REF] Waseda | The Structure of Non-crystalline Material: Liquid and Amorphous Solids[END_REF]:

( ) ( ) [ ] ( ) r d r q i - exp 1 r g ρ 1 q S 3 r r r • ∫ - + = (3.2)
In the methods of computer simulation, one considers a periodic system of cell side L. To be consistent with that, the r-range values for the calculated radial distribution functions ( ) r g is only significant up to 2 L . So, its truncation in the Fourier transformation leads to spurious oscillations in ( ) q S at small q-values. Such effects can be problematic if we wish to extrapolate the long-wavelength limit ( ) 0 S . It is well-known that this limit relates the atomic structure to the isothermal compressibility [START_REF] Egelstaff | An Introduction to the Liquid State[END_REF][START_REF] Hansen | Theory of Simple Liquids[END_REF]. So, instead of equation (3.2), the structure factor can be directly computed from canonical averages over the successive atomic configurations of the atomic positions { } α r r that are generated by MD simulation and over all q r vectors of the same magnitude. Explicitly we have:

( )

q q ρ ρ q S r r - = with [ ] ∑ = = N 1 α α q r . q i exp N 1 ρ r r r (3. 3)
Unfortunately, from equation (3.3) we observe that only the reciprocal lattice vectors

( ) L / π 2 n , n , n z y x
are accessible and meaningful [42]. So that, the smallest q-value is L 4π . This indicates that the density fluctuations cannot be examined over distances larger than 2 L/ .

We now consider succinctly the atomic transport properties without using more fundamental theories.

( ) ( ) ( ) ∑ = - = N 1 α 2 α α 2 0 r t r N 1 t ∆r r r (3.5)
The self-diffusion constant D can also be extracted from the recorded atomic velocities ( )

t v α r
at time t, by integral over the velocity autocorrelation function (VAF) ( ) t Z defined as follows [START_REF] Egelstaff | An Introduction to the Liquid State[END_REF][START_REF] Hansen | Theory of Simple Liquids[END_REF]:

( ) ( ) ( ) ∑ = ⋅ = N 1 α α α 0 v t v 3N 1 t Z r r (3.6)
So that we have the Green-Kubo like equation for the self-diffusion constant:

( )dt t Z D 0 ∫ ∞ = (3.7)
The equipartition theorem implicates that ( )

M T B k 0 Z =
where M, k B and T are the particle mass, the Boltzmann constant, and the temperature, respectively. Hence it is often convenient to define the normalized-VAF:

( ) ( ) ( ) 0 Z t Z t Z N = (3.8)
The velocity autocorrelation function on an equal footing with it corresponding spectral function ( )

ω Z almost provides essential information upon the diffusion process in liquid metal. It is obtained from Fourier transformation with respect to time:

( ) ( ) ( )dt t iω exp t Z ω Z - = ∫ +∞ ∞ - (3.9)
The integral over ω is simply M T k 2π B and the self-diffusion coefficient can be obtained as the zerofrequency limit:

( ) 0 Z D = (3.10)
Present standard molecular dynamics calculations using Verlet algorithm [START_REF] Allen | Computer Simulation of Liquids[END_REF] have been carried out for a (NVT)-ensemble with 4000 atoms enclosed in a cubic supercell with periodic boundary conditions.

The cubic cell size L is chosen so that it corresponds to the observed number density ρ at the given temperature (L=44.8 Å corresponding to ρ= 0.04458 atoms/Å 3 at T=463K for lithium, L=54.5Å corresponding to ρ =0.02486 atoms/Å 3 at T=378 K for sodium, and L=68.1 Å corresponding to ρ=0.0 1276 atoms/Å 3 at temperature T=338 K for potassium). The ionic pair potentials that we have used have been truncated at the cut-off radius close to L/2. In order to avoid any memory effect, we have carried out an annealing simulation by checking several MD-runs at different high temperatures (up to 5000K). ) is large enough to ensure accurate calculations of the static liquid structure and transport properties. We thus obtain reliable pair correlation function ( ) r g and from equation (3.2) accurate static structure factor ( ) q S . We sometimes use equation (3.3) in order to check the accuracy upon ( ) q S at long wavelengths. In this case, molecular dynamics calculations for a (NVT)-ensemble are performed with only 864 atoms. These calculations with 864 atoms are relatively accurate to study the diffusion properties.

IV. Results and discussion

A. The static ionic structure versus interatomic pair potentials

The accurate interatomic pair potentials for Li, Na and K liquid metals at temperature near the melting point (463 K, 378 K and 338 K respectively) are derived from the electronic structure calculations by using Shaw's nonlocal and energy-dependent OMP-model potential. The bare ion model is suitably screened by the dielectric function, taking account of its nonlocal nature and its energy-dependence, so that, among various forms of the local field-correction seen in the literature, we retain the one having a physically well-based background proposed by the Ichimaru and Utsumi. Hence, it appears clearly from our preliminary investigations that for alkali metals the ionic structure and the atomic transport properties depend moderately on the valence-exchange and correlation potential. We can draw the same conclusion with regard to the precise form of the depletion charge distribution. The latter has a weak impact on our results. As it, has been said earlier (see section II), we have considered it uniformly distributed within a spherical shell with radius C R On the other hand, the manner to deal with the "core shift" in order to obtain the well depths ( )

E A l
of the bare ion model potential is crucial (see section II). These one are calculated self-consistently (OMP-Self) or evaluated according to either the Animalu-Heine original approach (OMP-AH) or by one or the other of the two Ballentine-Gupta corrections to Animalu-Heine Scheme (OMP-BG1, OMP-BG2).

We note that the shape of the interatomic pair potential is quite sensitive to the magnitude of energydependent well depths ( )

E A l
. This statement for the liquid sodium case is depicted in fig. 1. From this figure, it is fair to say that both approaches denoted OMP-Self and OMP-AH, seem to be reasonably convergent. Comparatively, OMP-BG1-derived pair potential is slightly shifted toward the smallest inter-particle distances. However, this situation is not dramatic and all three approaches lead to calculated pair correlation curves ( ) r g (fig. 2) and to ones of the corresponding static structure factor ( ) q S (fig. 3) that are almost identical. One point deserves noticed is that these theoretical results successfully predict the experimental data [START_REF] Waseda | The Structure of Non-crystalline Material: Liquid and Amorphous Solids[END_REF]43,44]. On the other hand, OMP-BG2 calculations reported in fig. 1, yield pair potential that exhibits a marked shift of the first minimum to smaller interparticles distances. Consequently, it has an unrealistic much smaller diameter of the repulsive core. 12 So, the large discrepancies between predicted ionic structure from OMP-BG2 and the experiment are not surprising (figs. 2 and 3). To be compared, the OMP-Self-derived interatomic-pair potentials for Li and for K are plotted in fig. 4 together with the Na-one. Apart from a scaling factor, they are not significantly different. As one would expect, they show overall features that are typical of alkali metals. They consist of a strongly repulsive core at short distances plus an oscillatory tail at intermediate and long distances. These so-called "Friedel oscillations" with much damped amplitude at large r, are due to the logarithmic singularity of the dielectric function (2.8). It is worth to note that all potentials exhibit a deep negative minimum around the nearest neighbor typical distance: (Li :

3.16Å min R = , meV 74 min V - ≅ ; Na : 3.68Å min R = , meV 44 min V - ≅
; and for K :

4.52Å min R = , meV 40 min V - ≅
). The similarity of the potentials is reflected in the pair-correlation functions ( ) r g (for the needs of clarity, those of Li and K are not reported in this paper), whereas the respective static structure factors ( )

q S
are represented in figs. 5 and 6. Both functions have the main features of the hard-sphere-like structure. All curves also show an overall good agreement between our theoretical results and experimental ones [START_REF] Waseda | The Structure of Non-crystalline Material: Liquid and Amorphous Solids[END_REF][43][44][45].

B. Some properties of diffusion

Our results relative to the normalized VAF's for the three alkaline metals (Li, Na and K) in the liquid state are displayed in figs. 7. We may note at once that for all elements the main feature of the curves is characteristic of dense fluids, i.e. a damped-oscillatory behavior with negative regions. Negative regions are meaningful of thermal vibrating motion of an ion surrounded by others ions ("caging" effect). The picture is the one of a tagged particle which come into collision with near neighbors so that its velocity is, on average, in the reverse direction (backscattering). The magnitude of the oscillations so quickly decreases that only the first and the second oscillations both are apparent.

Besides, with the increase in temperature the damping is more important. It can even be so pronounced that the normalized VAF becomes almost monotonically decreasing without the first negative minimum. The example of potassium at 800 K illustrates this fact. In addition, one can notice that the magnitude of the first minimum is sensitive to the atomic mass. As regard the sequence (Li, Na, K), the depth is decreasing with the mass of diffusing atom. Now, let us consider the corresponding spectral density functions ( ) 13 diffusion constants, as we derived from the integration of VAF's (Equation 3.7), are reported in fig. 9 and summarized in Table II. Indeed, except some systematic deviations due to the computational conditions and that are upon the whole negligible, the D-values thus-obtained are consistent with those extracted from the Einstein relation (Equation 3.4) or founded upon the zero-frequency limit of spectral function ( ) 0 Z (Equation 3.10). The values of D calculated for Li, Na, and K at temperatures 463, 378, and 343 K, respectively are explicitly reported in table II and compared with experimental results [START_REF] Balucani | Dynamics of the Liquid State[END_REF]. In regard to the variation over a wide range of temperatures, it is convenient to fit the numerical results by considering Arrhenius' type behavior of the form ( )

      - = RT Q exp D T D 0 ; R
being the perfect gas constant. To compare our results to other calculations and to available measurements, we plot in fig. 9 ( )

D Ln
as a function in inverse temperature, so that the preexponential factor 0 D and the activation-like energy Q can be extracted (see table II). Broadly speaking, our results agree well with measurements [START_REF] Egelstaff | An Introduction to the Liquid State[END_REF][START_REF] Balucani | Dynamics of the Liquid State[END_REF] on the entire temperature-range that we have considered.

V. Summary

Through the theoretical results reported in this paper, we intended to understand and to describe both

(

  Unless explicitly stated otherwise, later on we use atomic units throughout:

  the nominal valence and 0 l is the highest value of l -angular momentum in the core ion model potential and of the angular momentum projection operator l P ˆ are:

  the exchange energy per electron is well known to be

FE.

  Nevertheless, one can note that in all respects Cowley assumption is arbitrary. Upon the whole, the previous procedures requiring the knowledge of the mean ionization energy (MIE) and of the binding energy per electron (BEE) yet has drawback. Unfortunately, there are no available data of BEE for any thermodynamic state. We believe that the values tabulated by Ese and Reissland[35] and those calculated by Cowley concern solely the solid state. Thus the evaluation of the Fermi energy on the absolute scale and of the core shift remains an intricate problem. To overcome this difficulty and following upon previous pioneer work ofTaut and Paasch [36], Hallers et al.[37] have devised a procedure in which OMP-model parameters

  consider the bare model potential for an ion implanted in a liquid metal. As it is previously emphasized, in the metal environment its parameters must be evaluated at the valenceenergy k E r shifted by k ∆ r namely at k k ∆ E r r -

9 III.

 9 Structural and atomic transport properties: molecular dynamics simulation

1 .fig. 2 . 3 .fig. 4 . 5 . 6 .

 123456 Figure captions

8 .

 8 fig. 9. Self-diffusion coefficients as a function of temperature calculated with the full-nonlocally screened OMP-model potential with self-consistent core shift (OMP-Self). The open and solid triangles correspond the experiments values for Na[START_REF] Meyer | [END_REF] and for K[49], respectively. (a) represent previous calculations with other potentials for Li[START_REF] Hsu | Proc. Nat. Acad. Sci[END_REF], Na[START_REF] Hsu | Proc. Nat. Acad. Sci[END_REF], and K[START_REF] Hsu | Proc. Nat. Acad. Sci[END_REF]. (b) correspond to our calculations. Note that Ln(D) is plotted versus 10 -3 / T

  we have considered that at first approximation, the "core shift"

										k ∆E r doesn't depend on
	energy level. Hence one need the Fermi energy level F E and the associated "core shift" F ∆E to be
	estimated in order to calculate the shifted Fermi energy	E	F *	=	E	F	-	∆E	F	. To do this, Animalu and
	Heine gave an explicit but crude approximations for F E and

F ∆E :

  We used at best this degree of freedom assuming that the depletion hole is uniformly distributed within a spherical shell whose radius C

			n d	and its Fourier transform dq n ,
	expression (2.13) for dq v	arises from the Poisson theorem. Unfortunately, model-potential theory
	does not provide	( ) r n d .

  To obtain accurate MD calculations, the equations of motion are integrated with MD-

			11
	time steps less than 1 fs. The number of such steps (	N steps =	36000
	F o r
	P
	e
	e r
	R
	e
	v i e
	w
		O n l
			y

Table II .

 II T(K)D (× 10 -4 cm 2 s -1 ) D 0 (× 10 -4 cm 2 s -1 ) Q (KJ /mol)

	Li	463	0.507	13.603	13.323
			0.61-0.68 a		12.442 b
	Na	378	0.436	10.442	10.182
			0.406-0.435 a		10.137 b
	K	343	0.398	10.719	8.762
			0.359-0.376 a		9.218 b

a Taken from Ref.

[START_REF] Balucani | Dynamics of the Liquid State[END_REF]

b Taken from Ref.

[START_REF] Egelstaff | An Introduction to the Liquid State[END_REF]

Page 32 of 32 http://mc.manuscriptcentral.com/pm-pml

Philosophical Magazine & Philosophical Magazine Letters

http://mc.manuscriptcentral.com/pm-pml Philosophical Magazine & Philosophical Magazine Letters