Fabrizio Cleri
email: fabrizio.cleri@univ-lille.fr

Fabio Landuzzi

Ralf Blossey

Mechanical evolution of DNA double-strand breaks in the nucleosome

Double strand breaks (DSB) in the DNA backbone are the most lethal type of defect induced in the cell nucleus by chemical and radiation treatments of cancer. However, little is known about the outcomes of damage in nucleosomal DNA, and on its effects on damage repair. We performed microsecond-long molecular dynamics computer simulations of nucleosomes including a DSB at various sites, to characterize the early stages of the evolution of this DNA lesion. The damaged structures are studied by the essential dynamics of DNA and histones, and compared to the intact nucleosome, thus exposing key features of the interactions. All DSB configurations tend to remain compact, with only the terminal bases interacting with histone proteins. Umbrella sampling calculations show that broken DNA ends at the DSB must overcome a free-energy barrier to detach from the nucleosome core. Finally, by calculating the covariant mechanical stress, we demonstrate that the coupled bending and torsional stress can force the DSB free ends to open up straight, thus making it accessible to damage signalling proteins.

Introduction

Double-strand breaks (DSB) in the double helix of the DNA molecule are defined as the cleavage of the phosphate-sugar backbone on both sides, the two cuts being comprised within 10 base pairs (bp) at most. Such an occurrence is only one among the many types of DNA lesions that a cell suffers at any time, [START_REF] Ward | DNA damage produced by ionizing radiation in mammalian cells: identities, mechanisms of formation, and reparability[END_REF][START_REF] Cadet | Oxidatively generated damage to the guanine moiety of DNA: mechanistic aspects and formation in cells[END_REF][START_REF] Roos | DNA damage and the balance between survival and death in cancer biology[END_REF] aside of single-strand breaks (SSB), base loss (AP site, removal of one purine or pyrimidine), base crosslinking or dimerization, various oxidative defects by reactive oxygen radicals. However, although being much less probable than most other types of lesions, and with relatively fast repair kinetics, [START_REF] Schipler | DNA double-strand-break complexity levels and their possible contributions to the probability for errorprone processing and repair pathway choice[END_REF] DSBs stand out as the most critical lesions to the DNA, since they ultimately lead to chromosome breakage and genome instability, cell mutation, or apoptosis. [START_REF] Iliakis | Mechanisms of DNA double strand break repair and chromosome aberration formation[END_REF] Moreover, several simpler lesions of the type mentioned above, can evolve into SSBs and DSBs upon further chemical processing, both during the repair process, and because of the interaction with other major nuclear proteins.

Because of their cytotoxic effectiveness, inducing DSBs in the DNA of malignant cells is one of the major objectives of chemo-and radiotherapy of cancer. Many powerful antitumor antibiotics, such as the enediyne C-1027, abstract hydrogen atoms at several C 0 sites in the backbone ribose, initiating an oxidative chain that leads firstly to SSBs (mainly at adenylate and thymidylate residues), and then to DSBs, typically cleaved at a distance of 1-2 bp. [START_REF] Xu | A single binding mode of activated enediyne C1027 generates two types of double-strand dna lesions[END_REF] High-energy radiation creates swarms of ionization products, both directly on the DNA structure and, most importantly, on the surrounding water molecules. [START_REF] Fielden | The early effects of radiation on DNA[END_REF] The free radicals produced in the process can attack the backbone and induce many different lesions, often clustered over short distances. Notably, DSBs are produced by ionizing radiation with a relatively high probability, and the terminations at 5 0 and 3 0 strand ends are typically more complex than for DSBs produced by enzymatic cutting, making their repair more complex and error prone. [START_REF] Schipler | DNA double-strand-break complexity levels and their possible contributions to the probability for errorprone processing and repair pathway choice[END_REF][START_REF] Wyman | DNA double-strand break repair: all's well that ends well[END_REF] The DNA of eukaryotic genomes is packaged into arrays of nucleosomes, which appear as cylindroid particles that make up the chromatin structure. About three quarters of the total nuclear DNA are included in the nucleosomes, with the remaining DNA acting as "linker", in a sort of beads-on-astring assembly. The nucleosome core particle consists of 147 bp of DNA tightly wrapped around a histone protein octamer containing two copies of four histone proteins (H2A, H2B, H3, and H4). The lysine-rich N-terminal tails of the histones extend from the protein core, making various contacts with the DNA minor groove (see Figure 1

below).

How the complex structural environment of chromatin is altered in the presence of DNA lesions is a longstanding question in the study of the cellular response to DNA damage. [START_REF] Kruhlak | Changes in chromatin structure and mobility in living cells at sites of DNA double-strand breaks[END_REF][START_REF] Bakkenist | Initiating cellular stress responses[END_REF] Incorporating a lesion within a nucleosome core particle introduces several features not present in naked (linker) DNA, which affect its reactivity. Firstly, wrapping the DNA around the octameric core introduces mechanical heterogeneity into the duplex, resulting in regions that are bent and/or in which base stacking is deformed. Secondly, the large number of lysine (Lys) and arginine (Arg) residues present in histone proteins (more than 22% on average, Lys being especially abundant in H2B, where it represents 16% of the total) can directly interact with the lesion; in particular, Lys side chains are directly involved in AP cleavage within nucleosomes, via Schiff base formation. [START_REF] Sczepanski | Rapid DNA-protein cross-linking and strand scission by an abasic site in a nucleosome core particle[END_REF] Furthermore, the whole nucleosome structure becomes less compact: histones at DSBs are susceptible to extraction in low salt, [START_REF] Xu | The p400 ATPase regulates nucleosome stability and chromatin ubiquitination during DNA repair[END_REF] implying a weaker interaction between DNA and histones at DSBs; and biophysical studies demonstrate that presence of DSBs lead to a localized chromatin expansion. [START_REF] Kruhlak | Changes in chromatin structure and mobility in living cells at sites of DNA double-strand breaks[END_REF] Computer simulations are increasingly demonstrating their utility, by allowing deeper analysis of experimental data, also in conditions where experiments are difficult to carry out. Molecular dynamics (MD) simulations of nucleic acids are today capable of following the system evolution over length and time scales approaching the real experimental set up. [START_REF] Pérez | Frontiers in molecular dynamics simulations of DNA[END_REF][START_REF] Beveridge | The ABCs of molecular dynamics simulations on B-DNA, circa 2012[END_REF] We recently completed a first MD study of the mechanical evolution of SSBs and DSBs in random-sequence DNA oligomers taken as representative of the exposed "linker" DNA between two nucleosomes in the chromatin fiber.. [START_REF] Landuzzi | Stability of radiation damaged DNA after multiple strand breaks[END_REF] We studied the mechanical response under tensile force of SSBs and DSBs with different spacing between the two strand cuts (or "nicks"). The results indicated that the absolute values of force necessary to break up a DSB-damaged, free DNA fragment can be very large, of the order of 100 pN, at elongations of ⇠20%. Such values of longitudinal stress and strain are unlikely to be observed in the normal dynamics of chromatin, nor during chromosome mitosis. Most importantly, however, that study demonstrated that thermal fluctuations are unable to provide the energy necessary to overcome the barrier to rupture, unless the two DSB cuts are separated by 2-3 base pairs at maximum.

In the present work we turn to investigating the mechanical evolution of DSBs in a nucleosome immediately after the backbone breaking event, by using very-large-scale MD simulations in the microsecond time scale. All-atom MD simulations of the nucleosome have started a few years ago, initially restricted to a 10-100 ns time scale, [START_REF] Ponomarev | Relaxation dynamics of nucleosomal DNA[END_REF][START_REF] Biswas | Role of histone tails in structural stability of the nucleosome[END_REF] and very recently extended to the microsecond time scale; [START_REF] Shaytan | Coupling between histone conformations and DNA geometry in nucleosomes on a microsecond timescale: atomistic insights into nucleosome functions[END_REF][START_REF] Pasi | Structure and dynamics of DNA loops on nucleosomes studied with atomistic, microsecond-scale molecular dynamics[END_REF] these works provided already a substantial description of many special features of DNA wrapped in a nucleosome, such as the effects of added torsion and bending, histone-DNA contacts, and much other. However, the structure and dynamics of DNA defects of any kind are yet unexplored, in the much wider context of the nucleosome. Here, we search for specific mechanical signatures induced by the DSB, by simulating microsecond-long trajectories of the entire system, embedded in a large box of water and neutralized with point ions (Figure 1). We use different analysis methods to characterize the mechanistic aspects of DSB structural evolution, at various positions in the nucleosome. Firstly, we look at the long-wavelength thermal fluctuations of the system, extracting the essential dynamics from the covariance matrix of the atomic displacements around the DSB region. Secondly, we determine the lability of the broken-DNA adhesion to the histone octamer, by force-pulling with the "umbrella sampling" method. Finally, we characterize the role of mechanical bending and torsional stress, in determining the evolution of the broken DNA ends at longer times.

Altogether, these analyses allow to trace a mechanical path of the evolution of individual strand breaks into a fully- developed DSB, up to the final fracture event, as well as suggesting the most probable late-stage mechanical evolution of the damaged nucleosome. We conclude that even the weakest DSBs can be resistant up to milliseconds against spontaneous disassembly by thermal forces, thanks to the strong DNAprotein interaction within the nucleosome; mechanical forces of some importance are needed to open up the DNA structure at the break site, in a manner that appears to imply also the intervention of external agents; once the DNA structure starts to be opened, its fate depends on the amplitude of the displacement, the DNA being able to fold back to its original configuration, or to straighten out from the nucleosome core, for a large enough initial opening; also in this case, external forces must help internal stress relaxation, to bring the DNA ends to a sufficiently wide opening. These crucial findings should have profound implications for the early stages of DNA damage detection and repair, for example implying that damage marker proteins (such as Ku70/80, which interacts strongly with broken DNA ends), should also be capable of exploiting complex mechanical actions, for the damaged DNA to be accessible to the repair agents subsequently recruited in cascade.

2/22

Results

DSB dynamics at different nucleosome positions

In our previous study on linker DNA fragments, [START_REF] Landuzzi | Stability of radiation damaged DNA after multiple strand breaks[END_REF] it was found that DSBs can be very stable against thermal fluctuations, unless the two cuts on the backbone are very closely spaced. In particular, we obtained an average bond lifetime of the order of 50 ns at T=350 K for the DSB with a single-bp A•• •T pair, and from these data we extrapolated room-temperature lifetimes of the order of hundreds of milliseconds for a DSB with 2-bp spacing, and up to several hours for a DSB with 3-bp spacing.

