Phuong H Nguyen

Philippe Derreumaux

A S-shaped Aβ42 Cross-β Hexamer Embedded into a Lipid Bilayer Reveals Membrane Disruption and Permeability

The interactions of amyloid oligomers with membranes are known to contribute to cellular toxicity. Numerous in vitro experimental studies reported on the insertion of oligomers of different sizes that can induce cell membrane disruption, extract lipids and form ion-permeable transmembrane pores. The current repertoire of amyloid-beta (Aβ) membrane-inserted folds that was subject to high-resolution structure NMR spectroscopy and computer simulations is devoid of any cross-β fibrillar structure. In this study, we explored the dynamics of a S-shaped Aβ42 cross-β hexamer model inserted into a lipid bilayer membrane by two atomistic molecular dynamics simulations. The initial model is characterized by the hydrophobic residues at the central hydrophobic core (residues 17-21, CHC) and the C-terminus (residues 30-42) embedded into the membrane. We observed major structural secondary, tertiary and quaternary rearrangements leading to two distinct species, hexamer and two trimers, accompanied by membrane disruption and water permeation. The simulations show that some configurations, but not the majority, have the CHC and C-terminus hydrophobic residues exposed to solvent. Overall, our computational results offer new perspectives to understand the relationship between Aβ42 assemblies and membrane permeability.

1.INTRODUCTION

Neurogenerative diseases are associated with the misfolding and aggregation of intrinsically disordered proteins which lack stable secondary structures in the bulk solution and at the surface of membranes. [START_REF] Dobson | Protein Misfolding, Evolution and Disease[END_REF][START_REF] Nasica-Labouze | Amyloid β Protein and Alzheimer's Disease: When Computer Simulations Complement Experimental Studies[END_REF][START_REF] Arai | Conformational Propensities of Intrinsically Disordered Proteins Influence the Mechanism of Binding and Folding[END_REF][START_REF] Nguyen | Structures of the Intrinsically Disordered Aβ, tau and αsynuclein Proteins in Aqueous Solution from Computer Simulations[END_REF] Numerous experimental studies explored the interactions of amyloid oligomers with membrane mimetic systems including liposomes, vesicles, micelles, monolayer and bilayer models and nanodiscs. They led to three membrane-induced toxicity models; the membrane destabilization model, the detergent (lipid extraction) model and the pore model, [START_REF] Quist | Amyloid Ion Channels: A Common Structural Link for Protein-misfolding Disease[END_REF][START_REF] Liu | Molecular Simulation Aspects of Amyloid Peptides at Membrane Interface[END_REF][START_REF] Nguyen | Amyloid Oligomers: A Joint Experimental/Computational Perspective on Alzheimer's Disease, Parkinson's Disease, Type II Diabetes, and Amyotrophic Lateral Sclerosis[END_REF][START_REF] Gonzalez-Garcia | Membrane Interactions and Toxicity by Misfolded Protein Oligomers[END_REF][START_REF] Berthelot | What Does Make an Amyloid Toxic: Morphology, Structure or Interaction with Membrane?[END_REF][START_REF] Ait-Bouziad | Discovery and Characterization of Stable and Toxic Tau/phospholipid Oligomeric Complexes[END_REF][START_REF] Lasagna-Reeves | The Formation of Tau Pore-like Structures is Prevalent and Cell Specific: Possible Implications for the Disease Phenotypes[END_REF][START_REF] Di Scala | Common Molecular Mechanism of Amyloid Pore Formation by Alzheimer's β-amyloid Peptide and α-synuclein[END_REF] which were demonstrated for amyloid-beta (Aβ) and tau oligomers [START_REF] Quist | Amyloid Ion Channels: A Common Structural Link for Protein-misfolding Disease[END_REF][START_REF] Lasagna-Reeves | The Formation of Tau Pore-like Structures is Prevalent and Cell Specific: Possible Implications for the Disease Phenotypes[END_REF][START_REF] Di Scala | Common Molecular Mechanism of Amyloid Pore Formation by Alzheimer's β-amyloid Peptide and α-synuclein[END_REF] contributing to Alzheimer's disease, and α-synuclein oligomers linked to Parkinson's disease. [START_REF] Di Scala | Common Molecular Mechanism of Amyloid Pore Formation by Alzheimer's β-amyloid Peptide and α-synuclein[END_REF] The "channel hypothesis" of Alzheimer's disease dates back from the electrophysiological measurements performed in 1993 by Arispe et al. [START_REF] Arispe | Giant Multilevel Cation channels Formed by Alzheimer Disease Amyloid Beta-protein [A beta P-(1-40)] in Bilayer Membranes[END_REF] Initially, they proposed the β-hairpin followed by a helix-turn-helix model of Aβ42 peptide to design three pore models, with the pore being formed either by the hairpins, middle helices or C-terminal helices. [START_REF] Durell | Theoretical Models of the Ion Channel Structure of Amyloid Beta-protein[END_REF] Recently, they designed βbarrels of six hexamers and pores formed by several hexamers spanning the bilayer. [START_REF] Shafrir | Models of Membrane-bound Alzheimer's Abeta Peptide Assemblies[END_REF] On the basis of SDS-PAGE, circular dichroism and electrical recording, specific Aβ42 pore forming β-barrel tetramers with an inner diameter of 0.7 nm in dodecylphosphocholine (DPC) micelle conditions were reported. [START_REF] Serra-Batiste | Aβ42 Assembles into Specific β-Barrel Pore-forming Oligomers in Membrane-mimicking Environments[END_REF] A tetrameric β-barrel, consisting of two distinct β-hairpin motifs and an asymmetric arrangement of eight antiparallel β-strands with an inner of diameter of 0.7 nm, was found by atomistic replica exchange dynamics simulations (REMD) to be destabilized for Aβ40 compared to its Aβ42 counterpart in a membrane model, explaining the absence of ionic currents of Aβ40 in planar lipid bilayers. [START_REF] Nguyen | Tetrameric Aβ40 and Aβ42 β-Barrel Structures by Extensive Atomistic Simulations. I. In a Bilayer Mimicking a Neuronal Membrane[END_REF] Additional REMD simulations showed that the lower AD propensity of aged rodents, characterized by mutations R5G, Y10F and H13R compared to human sequence, might be correlated to its tetrameric β-barrel stability in a model membrane. [START_REF] Ngo | Impact of the Rat R5G, Y10F, and H13R Mutations on Tetrameric Aβ42 β-Barrel in a Lipid Bilayer Membrane Model[END_REF] Another REMD simulation in a dipalmitoylphosphatidylcholine (DPPC) membrane showed that the wild-type Aβ42 sequence and its protective-disease A2T variant have similar propensity to form β-barrels. [START_REF] Ngo | Impact of A2T and D23N Mutations on Tetrameric Aβ42 Barrel within a Dipalmitoylphosphatidylcholine Lipid Bilayer Membrane by Replica Exchange Molecular Dynamics[END_REF] Other experimental studies showed that Aβ β-hairpin mimics form a triangular trimer unit that self-assembles into hexamers and dodecamers with annular pores, which are toxic toward neuronally derived cells. [START_REF] Kreutzer | Elucidating the Structures of Amyloid Oligomers with Macrocyclic β-Hairpin Peptides: Insights into Alzheimer's Disease and Other Amyloid Diseases[END_REF] Ion mobility-mass spectrometry revealed formation of hexameric βbarrels in a membrane-mimicking environment. [START_REF] Österlund | Native Ion Mobility-Mass Spectrometry Reveals the Formation of β-Barrel Shaped Amyloid-β Hexamers in a Membrane-Mimicking Environment[END_REF] Mass spectrometry and nuclear magnetic resonance (NMR) spectroscopy in DPC micelles at pH 9.0 provided the first high-resolution structure of Aβ42 tetramer and octamer and revealed edge conductivity pores as a mechanism for membrane damage. The tetrameric topology comprises a six stranded β-sheet core with the two external Aβ peptides forming β-hairpins. [START_REF] Ciudad | Aβ(1-42) Tetramer and Octamer Structures Reveal Edge Conductivity Pores as a Mechanism for Membrane Damage[END_REF] Overall, a wide range of concentric β-barrels is consistent with image averaged electron micrographs. [START_REF] Durell | The Amyloid Concentric β-Barrel Hypothesis: Models of Amyloid beta 42 Oligomers and Annular Protofibrils[END_REF] Inspired by many biological ion channels using α-helix bundles, simulations showed that trimer and tetramer of Aβ17-42 transmembrane α-helices are possible candidates to conduct Na + ions. [START_REF] Ngo | Probable Transmembrane Amyloid α-Helix Bundles Capable of Conducting Ca 2+ Ions[END_REF] Apart from the organization in α-helix bundles and β-barrels made of β-hairpins to build pores, Nussinov et al. proposed, on the basis of atomic force microscopy (AFM) imaging and modelling, that Aβ17-42 pentamer with U-shaped fibrillar state is able to insert into the lipid bilayer. [START_REF] Jang | Mechanisms for the Insertion of Toxic, Fibril-like β-Amyloid Oligomers into the Membrane[END_REF] The same strand-loop-strand fibril state was used to construct annular Aβ oligomers with a number of peptides varying between 12 to 36 and pore inner diameters varying between 2 and 4 nm. Molecular dynamics (MD) simulations in different lipid bilayers showed that pores made of 18 and 24 peptides with the C-terminal β-strands in contact with the lipid hydrophobic tails retained their morphologies. [START_REF] Jang | Disordered Amyloidogenic Peptides May Insert into the Membrane and Assemble into Common Cyclic Structural Motifs[END_REF] Of note, there is a structural convergence among diverse, toxic β-sheet ion channels, [START_REF] Jang | Structural Convergence among Diverse, Toxic Beta-sheet Ion Channels[END_REF] and the U-shaped conformation was also predicted by a basinhoping global optimization and an oligomer generation procedure for 6-and 8-mers of Aβ42 peptide in an implicit membrane environment. [START_REF] Strodel | Transmembrane Structures for Alzheimer's Aβ(1-42) Oligomers[END_REF] Overall, the repertoire of Aβ folds in membranes is diverse, and this is supported by REMD simulations which showed that Aβ11-40 trimer resulting from assembly of U-shaped fibril states or β-hairpins represents minimal seeds for the formation of amyloid fibrils, multiple β-rich aggregates and β-barrel pores in a membrane. [START_REF] Ngo | Stability of Aβ11-40 Trimers with Parallel and Antiparallel β-Sheet Organizations in a Membrane-Mimicking Environment by Replica Exchange Molecular Dynamics Simulation[END_REF] Membrane disruption and lipid extraction by amyloid oligomers were assessed by many experiments. They were demonstrated by cryo-electron tomography from the interactions of the amyloid fibril ends of β2-microglobulin and membrane surface. [START_REF] Milanesi | Direct Three-dimensional Visualization of Membrane Disruption by Amyloid Fibrils[END_REF] Nanoscale imaging experiments showed how Aβ42 fibrils displace the upper leaflet of the membrane as fibers laterally embed into the bilayer. [START_REF] Bode | Amyloid-β Oligomers have a Profound Detergent-like Effect on Lipid Membrane Bilayers by Atomic Force and Electron Microscopy[END_REF] Membrane disruption by incorporated-Aβ40 oligomers in 1-palmitoyl-2oleoyl-sn-glycero-3-phosphocholine (POPC) bilayers, [START_REF] Qiang | Time-Dependent Lipid Dynamics, Organization and Peptide-Lipid Interaction in Phospholipid Bilayers with Incorporated β-Amyloid Oligomers[END_REF] and incorporated-Aβ42 β-sheet tetramer and octamer in DPC micelles and lipid bilayer [START_REF] Ciudad | Aβ(1-42) Tetramer and Octamer Structures Reveal Edge Conductivity Pores as a Mechanism for Membrane Damage[END_REF] was also reported using solid-state and solution NMR spectroscopies, respectively. Melacini et al. evidenced Aβ40 fibrils inserted into small unilamellar vesicles and lipid bilayers. [START_REF] Ahmed | Atomic Resolution Map of the Soluble Amyloid beta Assembly Toxic Surfaces[END_REF] Combining wide-angle X-ray diffraction, solution NMR, dynamic light scattering, electron microscopy, ANS fluorescence and cell viability assays, they studied the toxic surfaces of Aβ40 that enable key interactions with the membrane. They found that toxicity scales with hydrophobic exposure of Aβ40 assemblies to solvent and their ability to interact with Aβ monomers and cell membranes. The central hydrophobic core (CHC, residues L17-A21) and the loop region (residues E22-K28) were found more accessible to monomer recognition in toxic Aβ assemblies. [START_REF] Ahmed | Atomic Resolution Map of the Soluble Amyloid beta Assembly Toxic Surfaces[END_REF] Many computational studies on the interactions of amyloid oligomers with membranes were reviewed in Ref. 7. Recent MD simulations explored the interactions of Aβ dimers at the surface or inserted into a membrane. Strodel et al. showed that Aβ42 dimer-membrane interactions do not change the lipid order parameters of the neuronal membrane. [START_REF] Fatafta | Amyloid-β Peptide Dimers Undergo a Random Coil to β-sheet Transition in the Aqueous Phase but not at the Neuronal Membrane[END_REF] Miller et al.

