Clément Calenge
email: clement.calenge@oncfs.gouv.fr

Matthieu Guillemain

Michel Gauthier-Clerc

Géraldine Simon

A new exploratory approach to the study of the spatio-temporal distribution of ring recoveries: the example of Teal

Keywords: product kernel algorithm, spatial distribution, time effect, Eurasian Teal, ring recoveries

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

Introduction

The ringing of migratory birds has long been used to study the the movements and spatial distribution of a bird population, e.g. to determine where individuals of a given breeding population wintered or, conversely, where the individuals of a known wintering population were breeding (e.g. [START_REF] Schüz | Atlas des Vogelzugs nach den Beringungsergebnissen bei palaearktischen[END_REF], see also the work carried out by HCC Mortensen, compiled by Jespersen and Veden Tåning 1950), or to define their flyways (e.g. [START_REF] Scott | Atlas of Anatidae populations in Africa and Western Eurasia[END_REF].

Purely spatial smoothing techniques have been used to explore the continuous variation in space of the probability of recoveries (Wand and Jones, 1995;[START_REF] Atkinson | Migratory movements of waterfowl: a web-based mapping tool[END_REF], but accounting for the time dimension in the analysis of such data may bring essential information on the temporal changes in the spatial distribution of the ringed animals. Recently, [START_REF] Keating | Modelling utilization distributions in space and time[END_REF] proposed an approach relying on the kernel method to smooth a distribution of points in both space and time. They developed the method to estimate the probability density distribution of presence of an animal monitored using GPS in both time and space, using relocation data. In the present paper, we test this method to assess if it could also be used valuably as an exploratory approach of the distribution of bird ring recoveries. We apply the technique to the very large Teal ring recovery dataset as an example, to assess if this can provide new insights into the migration patterns of this species throughout Europe. Our main aims were 1) to assess how many different routes these ducks use during spring and autumn migrations, 2) the extent and the rate at which they disperse from their ringing site during a given winter, 3) the consequences of known switches between wintering quarters from a winter to the next (abmigration) in terms of routes used during the migration journeys, and 4) if some areas are used markedly less than others by the species, at the continental scale.

Methods

Teal were caught at Tour du Valat in Camargue, southern France (43˚30' N, 04˚40'E) between January 1952 and February 1978 using standard dabbling duck funnel traps hidden in the vegetation [START_REF] Bub | Bird trapping and bird banding. A handbook for trapping methods all over the world[END_REF]. In total, 59,187 Teal were captured and ringed, of which most (96%, N=56,697) were so between November and March of each year, and constituted the dataset used in the present study. These ringing operations led to 9,114 recoveries (15% of the total number of ringed Teal), among which 4,652 (51%) were intra-annual recoveries, i.e. birds recovered before the month of November following the season of capture, and 4462 were inter-annual recoveries, i.e. birds recovered after this date. The different age and sex classes were not distinguished in the present exploratory analyses of this ring recovery dataset.

Each ring recovery is characterized by a set of spatial coordinates (X and Y), as well as a date of recovery. The whole set of recoveries therefore defines a point pattern in space and time. We aimed at identifying the changes of areas with the highest probability of recovery with time, which was expected to be a proxy for the migration corridor(s) of Teal. For this purpose, we used the product kernel algorithm (Wand and Jones, 1995) recently introduced in Ecology by [START_REF] Keating | Modelling utilization distributions in space and time[END_REF] to draw maps showing the distribution in space of this probability density at a given time t, and to combine the maps drawn for different times t into a movie to show the changes with time of this probability density. The method is described precisely in a document presented as Supplementary Electronic Material. The combination of maps into movies has the advantage, compared to maps, of easily providing the information relative to the speed of movements of animals.

All the estimations were carried out with the function kernelkc from the adehabitat package [START_REF] Calenge | The package adehabitat for the R software: a tool for the analysis of space and habitat use by animals[END_REF] for R software (R Development Core Team, 2008). Note that a subsample of the analysed data (slightly modified, to preserve copyright) is available as an example dataset of the adehabitat package (dataset "teal"), and that the analysis displayed here is reproduced in the "example" section of the function kernelkc help page, to allow the reader to reproduce the computations.

In our estimations of the probability density function (pdf) of recovery, we distinguished the first year following ringing of the bird from subsequent years. Indeed, because all animals were captured and ringed at the same point (i.e. Tour du Valat), it follows that the migration of the first year following ringing had a constrained departure. Furthermore, an earlier study with the same dataset suggested that, at least in terms of probability to switch from a supposed flyway to another (abmigration rate), the behaviour of the birds differed between the year of ringing and subsequent years (i.e. intra-versus inter-annual recoveries), but then did not change with time among inter-annual recoveries [START_REF] Guillemain | European flyway permeability and abmigration in Teal Anas crecca, an anlysis based on ringing recoveries[END_REF]. Therefore: (i) we estimated the pdf of recovery of ringed birds during the first year following the capture using the product kernel algorithm (i.e. intra-annual recoveries); (ii) we estimated this pdf in the years following the first year after the ringing (i.e. pooling all inter-annual recoveries). For the latter estimation, we considered the time as a circular variable.