Based on such results, we decided to use the most favorable DSB configuration in the present study, in order to increase the probability of eventually observing DNA break up. Therefore, we introduced in all models M1-M4 one single-bp DSB with a central A•• •T, which is the weakest bonded bp. We ran the MD simulations at the temperature of T=350 K, or about 77 C, in order to stimulate the thermal dynamics of the system, while remaining within a range of vibrational excitations that is still meaningful for the molecular force field used. MD trajectories were extended to ⇠1 µs for the M2-M4 models, and up to 1.8 µs for the M1, which displayed some potentially more interesting dynamic features. The reference model O with the intact nucleosome was simulated over a shorter trajectory of 500 ns. Shorter MD trajectories were also run at T=310 K for all models, for comparison.

We firstly present the results for the models M2-M4. For all three, we could not observe any substantial evolution of the DSB into a fully broken DNA, over the whole duration of the simulation, despite the relatively high temperature. While it cannot be excluded that such an event could be produced over longer times, this is an increase of more than a factor of 20 in lifetime compared to the free (linker) DNA; [START_REF] Landuzzi | Stability of radiation damaged DNA after multiple strand breaks[END_REF] once scaled down to 310 K, this translates in spontaneous DSB dissociation times at least in the 100-µs time scale, or longer, for the most favorable (i.e., least bound) DSB configuration, which represents therefore a lower bound for the dissociation time. Representative snapshots from the trajectories at the DSB sites are shown in Figure S 1.

The bonding configuration of the central base pair remains on average rather close to that of the pristine nucleosome, with the H-bonds providing a large fraction of the cohesive energy, and the mildly deformed stacking ensuring a substantial structure stability. An example can be observed in Figure 2, in which the time evolution of the H-bond lengths for the central A•• •T bp of model M3 are shown. The three bonds formed by the N1 (adenine), O1 and O2 (thymine) donors are indicated in red, blue and black, respectively. The relative strength of individual H-bonds in the A•• •T bp can be theoretically estimated [START_REF] Szatyłowicz | Characterizing the strength of individual hydrogen bonds in DNA base pairs[END_REF] to be about 10:4:1 for the N1:O1:O2. The last one is not usually accounted as a true H-bond, since it is very weak and with a length fluctuating around 2.8 Å. Indeed, the central N1 bond remains always in the range [1.8-2.1] Å RMS (note that the simulation is at high temperature); the side O1 is more dynamic than the corresponding bond in normal DNA, with an average length of 2.3 Å (⇠2.05 in normal DNA), and quite large RMS fluctuations due to the larger rotational freedom of the DSB about the central axis; the O2 length remains well beyond the definition of H-bond, fluctuating about an average of 3.2 Å. Overall, these interactions provide enough bonding to keep the DSB in place, even in this M3-dyad position that is the farthest from the histone protein core, among all the DSB configurations studied.

To characterize the dynamic motion of the DNA and of the closest protein residues around the DSB region, we performed for each model a study of the essential dynamics. [START_REF] Amadei | Essential dynamics of proteins[END_REF][START_REF] Narzi | Protein folding pathways revealed by essential dynamics sampling[END_REF] This method of analysis looks at a small subset of collective coordinates of the system, to extract the large-scale, anharmonic movements (bending, torsion, etc.) that dominate the global molecular dynamics.

We firstly perform the analysis for the regions surrounding each location M1-M4 in the pristine nucleosome, model O. Typically, the analysis is restricted to a length of about 7 DNA bp on each side of the DSB, plus the 15-20 histone residues in the closest neighborhood of the DSB. MD trajectories are sampled at a rate of 40 ps 1 . Such analysis of the undamaged system provides a spectrum of eigenvalues, from which we extract the first few significant ones, and an average reference configuration for each M1-M4 site. Then, we repeat the same analysis on each of the independent trajectories including a DSB at the M1-M4 positions, by using as reference molecular structure the corresponding average from model O, so as to highlight deviations from the normal DNA dynamics.

3/22
A key quantity providing information about the large-scale (or "long-wavelength") movements of the fragments implicated in the DSB comes from the study of the first few eigenvectors, and of their root-mean-squared fluctuation (RMSF) on a atom-by-atom basis. These new atomic variables capture the contribution of each group of atoms to the principal collective movements, as filtered out by the most important eigenvectors. For all the M1-M4 models, the first 4 eigenvectors are found to cover 65% of the weight, the 5-15 ones are responsible for another 20%, and all the remaining 3N-15 for the last ⇠15%. Such a distribution is less extreme for the O model, in which large-scale movements are quite more restricted, with the first 15 eigenvalues carrying about 55% of the total weight. The physical meaning of such principal eigenvectors can be appreciated, for example, with the plot of Figure S 3, in which the extreme configurations spanned by the large-scale motion of the first eigenvector, for the DNA fragments in models M1 and M3, are all simultaneously represented; the frames are colored from blue to red, the ordering showing how each atom's motion spans between the extremes of the eigenvector. It can be seen that the principal eigenvector for M3 describes quite homogeneous, local fluctuations of all DNA bases, with just a more evident oscillation along the stacking direction concentrated about the DSB; on the contrary, for the M1 this principal eigenvector describes a dramatic large-scale displacement of the central atoms making up the DSB, which tend to span ample areas across orthogonal planes, by turning about the backbone. This largely different behavior between M1 and the other models M2-M4 is discussed further in the following.

4/22
In Figure 3 we plot the RMSF for the first 4 eigenvectors of each DSB model; each plot compares the RMSF for the fragment of DNA surrounding the DSB (black lines), with the corresponding RMSF of the same fragment intact (red lines). For the M2-M4 models, it can be clearly seen that the RMSF of the DSB fragments is comparable to that of the same fragment in the reference model O; despite local quantitative variations, also of some importance between the various DNA bases, the black and red traces remain always close to each other, for each eigenvector, within a range of 0.1 in the arbitrary units of the RMSF. Moreover, the regions of the DSB and the base-pairs immediately adjacent (indicated by grey shaded areas) do not seem to display a peculiar or specific behavior, compared to the bp more distant from the DSB locations. Only the 1st and 3rd eigenvectors of M4 are somewhat outstanding compared to all the others, since they display an even distribution of displacements among all the bp. As it can be seen in the detailed eigenvector plots in Fig. S 4, this coordinated motion correspond to an ample twisting about the main axis, which exists both for the O and M4 model, therefore independently on the presence of the DSB. Notably, the M4 location is the most "straight" and less perturbed, compared to the rest of the strongly curved and bent DNA in the nucleosome core.

This same analysis for the RMSF of the histone residues closer to the DSB in each model, is shown in Figure 4. Also in this case, for the M2-M4 models it is hard to see a qualitative difference between the data for the intact fragments (red lines), and for the fragments with the DSB inserted (black lines). The lysine and arginine residues are overall more mobile than the others, as far as the 4 principal eigenvectors are concerned, describing a dynamic interaction with the DNA. However, with minor variations, this behavior is the same also in the The Schlitter entropy formula [START_REF] Schlitter | Estimation of absolute and relative entropies of macromolecules using the covariance matrix[END_REF] can be used to estimate an upper limit for the contribution to the free energy from the excess entropy due the presence of the DSB, as:

T DS DSB = T ⇣ hS MX i hS O i ⌘ (1)
with MX = M1,...M4, and h...i indicating the time average of the Schlitter entropy for each molecular fragment:

S = 1 2 k B ln ⇢ det  I + k B Te 2 h2 MC (2)
with C the covariance matrix of the atomic displacements, I the identity matrix and M the mass matrix, having respectively 1 and the atom masses on their diagonals, and 0 elsewhere. Table 1 reports the values for each DSB model, divided into DNA and histone contribution. The absolute DNA entropy S O from Eq 2 fluctuates about 18±0.5 kcal/mol/K for each bp, very homogeneously all along the most part of nucleosome, but increasing to 20 kcal/mol/K in the few terminal bps attaching to the straight segments. If the values of excess entropy of DNA are distributed to the 4 bases (green and red in Fig 1b) comprising the DSB, these correspond to an excess of 35 to 60% for the M2-M4 models, the excess per base being larger in the M4, in agreement with the somewhat larger mobility demonstrated in Fig 3 . On the other hand, the excess entropy for the histone residues selected for this analysis remains relatively small, for the three models M2-M4. Despite some difference in the total masses of the groups selected, even when expressed per unit mass instead of per-moles, the absolute entropy of the histones remains comparable, between the model O and the models including the DSB. This is a further confirmation of the relatively minor role played by histone dynamics in the M2-M4 models.

We now turn to describing the behavior of the DSB in the M1 model. Contrary to our expectations, this location in which the DSB is constrained between the histone core and the mobile H2B tail, and close to a DNA-protein contact, was the one to display the most interesting dynamics. The most evident change in the immediate environment of the DSB is the modification of the H2B tail, which can fold into very different interacting positions, starting from the outward extended conformation of the experimental crystallographic structure.

This behavior is shown in Figure 5, where the arrangement of the H2B tail is represented for three configurations, averaged over the respective MD trajectories: In this M1 model, the DSB is constantly enclosed between the two b -sheets of H3 and H4, which fluctuate about their equilibrium structure and interact with one side of the DSB, while the H2B tail experiences strong oscillations, coupling with the cut bases of the opposite DSB side. The time evolution of the four bases comprising the DSB (green-red colored in Fig 5 ) gives a qualitative appraisal of this strong interaction (Fig. S 2). Notably, the interacting portions of both the two b -sheets, and the H2B tail, include more than 60% of lysine and arginine residues, as expected given the strong electrostatic affinity of such amino acids for DNA (notably for G and T, [START_REF] Suzuki | A framework for the DNA-protein recognition code of the probe helix in transcription factors: the chemical and stereochemical rules[END_REF] ). The DNA ends at the DSB are clearly perturbed by such interactions, and it can no longer be said that the two broken backbones preserve a geometrical continuity, as it was instead observed for the M2-M4 models for the entire duration of the respective MD trajectories.