showed that inserted random/a-helix and U-shaped Aβ dimers disrupt a 1,2-dioleoyl-sn-glycero-3-phosphocholine (DOPC) membrane by locally decreasing the bilayer thickness. [START_REF] Press-Sandler | Molecular Insights into the Primary Nucleation of Polymorphic Amyloid β Dimers in DOPC Lipid Bilayer Membrane[END_REF] Li et al. investigated the effect of Aβ42 dodecamer and fibril on the surface of a neuronal membrane.

They showed that a dodecamer distorts the lipid membrane to a greater extent than fibril, and fibril does not change the membrane elastic properties. [START_REF] Nguyen | Amyloid β Dodecamer Disrupts the Neuronal Membrane More Strongly than the Mature Fibril: Understanding the Role of Oligomers in Neurotoxicity[END_REF] Strodel et al. also found membrane permeation induced by U-shaped fibrillar trimers and tetramers of human islet amyloid polypeptide in bilayers of anonic dipalmitoylphosphatidylglycerol (DPPG) lipids. [START_REF] Poojari | Membrane Permeation Induced by Aggregates of Human Islet Amyloid Polypeptides[END_REF] Using Aβ42 tetramers with a U-shaped topology spanning residues 17-42 inserted into the membrane and all N-terminal residues 1-16 located in the bulk solution at the top leaflet of the membrane, Strodel et al. reported on the translocation of a few water molecules through POPC, DPPC and palmitoyloleoylphosphatidylglycerol (POPG) bilayers using 500 ns MD simulations. Variation of bilayer thickness from 3.5 to 4.1 nm without any membrane destabilization was reported. [START_REF] Poojari | How the Amyloid-β Peptide and Membranes Affect Each Other: An Extensive Simulation Study[END_REF] The question we ask in this study is how a lipid bilayer is affected by placing a naked fibril cross section into the membrane and vice versa at an atomistic level. To this end, we explored the dynamics of a hexamer S-shaped Aβ42 fibrillar model inserted into a membrane (Figure 1) using atomistic MD simulations. The S-shaped fold, more disease-relevant than the U-shaped fold, displays a 3D structure with residues Q15-A42 forming a double horseshoe-like cross-βsheet structure, and residues D1-H14 partially ordered and residues A2-H6 in a β conformation. [START_REF] Wälti | Atomic-resolution Structure of a Disease-relevant Aβ(1-42) Amyloid Fibril[END_REF] We used this molecular structure rather than a prevalent Aβ40 fibril polymorph from Alzheimer's disease brain tissue [START_REF] Ghosh | Molecular Structure of a Prevalent Amyloidβ Fibril Polymorph from Alzheimer's Disease Brain Tissue[END_REF] because Aβ42 is known to be more toxic than Aβ40. Furthermore, the Sshape fold is shared by the amyloid fibrils of the human islet amyloid polypeptide (hIAPP). [START_REF] Gallardo | Fibril Structures of Diabetes-related Amylin Variants Reveal a Basis for Surface-templated Assembly[END_REF][START_REF] Röder | Cryo-EM Structure of Islet Amyloid Polypeptide Fibrils Reveals Similarities with Amyloid-β fibrils[END_REF]

RESULTS AND DISCUSSION

Impact of Membrane on Aβ42 Secondary Structure and Intermolecular

Interactions. We first report on the secondary structure averaged over all residues and along the amino acid sequence in the ss-NMR amyloid fibril and in the simulations. The β-strand/sheet content amounts to 53% in amyloid fibrils grown from AD brain with five β-regions spanning β1: A2-H6, β2: Q15-V18, β3: S26-K28, β4: A30-I32 and β5: V39-A42. [START_REF] Wälti | Atomic-resolution Structure of a Disease-relevant Aβ(1-42) Amyloid Fibril[END_REF] Averaged over the full trajectories, the (β-strand, helix, coil and turn) contents stabilize to (27.0, 0.0, 47.0 and 26.0%) in run 1 and (24.0, 1.0, 49.0 and 26.0%) in run 2. The β-strand, helix, coil and turn contents along the sequence averaged over the six chains are shown in Figure 2 using four consecutive timewindows. We see that the secondary structure along the sequence has almost converged between 500-2000 ns in each run. The membrane environment fully destabilizes the β1 and β3 strands (Figures 2A1-2B1) with a very high preference for coil/turn (Figures 2A3, 2A4, 2B3 and 2B4), modifies the β2 strand propensity with a 38% (run 1) and 74% (run 2) content, preserves the propensity of β4 strand (64% in run 1 and 75% in run 2) and reduces the β5 strand to residues 39 and 40. Run 1 is free of helix (Figure 2A2), while run2 displays 10% helix at residues 24-27 (Figure 2B2).

To explore the dynamics of the oligomer, we calculated the intermolecular side chain probability contact map between the trimer 1 and the trimer 2 as a function of time intervals (Figure 3). The ss-NMR fibril structure shows intermolecular side-chain M35-Q15 (region 1), and L17-M35 (region 2) contacts. [START_REF] Wälti | Atomic-resolution Structure of a Disease-relevant Aβ(1-42) Amyloid Fibril[END_REF] The MD-ensemble of conformations between 0 and 100 ns of both runs (Figure 3A-3E) reveal the two ss-NMR regions, but adds a few contact outliers between residues 34-36 (region 3). In run 1, the three contact regions are conserved in the time-interval 500-1500 ns (Figures 3B-3C), but only region 3 with the intermolecular M35-M35 non-native ss-NMR contact is present in the time-interval 1900-2000 ns (Figure 3D). This indicates that in run 1, the system remains a hexamer and has lost the symmetry of the ss-NMR fibril. In run 2, only the region 2 is conserved from 500 to 1500 ns (Figures 3F-3G), and the last 1900-2000 ns is free of any contacts (Figure 3H), indicating that the system has evolved to two separated trimers.

Overall, the results provide evidence of the significant change in the secondary and quaternary structures of Aβ oligomers embedded in a membrane model consisting of DOPC lipids.

The contribution of hexamer Aβ species to the development of AD pathology was reported by several studies. Bowers et al., using high resolution AFM to image oligomer populations, found that Aβ42 hexamers and dodecamers become the dominant oligomers after peptide solubilization, even at 1 μM concentrations and 5 min incubation times. [START_REF] Economou | Amyloid β-Protein Assembly and Alzheimer Disease: Dodecamers of Aβ42, but Not of Aβ40, Seed Fibril Formation[END_REF] Aβ42 pentamers/ hexamers were found to be the smallest detectable oligomers in the bulk solution using sedimentation velocity centrifugation, small angle neutron scattering and molecular modelling. [START_REF] Wolff | Aβ42 Pentamers/Hexamers are the Smallest Detectable Oligomers in Solution[END_REF] Hexamers in the bulk solution are largely disordered and random/coil. The presence of hexamers was also detected in transgenic mice brains exhibiting human Aβ pathology and in the cerebrospinal fluid (CSF) of AD patients. Furthermore, cell-derived hexameric Aβ seeded other human Aβ forms and enhanced amyloid deposition in vivo and neuronal toxicity in vitro. [START_REF] Vrancx | Mechanism of Cellular Formation and In Vivo Seeding Effects of Hexameric β-Amyloid Assemblies[END_REF] It is unknown how Aβ hexamers are structured in the interstitial space (ISS) and in the CSF.

The contribution of trimer Aβ species to the development of AD is also well documented.

Stable Aβ trimers were reported in pure buffer. [START_REF] Chen | Distinct Early Folding and Aggregation Properties of Alzheimer Amyloid-beta Peptides Abeta40 and Abeta42: Stable Trimer or Tetramer Formation by Abeta42[END_REF] Membrane-bound Aβ trimers correlated with toxicity towards cultured neurons, [START_REF] Jana | Membrane-bound Tetramer and Trimer Aβ Oligomeric Species Correlate with Toxicity towards Cultured Neurons[END_REF] and an Aβ trimer is likely the minimal oligomer size to insert into a membrane. [START_REF] Jang | Mechanisms for the Insertion of Toxic, Fibril-like β-Amyloid Oligomers into the Membrane[END_REF] Finally, Aβ trimers were found to induce progressive loss of hippocampal synapses. [START_REF] Shankar | Natural Oligomers of the Alzheimer Amyloid-beta Protein Induce Reversible Synapse Loss by Modulating an NMDA-type Glutamate Receptor-dependent Signaling Pathway[END_REF]

Impact of Membrane on Aβ42