Results and Discussion

The Teal migration

The two estimated pdfs are displayed in Figures 1, 2 and 3 (the movies are available as Supplementary electronic material). These pdfs clearly show distinct migration routes used by migratory Teal in spring and autumn (Figures 1 and2), which was also suggested by earlier analyses (e.g. Guillemain et al., 2009). Camargue Teal are considered as showing partial loop migration, using a route (South of the Alps) in spring and another one to some extent (North of the Alps) during the autumn, as other migratory birds do (reviews of loop migration studies in [START_REF] Newton | The migration ecology of birds[END_REF]. Interestingly, our visualization method indicates clearly that not all birds migrate via the South of the Alps in spring, even when only intra-annual recoveries are considered, i.e. only those birds which spent at least a part of the previous winter in Camargue are considered (Figure 1). In autumn, too, Teal apparently rely on two distinct routes, one via the Baltic states and the Netherlands, one via central Europe, to reach their winter quarter in Southwestern Europe (Figure 2). It is clear that the route North of the Alps is used to a greater extent in autumn than in spring (see comparison of Figures 1 and2).

Whether apparent loop migration in other bird species may also be the consequence of simultaneous use of multiple migratory routes remains to be tested and could be so with the product kernel algorithm.

A surprising result is that the spring migration of the population terminates earlier for intra-than for inter-annual recoveries (e.g. May 15th maps in Figure 1). On the other hand, autumn migration seems to occur over a long period of time, since some individuals are still in the Northern parts of Europe very late in the season (e.g. Figures 2 and3). This gives some indications regarding a previous hypothesis: analysing the Tour du Valat Teal ring dataset in an independent capture-recapture demographic manner, Devineau (2007) suggested that the Teal population should theoretically be declining (annual population growth rate λ being < 1). Conversely, the analysis of duck counts suggests the opposite both at the national [START_REF] Fouque | Distribution des effectifs hivernaux du fuligule milouin en France et tendances d'évolution sur les 16 derniers hivers[END_REF] and the international [START_REF] Delany | Waterbird population estimates Fourth Edition[END_REF] scales. Devineau (2007) attributed this pattern to the fact that the population of Teal wintering in Western Europe, though heavily hunted, was regularly renewed by individuals arriving late in the season from more Eastern and Northern countries (the population therefore being more "open" than thought in the past), which seems to be confirmed by the present maps. A nonnegligible number of Teal is still present there, e.g. around the Baltic Sea or in some areas of Central Europe even in mid-November.

Abmigration patterns are especially explicit if one compares the autumn migration for intra-and for inter-annual recoveries (Figure 2): from the early beginning of the autumn the northern migration route along the Baltic Coast and the Netherlands is used to a far greater extent for inter-than for intra-annual recoveries, leading to a greater use of the Atlantic and Northern coasts of France in mid-winter (Figure 3). Later on, this translates into a greater use of the same Northern migratory route during spring migration in inter-than in intra-annual recoveries (Figure 1). The high abmigration rate recorded from the very same dataset for Teal, with birds switching from the so-called "Mediterranean" to the so-called "North-Western European" flyways [START_REF] Guillemain | European flyway permeability and abmigration in Teal Anas crecca, an anlysis based on ringing recoveries[END_REF], can now be approached in a more dynamic fashion.

Abmigration in ducks is generally believed to result from the pairs forming in winter, the male following its mate to her breeding area (e.g. [START_REF] Alerstam | Bird Migration[END_REF][START_REF] Anderson | Population ecology of the mallard VI. The effect of exploitation on survival. United States Department of the Interior[END_REF]. The present visualization suggests that abmigration may also be the result of individuals switching to a new target wintering area while migrating in the autumn, not only later while pairing in winter: although relatively few Teal apparently disperse to the French Atlantic coast during the winter of ringing (intra-annual recoveries in Figure 1), and birds from different origins therefore have few opportunities to meet during this first winter pairing period, many birds use the Northern and the Atlantic coasts of France right from their first autumn after being ringed in the Camargue (cf. intra-annual recoveries in Figure 2), suggesting that the switch between wintering areas may occur before pairing.

The examination of the movies did not highlight any geographic area that is completely devoid of recoveries, i.e. avoided or flew-over by Teal during migrations, which accords with the acknowledged plasticity of this species in terms of habitat choice (e.g. [START_REF] Johnson | Green-winged Teal[END_REF]. However, the early spring maps clearly show that the recoveries do not progress gradually to the North-East over time along a clear front but, rather, that many distinct geographical areas get used by Teal simultaneously at the onset of the breeding season (e.g. April and May maps in Figure 1). Though this may reflect discontinuous human settlement distribution in these geographic areas (see below), this is also consistent with the idea that these ducks are not limited by their migratory speed in spring, nor have to wait until environmental conditions get favourable before colonizing their breeding areas (which supposedly occurs gradually with climate at least along a latitude gradient in spring). Rather, birds may try to get advantage from getting to or close to their breeding areas and start to nest as soon as possible, according to [START_REF] Elmberg | Early breeding teal Anas crecca use the best lakes and have the highest reproductive success[END_REF].

Hypotheses underlying the visualization of ring recoveries

It has to be kept in mind that the kernel estimation of the probability density to recover a Teal at a given place and a given date is not the probability density that a Teal is present at that given place. Although we were interested in the latter, our dataset only allowed to estimate the former, like ring recoveries always do. Actually, the former is equal to the product of: (i) the probability density that a Teal is present at this place and this date or a few days before this date; (ii) the probability density that a Teal dies at this place and this date or a few days before this date, given that it is present there; (iii) the probability that the dead Teal is found by someone at this date (which is relatively close to one when the bird is killed during hunting); and (iv) the probability that this person reports the recovery to the national ringing centre.

Although we need to assume that the first probability is adequately reflected by the data to draw reliable conclusions on the migration corridors of Teal in Europe, the other three are likely to vary in both space and time, a problem inherent to the use of ring recoveries in general, including in demographic studies [START_REF] Gauthier | Seasonal survival of Greater Snow Geese and effect of hunting under dependence in sighting probability[END_REF].

The probability that a Teal dies at a given place depends on many factors that vary in both space and time (greater mortality of the Teal during the hunting period as opposed to the breeding season due to the apparent higher effect of hunting versus natural mortality causes, [START_REF] Devineau | A comparison of greenwinged teal survival and harvest between Europe and North America[END_REF]. This may indeed have affected our results, but the fact that hunting seasons were generally very long in European countries during the period when these data were collected (much longer than today, e.g. see [START_REF] Ogilvie | A migration study of the Teal (Anas crecca) in Europe using ringing recoveries[END_REF] suggests that the results may not be too biased, at least at the national scale.

It remains of course that the hunting pressure may be higher in some regions than in others.

The probability that a Teal is found by someone also varies in space and time, and is of course higher during the autumn-winter hunting season than when birds die of natural causes in remote breeding areas with scarce human population. This explains of course why we may have failed to identify a clear migration corridor above the Baltic sea, since recovery probability is zero over seas. The probability that a shot game is killed but not retrieved by a hunter, termed the "crippling loss", has been estimated to be around 20% in ducks (Anderson and Burnham, 1990). However, we see little reasons to believe that the probability to retrieve a killed duck by a hunter varies markedly over the season, or differ strongly between geographic areas. We think that the influence of this non-uniform probability was proabbly smaller that in other studies, as Teal has traditionally been the second most harvested species in Europe after Mallard Anas platyrhynchos (e.g., [START_REF] Mooij | Protection and use of waterbirds in the European Union[END_REF], and because the species was hunted in most European countries (e.g., [START_REF] Lampio | Hunting seasons and methods in Europe 1969-70[END_REF]. In addition, ring recovery rate for Teal was far higher during our study period than it is nowadays (Guillemain et al. Unpub. data, see Robinson et al., 2009, for a similar 25% decline in Teal ring recovery rate per decade since the 1960s).

Finally, the probability to report a ring may also vary in space and time, depending on sociological factors characterizing local human populations (degree of information concerning ringing studies, national hunting rules, etc). As far as we know, ring report-ing rate for ducks has only been estimated in Mallard Anas platyrhynchos in North America, by means of reward ring studies [START_REF] Nichols | Band reporting rates for mallards with reward bands of different dollar values[END_REF]. Such studies have shown that ring reporting rate may well vary geographically [START_REF] Nichols | Geographic variation in band reporting rates for mallards based on reward banding[END_REF].

The extent to which differences in reporting rates between geographic areas (varying reporting rate across the year being less likely) biases the present results remains to be addressed. This remains a general question for the analysis of duck ring recoveries, also for demographic studies, and strongly calls for the development of reward ring studies in Europe (as also suggested by Devineau, 2007).

To conclude, the use of a method smoothing a point pattern in both space and time, such as the kernel product algorithm, allows to take into account the time-related nature of the data in their exploration. It also provides a smoother display of the data than simply movies of the occurrence of ring recoveries (as it "borrows" strength from the neighbouring recoveries in space and time for the estimation of the pdf at a given place and a given date). This gives a smoother aspect to the resulting maps and a more dynamic aspect to resulting movies, an essential feature for visual exploration of data [START_REF] Cleveland | Visualizing data[END_REF]. The kernel product algorithm appears as a promising method to the dynamic analysis of other bird ring recovery datasets. to the smallest area where the probability of recovering a ring, given that it is recovered at this date, is equal to 0.5. The dark grey areas correspond to the smallest area where this probability is equal to 0.7. The light grey areas correspond to the smallest area where this probability is equal to 0.9.

Zusammenfassung

  Ein neuer Forschungsansatz zur Untersuchung der räumlich-zeitlichen Verteilung von Ringwiederfunden am Beispiel von in der südfranzösischen Camargue beringten Krickenten (Anas crecca) Wir benutzen eine neue Methode namens "product kernel estimator", um die räumliche und zeitliche Verteilung von 9114 in der südfranzösischen Camargue beringten Krickenten (Anas crecca), die in ganz Europa wiedergefunden wurden, zu un-tersuchen. Diese Methode ermöglicht es, zeitliche Veränderungen der Wahrscheinlichkeitsdichte, einen beringten Vogel über ein das Untersuchungsgebiet abdeckende Raster wiederzufinden, zu untersuchen. Auf diese Weise haben wir zwei Zugkorridore für diese Art identifiziert, nördlich und südlich der Alpen, während frühere Analysen meist eine von der Camargue nach Süden führende Route vorgeschlagen haben, insbesondere im Frühling. Die nördliche Route scheint in größerem Maße im Herbst als im Frühling genutzt zu werden, was auf eine häufige Abmigration hindeutet (Tiere wechseln von einem Korridor zum anderen). Außerdem scheint die Wanderung der Population im Frühling schneller zu sein als im Herbst. Der gleichzeitige Wiederfund von Krickenten in vielen verschiedenen geographischen Bereichen im Frühling steht mit der Idee im Einklang, dass diese Enten nicht durch ihre Zuggeschwindigkeit eingeschränkt sind. Sie warten auch nicht auf günstige Umweltbedingungen, bevor sie ihre Brutgebiete besiedeln.

Fig. 1

 1 Fig. 1 Sample of maps for intra-(left column) and inter-annual recoveries (right column) of Teal ringed in the Camargue for the spring migration period. The black areas correspond tothe smallest area where the probability of recovering a ring, given that it is recovered at this date, is equal to 0.5. The dark grey areas correspond to the smallest area where this probability is equal to 0.7. The light grey areas correspond to the smallest area where this probability is equal to 0.9.

Fig. 2

 2 Fig. 2 Sample of maps for intra-(left column) and inter-annual recoveries (right column) of Teal ringed in the Camargue for the autumn migration period. The black areas correspond tothe smallest area where the probability of recovering a ring, given that it is recovered at this date, is equal to 0.5. The dark grey areas correspond to the smallest area where this probability is equal to 0.7. The light grey areas correspond to the smallest area where this probability is equal to 0.9.

Fig. 3

 3 Fig. 3 Ring recovery probability by mid-December for intra-(left column) and inter-annual recoveries (right column) of Teal initially ringed in the Camargue. The black areas correspond

Acknowledgements We are most grateful to Luc Hoffmann, Hubert Kowalski, Heinz Hafner, Alan Johnson, and the other people who ringed Teal at Tour du Valat for over 25 years. We would especially like to thank Marc Lutz, Paul Isenmann and the Centre de Recherche sur la Biologie des Populations d'Oiseaux (Muséum National d'Histoire Naturelle, Paris) for their help while computerizing the French Teal database. This work was partially funded by the French "Agence Nationale de la Recherche" (ANR) "Santé-Environnement/Santé au travail", the "European Union's Framework Program for Research and Technological Development (FP6)" and the "Région Provence-Alpes-Côte d'Azur".

Wand MP, Jones MC (1995) Kernel smoothing. Chapman and Hall/CRC, London


J.-P Champeaux

P Çarçabal

M Sence

P Moretto-Capelle

P Cafarelli

A simple "statistical" approach for fragmentation studies of doubly ionized cytosine, thymine and uracil bases

A simple statistical model describing the dissociation of molecular dications into correlated fragment pairs has been developed. This model is based on a combinatory approach in which all possible fragments are enumerated and is refined by taking into account the initial structure of the parent molecule, considering the number of chemical bonds to be broken to give rise to the fragments. We show how this model can be used as a tool to help interpreting experimental results of coincidence experiments. It shows that dissociation of doubly ionized molecules upon 100 keV proton irradiation is dominated by statistical processes but it also enables an easy identification of the dissociation products originating from non-statistical processes requiring further investigation, possibly conveying information on the radiation-molecule interaction itself.

Introduction

Over the last decade, ionization and fragmentation of molecules of biological interest such as DNA and RNA bases [1][2][3][4] , or more complex systems such as nucleoside or clusters of biomolecules [5,6] , submitted to high energy radiation (X or UV photons, ions, electrons), received a renewed interest. One of the motivations for such work is related to the advent of the wide-spread radio-therapy techniques, which prompted the need for a better understanding of the interaction between ionizing radiation and systems of biological relevance at the microscopic (molecular or atomic) scale. These studies focus on the initial physical stage (in the femto-or pico-second range) of the radiationmolecule interaction, and aim at identifying and characterizing its outcome, i.e. product ions and electrons that may lead, directly or indirectly, to tumor cell death.

The description of the ionization and fragmentation processes is a complex task. It has been investigated by Classical Trajectory Monte Carlo method [7,8] and by state-of-the-art quantum mechanical calculations from the description of the initial physical process of electron(s) removal (i.e ionization) from valence molecular orbitals [9] to the molecular relaxation leading to dissociation [10] . Another approach, based on the analysis of fragments size distributions, has also been adopted, bringing evidence for strong similarities between fragmentation of a variety of finite size systems like atomic nuclei, molecular clusters and fullerenes [11][START_REF] Bonasera | Proceedings Similarities and differences between atomic nuclei and clusters[END_REF][START_REF] Campi | [END_REF][14] . For instance, percolation theory has been applied on hydrogen clusters [15] and to the analysis of C 60 fragmentation in ion-C 60 collisions [11,16] , indicating a possible phase transition upon dissociation. In the case of the dissociation of biomolecules studied by mass spectrometric methods, rather complex mass spectra are observed and their interpretation can be challenging. Heinonen et al. [17] proposed a method based on a combinatorial approach in order to help the interpretation of molecular fragmentation.

We present here a simple statistical model based on fragments enumeration combined with bond breaking calculation. The mass spectra simulated with this model can be compared to experimental data and used as a tool to identify, among the many observed dissociation products of complex molecules, those resulting from specific (non-statistical) dissociation pathways. After describing the model, we show how it can be used to help the interpretation of our recent experimental results on the dissociation of cytosine, thymine and uracil (Figure 1), doubly ionized upon collisions with 100 keV protons.

Method

Experimental context

The model we have developed is a tool for the interpretation of experiments of multi-coincidence detection of fragments originating from gas phase bio-molecules doubly ionized by collisions with 100 keV protons. The experimental setup has been described in detail elsewhere [18] . In these experiments, coincidence detection relies on very high detection efficiency and very low detection rate. The density of the proton bunches impinging on the gaseous target molecules is kept very low, so that collision products are detected much less often than once per proton pulse cycle (typically, once every 100 pulses). Consequently, when one singly positive charged product (fragment or parent ion) is detected, it is considered to arise from a singly charged parent molecule. As a consequence, when two singly positive charged fragments are detected during one proton cycle, they are considered as a correlated fragment pair A + /B + originating from the same doubly ionized molecule, and not as two uncorrelated fragments issued from two singly charged molecules (low detection rate) or as two correlated fragments emitted from a triply ionized parent with a third ion product that eluded detection (high detection efficiency). Although we attempted to detect negatively charged products by operating the mass spectrometer in inverse mode, only positively charged products have been detected. Then, the correlated fragment pairs come from the following general reaction path:

M + p + M 2+ * + p + + 2e -A + * + B + * + {C} 0 * + p + + 2e -
The stars indicate that the fragments are not necessarily in their electronic ground state, and {C} refers to the neutral undetected fragments. This correlation treatment leads to experimental insights into the dissociation kinematics and in some cases, into the dissociation dynamics [18][19][20] .

Statistical model

The idea underlying our statistical approach is to consider that all possible pathways leading to any correlated A + /B + fragment pair can be explored by the system with an equal probability. This implies that all possible initial states can be populated with the same probability, independently of the amount of energy deposited in the molecule upon the collision process. We also assume that all rearrangements of the molecule or the fragments during dissociation are possible and equally probable. These assumptions imply that parent molecules have enough internal energy after double ionization to overcome all possible potential barriers in the dissociation and rearrangement pathways. Experiments performed at different collision energy (in the range 25-100 keV) support this assumptions as they revealed similar fragmentation mass spectra, independently of the incident protons energy. Therefore, specific dissociation schemes and associated reaction rates, that are often considered in lower energy processes such as photodissociation [21], are not taken into account here. In this approximation, every atom from the parent molecule can be combined to form all possible correlated A + /B + fragment pairs without any initial geometrical or charge density consideration. The only constrains imposed in the model are the initial number of atoms in the parent molecule, the singly charged state of the fragments and the valence of the atoms involved in the fragments. In this context, all dissociation pathways leading to two correlated mono-cation A i + and B j + , from a doubly charged parent molecule M 2+ , can be generalized into the two following cases:

• initial charge separation: M 2+ X + + X + {X-A i } 0 + A i + + {X-B j } 0 + B j + ; (a) • initial neutral(s) emission: M 2+ X 0 + X 2+ {X} 0 + {X-A i -B j } 0 + A i + + B j + ; (b)
where {α} 0 represents one or several neutral species with a total mass of α, and X and X stand for all species that can be possibly composed of a combination of atoms of M (M=X+X). For example in uracil C 4 N 2 O 2 H 4 (M=112 amu), a possible X/X pair is 1/111 with 4 possible distinct H/(U-H) pairs for the first step of the dissociation following (a) or (b).

The valence of each atom is used to verify that X + and X + for case (a), and X 2+ for case (b), are chemically possible. For example, in case (a), X=H 3 O + is acceptable while X=H 4 + is not. Similarly, X 0 = C 4 N 2 O 2 H 3 and X 2+ =H 2+ will not be treated in case (b).

For each A i + , we can compute and enumerate all the possible ways to produce A i + /B j + correlations following initial charge separation (a), or initial neutral emission (b), from the M 2+ parent molecule, by cycling over all possible X/X pairs and, for each X/X, all possible B j + (according to M Ai +M Bj M and atoms valence). Finally, by counting the number of pathways for each A i + /B j + pair we obtain, after normalization to unity, the probability distribution to produce a given A i + cation correlated with all possible B j + cations with respect to the initial equiprobability criteria postulated above. In the case of uracil with A + =H + , such enumeration yields a total number of dissociation pathways of nearly 1.4 million, the major one being the H + /28 + with more than 60000 ways giving rise to these correlated fragments, all redundancies being removed.

Structural extension of the statistical model

Since the statistical model assumes that all fragmentation and rearrangement pathways are equally probable, it is clear that it will largely overestimate some fragment pairs that are highly unlikely to occur for energetic, structural or physical reasons. For example, initial state separation (case (a)) is more probable than initial neutral evaporation (case (b)), simply because of coulombic explosion. Also, in spite of the very high degree of rearrangement required, the formation of fragments like e.g. H 3 + or C 3 + from uracil will have the same probability in this model as the formation of CH or OCNH fragments, made of concomitant atoms.

To circumvent the overestimation of fragmentation pathways involving unreasonably high degrees of rearrangement, one can consider locality criteria in the initial structure of the parent molecule by taking into account the number of chemical bonds that have to be broken to produce the fragments. To achieve this, a Monte Carlo algorithm, named CM(n), has been developed. This code identifies all fragments obtained by breaking a given maximum number of chemical bonds in the parent molecule.

For each event of the Monte Carlo procedure, the number of bonds to be broken is chosen randomly from 1 to n. The major difficulty is then to identify the connected parts constituting the fragments wich is achieved by a recursive procedure based on graph theory.Typically 10 4 events where generated to achieve convergence for cytosine, thymine and uracil, but the maximum number of events can be fixed by user. The number of bonds that are allowed to be broken also introduces, indirectly, the notion of excess energy available in the molecule after ionization. For instance, allowing a maximum number of 4 bond breakage (CM(4)) implicitly requires on average twice as much internal energy than allowing only 2 bond breakage (CM( 2)). This algorithm can be coupled to quantum chemistry codes to obtain the optimized geometries and energies of the fragments but this aspect is not developed here. Analog algorithms, using more accurate bond energy estimates, have been devoted to mass spectrometry analysis like the FiD code [17] .

Using fragments found by CM to restrict the A i and B j possibilities will drastically decrease the number of explored pathways following cases (a) or (b). The drawback of the geometrical extension is that it allows fragments made of initially bound atoms only, and rearrangement is forbidden. It has the advantage of removing highly improbable fragments but, since some rearrangement can physically occur before or during dissociation, this structural constrain may eliminate some physical pathways.

Results and discussion

Energetic considerations

In proton-molecule collision experiments, two major aspects of the interaction are not controlled: the impact parameter (distance between the incident proton and the target molecule) and the molecular orientation. Consequently, the amount of energy effectively transferred by the proton to the molecule, and where this energy is deposited, are not controlled. The only experimental information on the amount of deposited energy can be obtained from the degree of ionization of the target molecules. In the case of DNA/RNA bases, 90% of all events are due to single ionization of molecules (1 cation detected) and the remaining 10% are nearly entirely due to the double ionization (2 correlated cations detected). Higher degrees of ionization exist although they remain marginal (<< 1%).

To evaluate the order of magnitude of the available excess energy after ionization, we have calculated the vertical single, double and triple ionization potentials at the optimized geometry of the neutral parent in its electronic ground state, for the two lowest spin multiplicity of the ions, at the B3LYP/6-311G** level of theory with the Gaussian 98 software [22] . Results are given in Table 1. The minimum amount of energy required to produce dications lies around 25 eV, but the upper limit of the amount of energy deposited in the interaction is not easy to estimate due to the possible excitation of vibrational or electronic states in the ionized molecule. For 100 keV protons, the interaction is dominated by electronic processes and, using the Lindhart-Winter model, we have previously shown for adenine [2] that the number of events with an amount of deposited energy in the molecule greater than 100 eV falls drastically down to zero. Similar results have been obtained for all the bases treated here and the intensity of the high energy part of the spectrum follows an exponential law of the type I≈exp(-E/30), where E is expressed in eV. Nevertheless, since we detect some triple and higher ionization events, it means that the amount of excess energy in the doubly ionized molecules may reach several tens of eV. This relatively high excess of energy suggests that the systems may explore many pathways, including those involving potential energy barriers, as supposed by the initial assumptions of the statistical model.

Experimental fragmentation of dications

The 2D correlated time of flight spectrum obtained in 100keV proton collision with uracil is given in Figure 2. Characteristic masses are also indicated. Only one channel corresponding to two-body fragmentation and leading to the formation of cations of mass 43 amu (HNCO + ) and 69 amu (C 3 NOH 3 + ), is observed. The 2D TOF spectrum is dominated by all the other so-called 'correlation islands' corresponding to dissociation channels where two product ions and one or more neutrals are emitted. For cytosine, thymine and uracil, the most intense fragmentation channels following double ionization is the emission of a proton correlated with another cation. Table 2 summarizes the measured branching ratio of the overall proton emission upon double ionization of these molecules, with respect to the total number of recorded events.

Table2: Measured branching ratios (%) of H + emission from doubly ionized molecules. The number after the molecule name indicates the number of hydrogen atoms in the molecule.

Molecule

Branching Ratio cytosine (5) 20.4 ± 0.2 thymine (6) 24.4 ± 0.2 uracil (4) 27.7 ± 0.2

One can extract partial integrated spectra by 'slicing' the 2D spectrum along one selected mass to obtain the mass spectrum of the fragments correlated to the selected one. Such partial integrated spectra of the fragments correlated to the H + fragment (A i + =H + ) are given in Figure 3 for cytosine, thymine and uracil. These spectra exhibit very similar patterns. We can notice that the most intense channels correspond to B + fragments of mass between 20 and 50 amu and the signal intensity drastically decreases for heavier masses with a cutoff around 70 amu.

Calculated mass spectra for the proton emission channels from dications (A + =H + )

The simulated mass spectra of the fragments correlated to proton emission (H + /B j + ) obtained by applying our statistical model and its structural extension are compared to the experimental spectra in Figure 4. Qualitatively, the experimental spectra are relatively well reproduced by the statistical model. All the observed correlated fragments are found with a satisfactory intensity agreement especially for fragments of low masses, below 30 amu. At first sight, this suggests that most dissociation mechanisms involved in 100 keV proton induced dissociations are driven by statistics.

The cutoff in the H + correlated spectra observed around 70 amu is also well reproduced. It can be understood since the number of channels leading to an heavy fragment correlated to H + (e.g. (U-H) + /H + ) is necessarily lower than the number of possible ways to form H + correlated to intermediate mass fragments. The same conclusion is drawn for the H + /H + correlation, the number of ways to produce such correlation being very low.

In all cases, the statistical curve underestimates the number of events compared to the experiment for the low masses. This is expected since the statistical model overestimates the number of possible fragments by authorizing high degrees of molecular rearrangement, the global intensity of the distribution is then reduced by the normalization. This is slightly corrected if we impose locality criteria using Stat+CM(3) (red triangle in figure 4) that reduces the number of possible fragments. It is clear that for each group of mass peaks, the structural extension eliminates fragments at both ends of the groups. The main reasons for this is that the low mass end must be dominated by fragments resulting from pathways with more bond cleavages than limited by CM and the high mass end must be mostly due to fragments obtained after recombination, that are excluded by CM. Evidently, using this same locality criteria but increasing the number of possible bound breaks to 4 (Stat+CM(4)) in the molecule improves the agreement with the experimental results, as shown in the case of uracil where masses omitted by Stat+CM(3) appear, like fragments of mass 35 amu as well as 65-67 amu in uracil.

Statistical dissociation can be considered as random dissociation of 'hot' molecules where the energy is distributed over highly excited vibrational states of the system and not much is to be learned about the molecule-radiation interaction from such processes. It is then clear that attention must be paid to features of the experimental spectra that are not reproduced by the statistical model or by Stat+CM, as these features must relate to specific mechanisms, and possibly to the quantum nature of the interaction. In Figure 4, one can easily identify three regions where both the statistical model and its structural extension fail to reproduce the intensity of the observed spectra. These correspond to correlated fragment pairs at H + /12 + , H + /[35-40] + and H + /[50-55] + . The underestimation of the carbon peak (H + /12 + ) by the model is unexpected and cannot be interpreted on the basis of the models presented here, a contamination of the experiment due to the heating the oven can not be excluded. In the recent study of 100 keV proton dissociations of Halouracils [20] , it has been shown that fragments of mass 38-40 amu, correlated to the emission of the halogen atom (replacing the H atom in position 5 of Uracil) exhibit relative intensity trends upon halogen atom substitution, suggesting such specific dissociation pathways. The discrepancy between the experimental intensity and the number of pathways leading to these products, as predicted by the statistical approach presented here, may reveal quantum mechanisms increasing their probability.

Other correlated fragment pairs are predicted to exhibit significant intensities by the model but are not observed experimentally. This is the case for the H + /[55-60] + fragment pairs for all molecules studied here. This suggests that these fragments are likely to further dissociate, leading to secondary products in the experimental spectra, or are not formed at all.

Finally, the mass peaks for the H + /16 + and H + /17 + fragment pairs of uracil 2+ also show a striking disagreement with the intensities predicted by the models. These products (most likely, O + and OH + ) could be due to specific, highly efficient dissociation pathways from the uracil dication but we cannot exclude that they could also arise from residual water in the background vacuum of the chamber.

Another example that the statistical model can be used to point out peculiars feature in experimental spectra.

Fragmentation of dications: emission of fragment with mass 43 (A + =43 + )

The mass 43 amu fragment is another major component of the observed 2D correlated spectra [24] .

Figure 5 shows the comparison between the experimental partial integrated mass spectrum associated to mass 43 and the (Stat+CM(3)) calculation for uracil. The agreement between the experimental spectrum and the spectrum obtained from our model is qualitatively satisfactory in the mass range 10-60 amu except for fragments of mass 28 and 40 amu. Also the model completely underestimates the 43 + /69 + channel, which is the most intense correlation island of the 2D correlated mass spectrum (figure 2).

This indicates that these dissociation channels (28 + /43 + , 40 + /43 + , 43 + /69 + ) must be driven by molecular mechanisms with very high cross sections that considerably increase their probability, with respect to their statistical weight. In the case of 43 + /69 + , it can be understood since this dissociation channel corresponds to two-body dissociation that may occur on a very fast time scale, by simply breaking 2 chemical bonds, in a non-statistical regime.

Conclusion

A simple statistical approach based on the enumeration of the dissociation pathways, assuming equiprobability of all pathways has been developed to help interpreting the emission of correlated fragment pairs A + /B + observed in the dissociation of the pyrimidine DNA and RNA bases cytosine, thymine and uracil, after double ionization induced by 100 keV proton collisions.

The method gives a good agreement for low mass correlations and suggests that statistical dissociations are dominant for such high energy collision induced dissociations, where the amount of deposited energy is not controlled. Evidently, this approach cannot account properly for some particular fragmentation channels which may imply quantum effects that counter-balance their low statistical weight. For instance, in the case of the DNA/RNA bases investigated here, the fragment pair for which further investigation could be informative are : H + /12 + , H + /[35-40] + and H + /[50-60] + , as well as 28 + /43 + and [38 + -42 + /43 + ] for uracil.

The statistical approach can then be seen as a powerful tool to identify the products of non-statistical mechanisms that have to be investigated further, possibly conveying information on the radiationmolecule interaction itself. This can prove very helpful to guide and narrow down sophisticated and computationally expensive studies using ab-initio quantum chemistry or TDDFT approaches [10,25] .

Figure 1 :

 1 Figure 1: cytosine, thymine and uracil molecule

Figure 2 :

 2 Figure 2: 2D Correlated time of flight spectrum obtained from dissociation of Uracil dication. The grayed area highlights the 'slice' of the H + /B j + correlations events.

Figure 3 :

 3 Figure 3: Partial integrated mass spectra of fragments correlated to H + emission following dissociation of cytosine, thymine and uracil dication. The black triangles indicate the mass of the parent molecule.

Figure 4 :

 4 Figure 4: Comparison between calculated spectra associated to proton emission with experimental partial integrated spectra (black line) for cytosine, thymine and uracil. The pink spectra have been obtained with the statistical model alone, red triangles (resp, blue circles) are fragments distribution obtained by the structural extension with a maximum of 3 (resp. 4) bond breakage (Stat +CM(3), resp Stat+CM(4)). All spectra are normalized to unity. The black triangles indicate the mass of the parent molecule.

Figure 5 :

 5 Figure 5 : Comparison between the experimental partial integrated mass spectrum (black curve), the global stat distribution (pink curve), the Stat+CM(3) distribution (red triangles) and Stat+CM(4) (blue circles) of fragments correlated to the fragment cation of mass 43 in Uracil.

Table 1 :

 1 Vertical ionization potentials (in eV) for single (IP1), double (IP2) and triple (IP3) ionization of uracil, thymine and cytosine, calculated at the B3LYP/6-311G** level of theory. The numbers in parenthesis indicate the spin multiplicity of the ion.

					Italic values in the IP1 (2) column are
			experimental values from reference [23]
		IP1 (2)	IP1 (4)	IP2 (1)	IP2 (3)	IP3 (2)	IP3 (4)
	uracil	9.3 9.5	13.3	25.8	24.9	47.2	46.8
	thymine	8.9 9.1	12.9	24.4	24.0	45.2	45.2
	cytosine	8.6 8.9	12.2	24.2	23.0	44.5	44.2