By looking at the RMS fluctuation of the eigenvalues for the M1 model in Figs. 3,4, it can be seen that in this case the group of DNA bases adjacent to the strand breaks take up the majority of the weight, indicative of their participation in the ample fluctuations of the open DSB ends; also for the histones, it is readily apparent a more dramatic dynamics, especially by the lysine and arginine residues; finally, the values of excess entropies from Table I give a further confirmation of the peculiar large-scale dynamics of this DSB configuration. Because of these indications, we continued the MD trajectory up to 1.8 µs, but never observed a true mechanical destabilization of the DNA structure: the two DSB ends remain firmly in place, even if the two terminal bp on each end fluctuate quite wildly (see again

Free energy to detach broken DNA ends

As shown in the preceding Section, spontaneous dissociation of one or both DSB ends of a broken DNA from the nucleosome remains a difficult event, never observed in our simulations. DSB opening, and DNA detachment from the nucleosome are likely governed by a free energy barrier of adhesion, which even such a critical defect as a fully-cut DNA could not easily overcome simply by thermal fluctuations. The way to estimate the free-energy barrier in such a large and complex molecular system is to resort to controlled-force pulling, in order to impose the detachment, and then to use the intermediate structures along the reaction coordinate as starting points for the "umbrella" sampling of the potential of mean force. From the latter, the free energy barrier along the chosen reaction coordinate can be extracted.

As briefly described above, we used as reaction coordinate z the separation distance between the moving DSB end and the histone core surface. This was measured by taking the center of the DNA axis, at the average position of the C4 0 and P atoms of the last two bps, and projecting it on the closest histone surface atom, along the perpendicular to the superhelical axis.

Figure 6a shows the variation of z as a function of simulation time, at constant pulling force. It can be seen that the DNA broken end detaches from the histone surface in large steps (red segments), during which the internal energy builds up until some barrier is overcome; the final stage, indicated by the blue segment, is the complete detachment of the DSB end after t=1.1 ns, in which the free end is simply drifting at the constant speed of about 2 m/s (later on dropping to 1 m/s).

During the final stage of the pulling simulation, the DNA is forcefully unwrapped from the histone core, as it can be seen in Fig 6b . Here we show the distance from the core surface of three P atoms facing the histones, belonging to the bp 71-118 (contact site close to the DSB), 78-111 (middle site) and 82-107 (next contact site). The first contact site is detached in the interval t=1.-1.5 ns, as indicated by the black trace that follows the distance from the surface of of the P71 backbone phosphor. Then, under the continued pulling of the DSB end, also the P111 comes off, at t >3 ns (red trace); however, it may be noticed that this event is "cooperative", the P82 (blue trace) following the instantaneous opening of P111 at t=3.-3.4 ns, and then falling back into position, after which P111 is definitely "peeled off" the histone surface.

From this force-pulling simulation we can calculate the free energy profile of the barriers, which characterize the binding of the DNA end to the histone core surface. The potential of mean force (PMF, [START_REF] Irving | The statistical mechanical theory of transport processes. iv. the equations of hydrodynamics[END_REF] ) is a method to extract the free energy difference DA from a sequence of configurations, biased along a reaction coordinate that brings the system from a state a to a state b. In our case, the reaction coordinate is just the distance z defined above; the states a, b respectively represent the initial configuration at z =0, with the DSB end still attached to the histone surface, and the final configuration with the end detached, at z ⇠5 nm. The "umbrella sampling" technique, [START_REF] Torrie | Nonphysical sampling distributions in monte carlo free-energy estimation: Umbrella sampling[END_REF] is used to obtain the PMF at discrete values of z , and the discrete values of A(z ) between a and b are connected by the weighted-histogram method. [START_REF] Kumar | The weighted histogram analysis method for free-energy calculations on biomolecules[END_REF][START_REF] Hub | A free weighted histogram analysis implementation including robust error and autocorrelation estimates[END_REF] We extracted 100 configurations from the force-pulling simulation, spaced by 50 ps in the first 5 ns of the trajectory (corresponding to about 0.5 Å spacing along the reaction coordinate z =0 to z ⇠5 nm); each configuration was equilibrated for 2 ns at 310 K under constant-{NV T }, while biased with a harmonic "umbrella" potential of variable strength, progressively reduced to zero to obtain the unbiased limit. The force probability distribu- tion of the fluctuating DSB free-end at each value of z was reconstructed with the weighted-histogram analysis, and the free energy profile thereby extracted is shown in Figure 7.

From the noisy profile, a few features can be identified. The red circle defines the first barriers to the detachment of the DSB ends, corresponding to the red steps in Fig 6a ;  such barriers are quite small (<1 kcal/mol), and strongly depend on the choice of the point of application of the pulling force. The blue circle identifies the free-energy barrier for the detachment of the first contact at P71, about DA=1.8 kcal/mol or 3 k B T ; this does not represent a very large value, and should correspond to a ⇠5% Boltzmann probability of spontaneous detachment at T=310 K. The green circle roughly identifies the cooperative events leading to the detachment of P111, between 2 and 4 ns, with a sequence of DA overall not larger than 2-3 k B T . Further detachment events were not observed, with the above values of pulling force; in particular, P82 remained in place for >500 ns, even at larger deformations of the DSB free end, because of the H3 histone tail acting as a sort of brace that maintains the DNA firmly in place about that position. Much larger forces, or cooperative events of histone tail fluctuation, likely involving other nuclear proteins, seem to be necessary to pull the free DSB end further beyond the limits observed in the present simulations.

Internal stress relaxation and DSB structure

In the last part of our study, we turn our attention to the internal relaxation dynamics of the nucleosome including a broken DNA. To demonstrate what it is meant by "internal relaxation", we take two configurations along the final trajectory of the force pulling simulation of the M1 model described in Fig 6, C180 and C290, respectively extracted at times t=1.8 and 2.9 ns, well beyond the detachment stage that ends at 1.1 ns in the figure. Each of these two configurations is used as initial

8/22

structure for an MD simulation, and is then equilibrated and relaxed at 310 K and constant-{NV T }, without any external forces applied. The results of these two MD simulations are displayed in Figure S 5: starting from the two different initial conditions, after 40 ns the C180 tends to fold back into the initial M1 configuration, while C290 straightens out and increases its distance from the histone core. Notably, the C180 remains in a slightly open state, because of the free energy barrier to detachment that now has to be overcome in reverse. However, the important observation in both cases is that the folding back, or the straightening out, are driven entirely by the competition between the residual attraction between DNA and proteins (a "chemical" force), and the relaxation of internal constraints (mainly bending and torsion, therefore an "elastic" force).

The role of internal forces can be clearly understood by looking at the distribution of mechanical stress, which is a measure of the elastic energy accumulated by the bending and torsion of DNA while wrapping around the histones, and that is ready to be released if the structural constraints are softened, as it could be the case of a DSB cutting the DNA sequence. Recent developments led to alternative geometric derivations of the microscopic stress, [START_REF] Torres-Sanchez | Geometric derivation of the microscopic stress: A covariant central force decomposition[END_REF][START_REF] Tadmor | Modeling materials: continuum, atomistic and multiscale techniques[END_REF] based on the invariance of the free energy with respect to surface deformations, [START_REF] Doyle | Nonlinear elasticity[END_REF][START_REF] Mistura | The definition of the pressure tensor in the statistical mechanics of nonuniform classical fluids[END_REF] instead of the classical formulation based on invariance of momentum. [START_REF] Kirkwood | Statistical mechanics of fluid mixtures[END_REF][START_REF] Lutsko | Stress and elastic constants in anisotropic solids: molecular dynamics techniques[END_REF][START_REF] Cormier | Stress calculation in atomistic simulations of perfect and imperfect solids[END_REF][START_REF] Cleri | Representation of mechanical loads in molecular dynamics simulations[END_REF][START_REF] Subramaniyan | Continuum interpretation of virial stress in molecular simulations[END_REF][START_REF] Admal | A unified interpretation of stress in molecular systems[END_REF] ). The so-called CCFD scheme, [START_REF] Torres-Sanchez | Examining the mechanical equilibrium of microscopic stresses in molecular simulations[END_REF][START_REF] Torres-Sanchez | Geometric derivation of the microscopic stress: A covariant central force decomposition[END_REF] incorporated in a special-purpose version of the GROMACS code, ensures conservation of both linear and angular momentum under a generic stress-induced transformation.

The mechanical stress s (x) (a 3⇥3 tensor defined at any point x in space) is a meaningful way of representing the distribution of internal forces with respect to a given local direction vector. Once a DSB breaks the DNA backbone around the nucleosome, internal forces are going to be relaxed, and compete with the chemical (Van der Waals, electrostatic) forces from the interaction with the histone proteins. Looking back at Fig. S 5 for model M1, such a competition is very evident upon comparing the bottom configurations: in C185 the chemical forces overwhelm the internal stress, whereas in C290 the opposite holds, and the DNA ends up straightened out from the DSB site.

A tensor, such as the stress, can be meaningfully projected onto any direction vector, the choice of a particular projection being just a matter of convenience. In the present case, the "bent tube" structure of nucleosomal DNA makes it interesting to consider the stress projected onto its "tubular" surface. An intuitive way of looking at the mechanical stress as a "projected force" is through the surface traction vector, T(x) = s (x) ⌦ n. [START_REF] Ward | DNA damage produced by ionizing radiation in mammalian cells: identities, mechanisms of formation, and reparability[END_REF] The traction vector T(x) contains a great deal of information on the state of internal tension, compression, and torsion, of a complex structure like the DNA in the nucleosome. The portion of DNA wrapped around the histone core is forced to bend into nearly two full circles of diameter about 8 nm, a size much shorter than the persistence length of free DNA, x p '50 nm. Therefore, the DNA "tube" is here constrained in a geometry from which it should rather escape into a more straight structure, whenever possible, under the relaxation of internal forces. The state of tension and compression of a bent tube is described by a different projection of the traction vector, t(x) = T(x) • t, in this case along the unit vector t locally tangent to the continuous line sweeping the center of the tube.

Notably, a bent tube will experience a stretching (tensile) force in the half tube which lies outside the centerline with respect to the center of curvature, and a compressive force in the half lying inside the centerline, as shown in blue/orange in Figure 1c. The internal force is zero along the centerline itself, because of this called the "neutral axis" (also the helical axis of DNA). We computed the line tension t(x) all along the curved DNA pathlength, by averaging over slices of width 0.5 nm (see for example the white slice in Fig 1c), and by integrating separately over the inner and outer regions (orange and blue in the Figure ). Each slice is centered at the midpoint between the two P atoms of each bp, therefore adjacent slices have some overlap, to provide a smoother profile of the signal.

The twist stress is that part of the internal forces involved in the torsion about the central (neutral) axis of the tube. The DNA double helix is naturally twisted already in its normal B configuration; however, when it is bent in the nucleosome, the twist is necessarily modified with respect to the normal configuration. The twist component is obtained as well from the traction vector, as: w(x) = T(x) • (t ⇥ r), where r is the vector from the neutral axis to the point x, parallel to the local surface normal n (see again Fig 1c). The vector product between t and r defines a third vector threading like a spiral screw about the DNA "tube"; positive and negative values of w(x) indicate a rotational force (a torque) tending to over-or under-twist the DNA about its helical axis.