Intramolecular Side Chain Interactions. The intramolecular side chain contact maps were also compared. The ss-NMR fibril structure features interactions through the side chains L17, F19, F20 and V24 with the side chains A30, I32 and L34 stabilizing the packing between the β-sheets β2 and β3, and through the side chain I31 with the side chains V36, V39 and I41 stabilizing the β-sheets β4 and β5. [START_REF] Wälti | Atomic-resolution Structure of a Disease-relevant Aβ(1-42) Amyloid Fibril[END_REF] The conformations between 0 and 100 ns (Figures 4A and4E) reveal that, averaged over the two runs, the hydrophobic interactions stabilizing the C-terminus (residues 30-42) are still formed, but with lower lifetimes, contacts I31-V39 (probability of 75%), I31-I41 (probability of 36%) and I31-V36 (probability of 45%). The hydrophobic contacts stabilizing the interface β-sheets β2 and β3 are also partially recovered with three contacts A30-F19 (40%), I32-F19 (probability of 65%), A30-F20 (probability of 45%). In the three time-windows 500-1000, 1000-1500 and 1900-2000 ns, more non ss-NMR and a few ss-NMR contacts are formed (Figures 4B-4D and4F-4H). The C-terminus is still stabilized by the ss-NMR contacts I31-V39 and I31-V36 in the two runs. The interface between β2 and β3 is stabilized by the ss-NMR contacts V24-A30 and F20-A30 in run 2 and by no contacts with a probability > 30% in run 1. Overall, the intrapeptide contact maps provide evidence of the significant change in the tertiary structures of Aβ oligomers embedded in a membrane model consisting of DOPC lipids.

Insertion Depths of Aβ42

Residues into the Membrane. Next, we evaluated the time evolution of the insertion depths of residues into the membrane. There is a clear difference in the insertion depth profiles between the two runs (Figure 5). In the ensemble of conformations covering the 0-100 ns time interval of the two runs, residues D1-G9 and A21-G25 are exposed to solvent and all other residues are inserted into the membrane (Figures 5A and5C). The dynamics modifies the solvent-exposure of the residues, and the final 1900-2000 ns reveals distinct insert depths of the residues belonging to trimers 1 and 2. In run 1, residues D1-Y10(V12) and A21-G25 (A1-G29) remain exposed to solvent, the CHC is much less inserted into the membrane, residue L17 going from 0.9-1.4 nm (Figure 4A) to 1.8 nm (Figure 5B). In run 1 residues A30-A42 of trimer 2 are fully inserted into the membrane, while only residues A30-V40 of trimer 1 are fully inserted.

In run 2, shown in Figure 5D, all residues D1-G29, I41 and A42 of trimer 2 have insertion depths > 2.1 nm, and the C-terminus residues 30-40 are inserted into the membrane. In contrast in trimer 1 of run 2, the CHC has insertion depths < 2.1 nm and all residues D22-A42, except L34 and M35, are solvent-exposed. Figure S1 in Supplementary Information (SI) shows that the insertion depth of each residue, averaged over 1900-2000 ns in runs 1 and 2, does not change the results if we consider each chain.

The simulations also reveal that residues 1-10 in both trimers of run 1 and in trimer 1 of run 2, and residues 1-16 in trimer 2 of run 2 are solvent-exposed due to the complex network of interactions between the Aβ42 peptides, lipid bilayers and solvent. This result is consistent with many simulations of Aβ42 oligomers partially embedded in the membrane. [START_REF] Liu | Molecular Simulation Aspects of Amyloid Peptides at Membrane Interface[END_REF][START_REF] Nguyen | Amyloid Oligomers: A Joint Experimental/Computational Perspective on Alzheimer's Disease, Parkinson's Disease, Type II Diabetes, and Amyotrophic Lateral Sclerosis[END_REF][START_REF] Feng | Study of Structural Stability and Damaging Effect on Membrane for four Aβ42 Dimers[END_REF] The N-termini appear turn/coil, and we do not find any evidence of the formation of a β-hairpin in the Nterminus, as reported for Aβ42 dimers in DOPC bilayers. [START_REF] Press-Sandler | Molecular Insights into the Primary Nucleation of Polymorphic Amyloid β Dimers in DOPC Lipid Bilayer Membrane[END_REF] The high solvent accessibility of the Nterminus residues 1-16 may help understand the linking of liposomes through Aβ oligomers, as assessed by cryo-electron tomography. [START_REF] Tian | 3D-visualization of Amyloid-β Oligomer Interactions with Lipid Membranes by Cryo-electron Tomography[END_REF]

Membrane Deformation and Water

Permeability. An important question is whether the membrane deforms and disrupts during the dynamics. Figure 6 reports on the positional probability of all lipid atoms projected on the (y,z) plane where x is the direction of the fibril axis in run 1 as a function of time. We see insertion of lipid atoms within the bilayer at y distances varying between 6 and 12 nm and z distances between 4 and 8 nm with a probability of 0.2 in the time intervals 500-1000 ns, 1000-1500 ns and 1900-2000 ns in contrast to the CHARMM-GUI starting conformation. A very similar positional probability of lipids is observed in run 2 shown in Figure 7, albeit the positional probabilities are not superposable. Differences in the membrane dynamics between the two runs are observed in the curvature and in the thinning and thickening of the membrane. The z dimension of the membrane varies between the minimal and maximal values of 2.8 and 4.7 nm in run 1 and 2.5 and 5 nm in run 2 if we consider a lipid probability of 0.3. Overall, the membrane disruption is therefore much more pronounced than that observed in simulations of Aβ42 dimer [START_REF] Press-Sandler | Molecular Insights into the Primary Nucleation of Polymorphic Amyloid β Dimers in DOPC Lipid Bilayer Membrane[END_REF] and tetramer 38 in a membrane.

As a result of membrane deformation, we anticipated that waters molecules would insert into the lipid bilayer. This is supported by the positional probabilities of water molecules as a function of time intervals in particular at y distances between 8 and 12 nm and z distances between 4.5 and 7 nm in contrast to the CHARMM-GUI starting conformation. The ensemble is heterogeneous with no true dominant interpeptide and intrapeptide contacts. The MD structures of the oligomer at times 500, 1000, 1500 and 2000 ns in runs 1 and 2, shown in Figure S2, illustrate the wide variety of orientations and packings of the system from the starting structure.

The relationship between increased solvent exposure of hydrophobic residues, increased toxicity and disruption of lipid membranes was discussed for several amyloid protein systems, 51- on Aβ42 to reveal this relationship at atomistic level. In this study, we find a structural model (trimer 2 in run 2) disrupting lipid integrity that is fully consistent with the map of the soluble Aβ40 assembly toxic surface proposed by Melacini et al., namely the solvent-exposure of the CHC region and of the loop region (residues E22-K28) and shielding of residues A30-V40 from solvent. [START_REF] Ahmed | Atomic Resolution Map of the Soluble Amyloid beta Assembly Toxic Surfaces[END_REF] But all configurations disrupting membrane permeability do not satisfy these conditions as the CHC and C-terminus regions may be solvent-protected and solvent-exposed, respectively. Clearly exposing the C-terminus to solvent might be a template to promote fibril formation at the surface of the membrane.

It is important to discuss the initial membrane-embedded model with charged N-termini at both leaflets of the membrane. As assessed by various in vitro experiments, soluble oligomers in the extracellular space can adsorb and penetrate the cell membranes. [START_REF] Nguyen | Amyloid Oligomers: A Joint Experimental/Computational Perspective on Alzheimer's Disease, Parkinson's Disease, Type II Diabetes, and Amyotrophic Lateral Sclerosis[END_REF] β-barrel models suggested by Lal et al. have all N-terminus residues located in the bulk solution at the upper membrane leaflet, [START_REF] Jang | Mechanisms for the Insertion of Toxic, Fibril-like β-Amyloid Oligomers into the Membrane[END_REF][START_REF] Jang | Disordered Amyloidogenic Peptides May Insert into the Membrane and Assemble into Common Cyclic Structural Motifs[END_REF] and Aβ42 fibrils can displace the upper leaflet of the membrane. In this context, using cryo-electron tomography, Viles et al. showed that Aβ42 oligomers are able to embed within the upper leaflet of the bilayer. [START_REF] Tian | 3D-visualization of Amyloid-β Oligomer Interactions with Lipid Membranes by Cryo-electron Tomography[END_REF] Additionally, prior to Aβ release in the extracellular space after amyloid precursor cleavage, the residues 1-27 are in the bulk solution at the upper membrane leaflet and residues 28-42 in the interior of the membrane. [START_REF] Nguyen | Amyloid Oligomers: A Joint Experimental/Computational Perspective on Alzheimer's Disease, Parkinson's Disease, Type II Diabetes, and Amyotrophic Lateral Sclerosis[END_REF] There is however evidence that Aβ oligomerization begins intracellularly in cells derived from human brain, [START_REF] Walsh | The Oligomerization of Amyloid beta-protein Begins Intracellularly in Cells Derived from Human Brain[END_REF] and there is a generation of intracellular Aβ pool. [START_REF] Sannerud | Restricted Location of PSEN2/γ-Secretase Determines Substrate Specificity and Generates an Intracellular Aβ Pool[END_REF] In this case all N-terminus residues remain intracellular and all residues 17-42 need to cross the lower leaflet of the membrane. Finally, neutron reflection from tethered membranes showed that Aβ oligomers insert into the bilayer affecting both membrane leaflets, [START_REF] Valincius | Soluble Amyloid Beta-oligomers Affect Dielectric Membrane Properties by Bilayer Insertion and Domain Formation: implications for Cell Toxicity[END_REF] and NMR structures of incorporated-Aβ42 β-sheet tetramer and octamer showed N-terminal residues in both the upper and lower leaflets of the membrane, [START_REF] Ciudad | Aβ(1-42) Tetramer and Octamer Structures Reveal Edge Conductivity Pores as a Mechanism for Membrane Damage[END_REF] indicating that polar residues can cross the membrane. Overall, our initial model with charged N-terminus Aβ42 tails located in the bulk solution at the upper and lower leaflets of the membrane is supported experimentally and also corresponds to a mixture of both extracellular and intracellular aggregates, respectively.

CONCLUSIONS

The topic of Aβ embedding in membranes is of intense interest as the connection to Alzheimer's disease pathology is quite clear. In this study, we have determined the dynamics of a cross-β Aβ42 hexamer inserted into a neutral DOPC lipid bilayer. Using two extensive atomistic MD simulations, we showed significant changes in secondary, tertiary and quaternary structures of the cross-β model leading to two types of mainly random coil oligomers: hexamer and trimers.

Both species promote an increase of the hydrophobic surface of Aβ42 peptides exposed to solvent, insertion of water molecules and lipids into the bilayer and membrane destabilization.

Membrane-embedded Aβ species with increase in membrane conductance, which do not display Thioflavin T fluorescence, consist of a wide range of conformations and oligomer sizes, some forming perfect β-barrel pores [START_REF] Österlund | Native Ion Mobility-Mass Spectrometry Reveals the Formation of β-Barrel Shaped Amyloid-β Hexamers in a Membrane-Mimicking Environment[END_REF][START_REF] Jang | Mechanisms for the Insertion of Toxic, Fibril-like β-Amyloid Oligomers into the Membrane[END_REF][START_REF] Jang | Disordered Amyloidogenic Peptides May Insert into the Membrane and Assemble into Common Cyclic Structural Motifs[END_REF] or six stranded β-sheet and β-sandwich structures, [START_REF] Ciudad | Aβ(1-42) Tetramer and Octamer Structures Reveal Edge Conductivity Pores as a Mechanism for Membrane Damage[END_REF] and others adopting spherical species without any evidence of discrete channel or pore formation. [START_REF] Kayed | Permeabilization of Lipid Bilayers is a Common Conformation-dependent Activity of Soluble Amyloid Oligomers in Protein Misfolding Diseases[END_REF] Our computational work provides evidence of metastable hexamer and trimers characterized by the N-terminal residues located in the top and bottom leaflets of membrane with a minimal lifetime of several microseconds. These oligomers induce membrane permeability and water permeation that will facilitate the transport of hydrated ions from the intracellular to the extracellular space. These two assemblies, which are mainly coil/turn, are characterized by different solvent exposures of the residues. In the hexamer, we observe membrane insertion of the CHC region and residues 29-41. In the two separated trimers, we observe an increase of solvent exposure of the hydrophobic CHC and C-terminus residues.