In Figures 8a and8b we show for both configurations the tension profile t(x) along the helical axis of the DNA fragment right after the DSB, bp 72 to 167, i.e., the DSB end that was pulled out under force and subsequently relaxed; we neglect the first few bp immediately next to the DSB, too disordered for such a calculation. Two sets of data are shown in each panel, at the beginning of the relaxation (black lines), and after 40 ns (red lines); stress values are averaged over 100 frames with 10 ps spacing, in either case. In general, the terminal part of the DNA next to the DSB (indicated by a grey-shaded area in the panels) tends to lower values of both line tension and compression, for both configurations, compared to the rest of DNA beyond the dyad (bp 94-94). The data are rather noisy, such a noise being largely numerical (not arising from the thermal average), coming from the difficulty of properly identifying the local reference frame {n, t, b} at each point x along the fluctuating DNA helical axis (centerline). However, it may still be appreciated that for the C185 configuration, The extra twist stress (positive or negative) also contributes to the internal forces that are going to be relaxed, when the DSB cuts open the DNA, albeit to a much lesser extent, given its absolute value. In Figs.8c and 8d the profile of w(x) is shown, under the same conditions of the two panels above for the line-tension/compression. It can be noticed that, also for the twist stress, generally smaller values (1 MPa in modulus) are seen in the DSB tail comprised in the grey-shaded part. However, the large numerical noise does not allow in this case to draw a more firm conclusion, concerning the (likely minor) role of twist stress in the chemical vs. mechanical force competition in the two configurations.

Discussion

In the present work, we studied by very-large-scale molecular dynamics (MD) simulations the evolution under external force and temperature of double-strand breaks (DSB) in nucleosomal DNA. We collected and analyzed a large amount of raw data (more than 1.5 TBytes, and 5 million CPU hours on two large supercomputers), by running microsecond-long trajectories for 5 different, all-atom models of the experimental 1kx5 nucleosome structure. [START_REF] Davey | Solvent mediated interactions in the structure of the nucleosome core particle at 1.9 A resolution[END_REF] The basic model is made up of the canonical 8 histones, plus a 187-bp DNA comprising the 147 bp wrapped around the histone core and 20-bp terminations on each end, and embedded in large boxes of about 80-100,000 water molecules with Na + and Cl ions at 0.15 M physiological concentration. The pristine nucleosome configuration (model O) was modified, by inserting a DSB at four different positions in the DNA (models M1-M4), and the stability of the resulting structures was compared with model O nucleosome.

A general observation from the µs-long trajectories, is that 10/22 damaged DNA remains well attached to the nucleosome body, without qualitative differences compared to the intact DNA. Only the model M1, in which the DSB is tightly sandwiched between the histone H3 and the tail of histone H2B, displayed a dynamics substantially different from the corresponding region in model O, due to the increased interaction of the broken DSB ends with close-by histone residues; however, also this DSB configuration was stable over the entire observation time scale, which in this case was extended to 1.8 µs. In order to identify the free-energy barriers which maintain the broken DNA attached to the histone core, we carried out steered MD with a pulling force to "peel off" the free DSB end from the nucleosome; relatively small free-energy barriers of the order of 3 k B T were identified, which could allow spontaneous DSB end detachment at physiological temperatures, likely over longer time scales of hundreds of microseconds to milliseconds. Spontaneous unwrapping of DNA from the nucleosome core has been studied experimentally, [START_REF] Protacio | Coupled-enzymatic assays for the rate and mechanism of DNA site exposure in a nucleosome[END_REF][START_REF] Li | Rapid spontaneous accessibility of nucleosomal DNA[END_REF][START_REF] Tóth | Histone-and DNA-sequence-dependent stability of nucleosomes studied by single-pair FRET[END_REF] because of its relevance in gene regulation and DNA transcription; notably, such experiments were carried out on isolated nucleosomes, with a length of DNA just matching, or barely longer than needed to wrap the histone core (147 to ⇠180 bp). In such conditions, spontaneous detachment of the ends was indeed observed over the timescale of hundreds of milliseconds; simulations by coarse-grained MD methods roughly confirm such trends, [START_REF] Ettig | Dissecting DNA-histone interactions in the nucleosome by molecular dynamics simulations of DNA unwrapping[END_REF][START_REF] Voltz | Unwrapping of nucleosomal DNA ends: a multiscale molecular dynamics study[END_REF][START_REF] Kenzaki | Partial unwrapping and histone tail dynamics in nucleosome revealed by coarse-grained molecular simulations[END_REF] despite being strongly dependent on the empirical parametrization of each different force model. To such experiments it may be objected that the nucleosome constrained in the chromatin could have a rather different mechanics: our molecular stress calculations demonstrate that the circularly bent DNA has a strong internal driving force, from the relaxation of line tension and, to a lesser extent, of twist (torsional) stress. The reason may be found in the persistence length of the free DNA, which is much longer (⇠50 nm) than the average radius of curvature in the nucleosome (⇠8 nm), and pushes the DNA to regain the straight average conformation on that length scale; the fact that spontaneous fluctuations were observed [START_REF] Li | Rapid spontaneous accessibility of nucleosomal DNA[END_REF] both in presence of, and without binding proteins seems to support this view. In fact, our µs-long MD simulations were carried out with a soft restraining of the DNA linker (20 bp on each end), to simulate the effect of the background chromatin structure, and no fluctuations larger than thermal vibrations were detected for the terminal phosphors; on the other hand, the DSB free end, once extended beyond a distance of about 2.5 nm away from the histone core, tended to regain a straight conformation and detach completely, confirming the importance of stress relaxation as a main driving force in DNA unwrapping.

In conclusion, the results and potential implications of this study can be summarized by the following findings:

1. The closely-cut DSB remain relatively stable over long time scales, and display no sign of disassembly; interaction of DSB ends with histone surfaces and tails is a main factor in damaged-DNA dynamics. DSB configurations close to histone tails in fact display a more active internal dynamics, with a participation also from histone fragment fluctuations.

2. The free-energy barriers for detachment of DNA from histones are relatively low, of the order of a few k B T , implying that short sections of DNA could spontaneously unwrap over a time scale of >100 microseconds, from DSB broken ends, or from the linker sections at the nucleosome ends, as observed in some experiments. [START_REF] Protacio | Coupled-enzymatic assays for the rate and mechanism of DNA site exposure in a nucleosome[END_REF][START_REF] Li | Rapid spontaneous accessibility of nucleosomal DNA[END_REF][START_REF] Tóth | Histone-and DNA-sequence-dependent stability of nucleosomes studied by single-pair FRET[END_REF] At the same time, histone tails represent a major steric obstacle for unwrapping of larger DNA sections, notwithstanding the driving force from stress relaxation (3. below).

3. Fully-consistent molecular stress calculations on the DNA wrapped in a nucleosome revealed the existence of a strong internal driving force for straightening the circularly bent segments, from the relaxation of line-tension and torsional stress. This might be the main force leading to spontaneous unwrapping of DSB cut ends, as well as of nucleosome ends, opening the way to damage-signalling and repair proteins, and to remodelling factors, respectively. Notably, such proteins must be implicated in complex mechanical actions on the nucleosome, resulting from the competition between internal stresses and chemical forces, very likely both sequenceand position-dependent.

EXPERIMENTAL METHODS

Molecular structures of damaged nucleosome

We obtained the nucleosome molecular configuration from the RSCB Protein Database, entry 1kx5. [START_REF] Davey | Solvent mediated interactions in the structure of the nucleosome core particle at 1.9 A resolution[END_REF] This is an x-ray structure of the entire histone octamer with 147 DNA bp resolved at an average RMS of 1.94 Å, reconstituted from human nuclear extract expressed in E. coli; only 6 histone residues were unidentified in this experimental structure, with respect to the known histone sequences, therefore the model can be considered nearly complete. The 147 bp DNA is a palindromic sequence, chosen to maximize the degree of ordering and increase the x-ray spatial resolution. To obtain a model structure useful for our computer simulations, we removed all the crystallization water molecules and ions from the published structure, and added two DNA extensions of length 20 bp at each end of the nucleosomal DNA, with repeated sequence d(AGTC). [START_REF] Shaytan | Coupling between histone conformations and DNA geometry in nucleosomes on a microsecond timescale: atomistic insights into nucleosome functions[END_REF] DNA bases are numbered from 1 to 187 in each chain, one running clockwise and the other counter-clockwise, the dyad being located at basis 94 of each chain. This pristine nucleosome model without strand breaks is shown in Figure 1a, and will be labelled O in the foregoing.

DNA is wrapped left-handed about the histone core, making two nearly complete turns that join at the dyad symmetry point; the two DNA turns define two circles lying in two ideally parallel planes, with a superhelical symmetry axis perpendicular to the center of the circles (for a thorough discussion of nucleosome geometry and structure, see e.g. Ref. [START_REF] Mcginty | Nucleosome structure and function[END_REF] ). The relaxed DNA double helix makes a complete twist around its double-helical axis, about every 10.4 bp, defining a major and a minor groove; therefore, when turning around the histone core, the wrapped DNA makes 14 nearly full twists. Correspondingly, 14 contact points between DNA and proteins can be identified within the nucleosome structure, loosely situated at the minor groove locations facing inwards.

Based on these geometrical features, we defined 4 potentially interesting sites along this wrapped structure, where to place a DSB in a "mechanically significant" position, labelled 1 to 4 in Fig 1a . Correspondingly, we introduced a DSB at an inner contact 11/22 site (model M1); at an outer non-contact site (M2); at the dyad (M3); and at the entry point of the nucleosome (M4). To create a DSB at each such locations, we introduced 5 0 -OH and 3 0 -phosphate terminations at each end of the break, respectively between2 : bases C69-T68•••A120-G121 in M1; bases C73-A74•••T114-T113 in M2; bases A94-T95•••A94-T95 of both chains in M3; bases T22-A21•••T167-G168 in M4. In this way, the two backbone cuts of each DSB are spaced by 1 bp always comprising an A•••T pair (Fig 1b), which remains initially bonded by only its two hydrogen bonds, plus the stacking interactions on each intact side of the chain, while the other half of stacking is readily reduced, as soon as the MD relaxation starts.

The CHARMM-27 force field database [START_REF] Mackerell | All-atom empirical potential for molecular modeling and dynamics studies of proteins[END_REF][START_REF] Mackerell | Extending the treatment of backbone energetics in protein force fields: limitations of gas-phase quantum mechanics in reproducing protein conformational distributions in molecular dynamics simulations[END_REF] and its extension to treat nucleic acids [START_REF] Foloppe | All-atom empirical force field for nucleic acids: 1) parameter optimization based on small molecule and condensed phase macromolecular target data[END_REF][START_REF] Mackerell | All-atom empirical force field for nucleic acids: 2) application to molecular dynamics simulations of DNA and RNA in solution[END_REF] were used for the molecular bonding and nonbonding force parameters. Strict comparisons between CHARMM-27 and AMBER force fields [START_REF] Spiriti | DNA bending through large angles is aided by ionic screening[END_REF][START_REF] Reddy | DNA polymorphism: a comparison of force fields for nucleic acids[END_REF] ensure that the results of long-time, finite-temperature MD trajectories of nucleic acid fragments with largely different conformations are consistent, and able to correctly reproduce the key structural quantities (bond angles, hydrogen-bond structure, base tilt, twist, shuffle, etc.) compared to experimental data. However, in all cases great care must be taken by performing sufficiently long preparatory annealing cycles of the water and ion background, while keeping the nucleosome still, to obtain the right water density and allow a realistic arrangement of the counter-ions around the phosphate backbone, prior to starting the microsecond production runs.

Molecular dynamics simulations

For the molecular dynamics (MD) simulations we used the GRO-MACS 5.1 computer code. [START_REF] Berendsen | GROMACS: A message-passing parallel molecular dynamics implementation[END_REF][START_REF] Lindahl | GROMACS 3.0: a package for molecular simulation and trajectory analysis[END_REF] Nucleosome models O and M1-M4 were solvated in water box of size 14.5 or 18⇥19⇥10 nm [START_REF] Roos | DNA damage and the balance between survival and death in cancer biology[END_REF] with periodic boundary conditions in the three directions, containing about 82,600 or 110,500 TIP3P water molecules, plus 480 Na + and 250 Cl ions to ensure neutralization of the phosphate backbone charge, and a physiological salt concentration around 0.15 M. All the MD simulations were carried out at the temperature of 310 K and pressure of 1 atm, or 350 K and 50 atm for the thermal stability study. Because of the requirements of stress calculations (see below), we could not use standard Ewald-sum electrostatics but plain cut-off Coulomb forces. This is known to be at the origin of possible artifacts, therefore we adopted for both electrostatics and long-range non-bonding forces an unusually large cut off radius of 1.6 nm. The DNA terminal ends (linker) were restrained by soft harmonic constraints, allowing a fluctuation of ±5 Å, to represent embedding in the chromatin structure. We used rigid bonds for the water molecules, which allowed to push the time step to 2 fs for the thermal equilibration runs, and to 1 fs for the force-pulling simulations. Typical preparatory constant-{NPT } MD runs lasted between 10 and 20 ns; force-pulling simulations were carried out for 10 ns, and the subsequent force-free relaxation lasted up to 400 ns; thermal stability simulations at constant-{NV T } extended to ⇠1,000 ns for O and M2-M4, and up to 1,800 ns for the M1 model. Overall, the study used about 4.2 million hours of CPU time on 2048 IBM BlueGeneQ processors (IDRIS supercomputing center in Orsay), and about 800,000 hours on 896/1064 Broadwell Intel E5-2690 multi-core processors (CINES supercomputing center in Montpellier), with typical running times of 1.3 and 7 ns/hour on the IBM and Intel machine, respectively. About 1.5 Terabytes of raw data were accumulated over a period of 8 months, from March to October 2017, for subsequent post-processing. All-atom microsecond trajectories for the nucleosome (DNA+proteins, minus water and ions), stored in the compact GROMACS-xtc format at frame intervals of 40 ps, are freely available upon request to the authors.

Steered MD and umbrella sampling

Steered molecular dynamics (SMD) was performed on the fragments with the constant-force pull code available in GROMACS, only on the M1 model. In this case, we enlarged the water box to 18 nm in the x-direction, to allow possible outward extension of the broken DNA end, resulting in a system of 107,000 water molecules. Since the objective was to promote the detachment of one of the broken DSB ends from the nucleosome core, we applied a constant force parallel to the direction x and perpendicular to the superhelical axis, by means of a harmonic-spring fictitious potential attached to the C4 0 and P atoms of the last two base pairs at one DSB end. After some tests, the spring constants were set at 100 and 75 kJ mol 1 nm 2 , respectively for the two DNA strand ends farther and closer to the nucleosome surface. To provide a reaction force keeping the system in place, all the atoms of the H3 opposite to the DSB were retained by soft harmonic restraints, with a spring constant of 250 kJ mol 1 nm 2 . A pulling speed of 5 mm/s was used for most SMD simulations. Forces and displacements were recorded at intervals of 5-10 time steps. Umbrella sampling was performed by extracting 100 configurations spaced by 50 ps during the first 5 ns of the force pulling simulation; force bias was progressively reduced from 100 down to 10 kJ mol 1 nm 2 , to extract the zero-bias limit of the free-energy profile; the weighted-histogram analysis was used to interpolate and connect the data from discrete configurations.

Molecular stress calculation

We use the so-called covariant central-force decomposition scheme (CCFD, [START_REF] Torres-Sanchez | Examining the mechanical equilibrium of microscopic stresses in molecular simulations[END_REF][START_REF] Torres-Sanchez | Geometric derivation of the microscopic stress: A covariant central force decomposition[END_REF] ) for the intra-and intermolecular forces, which ensures conservation of linear and angular momentum of the molecular systems under very general conditions. The method is implemented in a special-purpose patch to GROMACS 4.6, which reads (all or part of) a MD trajectory for the selected subset of atoms for which stress is to be computed, and performs the entire analysis. Since the GROMACS-LS patch [START_REF] Torres-Sanchez | Examining the mechanical equilibrium of microscopic stresses in molecular simulations[END_REF][START_REF] Torres-Sanchez | Geometric derivation of the microscopic stress: A covariant central force decomposition[END_REF] constrains the code to run in serial rather than in parallel, care must be taken to define properly the subset of interest in order to avoid prohibitive computing times. We prepared simple scripts to extract the principal components of the stress, compare stress fields from different simulations, and write the outputs in the portable Gaussian-cube format for visualization. Comparison between stress fields from different MD runs poses an extra care, since the structures need to share exactly the same box size and center, to avoid numerical artefacts from the cancellation between large positive and negative values. According to the CCFD scheme, stress fields are calculated by GROMACS-LS on a continuous grid superposed on the molecular structure; however, stress components and individual force contributions (pair, angle, dihedral, etc.) can also be projected back on the atom sites by defining a conventional (but non unique) atomic volume.

Snapshots of MD simulations of the M2-M4 models of DSB, after 1 µs of dynamics atT=350 K.

Figure 1 .

 1 Figure 1. Nucleosome and DNA defect geometry. (a) Schematic of the experimental nucleosome structure 1kx5, with DNA in grey; in orange, the 20 bp added to the experimental structure on each end. Histone pairs H3, H4, H2A, H2B color-coded as shown below. Positions 1 to 4 indicate insertion sites of double-strand breaks (DSB). (b) Schematic structure of a DSB; arrows indicate the cut on each backbone; the green-red central base pair is an A•• •T for all models. (c) The circular bent-tube DNA geometry in the nucleosome. At any point along the neutral axis (centerline) a local reference frame is defined by three unit vectors: normal n (blue), tangent t (red), and the binormal b (green) directed to the center of curvature. The white slice represents a tube averaging zone for stress calculations.

Figure 2 .

 2 Figure 2. Hydrogen bond length at a DSB base pair. Plot of the H-bond length for the central A•• •T bp of the model M3, for a MD simulation of 950 ns at T=350 K. Red, blue and black traces relative to the respective arrows in the scheme above, with the A•• •T bp in ball-stick representation. The red arrow indicates the central H-bond with adenine N1, the blue arrow the side H-bond with thymine O1, and the black arrow the (very weak) H-bond with thymine O2. Horizontal dashed lines in the plot indicate the reference, T=350K H-bond distances of 1.9 (lower) and 2.1 Å (upper).

Figure 3 .

 3 Figure 3. Eigenvectors of the covariance matrix for DNA at at DSB. Plot of the RMS fluctuation of the 4 principal eigenvectors, for the DNA fragments including the DSB, of the four models M1-M4. For each eigenvector (respective origins shifted along the ordinate axis) the black line gives the atom-by-atom contribution of the fragment surrounding the DSB, while the red line gives the same quantity for the same fragment intact (from model O). On each abscissa, atoms are grouped by the base number, in two contiguous blocks divided by the dashed line, representing the two parallel strands; the grey shaded regions indicate the position of the DSB (the red-green bases of Fig 1b) for each model.

Figure 4 .

 4 Figure 4. Eigenvectors of the covariance matrix for histones at at DSB. Plot of the RMS fluctuation of the 4 principal eigenvectors, for the histone residues closest to the DSB in each of the four models M1-M4. For each eigenvector (respective origins shifted along the ordinate axis) the black line gives the atom-by-atom contribution of the fragment surrounding the DSB, while the red line gives the same quantity for the same fragment intact (from model O). On each abscissa, atoms are grouped by residues, with a spacing (also indicated by colored bars) proportional to the size (i.e., number of atoms) of each residue.

Figure 5 .

 5 Figure 5. Different conformations of histone H2B tail at a DSB. Representation of the time-averaged H2B histone configurations: in the reference O model (yellow ribbons), in the M1 model at T=310 K (cyan) and in the M1 model at T=350 K (blue); the tail of the histone can be seen taking different orientations in the former two, w/r to the last simulation. Two portions of the upper and lower turn of DNA about the histone core are represented in silver ribbons; the spheres on the back represent the H3/H4 residues closer to the DSB; the DSB bases are colored in red-green according to Fig 1b.

  absence of the DSB, therefore it reflects the usual affinity of such residues for the DNA bases. The M1 model, instead, is definitely different, as it was the case for the DNA analysis in Fig3 above, and it will be treated later in this Section.

22 Figure 6 .

 226 Figure 6. Force-induced detachment of a DSB end from nucleosome. (a) Plot of the reaction coordinate z as a function of the MD simulation time, for a pulling simulation at constant-force. The inset shows the definition of the z distance of the broken DNA end (red arrow), and the direction of the applied force vector (blue arrow), parallel to the x axis and perpendicular to the superhelical axis. (b) Detachment of DNA portions from the histone core surface, from 0 to 5 ns of MD simulation at constant force and 310K. The three traces (black=left ordinates; blue/red=right ordinates) represent the distance to the surface of the three P atoms indicated in the scheme on the upper-left corner, in which a quarter of turn of the DNA comprised between the dyad and the DSB is sketched. The yellow sphere is the geometric center of the nucleosome.

  Fig. S 3, and the motion indicated by blue arrows in Fig. S 2),while promoting a strong interaction with the protein surfaces.

Figure 7 .

 7 Figure 7. Free energy for DSB detachment. Zero-force extrapolated free-energy profile as a function of MD time, during the force-pulling simulation at T=310 K. The red, blue and green regions refer to the sequence of energy barriers for the DNA-histone detachment events, as described in the text.