Whether there is insertion of individual peptides which then form trimers and hexamers in the membrane or insertion of preformed trimers remains an open question, the insertion process depending in part on the physicochemical properties of the membrane. [START_REF] Jang | Mechanisms for the Insertion of Toxic, Fibril-like β-Amyloid Oligomers into the Membrane[END_REF][START_REF] Arce | Polymorphism of Amyloid β Peptide in Different Environments: Implications for Membrane Insertion and Pore Formation[END_REF][START_REF] Kleinschmidt | Folding of β-barrel Membrane Proteins in Lipid Bilayers -Unassisted and Assisted Folding and Insertion[END_REF] Homogeneous DOPC bilayers are a simple model to represent physiological cells. It is known that POPC-POPG bilayers, and brain total lipid extracts impact Aβ adsorption, and the macroscopic ion flux across Aβ channels in membranes, respectively. [START_REF] Yu | Molecular Interactions of Alzheimer Amyloidβ Oligomers with Neutral and Negatively Charged Lipid Bilayers[END_REF][START_REF] Lee | Amyloid β Ion Channels in a Membrane Comprising Brain Total Lipid Extracts[END_REF] Pure DOPC bilayers also ignore the contribution of cholesterol and GM1 lipids, [START_REF] Nguyen | Amyloid Oligomers: A Joint Experimental/Computational Perspective on Alzheimer's Disease, Parkinson's Disease, Type II Diabetes, and Amyotrophic Lateral Sclerosis[END_REF][START_REF] Nguyen | Tetrameric Aβ40 and Aβ42 β-Barrel Structures by Extensive Atomistic Simulations. I. In a Bilayer Mimicking a Neuronal Membrane[END_REF][START_REF] Fatafta | Amyloid-β Peptide Dimers Undergo a Random Coil to β-sheet Transition in the Aqueous Phase but not at the Neuronal Membrane[END_REF][START_REF] Nguyen | Self-Assembly of Amyloid-Beta (Aβ) Peptides from Solution to Near In Vivo Conditions[END_REF] and cells have extracellular matrices that result in very different behaviors compared to vesicles, liposomes and planar lipid bilayers in the absence of membrane protein receptors. [START_REF] Sakono | Amyloid Oligomers: Formation and Toxicity of Abeta Oligomers[END_REF] Going beyond DOPC bilayers and exploring computationally other β-rich aggregates in the presence of protein receptors, [START_REF] Sakono | Amyloid Oligomers: Formation and Toxicity of Abeta Oligomers[END_REF] and free phospholipids [START_REF] Ngo | Cholesterol Molecules Alter the Energy Landscape of Small Aβ1-42 Oligomers[END_REF][START_REF] Nguyen | Effect of Cholesterol Molecules on Aβ1-42 Wild-Type and Mutants Trimers[END_REF][START_REF] Chakravorty | Free Gangliosides Can Alter Amyloid-β Aggregation[END_REF] in the extracellular space could offer a deeper understanding of the relationship between amyloid oligomers and membrane permeability in order to develop new drugs. [START_REF] Errico | Quantitative Measurement of the Affinity of Toxic and Nontoxic Misfolded Protein Oligomers for Lipid Bilayers and of its Modulation by Lipid Composition and Trodusquemine[END_REF][START_REF] Doig | Why Is Research on Amyloid-β Failing to Give New Drugs for Alzheimer's Disease?[END_REF] In this context, new experimental methods of biomolecular processes in membrane spanning the microsecond time resolution and detecting populations of intermediates below the percent range would be very useful to validate MD-generated models.

MATERIAL AND METHODS

The initial structure of the hexamer Aβ42 fibril from solid state nuclear magnetic spectroscopy (ss-NMR) was taken from PDB entry 2NAO. [START_REF] Berman | The Protein Data Bank[END_REF] The fibril was inserted into the membrane of dimension (x =20 nm, y = 20 nm, z = 4.2 nm), and its orientation was selected to maximize the total number of hydrophobic CHC (L17-A21) and C-terminal (A30-A42) residues inserted into the membrane and the total number of polar/charged residues exposed to solvent. The bilayer, consisting of 1128 DOPC (neutral) lipids, was immersed in a box of 20 nm, 20 nm and 12.1 nm with about 99000 water molecules. Simulations were performed at pH 7. Aβ42 peptides have NH3 + and CO2 -termini, deprotonated Glu and Asp, protonated Arg and Lys, and neutral His with a protonated Nε atom. The system was neutralized by sodium ions leading to a total number of about 456000 atoms. CHARMM-GUI was used to set up the system. [START_REF] Lee | CHARMM-GUI Input Generator for NAMD, GROMACS, AMBER, OpenMM, and CHARMM/OpenMM Simulations Using the CHARMM36 Additive Force Field[END_REF] Figure 1 shows the initial structure of the hexamer generated by CHARMM-GUI with C2 symmetric molecules. Residues D1-E11 and E22-V24 are exposed to solvent, residues Q15-F20 and K28-A42 are embedded in the membrane, and residues V12-H14 and G25-N27 are in the subspace formed by head and glycerol lipid groups. It has to be noted the N-terminal residues of three chains are located in aqueous solution at the upper leaflet of the membrane (these three chains will be referred to as trimer 1) and the N-terminal residues of other three chains (referred to as trimer 2) are located in aqueous solution at the lower leaflet of the membrane.

The CHARMM36 force field was used for the Aβ42 peptides and the lipids, [START_REF] Pastor | Development of the CHARMM Force Field for Lipids[END_REF] and water molecules were modelled by the CHARMM-modified TIP3P model. [START_REF] Huang | CHARMM36m: An Improved Force Field for Folded and Intrinsically Disordered Proteins[END_REF] CHARMM36m, refined for intrinsically disordered proteins in the bulk solution, was not used because its relevance on transmembrane proteins remains to be determined. [START_REF] Huang | CHARMM36m: An Improved Force Field for Folded and Intrinsically Disordered Proteins[END_REF] Two MD simulations of 2 μs were performed at 303 K (run 1) and 310 K (run 2) in the NPT ensemble as a single simulation, regardless of its length, is not sufficient to show that the results are reliable or may be an artefact.

The choice of 310 K is also motivated by the fact that 2 μs at 303 K might be too short to see major oligomeric changes and two simulations at the same temperature starting from the same configurations with different starting velocities might not be adequate. Also, the use of two slightly different temperatures allows a better sampling of the configurational space. In this study, we used GROMACS (version 2019.4), [START_REF] Van Der Spoel | GROMACS: Fast, Flexible, and Free[END_REF] and a time step of 2 fs with the lengths of bonds with hydrogen atoms constrained using the LINCS algorithm. [START_REF] Hess | LINCS: a Linear Constraint Solver for Molecular Simulations[END_REF] Minimization and equilibration of 1 ns were performed with restraints on the backbone atoms of Aβ42 peptides. Temperature and pressure were controlled using the Noose-Hoover algorithm with a coupling constant of 1 ps, [START_REF] Evans | The Nose-Hoover Thermostat[END_REF] and the semi-isotropic Parrinello-Rahman algorithm with a coupling constant of 5 ps, [START_REF] Parrinello | Polymorphic Transitions in Single Crystals: A new Molecular Dynamics Method[END_REF] respectively. Electrostatic interactions were calculated using the particle mesh Ewald method with a cutoff of 1.2 nm, and Van der Waals interactions used forces switched to zero from 1.0 to 1.2 nm, as implemented in CHARMM-GUI. [START_REF] Lee | CHARMM-GUI Input Generator for NAMD, GROMACS, AMBER, OpenMM, and CHARMM/OpenMM Simulations Using the CHARMM36 Additive Force Field[END_REF] Prior to analysis, we emphasize that the protonation states we assigned assume canonical (solution) pKa values. While most of the polar/charged residues are initially exposed to solvent and notably the N-terminal residues and residues E22 and D23, K28 is embedded in the membrane. It is known that Lysine can shift its pKa down to 6 in hydrophobic environments, e.g., when fully inserted in a membrane. [START_REF] Panahi | Membrane Environment Modulates the pKa Values of Transmembrane Helices[END_REF] In our starting structure, as K28 partially feels the presence of the head and glycerol lipid groups, its pKa value is unknown. Assuming therefore that K28 is protonated is an approximation which might not be totally appropriate. In principle, we can solve this issue by using pH-replica exchange MD simulations, but this approach is not computationally feasible for a hexamer of 252 amino acids partially embedded in an atomistic lipid bilayer.

To characterize the conformational ensemble, secondary structures were determined using the DSSP method. [START_REF] Kabsch | Dictionary of Protein Secondary Structure: Pattern Recognition of Hydrogen-bonded and Geometrical Features[END_REF] Conformations were characterized by their intra and intermolecular (between trimer 1 and trimer 2) contact maps. A contact between two side chains was considered formed if there is at least one distance < 0.45 nm between all heavy atoms. We calculated the average depth of residues 1-42 into the lipid bilayer. One residue was considered inserted if the depth is < 2.1 nm. The insertion depth was calculated as the z component of the distance between the center of mass of each residue (including all backbone and side chain atoms) and the center of mass of all lipid atoms of the membrane. [START_REF] Nguyen | Molecular Dynamics Simulations of the Tau Amyloid Fibril Core Dimer at the Surface of a Lipid Bilayer Model: I. In Alzheimer's Disease[END_REF] Finally, we also calculated the positional probabilities of the lipids and water molecules with respect to the membrane surface.

ASSOCIATED CONTENT Supporting Information.

The Supporting Information is available free of charge at https://pubs.acs.org/doi/ Insertion depths of residues 1-42 belonging to each chain of the system in the time-interval 1900-2000 ns of run 1 and run 2, and snapshots of the system as a function of time in runs 1 and 2 are reported.

Figure 1 .

 1 Figure1. Initial CHARMM-GUI conformation of the system. View in the x direction parallel to the fibril axis. We show the ribbon structure of the hexamer, the phosphate atoms (tan) and the O21 and O22 atoms of the glycerol (black). Also shown are the Ca atom of V12 (yellow), Q15 (blue), V18 (orange), E22 (purple) and K28 (red).

Figure 2 .

 2 Figure2. Time-averaged secondary structure contents along the sequence using four timewindows of run 1 at 303 K (A1, A2, A3, A4) and run 2 at 310 K (B1, B2, B3, B4) using 0-500 ns (orange), 500-1000 ns (red), 1000-1500 ns (green) and 1500-2000 ns (blue).

Figure 3 .

 3 Figure 3. Intermolecular side chain contact map probabilities as a function of time intervals for the two runs, run 1 at 303 K and run 2 at 310 K.

Figure 4 .

 4 Figure 4. Intramolecular side chain contact map probabilities as a function of time intervals for the two runs (run 1 at 303 K and run 2 at 310 K) separated at least by four residues.