9 / 22 Figure 8 .

 9228 Figure 8. Stress profiles along DSB terminal ends. (a) Plot of the tension/compression force t(x) along part of the DNA model M1 configuration C185, at the beginning of the relaxation (black line) and after 40 ns (red line). Numbers along the upper and lower ordinates indicate the DNA bases in the two strands; the DSB is at the extreme left, bp 68-120; the grey rectangle is the main stress relaxation region, comprised between the DSB and the dyad (bp 94-94); color bars on the right refer to the two regions in Fig 1c. (b) Same as (a), for the M1 configuration C290. (c) Plot of the twist force w(x) along the DNA model M1 configuration C185. (d) Same as (c), for the M1 configuration C290.

3 4 5

 345 (a) Model M2, with DSB at the outer non-contact site, showing the central A74•• •T114 bp still well bonded. Grey spheres represent a portion of the H3 histone flanking the defect. (b) Model M3, with DSB at the dyad. Grey spheres represent a portion of the H3 tail. (c) Model M4, with DSB at the entry point of nucleosomal DNA. Grey spheres represent a portion of the H3 tail close to the break, which has folded into a double a-helix. . . . . . . . 18 2 Time-sequence of the evolution of the DNA-histone contacts for the DSB at the M1 position. The groups in grey VdW-spheres are the Leu82-Arg83-Phe84-Gln85 of H3, Lys77-Arg78-Lys79-Thr80 of H4 (left side); and Lys9-Gly10-Ser11-Lys12-Lys13, Lys24-Lys25-Arg26 of H2B (right side). DNA bp around the DSB are colored red-green according to the scheme of Fig.1b. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 Configurations associated with first eigenvector of the covariance matrix. On the left, representative configurations of the DNA fragments close to the DSB, in models M3 (above) and M1 (below). On the right, simultaneous plot of the configurations spanned by the principal motions associated with the first eigenvector, for each model. DNA fragments are aligned with their main axis vertical, the DSB being at the center. The superimposed frames are colored from blue to red, the ordering reflects a virtual motion spanning between the eigenvector extremes. A long stick spanning between the two colors identifies a large motion of the corresponding atom; a shorter stick identifies a local oscillation, of smaller amplitude. . . . . . . . . . . . . . . . 20 Configurations associated with eigenvector 1-4 of DSB model 4. Simultaneous plot of the configurations spanned by the principal motions associated with the eigenvectors 1-4, for the DNA fragment close to the DSB M4 (represented in the central panel). The superimposed frames are colored from blue to red, the ordering reflects a virtual motion spanning between the eigenvector extremes. Also in this case, a long stick spanning between the two colors identifies a large motion of the corresponding atom; a shorter stick identifies a local oscillation, of smaller amplitude. It can be readily appreciated that eigenvectors 1 and 3 correspond to a coordinated, twisting motion of the entire fragment, while eigenvectors 2 and 4 correspond to smaller and less cooperative deformations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 Stress relaxation in a detached DSB. Two configurations extracted from the force-pulling MD simulation of DSB at position M1, at t=1.8 ns (C180, left, red ribbons) and t=2.9 ns (C290, right, blue ribbons). The pull force was applied only at the C 0 -P atoms of the 2 last bp on the upper end of the DSB. The terminal portions of the pulled DSB end are highlighted as a thicker tube, for clarity. Row above: the two configurations at the start of the relaxation. Row below: the two configurations after 40 ns of MD equilibration/relaxation at T=310 K without any external forces applied. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

Table 1 .

 1 Excess entropy of DSB fragments. Upper limit of the excess entropy contribution T DS to the free energy at T=310 K, estimated from the Schlitter formula, Eq 2, for each molecular fragment in the different DSB models.

	DSB configuration DNA (kcal/mol) histones (kcal/mol)
	M1	89.7	29.8
	M2	25.6	6.9
	M3	38.1	2.7
	M4	47.9	5.3

The symbol 0 ⌦ 0 indicates the tensor product between the stress and the vector n, in practice the matrix product between the 3⇥3 matrix of the stress at each point x, and the 3-component vector locally perpendicular to the surface at x (see the local reference frame {n, t, b} in Fig1c).

The -symbol indicates the break site along each backbone, the ••• indicate the central interacting base pair.

Acknowledgments

Useful discussions about definitions of continuum stress with S. Giordano (IEMN Lille) are acknowledged. We gratefully thank partial funding from the SIRIC OncoLille under project "ModCel" for 2015 and 2016, and a three-year doctoral grant to F. L. co-funded by the

12/22

President of Lille-I University and the Government of the Region Nord-Pas de Calais (now Hauts-de-France). Computer resources provided by the CINES and IDRIS French Supercomputing Centres, under grants a0020707225 and x2016/077225, and by generous extensions thereof.

Author contributions

Declaration of interests

The authors declare no competing interests.

Supplemental data

Supplemental Data include five figures and can be found with this article online at http://www.cell.com/structure/supplemental/SXXXX.

List of Figures


Aurélien Fabas
email: aurelien.fabas@outlook.com

Rouaix-Vande Put

Sébastien Doublet

Didier Domergue

Mehdi Salem

Daniel Monceau

Metal dusting corrosion of austenitic alloys at low and high pressure with the effects of Cr, Al, Nb and Cu

Keywords: A, Superalloys C, High temperature corrosion C, Carburization C, Oxidation

Commercial and model alloys were exposed at 570 °C to metal dusting conditions at 1 and 21 bar in CO-H 2 -H 2 O and CO-H 2 -CO 2 -CH 4 -H 2 O atmospheres, respectively, with similar a c and P O2 . γ" precipitation beneath the surface of 625 superalloy was due to Nb diffusion toward the Cr-depleted zone formed during oxidation. 693 superalloy was more carburised than 690 superalloy. α-Al 2 O 3 formation was believed to induce cracks through the oxide scale, resulting in the alloy carburisation. Cu addition to a FeNiCr model alloy lead to the formation of an almost continuous Cu layer at the oxide/metal interface, greatly improving the metal dusting resistance.

Introduction

Metal dusting is a catastrophic type of corrosion affecting metallic materials in highly carburising atmospheres, characterised by a carbon activity (a c ) greater than unity (a c ≫ 1) [START_REF] Hochman | Metal deterioration in carbon monoxide and hydrocarbons at elevated temperatures[END_REF][START_REF] Hochman | Catastrophic deterioration of high-temperature alloys in carbonaceous atmospheres[END_REF]. This particular corrosion phenomenon is encountered in several industrial processes, such as steam methane reforming [START_REF] Holland | Metal dusting failures in methane reforming plant[END_REF], ammonia [START_REF] Eberle | An industrial experience of attack on metals by synthesis gas from methane-oxygen combustion[END_REF] and methanol production [START_REF] Szakálos | Metal dusting on an alumina forming Nibase alloy[END_REF] and coal gasification [START_REF] Perkins | Metal-dusting corrosion in coal-gasification environments[END_REF] plants. The degradation processes and the degradation rates are highly dependent on the alloy composition [7 9] and microstructure [START_REF] Grabke | Resistance of 9-20%Cr-steels against metal dusting[END_REF][START_REF] Grabke | Effects of grain size, cold working, and surface finish on the metal-dusting resistance of steels[END_REF], and also on the environment [6,12 14]. Many studies were performed to understand corrosion mechanisms, mainly at atmospheric pressure [5,7,15 19]. Some studies were also performed at high pressures [START_REF] Szakálos | Metal dusting on an alumina forming Nibase alloy[END_REF][START_REF] Albertsen | A model for hightemperature pitting corrosion in nickel-based alloys involving internal precipitation of carbides, oxides, and graphite[END_REF]. Only a few publications dealt with a comparison between tests performed in both pressure domains [5, 21 24]. The study of Rouaix Vande Put et al. [START_REF] Rouaix-Vande Put | Performance of chromiaand alumina-forming Fe-and Ni-base alloys exposed to metal dusting environments: the effect of water vapor and temperature[END_REF] was carried out at constant a c and increasing total pressure with different oxygen partial pressures (P O2 ). Nishiyama et al. studied the influence of carbon and oxygen activities at low [START_REF] Nishiyama | Effect of syngas compositionon metal dusting of alloy 800H in simulated reforming gas atmospheres[END_REF] and high [START_REF] Nishiyama | Metal dusting behaviour of alloy 800H in laboratory carbonaceous environments under high pressure[END_REF] pressure. However, to assess the effect of the total pressure, new experiments at constant a c and P O2 are necessary since these two parameters are recognised as important in the metal dusting literature. The resistance to metal dusting of several alloys in atmospheric pressure and high pressure environments at close carbon and close oxygen activities is reported in the following. A special attention is also given to the effects of Cr, Al, Cu and Nb. A discussion about the relevance of parameters other than the carbon activity, based on the current experiments, is already available [START_REF] Put | Relevance of other parameters than carbon activity in defining the severity of a metal dusting environment[END_REF].

Materials and experiments

Materials

Two metal dusting experiments were performed at 570 °C. The first one was carried out at atmospheric pressure (AP) in a 47.25%CO 47.25%H 2 5.5%H 2 O gas mixture with a gas velocity of 1.8 × 10 -5 m/ s. The second test was operated at high pressure (HP) 21 bar in a 49.1%H 2 12.8%CO 3.1%CO 2 1.6%CH 4 33.4%H 2 O gas mixture with a gas velocity of 2.8 × 10 -3 m/s.

According to the syngas reaction, carbon activities of 33 and 32 were calculated for the AP and HP experiments, respectively. Oxygen partial pressures of 2.1 × 10 -27 bar and 7.1 × 10 -26 bar were esti mated for the AP and HP corrosion tests respectively, considering the water decomposition reaction. A more precise description of the AP and HP rigs and experiments is reported in references [START_REF] Fabas | Mechanism of metal dusting corrosion by pitting of a chromia-forming alloy at atmospheric pressure and low gas velocity[END_REF][START_REF] Fabas | Modelling of the kinetics of pitting corrosion by metal dusting[END_REF] respectively.

Five commercial alloys (800HT, HR120, 625, 690 and 693) were tested in both rigs. Samples assigned to the AP and HP tests were made from different batches, except for the HR120 alloy. Two model alloys fabricated by Spark Plasma Sintering (SPS) were also tested under HP. A ternary FeNiCr alloy was made based on the composition of HR120 alloy, and a quaternary FeNiCrCu alloy was prepared substituting Ni by Cu. Model alloys were sintered from powders of the pure elements. Because of the initial oxidation of the Fe raw powder, the resulting microstructure contained a few volume percent of internal chromia. The chemical compositions of the alloys are reported in Table 1. All samples were ground to a P600 grit. The edges were chamfered to limit any edge effect. The samples for the AP experiment were 20 × 10 mm 2 rectangles and 2 mm thick, except HR120 and 693 samples which were 1 mm thick. Samples tested at HP were disks with a diameter ranging between 10 and 14 mm and a thickness comprised between 1 and 2.6 mm. The central hole was about 1.5 mm wide. Samples made of the same material were placed upstream and downstream to ensure that the atmosphere conditions were constant along the reaction chamber in both corrosion rigs.

Sample characterisation

For both experiments, the test was interrupted about every 500 h. At every interruption samples were removed and cleaned in ethanol. Samples were systematically observed using an optical microscope in order to follow the evolution of the corrosion damage on the surface. The samples were also weighed to determine the mass change during the test. After final removal the specimens were characterised by X ray diffraction (XRD) on a Seifert 3000TT apparatus with a copper anti cathode (λ = 1.54056 Å), with a low angle of incidence of 2°, Raman spectroscopy analyses were performed on a 64000 Jobin Yvon spectro meter equipped with a confocal microscope (λ = 514 nm). Scanning Electron Microscopy (SEM) using Secondary Electron (SE) and Backscattered Electron (BSE) modes were carried out on a SEM/FIB (FEI HELIOS 600i), coupled with an Energy dispersive X ray (EDX) spectroscopy analyser (from Aztec Advanced) equipped with a SSD detector. Analyses were performed using a 5 keV accelerating voltage. The two above techniques were used on surface and cross sections of the samples. Cross sections were prepared cutting the sample after having protected the surface with a TEM resin and grounding and polishing to obtain a mirror polished surface. The experimental procedure for AP and HP tests is fully detailed in references [START_REF] Fabas | Mechanism of metal dusting corrosion by pitting of a chromia-forming alloy at atmospheric pressure and low gas velocity[END_REF][START_REF] Fabas | Modelling of the kinetics of pitting corrosion by metal dusting[END_REF] respectively.

Results

Overview

The net mass changes of the alloys during the AP and HP experiments are shown on Fig. 1. The Fe rich commercial alloys (800HT, HR120) presented large mass losses in both corrosion rigs while the Ni based superalloys (625, 690, 693) showed no significant mass change. The SPS FeNiCr samples were severely attacked, unlike the Cu enriched ones. The present paper focuses on the resistance of the three commercial superalloys and of the FeNiCr + Cu model alloy fabricated by SPS. The oxidation and the metal dusting degradation of 800HT and HR120 alloys will be studied in details in forthcoming articles see also references [START_REF] Put | Relevance of other parameters than carbon activity in defining the severity of a metal dusting environment[END_REF][START_REF] Fabas | Mechanism of metal dusting corrosion by pitting of a chromia-forming alloy at atmospheric pressure and low gas velocity[END_REF]. All samples presented at first a mass gain, not visible due to the scale of the figure. This mass gain can be linked to the oxidation and carburation of the samples, i.e. the incubation period before the large mass loss.

Macrophotographs of superalloy specimens exposed 6000 h in AP conditions and 7966 h under HP conditions are presented on Fig. 2. For these materials no sign of degradation was visible. However, the oxide scale formed under HP seemed homogeneous, while variations in colours were visible on the AP specimens. Colour variations were localised where the samples rested on the sample holder. Thus, this is expected to be due to the modification of the gas flow in these areas.

As each material showed a similar behaviour in both conditions and for all exposure durations, only the superalloys exposed 7966 h in the HP environment will be presented.

Superalloys

Surface analyses

SEM observations showed that the oxide scale formed on the 625 (Fig. 3a) and 690 superalloys (Fig. 3b) was continuous, although having an inhomogeneous thickness. The oxide scale formed on the 693 superalloy exhibited Al and C rich areas (Fig. 3c) and locally a needle like morphology (Fig. 3d). Chromia and Fe Cr spinel were identified as the main constituents of the oxide scale for every material, based on XRD and Raman spectroscopy analyses. The fluorescence of Cr 3+ in α Al 2 O 3 i.e. the presence of α Al 2 O 3 was hardly detected during Raman analyses on the 693 sample exposed 1987 h under HP. Its detection became easier after 7966 h in the same conditions, and even more after 14,177 h of exposure.

Cross sections

The cross section of the 625 sample exposed 7966 h at HP showed needle like precipitates (Fig. 4a andb). After 14,177 h of exposure, these precipitates can form locally a continuous layer at the alloy oxide interface. EDS maps obtained after FIB milling revealed that these needles contained Nb. They were finally identified as being made of the γ" Ni 3 Nb phase, previously detected by XRD. These needles precipi tated at or close to the Cr depleted zone that stemmed from the oxide scale formation (Fig. 4b and Fig. 5). Other precipitates containing Cr or Si were also visible around the γ" needles. According to the work of Garcia Fresnillo et al. [START_REF] Garcia-Fresnillo | Oxidation behaviour and microstructural stability of alloy 625 during long-term exposure in steam[END_REF] on the oxidation of 625 superalloy at 700 °C and 800 °C, these precipitates could be Cr carbides and silica. Al rich oxide precipitates were also visible below the oxide scale. These ones are expected to be internal alumina, as observed earlier for the same alloy during oxidation tests at 760 °C [START_REF] Deodeshmukh | Long-Term Oxidation and Chromia Vaporization Behavior of High-Temperature Foil Alloys in Water vapor Containig Environment[END_REF], 900 °C and 1000 °C [START_REF] Chyrkin | Sub-scale depletion and enrichment processes during high temperature oxidation of the nickel base alloy 625 in the temperature range 900-1000[END_REF].

The cross sections of 690 and 693 superalloys exposed 7966 h under metal dusting conditions at HP showed a carburisation for both alloys (Fig. 4c andd). However, this carburisation is more severe for 693 alloy. In addition, the carburised depth of 693 increased with the exposure time, while it remained constant for 690. EDS maps performed on a 693 cross section obtained by FIB milling showed that the oxide scale was mainly composed of Al and Cr (Fig. 6). No Cr depletion was visible beneath the oxide scale. The Ni and Al enrichments close to the Cr carbides were due to the local Cr consumption. It has to be noticed that the O detection at the location of Cr carbides is an artefact due to the overlap between the K band of O and the L band of Cr. Indeed, if the observed precipitates were Cr rich oxides, alumina would have formed, which was not the case.

SPS alloys

Surface analyses

As explained previously, the SPS FeNiCr alloy was deeply attacked, unlike SPS FeNiCr + Cu samples (Fig. 1). The degradation of the ternary specimens took place by pitting, while the Cu enriched alloy showed no macroscopic evidence of degradation (Fig. 7). The attack of the ternary material occurred preferentially at the edges of the sample. Its location could be explained by damage of the oxide scale during temperature change as mechanical stress is more important at the edges or by chemical failure due to the larger Cr depletion at the sample edges. It could also be explained by the manufacturing process. Indeed, samples were taken in the transverse direction of a sintered bar fabricated in a graphite mold. The carbon could have diffused close the metal powder mold interface, leading to a pre carburisation of the external parts of the sample. To understand the resistance to metal dusting corrosion of the Cu rich material, its cross section was char acterised and compared to its initial microstructure.

Cross sections

Before exposure, the Cu rich alloy showed the precipitation of a Cu rich phase at the grain boundaries (Fig. 8). In addition, much finer Cu rich precipitates were visible inside the grains and grain boundaries were surrounded by a precipitate free zone. These precipitates are believed to be almost pure FCC ε Cu. It is proposed that the intergranular precipitation took place during sintering, while the intragranular precipitates formed during cooling.

The cross section characterisation of the sample exposed 7966 h under the HP conditions revealed the formation of a Cr rich oxide scale (Fig. 9). This oxide scale was identified as chromia by both XRD and Raman spectroscopy. Its formation led to a Cr depletion. EDS analyses measured between 10 to 15 at.% of Cr within the two micrometres beneath the alloy oxide interface. Beneath the oxide scale, i.e. in the Cr depleted alloy, the number of Cu rich precipitates was lower than in the bulk, and the precipitates were larger, forming an almost continuous layer. The fine secondary precipitates observed before exposure are not present in the Cr depleted zone of the corroded alloy.

Discussion

Ni based superalloys

The precipitation of a Nb rich phase has already been reported in previous studies on 625 alloy under oxidising conditions above 700 °C [ 28 31]. Between 900 °C and 1000 °C, such precipitation was accom panied by the appearance of a precipitate free zone [START_REF] Chyrkin | Sub-scale depletion and enrichment processes during high temperature oxidation of the nickel base alloy 625 in the temperature range 900-1000[END_REF]. The authors explained the formation of this particular microstructure as the result of the Nb activity decrease due to the Cr depletion that stemmed from oxidation. The lower Nb activity in the Cr depleted zone induced the outward Nb diffusion toward the alloy oxide interface where it precipitated, forming a δ Ni 3 Nb layer. Such uphill diffusion induced also a depletion in Nb below the Cr depleted zone, leading to the dissolution of the δ Ni 3 Nb precipitates, hence the formation of a precipitate free zone. The experiments performed between 700 °C and 800 °C did not present such precipitate free zone [START_REF] Garcia-Fresnillo | Oxidation behaviour and microstructural stability of alloy 625 during long-term exposure in steam[END_REF][START_REF] Deodeshmukh | Long-Term Oxidation and Chromia Vaporization Behavior of High-Temperature Foil Alloys in Water vapor Containig Environment[END_REF][START_REF] Gorman | Microstructural analysis of IN617 and IN625 oxidised in the presence of steam for use in ultra-supercritical power plant[END_REF]. However, the mechanism proposed by Chyrkin et al. [START_REF] Chyrkin | Sub-scale depletion and enrichment processes during high temperature oxidation of the nickel base alloy 625 in the temperature range 900-1000[END_REF] can be adapted to explain the current observations. The detection of γ" Ni 3 Nb in the present study, instead of the δ phase, can be justified by the experi mental conditions. The TTT diagram of 625 alloy [START_REF] Schnabel | The precipitation and recrystallization behavior of the nickel base alloy INCONEL alloy 625[END_REF] shows that the γ" Ni 3 Nb phase is more likely to form at 570 °C than the δ phase, the latter being expected for higher temperatures [28 31]. The presence of γ" Ni 3 Nb needles instead of a continuous layer can also be explained by the experimental conditions. Indeed, the characteristic diffusion length of Nb in pure Ni is at least one order of magnitude below the ones obtained for the conditions leading to a continuous layer formation (Table 2). The needle shape can be explained by the lattice misfit between the γ matrix (austenitic) and the γ" precipitates [START_REF] Paulonis | Precipitation in nickel-basealloy 718[END_REF].

The reason for the great resistance to metal dusting corrosion of 625 alloy is then, as expected for high alloyed Ni based materials, its low carbon permeability [START_REF] Zhang | Coking and dusting of Fe-Ni alloys in CO-H2-H2O gas mixtures[END_REF] and its high Cr content which allows the formation and healing of a protective oxide scale [START_REF] Klöwer | Metal dusting of nickel-base alloys[END_REF]. In addition, its high Mo and Nb contents contribute to enhance its corrosion resistance, trapping carbon [START_REF] Strauß | Role of alloying elements in steels on metal dusting[END_REF] and letting more Cr available to form the protective oxide scale. The effect of the Nb is believed to be enhanced by its uphill diffusion.