Figure 5 .

 5 Figure 5. Insertion depths of the residues belonging to trimer 1 and trimer 2 in two time-intervals, 0-100 ns and 1900-2000 ns. (A) and (B) refer to run 1 at 303 K, and (C) and (D) to run 2 at 310 K.

Figure 6 .

 6 Figure 6. Run 1 (303 K) simulation results of water (A) and lipid (B) positional probabilities projected on the (y,z) plane where x is fibril axis. We show the CHARMM-GUI values, free of equilibration, and the averaged values using three time-intervals.

Figure 7 .

 7 Figure 7. Run 2 (310 K) simulation results of water and lipid positional probabilities projected on the (y,z) plane where x is fibril axis. We show the CHARMM-GUI values, free of equilibration, and the averaged values using three time-intervals.

Notes

The authors declare no competing financial interest. ACKNOWLEDGMENTS. This work was supported by grants from the French IDRIS and CINES computer centers, and the "Initiative d'Excellence" program from the French State, Grant "DYNAMO", ANR-11-LABX-0011.


Student Member, IEEE Laurane Charrier

Member, IEEE Yajing Yan

Senior Member, IEEE Emmanuel Trouvé

Elise Colin Koeniguer

Jérémie Mouginot

Romain Millan

Fusion of multi-temporal multi-sensor velocities using temporal closure of fractions of displacements

Keywords: velocity time series, multi-sensor fusion, temporal inversion, glacier, intra-annual variations

Numerous glacier velocity observations, derived from spaceborne imagery, are available online, but it remains difficult to analyze them because they are measured with different temporal baselines, by various sensors. In this study, we propose a novel formulation of the temporal closure to fuse multi-temporal multi-sensor velocity observations without prior information on the displacement behavior and the data uncertainty. We establish a system of linear equations between combinations of displacement observations and fractions of estimated displacements. The proposed approach provides a velocity time-series with a regular and optimal temporal sampling, the latter representing a compromise between the temporal resolution and the signal-tonoise ratio. The proposed approach is first evaluated on synthetic datasets and second on Sentinel-2 and Venµs velocity observations over the Fox glacier in New Zealand. The results show the intraannual variability of Fox glacier surface velocity with a reduced uncertainty and complete temporal coverage.

I. INTRODUCTION

W HILE numerous glacier velocity observations are avail- able online [START_REF] Gardner | Increased west antarctic and unchanged east antarctic ice discharge over the last 7 years[END_REF]- [START_REF] Millan | Mapping surface flow velocity of glaciers at regional scale using a multiple sensors approach[END_REF] with a temporal resolution up to 2 days and a spatial sampling up to 50 m [START_REF] Millan | Mapping surface flow velocity of glaciers at regional scale using a multiple sensors approach[END_REF], [START_REF] Millan | Ice velocity and thickness of the world's glaciers[END_REF], intraannual variations of glacier velocity remain poorly understood at regional scales, especially in mountain areas. This is due to the difficulty to analyze scene-pair velocity observations since they span different temporal baselines [START_REF] Charrier | Extraction of velocity time series with an optimal temporal sampling from displacement observation networks[END_REF], are based on images from different sensors, and still contain noise, artifacts, and missing data when outliers have been filtered out. Therefore, a question arises: how to exploit the amount of available multisensor and multi-temporal glacier velocities in order to study intra-annual velocity variations with a reasonable uncertainty at a global scale?

A way to fuse multi-temporal velocity observations is to use a Temporal Inversion (TI) based on the temporal closure of the displacement observation networks, inspired from the Small BAseline Subset (SBAS) in the InSAR community. Currently, most of the TI approaches consider mono-sensor datasets [START_REF] Charrier | Extraction of velocity time series with an optimal temporal sampling from displacement observation networks[END_REF]- [START_REF] Lacroix | Self-entrainment motion of a slow-moving landslide inferred from landsat-8 time series[END_REF] or apply TI separately to each sensor datasets [START_REF] Ding | Quantifying the spatio-temporal patterns of dune migration near minqin oasis in northwestern china with time series of landsat-8 and sentinel-2 observations[END_REF], [START_REF] Dille | When image correlation is needed: Unravelling the complex dynamics of a slow-moving landslide in the tropics with dense radar and optical time series[END_REF]. A few studies imply multi-sensor datasets inside the same system of linear equations [START_REF] Bontemps | Inversion of deformation fields time-series from optical images, and application to the long term kinematics of slow-moving landslides in peru[END_REF], [START_REF] Pepe | The use of c-/x-band time-gapped sar data and geotechnical models for the study of shanghai's ocean-reclaimed lands through the sbas-dinsar technique[END_REF] whereas it could improve the results by adding more redundancy as stressed L. Charrier, Y. Yan, Emmanuel Trouvé are with LISTIC, Université Savoie Mont Blanc, Annecy, France. E. Colin Koeniguer and L. Charrier are with DTIS, ONERA, Paris, France. J. Mouginot is with IGE, Université Grenoble Alpes, Grenoble, France. R. Millan is with Department of Geosciences and Natural Resource Management, University of Copenhagen, Copenhagen K, Denmark.

by [START_REF] Dille | When image correlation is needed: Unravelling the complex dynamics of a slow-moving landslide in the tropics with dense radar and optical time series[END_REF]. In [START_REF] Bontemps | Inversion of deformation fields time-series from optical images, and application to the long term kinematics of slow-moving landslides in peru[END_REF], [START_REF] Pepe | The use of c-/x-band time-gapped sar data and geotechnical models for the study of shanghai's ocean-reclaimed lands through the sbas-dinsar technique[END_REF] irregular Leap Frog (LF) time series (time series of displacements between consecutive dates) are inverted. This kind of time series are adapted to reconstruct a Cumulative Displacement (CD) time series afterward. While CD time series are useful to study slow-moving landslides or dunes [START_REF] Lacroix | Self-entrainment motion of a slow-moving landslide inferred from landsat-8 time series[END_REF], [START_REF] Ding | Quantifying the spatio-temporal patterns of dune migration near minqin oasis in northwestern china with time series of landsat-8 and sentinel-2 observations[END_REF], [START_REF] Bontemps | Inversion of deformation fields time-series from optical images, and application to the long term kinematics of slow-moving landslides in peru[END_REF], it is more adapted to convert the LF displacement time series into LF velocity time series for fast-moving targets with strong seasonality patterns, such as glaciers [START_REF] Altena | Extracting recent short-term glacier velocity evolution over southern alaska and the yukon from a large collection of landsat data[END_REF]. Because velocities correspond to the average of instantaneous velocity over the temporal baseline, the temporal sampling of the velocity series should be regular (all the LF velocities should span the same temporal baseline). Hence, [START_REF] Charrier | Extraction of velocity time series with an optimal temporal sampling from displacement observation networks[END_REF] proposed to include combinations of velocity observations in the TI to obtain a Regular LF (RLF) with a given temporal sampling. This approach, which will be referred as Temporal Inversion using COmbinations of displacements (TICO), allows obtaining a RLF with a temporal sampling matching a compromise between temporal resolution and signal-overnoise ratio.

This study is an extension of [START_REF] Charrier | Extraction of velocity time series with an optimal temporal sampling from displacement observation networks[END_REF] in order to obtain a RLF with an optimal temporal sampling by taking advantage of every multi-sensor observation, even if outliers have been removed beforehand. For multi-sensor datasets, each sensor does not acquire images on the same temporal grid: hence, there is less redundancy and more disconnected subsets of observed displacements. This is even more true when outliers have been removed beforehand. Therefore, we introduce fractions of displacement in the TICO formulation proposed in [START_REF] Charrier | Extraction of velocity time series with an optimal temporal sampling from displacement observation networks[END_REF]. This extension is called Temporal Inversion using COmbinations of Fractions of displacements (TICOF).

Besides, evaluation of observation uncertainties remain crucial to invert velocity time series [START_REF] Charrier | Extraction of velocity time series with an optimal temporal sampling from displacement observation networks[END_REF], [START_REF] Altena | Extracting recent short-term glacier velocity evolution over southern alaska and the yukon from a large collection of landsat data[END_REF], [START_REF] Bontemps | Inversion of deformation fields time-series from optical images, and application to the long term kinematics of slow-moving landslides in peru[END_REF]. The shape of the similarity function used in image matching appear to be a robust estimation of precision for both off-and onglacier pixels [START_REF] Charrier | Extraction of velocity time series with an optimal temporal sampling from displacement observation networks[END_REF], [START_REF] Bontemps | Inversion of deformation fields time-series from optical images, and application to the long term kinematics of slow-moving landslides in peru[END_REF], [START_REF] Burgess | Flow velocities of alaskan glaciers[END_REF], [START_REF] Altena | Correlation dispersion as a measure to better estimate uncertainty of remotely sensed glacier displacements[END_REF]. However, most of the online datasets do not provide the full similarity function values [START_REF] Gardner | Increased west antarctic and unchanged east antarctic ice discharge over the last 7 years[END_REF]- [START_REF] Millan | Mapping surface flow velocity of glaciers at regional scale using a multiple sensors approach[END_REF]. Hence, there is a need to get a robust estimation of uncertainty without any information on surface correlation values.

Therefore, this study has two objectives: 1) to take into account that velocity observations have different uncertainties which are not necessarily provided by the authors 2) to propose a new system to invert a RLF time series with an optimal temporal sampling from multi-sensor datasets, which can present missing data. We apply the proposed TICOF approach to Sentinel-2 and Venµs velocity observations over the Fox glacier located in the Southern Alps of New Zealand.

II. STUDY AREA AND DATA

The study focuses on the Fox glacier, in the Southern Alps of New Zealand. This glacier is a fast-flowing glacier with a strong seasonality [START_REF] Millan | Mapping surface flow velocity of glaciers at regional scale using a multiple sensors approach[END_REF], [START_REF] Purdie | Seasonal variation in ablation and surface velocity on a temperate maritime glacier: Fox glacier, new zealand[END_REF], situated under a maritime climate (Fig. S1). The considered velocity observations are obtained from Sentinel-2 and Venµs optical images. Sentinel-2 has a repeat cycle of 5 days and Venµs 2 days. Displacement offsets have been computed by [START_REF] Millan | Mapping surface flow velocity of glaciers at regional scale using a multiple sensors approach[END_REF] using a modified version of the cross-correlation algorithm AMPCOR from ROI PAC [START_REF] Millan | Mapping surface flow velocity of glaciers at regional scale using a multiple sensors approach[END_REF], [START_REF] Millan | Ice velocity and thickness of the world's glaciers[END_REF]. Outliers have been removed by rejecting offsets that deviate more than three pixels from the median offset computed over a spatial window of 9 × 9 pixels [START_REF] Millan | Mapping surface flow velocity of glaciers at regional scale using a multiple sensors approach[END_REF]. Velocity observations span temporal baselines from 5 to 100 days and from 330 to 400 days. The spatial sampling of both Sentinel-2 and Venµs velocity maps is 50 m. The central dates of Sentinel-2 velocity observations range from July 21, 2016 to October 16, 2018, whereas those of Venµs velocity observations cover the period of January 26, 2018 to May 21, 2019. The mean percentage of remaining values after outliers removal is 26 % for Sentinel-2 and 14 % for Venµs as illustrated in Fig. S2.

III. METHOD A. Data quality indicators

In many cases, uncertainties are not provided with velocity observations. In order to quantify the quality of velocity observations without any information given by the authors, two criteria are defined. The first one relies on the spatiotemporal coherence of velocity directions, the second on the spatio-temporal coherence of respectively x-and y-velocity component values [START_REF] Lüttig | A combined approach for filtering ice surface velocity fields derived from remote sensing methods[END_REF]:

1) The cosine of the angle between each displacement vector and the spatio-temporal median [START_REF] Burgess | Flow velocities of alaskan glaciers[END_REF] vector, called Median Angle (MA) here, is defined as:

cos(⃗ u i,j,t , ⃗ v w i,j ) = ⃗ u i,j,t • ⃗ v w i,j ∥⃗ u i,j,t ∥ ∥⃗ v w i,j ∥ (1)
where ⃗ u i,j,t is the displacement vector at the pixel (i, j) and time t and ⃗ v w i,j is the median vector in time over a spatial window w, here set to 3×3 pixels, centered on the pixel (i, j). If the two vectors are co-linear, as expected, this criterion will be equal to 1. The negative values are set to 0 to penalize vectors which have a direction from 90 • to 180 • away from the median direction.

2) The modified z-score, a standardized and robust score, measuring the outlier strength [START_REF] Maronna | Robust statistics: theory and methods (with R)[END_REF], is computed respectively for the x and y component of the velocity vector as:

MZscore = x i,j,t -xw i,j 1.483 MAD (2)
where x i,j,t stands for the x or y component of the velocity at the pixel (i, j) and time t, xw i,j is the median of the x or y component of the velocity in time over a spatial window w, here set to 3×3 pixels, centered on the pixel (i, j). MAD is the Median Absolute Deviation defined as: median|x i,j,t -xw i,j |. Note that the MZscore considers a robust estimate of the dispersion: 1.483 MAD, unlike the standard deviation which is not robust to outliers [START_REF] Gardner | Increased west antarctic and unchanged east antarctic ice discharge over the last 7 years[END_REF], [START_REF] Millan | Mapping surface flow velocity of glaciers at regional scale using a multiple sensors approach[END_REF].

The final confidence value of the two velocity components is the same for a given pixel and a couple of dates. It is computed as the multiplication of the three criteria scaled between 0 and 1: the scaled M Zscore for x and y components of the velocity and the MA. This confidence value will be used in the inversion as an initial weight since it indicates data quality. The classical TI formulation is AX = Y , where Y is a vector of dimension n containing the displacement observations d ti,tj between dates t i and t j , X is a vector of dimension p representing the LF displacement estimation dτ k ,τ k+∆τ between dates τ k and τ k+∆τ . ∆τ days is the temporal sampling which is identical for every LF displacement. A is the design matrix of dimension n×p which links the vectors X and Y . As shown in [START_REF] Charrier | Extraction of velocity time series with an optimal temporal sampling from displacement observation networks[END_REF], the classical TI does not enable the user to obtain a LF time series using all available displacement observations if the temporal sampling ∆τ is larger than the minimal time interval between two observed dates. For example, in Fig. 1, the minimal time interval between two observed dates is 1 day. For a temporal sampling of ∆τ = 5 days, d t0,t2 , d t2,t4 , d t2,t10 and d t6,t10 (in orange) cannot be written as any combinations of dτ k ,τ k+∆τ .

B. Temporal Inversion using COmbinations of Fractions of displacement (TICOF)

However, it can be established that dt0,t5 + dt5,t10 = d t0,t2 + d t2,t10 . Therefore, [START_REF] Charrier | Extraction of velocity time series with an optimal temporal sampling from displacement observation networks[END_REF] proposed a novel TI formulation, namely Temporal Inversion using COmbinations of displacement observations (TICO), to include combinations of displacement observations. The TICO is formulated as BY = A ′ X where B is a design matrix of dimension n ′ × n and A ′ is the design matrix of dimension n ′ ×p, with n ′ the number of displacement combinations as specified in [START_REF] Charrier | Extraction of velocity time series with an optimal temporal sampling from displacement observation networks[END_REF].

The TICO formulation is adapted when a high amount of observed dates in Y match the dates τ k or τ k+∆τ because it makes possible to link a combination of dτ k ,τ k+∆τ with a combination of d ti,tj . However, when the displacement observations are based on different sensors, the observed dates may not be on the same temporal grid. The dates where acquisitions of both satellites exist depend on the common multiples of the repeat cycles of satellites 1 and 2, here ∆t 1 and ∆t 2 . There are even less common acquisition dates between sensors when missing data exist due to outliers removal.

Therefore, we propose an extension of TICO [START_REF] Charrier | Extraction of velocity time series with an optimal temporal sampling from displacement observation networks[END_REF] for datasets containing multi-sensor displacements and missing data. When a displacement observation can not be included in the TICO formulation, a linear relation between d ti,tj and a fraction of displacements dτ k ,τ k+∆τ is established, assuming linear displacement in the considered temporal interval. For example, in Fig. 1,2 5 dt0,t5 = d t2,t4 . In a more general form, d ti,tj is written as:

d ti,tj = k w k dτ k ,τ k+∆τ (3)
with

w k =                      0, if τ k+∆τ < t i 0, if τ k > t j 1, if τ k ≥ t i ∧ τ k+∆τ ≤ t j τ k+∆τ -ti τ k+∆τ -τ k , if τ k < t i < τ k+∆τ ∧ τ k+∆τ < t j tj -τ k τ k+∆τ -τ k , if τ k < t j < τ k+∆τ ∧ t i < τ k tj -ti τ k+∆τ -τ k , if t i > τ k ∧ t j < τ k+∆τ
The TICOF formulation is: BY = A ′ X where A ′ is filled, in consequence, with 0, 1 and/or fractions (equation S.1). Since the velocity observations have different but not necessarily known uncertainties, the inversion of BY = A ′ X is obtained by an Iterative Reweighted Least Square (IRLS) using the Tukey's biweight function and a Tikhonov regularization [START_REF] Charrier | Extraction of velocity time series with an optimal temporal sampling from displacement observation networks[END_REF]. The weight at the first iteration corresponds to the confidence indicators presented in section III-A.

The use of fractions of displacements is also found in [START_REF] Samsonov | Measuring the state and temporal evolution of glaciers in alaska and yukon using syntheticaperture-radar-derived (sar-derived) 3d time series of glacier surface flow[END_REF] to adjust the first and last date of descending and ascending displacements. This idea assumes a constant velocity over the relevant temporal sampling. For large temporal sampling and/or non-linear displacement behavior, a bias may be introduced in the final estimation of the LF time series. To our knowledge, this bias has never been discussed in the literature.

IV. SYNTHETIC DATASET

To assess the bias introduced by the use of fraction in the case of non-linear displacements (e.g. seasonal variation), synthetic simulations are performed. The synthetic instantaneous velocity is taken as: v(t) = a + b sin( 2π T t) + c cos( 2π T t) with T = 365.25 as in [START_REF] Greene | Detecting seasonal ice dynamics in satellite images[END_REF]. N 1 displacements based on Sentinel-2 repeat cycle and N 2 displacements based on Venµs repeat cycle are simulated. To make sure that the coefficients a, b and c represent well the data, instead of an arbitrary choice, these coefficients are estimated by an IRLS inversion by adding a regularization term corresponding to a displacement model with the coefficients a, b and c as parameters as performed in [START_REF] López-Quiroz | Time series analysis of mexico city subsidence constrained by radar interferometry[END_REF]. The corresponding system of equations is given in equation S.2. The system is solved for Sentinel-2 data on the point represented in blue in Fig S1 . The coefficients are found to be: a = -0.49, b = -0.0788 and c = 0.018.

In Fig 2, LF time series obtained from a moving average of ∆τ days (Fig. 2b), the traditional TI (Fig. 2c) and the TICOF formulation (Fig. 2d) are presented. The moving average corresponds to the average of all the velocities which overlap each period of ∆τ days. The results are represented for a temporal sampling of 20 days. Because velocities are the average of the instantaneous velocity over the temporal baseline, long temporal baseline velocities tend to pull the moving average results toward the annual averaged velocity, here equal to a * 365.25 = -178.85 m/y. In addition, because no displacement observations can be included in the traditional TI when only Venµs data are available, the LF displacements in Fig. 2d) are equal to 0 after October 2018.

Moreover, the Root Mean Square Error (RMSE) between the synthesised and the estimated velocities is more than 80 m/y lower with the TICOF formulation than with the TI or TICO formulation, and about 20 m/y with TICOF than with the moving average, for different temporal sampling (Fig. 3). However, the mean residual of the inversion (the mean difference between AX and Y or A ′ X and BY ) characterizes the bias introduced by the assumption of a constant velocity over the considered temporal sampling. For a synthetic dataset without noise, the residual of the inversion should be null. This is the case for the traditional TI or TICO (Fig. S3). On the contrary, the increase of the TICOF residual which ranges from 0 to 3 m/y, for 0 to 120 days temporal sampling, is caused by the bias introduced by the fractions. This bias is around 25 times lower than the benefit of using TICOF rather than TICO, in case of multi-sensor datasets (Fig. 3). The use of TICOF is advantageous: 1) for multi-sensor datasets and/or 2) datasets with a lot of missing data and/or 3) when the temporal sampling is not a multiple of the satellite repeat cycle (Fig. 3 andS4). In the absence of any of these conditions, Fig. S4 illustrates that the RMSE from TICO is up to 5 m/y better than the one from TICOF. This error also corresponds to the bias introduced by the fractions.

V. REAL DATASET

A. Uncertainty of the results for different temporal samplings

To analyze the uncertainty of the LF time series according to their temporal sampling, the RMSE is computed over the study area delimited by the orange square in Fig. S1, assuming that the true velocity magnitude is null over the stable areas. RMSE is performed over ice-free areas according to the Randolph Glacier Inventory (RGI) V6.0.

Fig. 4 reveals that the RMSE is lower after inversion than before inversion for a temporal sampling ranging from 5 to 20 days (the RMSE is 68% lower for a 5-days temporal sampling to 13 % lower for a 20 days temporal sampling) and remains slightly higher with larger temporal sampling (around 12 m/y higher). According to the RMSE, the performance of the TICOF inversion is degraded with a temporal sampling larger than 20 days. It may be caused by the validity of the assumption beyond the use of fractions, which increases with the temporal sampling. Hence, the use of TICOF inversion appears beneficial for temporal samplings, here, lower than 20 days. From 5 to 20 days, the lowest RMSE is reached for a temporal sampling of 20 days, hence 20 days is selected to match a compromise between uncertainty and temporal resolution. The velocities are about 0 using the TI formulation after October 2018 because the velocity observation based on Venµs cannot be included inside the system. Fig. 3. RMSE between the synthetic multi-sensor velocities and the estimated velocities over different temporal baselines using the traditional TI formulation (in red), the TICOF formulation (in blue), the TICO formulation (in green) and the moving average (in gold).

B. Spatio-temporal analysis of the results