The severe carburisation of 693 alloy could be linked to its higher Al content compared to the 690 alloy. The increasing Cr 3+ fluorescence signal of α Al 2 O 3 with exposure duration, coupled to the needle like morphology observed for the oxide on the 693 specimen exposed 7966 h under HP, could be explained by the formation of transient alumina that turned into α Al 2 O 3 in time. Such allotropic transforma tion is accompanied by a large molar volume contraction [START_REF] Lipkin | Stress development in alumina scales formes upon oxidation of (111) NiAl single crystals[END_REF] which may induce cracks through the oxide scale, as often observed during the oxidation of nickel aluminides at higher temperatures, e.g. by Tolpygo et al. [START_REF] Tolpygo | Microstructural study of the theta-alpha transformation in alumina scales formed on nickel-aluminides[END_REF]. Such cracks could constitute a direct access of the corrosive atmosphere to the metallic alloy, leading to the observed severe carburisation. Nevertheless, no transient alumina could be detected by Raman spectroscopy and XRD, but this may be due to the fact that the scale was very thin. EDS mapping on Fig. 6 did not discriminate the presence of Al in the oxide scale, the latter being too thin compared to the interaction volume of the EDS analysis. As there is no clear evidence of transient alumina formation, it is also possible that α alumina formed directly from a Cr and Al rich oxide. Indeed, XRD coupled to EDS and Raman analysis suggested that the corundum phase was (Cr,Al) 2 O 3 . A second proposal can be made to explain the high level of carburisation of 693 alloy. It is known that Al enrichment over the miscibility limit of the (Cr,Al) 2 O 3 leads to a spinodal decomposition into chromia and α alumina [START_REF] Schultz | Separation of phases by spinodal decomposition in the systems Al2O3-Cr2O3 and Al2O3-Cr2O3-Fe2O3[END_REF] and is accompanied by a molar volume contraction. Then this could also be responsible for the crack formation within the oxide scale, explaining the severe carburisation of the alloy. This phenomenon would not occur on 690 because of its lower Al content. The slight carburisation of this material is supposed to stem from the early stages of exposure, when the initial competition between carburisation and oxidation took place. The resistance to metal dusting corrosion of both 690 and 693 alloys, despite the detrimental effect of Al for the latter, is explained by their high Cr content, allowing the formation and the healing of their protective chromia scale.

SPS fabricated alloys

The formation of Cu rich precipitates at the alloy oxide interface of the Cu rich specimen exposed 7966 h to the HP environment, forming an almost continuous FCC Cu rich layer, could be explained as follows. As Cr content diminished as it was consumed by the oxide growth, Cu concentration increased and the alloy composition entered in the Fig. 6. SEM (SE) micrograph of the FIB cross-section and its associated elementary maps of 693 superalloy after 7966 h of exposure at HP. miscibility gap. This lead to the precipitation of a Cu rich phase [START_REF] Zhang | Effect of copper on metal dusting of austenitic stainless steels[END_REF]. This precipitation was also enhanced by the fact that interfaces are preferential sites for precipitation. As the oxidation of the sample continues, the Cr depletion increases, favouring the growth of the Cu precipitates and forming the observed microstructure at the alloy oxide interface. During cooling, it can be supposed that the growth of the precipitates was favoured to the detriment of the precipitation of secondary precipitates, explaining their absence in the Cr depleted zone. However, more experiments and characterisations are necessary to firmly confirm the proposed mechanism for FCC Cu precipitation. Cu addition reduces carbon formation catalysis, as Cu possesses filled d orbitals [START_REF] Andreoni | Binding and dissociation of CO on transition-metal surfaces[END_REF][START_REF] Nishiyama | Improving metal dusting resistance of transition-metals and Ni-Cu alloys[END_REF]. It diminishes also graphite nucleation [START_REF] Zhang | Alloying with copper to reduce metal dusting of nickel[END_REF] and carbon solubility [START_REF] Zhang | Alloying with copper to reduce metal dusting of nickel[END_REF], improving metal dusting resistance. The resistance to metal dusting corrosion of the Cu rich samples can then be explained by their high Cu content which reduced the carbon transfer into the material, favouring the formation and healing of the protective oxide scale. The formation of the Cu precipitates at the alloy oxide interface, where they could form an almost continuous layer for longer exposure times, is expected to enhance its effect. The influence of this precipita tion on the oxide adherence and thermodynamical stability, but also on the uphill Cr diffusion, have to be taken in consideration. Cu precipita tion could indeed be beneficial for metal dusting resistance but detrimental for the oxidation resistance. It has to be noticed that it was reported in the literature that Cu addition to commercial Fe Ni Cr alloys above the Cu solubility limit had a detrimental effect on metal dusting resistance [START_REF] Zhang | Effect of copper on metal dusting of austenitic stainless steels[END_REF]. The authors explained this behaviour by the precipitation of the Cu rich phase, whose interface with the matrix was considered as a preferential graphite germination site. This detrimental effect was not observed in the current study, despite a pronounced Cu rich phase precipitation. This difference could be explained by the use of model alloys, avoiding the detrimental effect of minor elements contained in the commercial alloys used by Zhang and Young [START_REF] Zhang | Effect of copper on metal dusting of austenitic stainless steels[END_REF]. The difference could also be explained by a different surface preparation, as the samples of Zhang and Young were polished to a 3 μm finish while the specimens used in the present study were ground. Indeed, it is well known that grinding enhances the resistance to metal dusting of alloys which form a protective oxide scale [START_REF] Grabke | Effects of grain size, cold working, and surface finish on the metal-dusting resistance of steels[END_REF].

Conclusions

The resistance to metal dusting of three commercial Ni based superalloys and of a Cu rich FeNiCr model alloy was investigated under high pressure for durations up to 14,000 h. 625 alloy showed great resistance to metal dusting due to its Cr rich composition, allowing the formation of a protective oxide scale. Cr depletion resulting from its selective oxidation lead to a Nb uphill diffusion and γ" Ni 3 Nb needles formed below the surface, beneficial for the resistance to metal dusting.

Despite a more or less severe carburisation regarding the material, 690 and 693 alloys also showed no evidence of metal dusting degrada tion because of their high Cr content. The slight carburisation observed on 690 alloy is supposed to come from the early stages of exposure, during the initial competition between oxidation and carburisation. The severe carburisation of 693 alloy stemmed from a direct access of the atmosphere to the metallic material through cracks in the oxide scale. These cracks could be explained by the molar volume contraction induced either by a transient alumina formation that turned into α Al 2 O 3 or by a Cr and Al rich oxide decomposition into chromia and α Al 2 O 3 . Al is thought to be detrimental in the studied conditions.

The FeNiCrCu quaternary alloy resistance to metal dusting is explained by its Cu content, as this element decreases carbon transfer and graphite nucleation rates. The beneficial influence of Cu might be enhanced by its precipitation at the alloy oxide interface. The influence of the precipitation on the oxidation of the material should be investigated as it might influence the thermodynamics and the kinetics of the oxidation due to the localisation of the precipitates.

Fig. 1 .

 1 Fig. 1. Net mass change of the commercial and model alloys during metal dusting tests at (a) AP and (b) HP. (c) Higher magnification of the HP mass changes up to 7000 h.

Fig. 2 .

 2 Fig. 2. Optical macrophotographs of the superalloys 625 (a and c), 690 (b and d) and 693 (c and f) exposed 6000 h at AP (a-c) and 7966 h at HP (d-f).

Fig. 3 .

 3 Fig. 3. SEM (BSE) images of the surface of the superalloys (a) 625, (b) 690 and (c and d) 693 after 7966 h of exposure at HP.

Fig. 4 .Fig. 5 .

 45 Fig. 4. SEM (BSE) images of the FIB cross section of the superalloys (a) 625 before exposure and (b) 625, (c) 690 and (d) 693 after 7966 h of exposure at HP.

Fig. 7 .

 7 Fig. 7. Optical macrophotographs of the alloys (a) SPS FeNiCr after 6071 h of exposure at HP and (b) SPS FeNiCr + Cu after 7966 h of exposure at HP.

Fig. 8 .

 8 Fig. 8. SEM (SE) micrograph of the FIB cross-section and its associated elementary maps of the SPS FeNiCr + Cu alloy before exposure.

Fig. 9 .

 9 Fig. 9. SEM (SE) micrograph of the FIB cross-section and its associated elementary maps of the SPS FeNiCr + Cu alloy after 7966 h of exposure at HP.

Table 1

 1 Alloy composition (at%).

	Alloy		Ni	Fe	Cr	Al	Co	Mo	Nb	Mn	Si	Ti	Cu	C	V	Technique
	800HT	AP	31.3	43.3	20.6	0.9	0.07	0.16	0.01	1.0	1.2	0.7	0.2	0.25	0.08	OES a
		HP	28.4	45.3	22.4	1.3	-	-	-	0.8	1.0	0.6	0.2	-	-	EDX
	HR120		38.2	33.3	25.6	0.2	0.02	0.13	0.4	0.8	0.8	0.02	0.07	0.33	0.08	OES
	625	AP	65.7	4.4	22.3	0.5	0.04	4.2	1.7	0.1	0.4	0.4	0.01	0.17	0.03	OES
		HP	64.1	1.0	25.3	0.4	0.2	5.5	2.3	0.5	0.8	0.3	-	-	-	Certificate
	690	AP	60.1	9.4	28.4	0.6	0.02	0.02	-	0.3	0.6	0.4	0.01	0.09	0.04	OES
		HP	59.2	8.5	30.7	-	-	-	-	0.1	0.6	-	-	0.9	-	Certificate
	693	AP	58.0	3.0	30.8	7	-	-	-	0.3	0.1	0.5	-	-	-	EDX
		HP	59.0	3.2	30.4	6.2	-	-	0.4	0.2	0.1	0.5	0.03	0.08	-	Certificate
	FeNiCr	HP	35.3	40.8	23.9	-	-	-	-	-	-	-	-	-	-	EDX
	FeNiCr + Cu	HP	29.2	36.0	26.5	-	-	-	-	-	-	-	8.3	-	-	EDX
	a Optical Emissions Spectrometer.

Acknowledgements

This work has been supported by the French national research agency through the project ANR SCAPAC 11 RNMP 0016 in partnership between Air Liquide, Veolia, Sedis, University of Nancy, CIRIMAT laboratory. The authors acknowledge Sabine Leroux of the ICA laboratory for daily monitoring of the atmospheric pressure experi ment.