The spatio-temporal evolution of the velocity magnitude is computed over a longitudinal profile defined in Fig. S1. The results from the TICOF formulation using a 20-days temporal sampling are compared with the results from a 20-days moving average. It can be seen in Fig. 5 that there are no available velocity observations from April 2017 to November 2017 at the center part of the glacier tongue.

After inversion, there are visually fewer outliers (e.g. there is no more velocity with a zero value at 3.8 km). Moreover, even if TICOF is pixel-wise, the LF time series seems to be spatially coherent: the velocity field is smooth through the period. A seasonal trend can also be observed: the velocities increase between the two periods September to November 2017 and April to June 2018 with a median increase of 40% along the longitudinal profile. From April/June 2018 to September/November 2018, the velocities decrease with a median decrease of 30%. By using all the available velocity observations, TICOF recovers a seasonal trend when TI or TICO fails because only few velocity observations could have been included in TI or TICO in this case (Fig. S5).

VI. CONCLUSION

The proposed approach enables the fusion of glacier velocity observations from different temporal baselines and sensors without prior information on the glacier dynamics and the data uncertainties. It relies on a new formulation of the Temporal Inversion (TI) based on the temporal closure of the displacement observations network. This new formulation establishes linear equations between combinations of displacement observations and fractions of LF displacements and is named TICOF (Temporal Inversion using COmbinations and Fractions of displacements). The inversion is solved using an Iterative Reweighted Least Square inversion where the initial weights are based on a combination of two statistical criteria: the Median Angle and the modified z-score. The use of displacement fractions assumes a constant velocity over the considered time interval. The bias caused by the validity of this assumption is evaluated on synthetic data. For a sinusoidal signal with a period of 1 year, the bias introduced is below 0.1 m/y for a temporal sampling of 5 to 20 days. This is compared to the improvement in velocity uncertainty obtained after TICOF. For time sampling less than 20 days, the use of TICOF is still beneficial (the uncertainty is reduced by 13%). The 20-day LF velocity time series makes it possible to highlight the strong seasonality of the Fox glacier with lower uncertainty, fewer outliers, and no more gaps. The fluctuation of Fox surface velocities is on average 35 % between spring and autumn from 2017 to 2018. Since this approach does not require prior information, it can be applied to any kind of available velocity datasets computed by any research team. S4. RMSE between the synthetic mono-sensor velocities and the estimated velocities over different temporal samplings using the traditional TI formulation (in red), the TICOF formulation (in blue), the TICO formulation (in green) and the moving average (in gold). The vertical grey lines correspond to the temporal samplings which are a multiple of the satellite repeat cycle, here 2 days. For these temporal sampling, TICO is more advantageous than TICOF. For other temporal sampling, TICOF performs better. Note that, for some temporal sampling, it is not possible to include any observations inside the TI system, therefore there is no corresponding value of RMSE.

The system of equations corresponding to the TICOF is: where λ corresponds to the damping factor of the regularization. ∆τ is the temporal sampling of the LF time series [ dτ 1 ,τ 1+∆τ , ..., dτ 1+(n-1)∆τ ,τ 1+n∆τ ].

     1 0 0 . . . 0 0 0 1 1 . . . 0 0 . . . . . . . . . . . . . . . 0 0 0 . . . 1 0 0 0 0 . . . 0 1               dτ 1 ,τ 1+2∆τ dτ 1 ,τ 1+ 3 2 ∆τ dτ 1+ 3 2 ∆τ ,τ 1+3∆τ . . . dτ 1+(n-2)∆τ ,τ 1+(n-1 2 )∆τ dτ 1+(n-1)∆τ ,τ 1+n∆τ          =                1 1 0 . . . 0
[dτ 1 ,τ 1+2∆τ , ..., dτ 1+(n-1)∆τ ,τ 1+n∆τ ] stands for the displacement observations.

The system of equations including a model based regularization used to determine simulation parameters in section IV is:

        1

Fig. 1 .

 1 Fig. 1. Illustration of a) the traditional TI formulation b) the proposed TICOF formulation. The displacement observations from satellite 1 is in blue while the displacement observations from satellite 2 is in orange.

Fig. 2 .

 2 Fig. 2. a) Synthetic velocities spanning a given temporal baseline are represented by blue and orange bars, for velocities based on Sentinel-2 and Venµs images respectively. These velocities correspond to the average of the instantaneous velocity (represented by a black curve) over the considered temporal baseline. Hence, small temporal baseline velocities are closer to the instantaneous velocity than long temporal baseline velocities. The dots represent the central date of each velocity. RLF time series obtained using the b) moving average c) traditional TI formulation d) TICOF formulation, with a 20-day temporal sampling.The velocities are about 0 using the TI formulation after October 2018 because the velocity observation based on Venµs cannot be included inside the system.

Fig. 4 .

 4 Fig. 4. RMSE over stable ground for different temporal samplings before (in gold) and after TICOF inversion (in blue).

Fig. 5 .

 5 Fig. 5. Spatio-temporal evolution of velocity magnitude along a longitudinal profile defined in Fig. S1. The left-hand figure correspond the results from a 20-days moving average, the results on the right corresponds to the results of the TICOF inversion.

Fig. S1 .

 S1 Fig. S1. Fox glacier in the Southern Alps of New Zealand. The red dotes are spaced 1 km apart. The background image is an optical image from the Google Earth collection.

Fig. S2 .

 S2 Fig. S2. Temporal baselines of the velocity observations according to their central date for Sentinel-2 velocities in blue and Venµs velocities in orange. The colormap which ranges from light to dark blue or orange indicates the percentage of non-masked data for each temporal baseline and central date.

Fig. S3 .

 S3 Fig. S3. Residual of the inversion on synthetic data with the traditional TI formulation (in red) and the TICOF formulation (in blue).

Fig.

  Fig.S4. RMSE between the synthetic mono-sensor velocities and the estimated velocities over different temporal samplings using the traditional TI formulation (in red), the TICOF formulation (in blue), the TICO formulation (in green) and the moving average (in gold). The vertical grey lines correspond to the temporal samplings which are a multiple of the satellite repeat cycle, here 2 days. For these temporal sampling, TICO is more advantageous than TICOF. For other temporal sampling, TICOF performs better. Note that, for some temporal sampling, it is not possible to include any observations inside the TI system, therefore there is no corresponding value of RMSE.

  . . . 0 -[τ1+∆τ -τ1] -T 2π [sin( 2π T τ1+∆τ ) -sin( 2π T τ1)] -T 2π [cos( 2π T τ1+∆τ ) -cos( 2π T τ1)] . . . λ -[τ1+n∆τ -τ 1+(n-1)∆τ ] -T 2π [sin( 2π T τ1+n∆τ ) -sin( 2π T τ 1+(n-1)∆τ )] -T 2π [cos( 2π T τ1+n∆τ ) -cos( 2π T τ 1+(n-1)∆τ )] 1where T is the period of the sinusoidal signal. a, b, c, d are the coefficients of the model.

Fig. S5 .

 S5 Fig.S5. Spatio-temporal evolution of velocity magnitude along a longitudinal profile defined in Fig.S1. The results are from the TI (left) and TICO (right) inversion, with a temporal sampling of 20 days. Because the velocity observations are multi-sensor, with missing data due to outlier removal, all the velocity observations cannot be included inside the system. Therefore, the estimations are either not possible (from 2.4 to 5.1 km with TI) or mainly constrained by the Tikhonov regularization. This results in very smooth variations, incoherent with the observations.

Acknowledgment: We acknowledge CNES (French Spatial Agency) and ISA (Israel Space Agency) for the VENµs observations acquired through the VENµs Scientific Mission, and ESA (European Space Agency) for the Sentinel-2 observations. JM and RM acknowledge support from the CNES.


Chupin Olivier
email: olivier.chupin@ifsttar.fr

Chabot Armelle

Bodin Didier
email: didier.bodin@arrb.com.au

Piau Jean-Michel

-structure-application-case-thick-flexible-pavement-olivierarmelle

Keywords: thick flexible pavement, viscoelasticity, equivalent elastic system

Bituminous pavements are known to exhibit a thermo-viscoelastic behavior which is strongly dependent on temperature and load velocity. However, for many applications (pavement design, FWG analysis…) it is more convenient to deal with elastic calculations which facilitate parametric studies. But to be accurate such an approximation requires rules defining the right choice of elastic data sets. In this context, this paper presents a method to derive an equivalent elastic system to a multilayer viscoelastic pavement under given conditions of temperature and subject to loads moving at constant speed. Thus, the tool proposed hereafter should ease the elaboration of such rules. This paper explains the proposed method and illustrates its use on the case of a thick flexible pavement. By the way, we reexamine the "10Hz rule" of the French pavement design method within this framework. We recall that this rule assumes that the equivalent modulus of a bituminous layer is approximately equal to the norm of its complex modulus computed at 10Hz and at the temperature of the layer.

INTRODUCTION

It is well known that the resilient behavior of asphalt concrete under traffic load is viscoelastic and highly sensitive to temperature, these features being inherited from those of bitumen binders. To take account of these features numerical programs such as ViscoRoute© 2.0 were developed more or less recently [START_REF] Hopman | A Viscoelastic Multilayer Computer Program[END_REF][START_REF] Siddharthan | Pavement strain from moving dynamic 3D load distribution[END_REF][START_REF] Duhamel | ViscoRoute: visco-elastic modeling for asphalt pavements[END_REF][START_REF] Chabot | Viscoroute 2.0: a tool for the simulation of moving load effects on asphalt pavement[END_REF][START_REF] Chupin | Influence of sliding interfaces on the response of a layered viscoelastic medium under a moving load[END_REF]. The high accuracy potential of these tools was demonstrated by way of comparison with experimental measurements stemming from instrumented pavements under simulated traffic, when varying amplitude and velocity of the loads, the pavement geometry or the profiles of temperature. Nevertheless for many applications such as pavement design or survey and back analysis, the use of elastic-based programs still remains at first step a very easy, fast and useful mean of modeling, which is also capable of providing an accurate quantitative description of pavement behavior. However the main condition for this is to make the right choice of the "equivalent elastic moduli" to be affected to the viscoelastic layers of the structure. In this context, the present paper aims at developing a methodology based on viscoelastic simulations and back calculation for determining these moduli. In addition, the results obtained from the method developed are eventually compared to the choice of modulus recommended by the French design method [START_REF] Corte | Design of pavement structures: the French technical guide[END_REF].

CALCULATION METHOD

The calculation method relies on the use of software ViscoRoute© 2.0 dedicated to the computation of the dynamic response (3D) under moving loads of pavements incorporating viscoelastic material layers. The software program is based on a semi-analytical solution technique described in other references to which the reader is referred for a comprehensive description [START_REF] Duhamel | ViscoRoute: visco-elastic modeling for asphalt pavements[END_REF][START_REF] Chabot | Viscoroute 2.0: a tool for the simulation of moving load effects on asphalt pavement[END_REF][START_REF] Chupin | Influence of sliding interfaces on the response of a layered viscoelastic medium under a moving load[END_REF]. The viscoelatic model is considered as the reference to simulate the pavement response at a given temperature and loading speed. The same tool is then used to calculate the response of the associated elastic structure of same geometry and subject to similar loading. The Young moduli of the different layers will be determined in order to obtain, in "some sense", comparable responses for the two structures. For this purpose, we define a mathematical function (E) to minimize which expresses the "distance" between the two models:

( ) ( ) ( ) ( ) 2 ,max ,max 1 ,max , , , , , ve e n
i i i i v e i i i c V z c z c V z    =   - =        E E (1) max v e
c and max e c denote the selected quantities (e.g. deflection, strain, etc.) that we wish to be close in the viscoelastic and the elastic cases, respectively. The "equivalence criterion"  is defined using n quantities which may be calculated at different depths z in the structure. max v e c represents the reference computed using ViscoRoute© 2.0 in which asphalt layers are modeled through the Huet-Sayegh model [START_REF] Huet | Etude par une méthode d'impédance du comportement viscoélastique des matériaux hydrocarbonés[END_REF][START_REF] Huet | Coupled size and boundary-condition effects in viscoelastic heterogeneous and composite bodies[END_REF][START_REF] Sayegh | Contribution à l'étude des propriétés viscoélastiques des bitumes purs et des bétons bitumineux[END_REF] (Figure 1). max v e c is thus a function of speed V and temperature  so the minimization technique is run for given values of these parameters. max e c depends on E which is an array of size nve (number of viscoelastic layers) composed of the unknown elastic moduli.

E 0 E  -E 0 ( , t k ) (t h ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )
* 0 0 2 0 1 2 1 exp k h E E E E j j A A A             - - - = + + + = + + FIGURE 1
Schematic of the Huet-Sayegh rheological model and its complex modulus as a function of pulsation and temperature (h>k>0; j 2 =-1)

The minimization is performed according to a gradient method with respect to E for which an increment of modulus, dE, is computed at each iteration as follows:

( ) ( ) ( )

2 with       =    d = - E E E E E E
(2) The minimization process is stopped when

( )     E
.  and  are small positive quantities.

For pavement design applications, the important issue is to obtain comparable true and equivalent critical strains or stresses governing the material damage criteria. Thus here we define  as a function of the maximal strain and stress at the bottom of the asphalt layers, considered as the most relevant quantities for asphalt fatigue prediction. Note that this function can be extended to other quantities. For instance, if the application is focused on deflections or curvature then the whole deflection bowl could be used in the  function.

STUDY OF A THICK FLEXIBLE PAVEMENT

The calculation method is applied to the study of a typical thick flexible pavement modeled with either two or three viscoelastic layers. In the first case, only two unknown moduli must be determined. This case will be used to illustrate the shape of the  function and evaluate the calculation method. On the other hand, the case with three viscoelastic layers aims at illustrating a complicated situation involving thermal gradient in the structure. For both cases, we also examine the validity of the "10 Hz rule" used in the French pavement design method [START_REF] Corte | Design of pavement structures: the French technical guide[END_REF] to assess the elastic stiffness of bituminous layers (see further).

FIGURE 2 Cross-section of the pavement structure considered A schematic representation of the pavement structure is shown in Figure 2. It is composed of viscoelastic layers supported by a semi-infinite elastic foundation. For the analysis of the two-viscoelastic-layer case we consider that the two GB3 layers are melded. In this section, the pavement is subjected to the French standard dual-wheel moving at constant speed V=72km/h=20m/s. The weight of each wheel (32,5kN) is assumed to be uniformly distributed at the surface of the pavement over a disk of radius 0.125m (Figure 3).

Choice of  and Validation of the Calculation Method

In the present case, function  depends on two elastic moduli and its minimum can be graphically determined in a simple way. To validate the calculation method these moduli can be compared to those obtained by minimization. To evaluate the impact of the selected function  on the back-calculated moduli, two functions composed of three quantities (n=3) are used below.

The first function 1 includes T and T computed at the bottom (z=0.33m) of the pavement as well as T computed at the top (z=0.04m) (see Figure 2). T denotes the maximum between the absolute values of the longitudinal strains xx and yy; T is the equivalent for stress.

In the second function denoted 2, T in z=0.04m is replaced by zz calculated at the top of the soil foundation (z=0.35m). Figure 4 shows a plot of the logarithm (log10) of these two functions.

In this figure, E1 and E2 are respectively the values of the Young modulus of layer 1 and 2 attributed for elastic structure. Besides the viscoelastic computation is performed for V=20m/s and =30°C (temperature is uniform in both layers).

These two functions exhibit a single minimum which is more apparent for 1 in contrast with 2. Consequently, function 1 is probably more conducive to determining numerically the equivalent moduli by means of gradient descent. Note that depending on its definition, function  may have several local minima and the choice of the quantities included in  should be made with caution. For instance, 2 has not a well defined minimum probably because the value of zz in z=0.35m is more sensitive to the mechanical properties of the soil foundation than to E1 and E2. In the rest of the paper, we use function 1 (=1). Under the present conditions, this function leads to values of the equivalent moduli obtained by minimization of 1625 and 2712 MPa for E1 and E2, respectively. These values match those of the minimum shown in Figure 4 thus validating the calculation method.

Detailed Comparison Between the Viscoelastic and Elastic Equivalent Solutions (=30°C)

The viscoelastic response of the pavement is compared to that obtained in the elastic case using the equivalent moduli computed before (i.e. for V=20m/s and =30°C). For that purpose, Figure 5 shows horizontal profiles in the x-or y-direction (at a given depth z) of the mechanical fields used to define the equivalence criterion, i.e. to calculate function . We observe that these profiles provide as expected a good estimate of the maximal values T and T. Nonetheless, the shapes of the profiles are somehow different since creep/relaxation inherent to viscoelasticity leads to non-symmetric profiles in the moving direction and to a shifting of the maximum positions. For these fields, in addition to some differences in shape between the elastic and the viscoelastic response non negligible differences are also observed at peak values. This could be expected since no constraint is imposed to match these components of the viscoelastic response.

Nevertheless, the maximum of these fields remain lower than T and T and in the present example this difference would not impact design based on maximum strain or stress. Globally, the elastic calculation leads to quite acceptable results.

Equivalent Moduli as a Function of Temperature (V=20m/s) -The "10 Hz Rule"

The variation of the equivalent moduli with respect to temperature is now studied. The minimization is run for different temperatures which are assumed to be homogeneous with the same value in both viscoelastic layers. The temperature is varied successively from 0 to 40°C. Figure 7 shows the results obtained for the BB and the GB3 layers. As expected the equivalent moduli decrease with temperature (V=20m/s), i.e. as the viscoelastic materials soften. At =30°C the equivalent moduli are those already computed in the previous section (E1=1625MPa and E2=2712MPa). For the different temperatures, the equivalent moduli are also compared to the magnitude of the complex modulus of the corresponding layer computed at temperature  and frequency f=/2=10Hz, which is the rule applied in the design of pavements in France to account for variations of temperature. The values of the complex modulus for each layer are plotted using circle markers in Figure 7. These values are in good agreement with the equivalent back-calculated moduli indicating that the rule under consideration performs well in the present case. We consider in this section a three-viscoelastic-layer structure subjected to a tridem loading configuration. A thermal gradient is imposed in the asphalt concrete part of the structure. The temperature in the viscoelastic materials is assumed uniform per layer and equal to 35, 25 and 15°C for layer 1 (BB), 2 (GB3) and 3 (GB3), respectively. The same function  as previously is used in the minimization. Figure 8 shows profiles of the selected fields used to compute function . The comparison between the viscoelastic and the elastic equivalent 1029 4461 9375 responses leads to similar conclusions as evoked for the two-viscoelastic-layer structure, i.e. a good match of the maximum values with some differences in the shape of the profiles. However, the back calculated moduli obtained and summarized in Table 2 are much lower than the current design moduli stemming from the "10Hz rule". In terms of mechanical response shown in Figure 8, this yields underestimation of T in z=0.04m and 0.33m, and overestimation of T in z=0.33m when considering the design rule. The application of this rule is thus not well adapted to the present case whereas the elastic calculation still performs well provided that the equivalent elastic moduli are correctly chosen. The differences observed could be attributed to a change of the appropriate loading frequency (load pulse duration) to be selected at different depth [START_REF] Denneman | Improved Design Procedures for Asphalt Pavements: Pavement 1 Temperature and Load Frequency Estimation[END_REF] for the thick asphalt pavements. As a remark, the loading frequencies corresponding to the back calculated moduli (inversion of the Huet-sayegh law) are 7.6, 2.0 and 1.5 Hz for the layers 1, 2 and 3, respectively.

CONCLUSION

The aim of this paper was to present a method intended to help the setup of rules for defining equivalence between viscoelastic and elastic calculations of pavements, depending on the context of use and of a desired accuracy. In this framework, the proposed method was used to compute the equivalent moduli to be affected to bituminous layers of a thick flexible pavement subject to moving loads. This method has proven to provide accurate estimate of the equivalent moduli for various conditions of loading and temperature. The back calculated moduli obtained using the developed tool have also been compared to those recommended by the French pavement design method. Depending on the situation under consideration, the difference between these moduli was found more or less important. Nevertheless, to conclude on this aspect it would be necessary to define an acceptable difference after integrating weighting of the results to account for traffic, weather conditions, etc. which are elements considered in the design method. Finally, the type of approach presented here could be used to extrapolate the pavement design method to other contexts as done already for airfield pavements in France for which the equivalent moduli are determined as a function of temperature, load speed and thickness of the layers.

4 1 .FIGURE 3

 13 FIGURE 3 Top views of the loading configurations applied to the pavement 3.1 Case of a Two-Viscoelastic-Layer Structure

FIGURE 4

 4 FIGURE 4 Plot of function  (E1,E2) and  (E1,E2) for the two-viscoelastic-layer pavement

FIGURE 5 FIGURE 6 and viscoelastic computations 5 6 7

 566 FIGURE 5 Profiles in the x-or y-direction of the fields used in the minimization process: 1 comparison between elastic and viscoelastic computations 2 3

FIGURE 7 3 . 2

 732 FIGURE 7 Variation of the equivalent moduli with respect to temperature and comparison with the norm of the complex modulus computed for  and 10Hz 3.2 Case of a Three-Viscoelastic-Layer Structure

FIGURE 8

 8 FIGURE 8 Profiles versus x of the fields used in the minimization process of the threeviscoelastic-layer structure: comparison between viscoelastic, elastic equivalent, and "10 Hz rule" computations

Table 1

 1 gives the values of the Huet-Sayegh parameters of the different viscoelastic layers.

	BB: viscoelastic	0.08m
	GB3: viscoelastic	0.13m
	GB3: viscoelastic	0.13m
	Soil : Esoil = 50MPa, =0.35	
	z

TABLE 1 Values of the Huet-Sayegh parameters

 1

		E0 (MPa)	E (MPa)		k	h	A0	A1	A2
	GB3	11	28000	2.0	0.18	0.5	2.8	-0.4	1.7E-3
	BB	11	18000	2.0	0.18	0.45	2.3	-0.

TABLE 2 Back calculated moduli for the different layers and design moduli from the "10Hz rule"

 2

		Layer 1 -BB (35°C) Layer 2 -GB3 (25°C) Layer 3 -GB3 (15°C)
	Eeq (Mpa)	929	2695	6379
	|E(theta,10Hz)| (MPa)