Charlotte Bouquiaux

Claire Tonnelé

Frédéric Castet

Benoît Champagne
email: benoit.champagne@unamur.be

Second-Order Nonlinear Optical Properties of an Amphiphilic Dye Embedded in a Lipid Bilayer. A Combined Molecular Dynamics -Quantum Chemistry Study

The second harmonic signal of an amphiphilic dye embedded in a lipid bilayer has been calculated by combining molecular dynamics simulations and quantum chemistry calculations based on density functional theory. This computational approach provides insight into the morphology and dynamics of the fully hydrated biological system, as well on the relationships linking the geometry and the environment of the dye to the amplitude of its second-order nonlinear optical response. The results point out a significant enhancement of the dynamic first hyperpolarizability of the dye induced by its interaction with the membrane, and highlight the relative importance of dynamical, steric, and electrostatic effects. This computational scheme is thus particularly relevant for rationalizing the nonlinear optical contrasts revealed by second harmonic imaging microscopy of exogenous dyes embedded in biological media.

Introduction

The cell membrane is at the heart of many biological processes: without membranes, life as we know it would not be possible. A membrane is a selective semi-permeable barrier separating cellular compartments that acts as a transporter. [START_REF] Goodman | Medical Cell Biology[END_REF] Furthermore, transmembrane ionic gradients (mainly involving Na + , K + and Cl -ions) create a membrane potential that in turn provides an electrical driving force for many biological mechanisms. [START_REF] Goodman | Medical Cell Biology[END_REF][START_REF] Friedman | Membrane-Ion Interactions[END_REF] The fundamental functions of membranes are dictated by their compositions: 1 lipids represent approximatively 50% in mass of most animal biomembranes, the rest being mostly proteins. [START_REF] Alberts | The Lipid Bilayer from Molecular Biology of the Cell[END_REF] Amongst lipids, there is an enormous diversity of structures owing to the variety of fatty acids, their linkage positions, and their headgroups. [START_REF] Harayama | Understanding the Diversity of Membrane Lipid Composition[END_REF] Interestingly, changes in lipid composition have been reported in numerous diseases. [START_REF] Goodman | Medical Cell Biology[END_REF][START_REF] Harayama | Understanding the Diversity of Membrane Lipid Composition[END_REF] The structure of lipid bilayers is described as fluid and flexible due to the possible motions of lipid molecules. These motions include rotations around the lipid long axes, "wiggles" of the fatty acids chains, translations within one layer, and motions to the other layer. [START_REF] Goodman | Medical Cell Biology[END_REF][START_REF] Alberts | The Lipid Bilayer from Molecular Biology of the Cell[END_REF] As a matter of fact, the investigation of the structure and dynamics of a membrane is of particular importance. However, most biomolecules possess only few natural moieties with exploitable optical properties [START_REF] Thomas | Optical Imaging Probes for Biomolecules: An Introductory Perspective[END_REF] -with the exception of structural proteins such as collagen [START_REF] Freund | Optical Second-Harmonic Microscopy, Crossedbeam Summation, and Small-Angle Scattering in Rat-tail Tendon[END_REF][START_REF] De Wergifosse | How the Second-Order Nonlinear Optical Response of the Collagen Triple Helix Appears: A Theoretical Investigation[END_REF] and cellulose [START_REF] Brown | Macromolecular Structure of Cellulose Studied by Second-Harmonic Generation Imaging Microscopy[END_REF] -so that exogenous dyes are often employed to improve the contrast in tissues for being detected by commercially available microscopes. Among the linear and nonlinear optical phenomena operating in microscopy, fluorescence and Second Harmonic Generation (SHG) are two main techniques whose combined use provides complementary information. [START_REF] Moreaux | Membrane Imaging by Simultaneous Second-Harmonic Generation and Two-Photon Microscopy[END_REF][START_REF] Zipfel | Live Tissue Intrinsic Emission Microscopy using Multiphoton-Excited Native Fluorescence and Second Harmonic Generation[END_REF][START_REF] Millard | Second Harmonic Imaging of Exocytosis at Fertilization[END_REF] These optical processes differ by the nature of the initial and final molecular electronic states, and depend on the illumination wavelengths and intensities. AminoN-aphtylEthenylPyridinium (ANEP) dyes, like di-8-ANEPPS (Figure 1) developped by Loew and colleagues, [START_REF] Fluhler | Membrane Binding, and Potentiometric Responses of New Charge Shift Probes[END_REF][START_REF] Loew | Potentiometric Dyes: Imaging Electrical Activity of Cell Membranes[END_REF] are often used in this context as labelling agents. These push-pull compounds are not only fluorescent probes, but also exhibit SHG responses. SHG is a second-order nonlinear optical (NLO) process, first demonstrated by Franken et al. [START_REF] Franken | Generation of Optical Harmonics[END_REF] in crystalline quartz in 1961. It was then developed by Shen, [START_REF] Shen | Surface Properties Probed by Second-Harmonic and Sum-Frequency Generation[END_REF] Eisenthal, [START_REF] Eisenthal | Liquid Interfaces Probed by Second-Harmonic and Sum-Frequency Spectroscopy[END_REF] and others for the study of surfaces and artificial monolayers. The use of SHG in microscopy was proposed by Sheppard et al. [START_REF] Sheppard | The Scanning Optical Microscope[END_REF] in 1977. Since then, there has been an expansion in the use of Second Harmonic Imaging Microscopy (SHIM), partially thanks to advances in laser technology. [START_REF] Campagnola | High-Resolution Nonlinear Optical Imaging of Live Cells by Second Harmonic Generation[END_REF] SHG presents several advantages over fluorescence: 19 i) it is intrinsically specific to interfacial regions, making it an ideal technique to study the biophysics of model membranes. Indeed, the observed signal is only due to non-centrosymmetrical molecules and molecular assemblies, whereas in fluorescence, a background is usually observed arising also from isotropic regions; ii) the wavelength of 1064 nm commonly used in SHG measurements should prevent photodamage and photobleaching; iii) the SHG signal retains the phase and directional data after interaction with the sample. Information about the molecular organization of the chromophore can therefore be extracted from SHG imaging data, by taking advantage of the SHG polarization anisotropy; 20 and iv) finally, the SHG intensity can be modulated by an applied electric field, [START_REF] Bouevitch | Probing Membrane Potential with Nonlinear Optics[END_REF] and therefore be used as an efficient probe of the membrane potential. SHG responses of biological systems have been widely studied throughout the literature. [START_REF] Brown | Macromolecular Structure of Cellulose Studied by Second-Harmonic Generation Imaging Microscopy[END_REF][START_REF] Campagnola | High-Resolution Nonlinear Optical Imaging of Live Cells by Second Harmonic Generation[END_REF][START_REF] Campagnola | Three-Dimensional High-Resolution Second-Harmonic Generation Imaging of Endogenous Structural Proteins in Biological Tissues[END_REF][START_REF] Moreaux | Membrane Imaging by Second-Harmonic Generation Microscopy[END_REF][START_REF] Moreaux | Coherent Scattering in Multi-Harmonic Light Microscopy[END_REF][START_REF] Campagnola | Second-Harmonic Imaging Microscopy for Vizualizing Biomolecular Arrays in Cells, Tissues and Organisms[END_REF][START_REF] Mohler | Second Harmonic Generation Imaging of Endogenous Structural Proteins[END_REF][START_REF] Williams | Interpreting Second-Harmonic Generation Images of Collagen I Fibrils[END_REF][START_REF] Millard | Second Harmonic Imaging of Exocytosis at Fertilization[END_REF][START_REF] Asselberghs | Second-Harmonic Generation in GFP-like Proteins[END_REF][START_REF] Deniset-Besseau | Measurement of the Second-Order Hyperpolarizability of the Collagen Triple Helix and Determination of Its Physical Origin[END_REF][START_REF] De Meulenaere | Improving the Second-Order Nonlinear Optical Response of Fluorescent Proteins: The Symmetry Argument[END_REF][START_REF] Timr | Nonlinear Optical Properties of Fluorescent Dyes Allow for Accurate Determination of Their Molecular Orientations in Phospholipid Membranes[END_REF][START_REF] Okur | Chemistry of Lipid Membranes from Models to Living Systems: A Perspective of Hydration, Surface Potential, Curvature, Confinement and Heterogeneity[END_REF] In particular, a lot of works focused on the study of the transmembrane potential. [START_REF] Bouevitch | Probing Membrane Potential with Nonlinear Optics[END_REF][START_REF] Peleg | Nonlinear Optical Measurement of Membrane Potential Around Signle Molecules at Selected Cellular Sites[END_REF][START_REF] Millard | Sensitivity of Second Harmonic Generation from Styryl Dyes to Transmembrane Potential[END_REF][START_REF] Dombeck | Optical Recording of Fast Neuronal Membrane Potential Transients in Acute Mammalian Brain Slices by Second-Harmonic Generation Microscopy[END_REF][START_REF] Nuriya | Imaging Membrane Potential in Dentritic Spines[END_REF][START_REF] Nuriya | Membrane Potential Dynamics of Axons in Cultured Hippocampal Neurons by Second-Harmonic-Generation Imaging[END_REF] However, although several computational approaches allowing to predict and interpret macroscopic SHG responses of solutions, [START_REF] Benassi | General Strategy for Computing Nonlinear Optical Properties of Large Neutral and Cationic Organic Chromophores in Solution[END_REF][START_REF] Pielak | Dynamical Behavior and Second Harmonic Generation Responses in Acido-Triggered Molecular Switches[END_REF] periodic solids, [START_REF] Reis | Calculation of Macroscopic, First-, Second-, and Third-Order Optical Susceptibilities for the Urea Crystal[END_REF][START_REF] Kirtman | Major Intermolecular Effects on Nonlinear Electrical Response in a Hexatriene Model of Solid State Polyacetylene[END_REF][START_REF] Champagne | Advances in Chemical Physics[END_REF][START_REF] Seidler | Investigation of the Linear and Second-Order Nonlinear Optical Properties of Molecular Crystals Within the Local Field Theory[END_REF][START_REF] Rérat | Computation of Second Harmonic Generation for Crystalline Urea and KDP. An ab Initio Approach through the Coupled Perturbed Hartree-Fock/Kohn-Sham Scheme[END_REF][START_REF] Seidler | QTAIM-Based Scheme for Describing the Linear and Nonlinear Optical Susceptibilities of Molecular Crystals Composed of Molecules with Complex Shapes[END_REF] , interfaces, [START_REF] Licari | Exploring Fluorescent Dyes at Biomimetic Interfaces with Second Harmonic Generation and Molecular Dynamics[END_REF] and self-assembled monolayers [START_REF] Tonnelé | Nonlinear optical responses of self-assembled monolayers functionalized with indolino-oxazolidine photoswitches[END_REF][START_REF] Tonnelé | Nonlinear Optical Contrast in Azobenzene-Based Self-Assembled Monolayers[END_REF] have so far been developed, much less has been done for complex dynamical systems like dyed lipid bilayers. Few theoretical chemistry studies have also been performed to assess their SHG responses but either they truncate the chromophore complex surrounding, they do not account for geometry relaxation effects due to this surrounding, they rely on low-level methods to evaluate the β tensor components, or even they do not account for the dynamics of the system. [START_REF] De Meulenaere | Improving the Second-Order Nonlinear Optical Response of Fluorescent Proteins: The Symmetry Argument[END_REF][START_REF] Rusu | Monitoring Biological Menbrane-Potential Changes: A CI QM/MM Study[END_REF][START_REF] Murugan | Association Dynamics and Linear and Nonlinear Optical Properties of an N-Acetylaladanamide Probe in a POPC Membrane[END_REF][START_REF] Osella | Investigation into Biological Environments through (Non)linear Optics: A Multiscale Study of Laurdan Derivatives[END_REF][START_REF] De Wergifosse | ONIOM Investigation of the Second-Order Nonlinear Optical Responses of Fluorescent Proteins[END_REF] On this basis, our global objective is to better understand the SHG responses of chromophores in lipid bilayers and to develop appropriate computational methods to unravel SHIM signatures. We focus in this report on the elaboration of a two-step simulation procedure and its application for a model system composed of a di-8-ANEPPS probe inserted into a hydrated 1,2-DiPalmitoylPhosphatidylcholine (DPPC) membrane. Molecular Dynamics (MD) simulations are first performed to investigate the morphology of the ANEPPS@DPPC system and the dynamics of the chromophore inserted into the membrane. Then, the first hyperpolarizability (β) responses of the di-8-ANEPPS probe are evaluated using density functional theory (DFT) for selected snapshots extracted from the MD simulations. This computational approach provides insights into the relationships linking the geometry of the chromophore to the amplitude of its SHG responses and the related visible absorption spectra, into the impact of structural dynamics, as well as into the effects of the lipid bilayer environment.

Computational and Theoretical Methods

Molecular Dynamics simulations

MD simulations of the ANEPPS@DPPC system were carried out using NAMD2.12, 53 while VMD1.9.4 [START_REF] Humphrey | Visual Molecular Dynamics[END_REF] was used for data visualization and analysis. The system consists of one di-8-ANEPPS dye molecule inserted into a lipid bilayer formed by 125 DPPC molecules (61 and 64 in the leaflet with and without the chromophore, respectively) and 3840 water molecules (see After the starting structure was generated using Packmol, [START_REF] Martinez | A Package for Building Initial Configurations for Molecular Dynamics Simulations[END_REF] the system was modeled following a procedure similar to those commonly used in the literature, and described hereafter. [START_REF] Reeve | Dyes for Biological Second Harmonic Generation Imaging[END_REF][START_REF] Rusu | Monitoring Biological Menbrane-Potential Changes: A CI QM/MM Study[END_REF][START_REF] Van Der Ploeg | Molecular Dynamics Simulation of a Bilayer Membrane[END_REF][START_REF] Feller | Molecular Dynamics Simulation of Unsaturated Lipid Bilayers at Low Hydration: Parameterization and Comparison with Diffraction Studies[END_REF][START_REF] Marrink | Adhesion Forces of Lipids in a Phospholipid Membrane Studied by Molecular Dynamics Simulations[END_REF][START_REF] Forest | Membrane Simulations: Bigger and Better?[END_REF][START_REF] Mashl | Molecular Simulation of Dioleoylphosphatidylcholine Lipid Bilayers at Differing Levels of Hydration[END_REF][START_REF] Mukhopadhyay | Molecular Dynamics Simulation of a Palmitoyl-Oleoyl Phosphatidylserine Bilayer with Na + Counterions and NaCl[END_REF][START_REF] Sankararamakrishnan | Surface Tension Parametrization in Molecular Dynamics Simulations of a Phosholipid-Bilayer Membrane: Calibration and Effects[END_REF][START_REF] Benz | Experimental Validation of Molecular Dynamics Simulations of Lipid Bilayers: A New Approach[END_REF][START_REF] Barzoukas | A Two-State Description of (hyper) Polarizabilities of Push-Pull Molecules Based on a Two-Form Model[END_REF][START_REF] Klauda | Update of the CHARMM All-Atom Additive Force Field for Lipids: Validation on Six Lipid Types[END_REF][START_REF] Robinson | Di-8-ANEPPS Emission Spectra in Phospholipid/Cholesterol Membranes: A Theoretical Study[END_REF][START_REF] Jämbeck | Derivation and Systematic Validation of a Refined All-Atom Force Field for Phosphatidylcholine lipids[END_REF][START_REF] Jurkiewicz | Structure, Dynamics, and Hydration of POPC/POPS Bilayers Suspended in NaCl, KCl, and CsCl Solutions[END_REF][START_REF] Mori | Molecular Dyanmics Simulations of Biological Membranes and Membranes Proteins Using Enhanced Conformational Sampling Algorithms[END_REF][START_REF] Lyubartsev | Force Field Development for Lipid Membrane Simulations[END_REF][START_REF] Macchiagodena | Computational Study of DPAP Molecular Rotor in Various Environments: from Force Field Development to Molecular Dynamics Simulations and Spectroscopic Calculations[END_REF] The partial charges were calculated for the chromophore, the water, and the lipid molecules using the electrostatic potential (ESP) model [START_REF] Momany | Determination of Partial Atomic Charges from Ab Initio Molecular Electrostatic Potentials. Application to Formamide, Methanol, and Formic Acid[END_REF][START_REF] Singh | An Approach to Computing Electrostatic Charges for Molecules[END_REF] at the DFT/M06-2X level with the cc-pVDZ basis set, as implemented in the Gaussian16 package. [START_REF] Frisch | Gaussian 16 Revision A.03[END_REF] Periodic boundary conditions were applied in all three directions, so that the simulation is actually that of a multilamellar system. After an initial minimization of a low-density sample (with a starting box side of 80 Å), the ANEPPS@DPPC system was equilibrated for 140 ns in the NPT ensemble (P = 1 atm and T = 315.15 K 57 ). This temperature was selected to ensure that the bilayer adopts the liquid crystalline (L α ) phase, the L β to L α phase transition temperature of DPPC being 315 K (see Section 1.3 in Supporting Information (SI)). A 240 ns production with a time step of 1 fs was then performed under the same conditions. The pressure and temperature were maintained using the Langevin piston [START_REF] Feller | Constant Pressure Molecular Dynamics Simulation: The Langevin Piston Method[END_REF] and Langevin temperature bath, respectively. Lennard-Jones interactions were truncated at 12 Å using a switching function starting at 10 Å, while long-range Coulombic interactions were treated using the particle-mesh Ewald (PME) technique. [START_REF] Piana | Evaluating the Effects of Cutoffs and Treatment of Long-range Electrostatics in Protein Folding Simulations[END_REF] Given the high sensitivity of the β responses to the molecular geometry, a modified version of the General AMBER Force Field (GAFF) [START_REF] Wang | Development and Testing of a General Amber Force Field[END_REF] was prepared for ensuring an accurate description of the equilibrium geometry and torsional degrees of freedom of the di-8-ANEPPS molecule. As detailed in the Supporting Information, atom types and force constants were finely tuned in the Force Field (FF) in order to accurately reproduce the bond lengths of the DFT equilibrium structure, as well as the associated bond length alternation (BLA), which is prerequisite for reliably describing the NLO responses of push-pull π-conjugated compounds.

Relaxed potential energy scans (PESs) were computed at the M06-2X/cc-pVDZ level on a reduced fragment of di-8-ANEPPS, and fitted following the procedure described in Ref. 81.

Lipid molecules of the DPPC bilayer were described using the CHARMM 36 [START_REF] Klauda | Update of the CHARMM All-Atom Additive Force Field for Lipids: Validation on Six Lipid Types[END_REF][START_REF] Pastor | Development of the CHARMM Force Field for Lipids[END_REF] (C36) FF, while the SPC/Fw 83 FF was used for water molecules. The quality of the solvent description was verified in terms of density (1.06 g cm -3 against 1.00 g cm -3 for the experimental value) as well as radial distribution function (see Section 1.2 in SI).

Calculation of the first hyperpolarizability

Calculations of the static and dynamic NLO responses were performed at the TDDFT/M06-2X/6-311+G* level using the Gaussian16 77 package. Frequency-dependent calculations were carried out using an incident wavelength of 1064 nm. We focused in this work on the hyper-Rayleigh scattering (HRS) response, β HRS , and the associated depolarization ratio (DR), which are associated respectively to the second harmonic light intensity and to the symmetry of the molecular scatterer. Full expressions of these terms can be found in Ref. 84. In the case of the di-8-ANEPPS molecule embedded in the membrane, we also analyzed the diagonal β ZZZ component of the first hyperpolarizability tensor, with Z the direction normal to the bilayer.

The NLO response of the isolated di-8-ANEPPS molecule (i.e. non interacting with the membrane), was computed in gas phase and in water using the integral equation formalism (IEF) of the polarizable continuum model (PCM). [START_REF] Mennucci | Evaluation of Solvent Effects in Isotropic and Anisotropic Dielectrics and in Ionic Solutions with a Unified Integral Equation Method: Theoretical Bases, Computational Implementation, and Numerical Applications[END_REF][START_REF] Tomasi | Quantum Mechanical Continuum Solvation Models[END_REF][START_REF] Miertus | Electrostatic Interactions of a Solute with a Continuum. A direct Utilization of Ab Initio Molecular Potentials for the Prevision of Solvent Effects[END_REF] To evaluate the NLO response of the di-8-ANEPPS molecule inserted into the membrane while accounting for the dynamical fluctuations of its geometry, 120 structures containing both the chromophore and the lipid bilayer were extracted from uncorrelated snapshots of the MD production run (every 2 ns).

M06-2X/6-311+G* calculations were subsequently carried out on the 120 structures of the di-8-ANEPPS molecule, each being surrounded by ESP [START_REF] Momany | Determination of Partial Atomic Charges from Ab Initio Molecular Electrostatic Potentials. Application to Formamide, Methanol, and Formic Acid[END_REF][START_REF] Singh | An Approach to Computing Electrostatic Charges for Molecules[END_REF] charges of the neighboring lipids in order to account for electrostatic environment effects. Each lipid possessing at least one atom in a 5 Å radius defined around any atom of di-8-ANEPPS was considered entirely. The water solvent surrounding the ANEPPS@DPPC system was accounted for using IEFPCM.

In addition, to relate the nonlinear to the linear optical properties, the vertical excitation energies and transition dipoles of the low-lying excited states of the di-8-ANEPPS molecule were also calculated at the TDDFT/M06-2X/6-311+G* level. The computed first hyperpolarizabilities were rationalized by relying on the two-state approximation (TSA), [START_REF] Oudar | Hyperpolarizabilities of the Nitroanilines and their Relations to the Excited State Dipole Moment[END_REF] which assumes that only one electronic excited state (here, generally the S 2 state) contributes to the sum-over-state expansion of the second-order NLO response:

β T SA ∝ µ 2 ∆µ ∆E 2 × ∆E 4 [∆E 2 -( ω) 2 ][∆E 2 -(2 ω) 2 ] ( 1
)
where ω is the energy of the incident photons, ∆E is the S 0 → S 2 excitation energy,

µ = || S 0 | - → µ |S 2 ||
is the norm of the associated transition dipole, and

∆µ = || - → µ S 2 -- → µ S 0 ||
is the norm of the dipole moment variation between the S 0 and S 2 electronic states.

3 Results and discussion

Morphology of the ANEPPS@DPPC system

Several structural parameters have been evaluated to characterize the membrane, namely the average area per lipid A , the membrane thickness, the tilt angle of DPPC (including the lipid chains and polar head), and the hydrocarbon chain order parameter. Full details on these parameters and their calculations are provided in SI. An average area per lipid of 56.3 ± 1.6 Å 2 is obtained for the leaflet without the chromophore, which is in the range predicted experimentally. The membrane thickness, measured as the average distance between the phosphorous atoms in the two bilayer leaflets, amounts to 39.3 ± 0.3 Å 2 , in good agreement with previously reported data (see SI for details). [START_REF] Lewis | Lipid Bilayer Thickness Varies Linearly with Acyl Chain Length in Fluid Phosphatidylcholine Vesicle[END_REF] Given the nature of its constituent moieties and overall flexibility of the DPPC molecules, three tilt angles were defined to get full insight into their orientation, namely θ P -N defined between the phosphate group and the nitrogen atom of the choline group in the polar head, and two θ C-C angles giving the orientation of the lipid alkyl chains (SI, Figures S7a andS7b).

For θ P -N , an average value of 76.0 ± 3.1°is obtained (SI, Figure S8a). This large angle allows for favorable interactions between the glycerol moieties and the water layer. The magnitude of the tilt angle is indeed controlled by the balance between the need for efficient hydrocarbon chain packing and the hydrophilic character of the polar head group. [START_REF] Chapman | Physical Studies of Phospholipids. VI. Thermotropic and Lyotropic Mesomorphism of some 1,2-diacyl-phosphatidylcholines (Lecithins)[END_REF] In relatively good agreement with experimental measurements reported in the literature, [START_REF] Levine | Physical Studies of Membrane Structure[END_REF] an average value 27.0 ± 3.1°is computed for θ C-C (see SI, Figure S8b), indicating a preferential orientation along the normal to the surface towards the interior of the membrane. Our description of the membrane is further validated by an hydrocarbon chain order parameter of 0.225 ± 0.003 , fully consistent with the experimental value of 0.20 ± 0.02, which is characteristic of the L α phase. [START_REF] Seelig | Lipid Conformation in Model Membranes and Biological Membranes[END_REF] The di-8-ANEPPS chromophore tends to orient preferentially along the normal to the membrane interface, with its polar head group interacting with the hydrophilic part of the lipid while alkyl tails are pointing towards the hydrophobic interior of the membrane. A more detailed picture is provided by the analysis of the tilt angle θ N -N of the long molecular axis with respect to the normal of the membrane plane (Figure 3a). A value of 90°indicates a parallel orientation of the chromophore with respect to the membrane surface, while values with respect to the bilayer normal. A very similar distribution is obtained on the 120 frames used for the subsequent evaluation of the NLO response, with an average value of 18.9 ± 10.7°, confirming that it constitutes a representative set of structures.

Mass density profiles of the relevant components of the system have also been calculated

to get insight into their relative position, and are displayed in Figure 4: Mass density profiles along the Z axis for the different components of the ANEPPS@DPPC system: DPPC (red) and its N 28 (orange) and P 23 (yellow) atoms, Chromophore (CHR, dark green) and its S 23 (light green), N 18 (yellowish green), and N 19 (dark blue) atoms, as well as water molecules (light blue). However, it is notably amplified by a few extreme β ZZZ values ranging between 1500 10 3 a.u. and 3000 10 3 a.u.. When excluding these structures from the statistics, the average dynamic β ZZZ and its standard deviation reduces to 444 10 3 ± 217 10 3 a.u., which is still seven times larger than the static ones. These larger β ZZZ responses are the consequence of resonant or near-resonant conditions, as confirmed by TDDFT/M06-2X/6-311+G* calculations of the excitation energies. Structures exhibiting intense NLO responses all display excitation wavelengths close to 532 nm, which corresponds to the second harmonic wavelength of a 1064 nm incident beam. Still, they are only a few, as shown by their distribution and the resulting visible absorption spectrum (Figure 6). Note that, although technological applications might take advantage of resonance-enhanced NLO responses, excitation-induced electron density redistributions can lead to photo-damages and should be avoided for biological applications. A linear correlation was also unravelled between the β 1064nm ZZZ /β ∞ ZZZ ratio and the first dipole-allowed excitation energy ∆E of the chromophore, or more specifically the dispersion factor at the second harmonic frequency, defined as 1/(∆E -2 ω) with ω the energy of the incident photons (Figure 7). This strong correlation demonstrates that the two state model describes most of the frequency dispersion effects of β ZZZ .

NLO response of the chromophore within the lipid bilayer

β 1064nm ZZZ /β ∞ ZZZ 1/(∆E -2hω) (eV -1 )
Figure 7: Evolution of the β 1064nm ZZZ /β ∞ ZZZ ratio as a function of the excitation energy (∆E, eV) of the chromophore. For readability, few resonant frames were not included. The β and ∆E values were calculated at the IEFPCM/TDDFT/M06-2X/6-311+G* level for the chromophore surrounded by the lipid charges and implicit solvation. ω is the energy of the incident photons.

Structure-property relationships

In this section, we decipher the relationships between the structure of the chromophore and its NLO response. By definition, the normal β ZZZ component discussed above is expected to depend on the orientation of the chromophore within the bilayer, which is characterized by the θ N -N angle (Figure 3). The correlation between β ZZZ and θ N -N is evidenced in Figure 8a, where the third power of cos(θ N -N ) is used as abscissa because β is a rank-3 tensor. The correlation plot confirms that small β ZZZ values are associated with large θ N -N angles.

To establish structure-property relationships that are independent of the orientation of the chromophore within the lipid bilayer (or at least to reduce the orientational effects), the focus switches now from β ZZZ to β HRS , which accounts for averages over molecular orientations. Figure 8b shows the correlation between β ZZZ /cos 3 (θ N -N ), which approximates a hypothetical tilt-independent normal β tensor component, and the β HRS values. Overall, these two quantities are strongly correlated (with a correlation coefficient R 2 = 0.95), which indicates the adequacy of considering β HRS for further analyses. Like for β ZZZ , the time evolutions for the static and dynamic β HRS values given in Figure 9 also exhibit strong variations while their cumulative averages rapidly converge. The average static β HRS amounts to 35 10 3 ± 7 10 3 a.u. whereas at 1064 nm it is more than three times larger and attains 125 10 3 ± 40 10 3 a.u.. Then, the average static and dynamic DR amount to 4.88 ± 0.02 and 4.99 ± 0.01, respectively, and are typical of one-dimensional π-conjugated push-pull molecules. The small standard deviation on DR indicates that it is virtually unaffected by the dynamical nature of the system, which can be correlated to the small departures of the di-8-ANEPPS geometry from planarity (vide infra) along the dynamics. Then, we assess how the dynamical fluctuations in the geometry of the chromophore affect the β HRS response, focusing in particular on two structural parameters, namely the BLA and the torsional angles θ 1-3 , as defined in Figure 1. Average values and standard deviations are reported in Table 1 while individual values are reported in SI (Figures S13a, S13b, S13c, S13d). Our simulation shows a rather small BLA value of 0.06 ± 0.03 Å, consistent with the DFT geometry used for the FF parametrization, and associated with a strong π-electron delocalization between the donor and acceptor groups. However, about 20 % of the extracted structures exhibit a value outside the "one standard deviation" range (S13a). While some of them can be linked to extreme values of β HRS , there is no global correlation between the two quantities, as shown in Figure S14a. Moreover, along the simulation, the chromophore deviates from planarity with average out-of-plane distortions of about 10°. Here also extreme values of the torsional angles are observed, but no correlation is found between β HRS and any of the investigated dihedrals (see Figure S14b, S14c, and S14d). The lack of correlation between β HRS and the selected structural parameters indicates that the variations in the NLO response cannot be ascribed to the dynamical fluctuations of a single parameter, but result from global distortions of the molecular structure. Note also that the fluctuations of β values can also be induced by variations of the lipid environment (i.e. of the position of the ESP charges surrounding the chromophore), which makes the analysis even more complex.

Environment effects

In this part we analyze the influence of the environment on the NLO responses of the di-8-ANEPPS molecule. In Table 2, the β HRS and DR values computed for the isolated chromophore or within an implicit water solvent using the DFT or GAFF equilibrium geometries are compared with those calculated on structures extracted from MD simulations. The latter account for the structural dynamics of the system and include the electrostatic environment due to the lipid bilayer. Environmental effects are first highlighted by the large increase of the static β HRS and the small decrease of the dynamic β HRS when going from the gas phase to the implicit aqueous phase (19 10 3 a.u. versus 32 10 3 a.u. and 86 10 3 a.u. versus 80 10 3 a.u., respectively). This is attributed to the much larger static dielectric constant of water, in comparison to its optical one, which enhances the push-pull character of the chromophore.

Then, moving to the bilayer environment leads to a substantial increase of the dynamic β HRS , from 97 10 3 to 255 10 3 ± 198 10 3 a.u., which highlights both the specific effects of the environment and of the dynamical behavior. Comparison with additional gas phase β HRS calculations on the MD structures (Figure S15 in SI) demonstrates that i) in the static limit β HRS is almost not impacted by structural dynamics (β HRS = 19 10 3 ± 2 10 3 a.u. from MD versus 20 10 3 a.u. for the equilibrium structure) but its variations originate from the polarization due to the surrounding (β HRS = 32 10 3 ± 5 10 3 a.u. in the bilayer versus 19 10 3 ± 2 10 3 a.u. for the isolated chromophore) while ii) at λ = 1064 nm, the dynamical behavior plays also a role, which is attributed to the variations of the excitation energies and the presence of near-resonance situations. In fact, when only accounting for the dynamical structural behavior of the chromophore in the bilayer, which anyway is constrained by the DPPC molecules of the bilayer, the dynamic β HRS increases from 97 10 3 to 170 10 3 ± 88 10 3 a.u.. Then, when taking into account the effects of the surrounding on the β HRS calculations, the latter further increases from 170 10 3 ± 88 10 3 a.u. to 255 10 3 ± 198 10 3 a.u.. This second effect is consistent with the modifications of the visible absorption spectrum of the chromophore, which gets broader and exhibits a vibronic-like structure when accounting for the bilayer surrounding (Figure 6). In particular, the latter is responsible for an increased absorption intensity close to the SHG wavelength. Note that the corresponding excitation energy of the chromophore in its rpGAFF ground state equilibrium structure is even larger and attains 2.79 eV, in comparison to the maximum absorption energies of 2.55 and 2.60 eV for the MD structures with and without the bilayer surrounding, respectively. Finally, all DR values are very close to 5, the typical value of one-dimensional push-pull π-conjugated compounds but still, accounting for the bilayer environment reduces strongly the standard deviation (see also Figure S16 in SI).

The striking increase of β HRS after inclusion of the molecular probe into the lipid bilayer is particularly relevant for second harmonic imaging microscopy because it demonstrates that the surrounding has no detrimental effect on the SHG signal but, on the contrary, it enhances it.

Conclusions and outlook

A combined Molecular Dynamics (MD) -Quantum Chemistry (QC) method has been elaborated to investigate the embedding effects on the first hyperpolarizability of an amphiphilic dye in a lipid bilayer. The system consists of di-8-ANEPPS, a typical AminoN-aphtylEthenylPyridinium dye in a hydrated 1,2-DiPalmitoylPhosphatidylcholine (DPPC) membrane. The method associates MD simulations where classical force fields are employed with time-dependent density functional theory (TDDFT) calculations to evaluate the first hyperpolarizability of the dye in the charge electric field due to the surrounding lipids. Still, the force field of the dye had to be re-parameterized with respect to DFT geometries in order to account for the strong dependence of the first hyperpolarizability on the molecular geometry. The results highlight i) the huge increase of the first hyperpolarizability of the dye when inserted into the bilayer lipids, ii) the double origin of this enhancement resulting both from the dynamical structural flexibility of the chromophore and from the polarization of the bilayer surrounding, iii) the orientational motions of the dye in the bilayer, of which the backbone axis is tilted by about 20 ± 10°with respect to the bilayer normal, iv) the relationship between this tilt angle and β ZZZ , the diagonal β tensor component oriented along the bilayer normal, v) the one-dimensional NLOphore character of the dye, evidenced by depolarization ratios close to 5 and the validity of the two-state model, as well as vi) the difficulty to relate the β ZZZ and β HRS values to single geometrical parameters of the dye like the BLA and the torsion angles. These huge surrounding effects, which are beneficial for the SHIM technique, call for further investigating the role of the nature of the lipid bilayer on the second-order NLO responses as well as on the variations of these responses as a function of external perturbations, like those created by membrane potentials. In particular, future works will focus on investigating the correlations between the SHG responses and the relative orientations of the chromophore and of the lipid chains. So, focusing on more and more complex systems will finally enable us to closely model real bilayers and to perform comparisons with experimental data.

nonlinear optical properties (structure-β relationships, gas phase β HRS for the MD structures, and DR values for the MD structures). This material is available free of charge via the Internet at http://pubs.acs.org/. This material is available free of charge via the Internet at http://pubs.acs.org/.

Figure 1 :

 1 Figure 1: Structure of di-8-ANEPPS (top) and DPPC (bottom) with key atoms labeling, segment used to define the bond length alternation (BLA = 1 2 [d C 5 -C 6 -2d C 6 -C 7 + d C 7 -C 8 ]) of di-8-ANEPPS (red), and atoms used to define its torsion angles [θ 1 = θ(C 1 -C 5 -C 6 -C 7 ), θ 2 = θ(C 6 -C 7 -C 8 -C 13 ), and θ 3 = θ(C 16 -C 15 -N 19 -C 26 )].

Figure 2 )

 2 Figure 2), providing a 3840/125 = 30.72 water/lipid number ratio. This ratio corresponds to a weight content c = m H 2 O /(m H 2 O + m DP P C ) = 0.42. Note that experimental DPPC bilayers in the L α phase are fully hydrated for c = 0.36 55,56 or c = 0.40.57

Figure 2 :

 2 Figure 2: Simulation box containing 125 DPPC, the di-8-ANEPPS chromophore, and 3840 water molecules.

  of 0°and 180°correspond to perpendicular arrangements. The distribution of θ N -N (Figure3b) reflects the dynamical fluctuations, with values spanning a range of about 50°. The average tilt angle over the whole production run amounts to 19.3 ± 10.6°, confirming a quasi vertical orientation of the probe in the lipid bilayer. Alternative analyses of the chromophore tilt angle were performed by considering cos(θ N -N ) instead of θ N -N as well as their violin plots with interquartile ranges (FiguresS9 and S10) and they all conclude about the tilt of the chromophore. This global orientation of the di-8-ANEPPS molecule within the bilayer seems to be dictated by the lipid alkyl chains, which are also globally tilted (by 27.0 ± 3.1 °)

Figure 4 .Figure 3 :

 43 Figure 3: (a) Definition of the θ N -N tilt angle between the vector connecting the N 18 and N 19 atoms and the normal to the bilayer interface (Z axis) and (b) distribution of θ N -N (°) over the whole MD run.

Figure 5 65 10 3 ±

 5653 Figure 5 presents the time evolution of the static and dynamic β ZZZ hyperpolarizability components, as well as of their cumulative averages. Even though individual values vary strongly over time, the cumulative averages rapidly converge. The average static β ZZZ amounts to 65 10 3 ± 15 10 3 a.u., while frequency dispersion effects enhance the NLO response by one order of magnitude, with β ZZZ = 510 10 3 ± 418 10 3 a.u.. This large standard deviation highlights the influence of the dynamical structural changes on the dynamic first hyperpolarizability.

Figure 5 :

 5 Figure 5: Time evolution of β ZZZ (10 3 a.u.) [cross] and its cumulative average [line] for the chromophore surrounded by lipid charges and implicit solvent molecules, as calculated at the IEFPCM/TDDFT/M06-2X/6-311+G* level. (a) λ = ∞ and (b) λ = 1064 nm.

Figure 6 :

 6 Figure6: Distribution of the lowest-energy dipole-allowed excitation energies and oscillator strengths of the chromophore surrounded by lipid charges and implicit solvent molecules and visible absorption spectrum as obtained from their convolution with Gaussian functions having FWHM of 0.05 eV, in comparison to their analogs obtained without accounting for the effects of the surrounding (gas) while keeping the same geometries. The calculations were performed at the TDDFT/M06-2X/6-311+G* level. The vertical dashed line corresponds to the energy of the SHG photons.

Figure 8 :

 8 Figure 8: (a) Static β ZZZ (10 3 a.u.) as a function of cos 3 (θ N -N ) and (b) relationship between β ZZZ /cos 3 (θ N -N ) and β HRS (both in (10 3 a.u.).

Figure 9 :

 9 Figure 9: Time evolution of β HRS (10 3 a.u.) [cross] and its cumulative average [line] for the chromophore surrounded by lipid charges and implicit solvent molecules, as calculated at the IEFPCM/TDDFT/M06-2X/6-311+G* level. (a) λ = ∞ and (b) λ = 1064 nm.

Table 1 :

 1 Average values and standard deviations of the BLA (Å) and of the dihedral angles (°) of di-8-ANEPPS calculated over the 120 snapshots extracted from the MD simulation of the ANEPPS@DPPC system.

	BLA	θ 1	θ 2	θ 3
	0.06 ± 0.03 171.4 ± 7.4 170.8 ± 7.1 170.4 ± 6.2

Table 2 :

 2 Static and dynamic (λ = 1064nm) β HRS (10 3 a.u.) and DR of di-8-ANEPPS, calculated at the TDDFT/M06-2X/6-311+G* level for different structures, either corresponding to (ground state) equilibrium geometries or obtained from MD simulations. In the former case, the β calculations were performed for the isolated molecule or by accounting for solvent (water) effects using the IEFPCM scheme while in the latter case, they were performed with or without accounting for the bilayer environment. Level of approximation for optimizing the geometries b Treatment of the surrounding in β HRS calculations

		Structures(s)	Geometries a	Surrounding b	β HRS	static	DR	β HRS	dynamic	DR
		Equilibrium	DFT	Isolated	19	4.88	86		4.98
		Equilibrium	DFT	IEFPCM	32	4.83	80		4.97
		Equilibrium Standard GAFF	Isolated	9	4.73	31		4.88
		Equilibrium	rpGAFF	Isolated	20	4.87	97		4.98
	a	MD MD	rpGAFF rpGAFF	Isolated Bilayer	19 ± 2 4.77 ± 0.10 170 ± 88 4.89 ± 0.33 32 ± 5 4.88 ± 0.02 255 ± 198 4.99 ± 0.01

Acknowledgement

The calculations were performed on the computers of the Consortium des Équipements de Calcul Intensif (CÉCI, http://www.ceci-hpc.be) and particularly those of the Technological Platform of High-Performance Computing, for which the authors gratefully acknowledge the financial support of the FNRS-FRFC, of the Walloon Region, and of the University of Namur (Conventions No. 2.5020.11, GEQ U.G006.15, 1610468, and RW/GEQ2016), as well as on zenobe, the Tier-1 facility of the Walloon Region (Convention 1117545).

Supporting Information Available

Details on i) the molecular dynamics simulation (FF reparameterization for di-8-ANEPPS, assessment of the validity of the FF for water and the DPPC lipid bilayer) and on ii) the

Graphical TOC Entry


Ilja Van Meerveld

Eric Sauquet

Francesc Gallart

Catherine Sefton

Jan Seibert

Kevin Bishop

Aqua Temporaria Incognita

. Their effective management and protection, therefore, requires an understanding of both natural and artificial causes of intermittence.

.

They include intermittent streams that flow seasonally, ephemeral streams that only flow in response to rainfall or snowmelt events, and episodic streams that contain flowing water only during

extreme rainfall events [START_REF] Buttle | An Overview of Temporary Stream Hydrology in Canada[END_REF]. The terms used for these non-perennial streams vary [START_REF] Busch | What's in a Name? Patterns, Trends, and Suggestions for Defining Non-Perennial Rivers and Streams[END_REF] and more classes can be defined to describe the occurrence of pools with standing water [START_REF] Gallart | TREHS: An open-access software tool for investigating and evaluating temporary river regimes as a first step for their ecological status assessment[END_REF]. Here we use the term temporary stream to refer to all nonperennial streams, but we acknowledge that the exact naming and definition of these water bodies can have important implications for their legal protection [START_REF] Caruso | Science and policy integration issues for stream and wetland jurisdictional determinations in a semi-arid region of the western U[END_REF][START_REF] Nikolaidis | Towards sustainable management of Mediterranean river basins: policy recommendations on management aspects of temporary streams[END_REF][START_REF] Magand | Intermittent rivers and ephemeral streams: what water managers need to know[END_REF]. In Switzerland, for instance, streams that flow on average less than 347 days per year (over a 10 year period) are considered non-permanent streams and regulations for permanent streams, such as requirements for permits to discharge or withdraw water, might not apply. [START_REF] Bishop | Aqua Incognita: the unknown headwaters[END_REF] argued that the majority of the total stream length needs to be represented as a blank space on assessment maps because they are not part of any systematic documentation.

Temporary streams in headwater catchments are rarely included in assessments and often not even shown on maps. [START_REF] Levick | The Ecological and Hydrological Significance of Ephemeral and Intermittent Streams in the Arid and Semi-arid American Southwest In[END_REF] reported that 59% of the streams in the U.S. (excluding Alaska) are temporary, but their survey was based on 1:100,000 scale topographic maps and did not include stream segments shorter than 1.6 km, such that it excluded all temporary headwater streams.

Analyses based on the 1:24,000 scale National Hydrography Dataset (NHDPlus) suggested that 42% of all stream segments in the upper Colorado river are first-order streams with intermittent flow and that temporary streams make up 73% of the total stream length [START_REF] Caruso | Biophysical-Regulatory Classification and Profiling of Streams Across Management Units and Ecoregions1[END_REF]. All studies that have actually mapped temporary streams in the field have shown that they are far more prevalent than indicated by the dashed blue lines on maps [START_REF] Hansen | Identifying stream types and management implications[END_REF][START_REF] Fritz | Comparing the Extent and Permanence of Headwater Streams From Two Field Surveys to Values From Hydrographic Databases and Maps[END_REF]. For example, the Swiss national topographic map shows 0.68 km of streams in the 13 ha upper Studibach catchment but repeated field mapping has shown that there are at least 3.77 km of streams, of which 2.66 km (71%) did not have flowing water during the dry summer of 2018 [START_REF] Van Meerveld | Expansion and contraction of the flowing stream network alter hillslope flowpath lengths and the shape of the travel time distribution[END_REF]. Similarly, field mapping of stream heads during wet conditions in the 68 km 2 Krycklan catchment in northern Sweden showed that 76% of the fully expanded network was missing on the official map (Ågren et al., 2015). A lack of knowledge about the location and extent of temporary streams hampers their protection [START_REF] Caruso | Science and policy integration issues for stream and wetland jurisdictional determinations in a semi-arid region of the western U[END_REF][START_REF] Caruso | Biophysical-Regulatory Classification and Profiling of Streams Across Management Units and Ecoregions1[END_REF].

Temporary streams have high biodiversity and are home to many endemic species [START_REF] Stanley | Ecosystem expansion and contraction in streams[END_REF][START_REF] Meyer | The contribution of headwater streams to biodiversity in river networks[END_REF][START_REF] Stubbington | Temporary streams in temperate zones: recognizing, monitoring and restoring transitional aquatic-terrestrial ecosystems[END_REF]; the dry riverbed is an egg bank for aquatic invertebrates and seed bank for aquatic plants [START_REF] Brock | Drought and aquatic community resilience: the role of eggs and seeds in sediments of temporary wetlands[END_REF][START_REF] Steward | When the river runs dry: human and ecological values of dry riverbeds[END_REF]. The onset and cessation of flow significantly affects the species assemblage [START_REF] Pařil | Short-term streambed drying events alter amphipod population structure in a central European stream[END_REF]Sarremejane et al., in press). Connectivity of previously disconnected stream segments increases streamflow (e.g., [START_REF] Godsey | Dynamic, discontinuous stream networks: hydrologically driven variations in active drainage density, flowing channels and stream order[END_REF][START_REF] Jensen | Headwater stream length dynamics across four physiographic provinces of the Appalachian Highlands[END_REF][START_REF] Pate | Streamflow permanence in headwater streams across four geomorphic provinces in Northern California[END_REF]. Sediment and organic material that has collected in the dry river bed is flushed during the onset of flow, leading to high sediment and nutrient fluxes [START_REF] Hladyz | Flows and hypoxic blackwater events in managed ephemeral river channels[END_REF][START_REF] Fortesa | Analysing hydrological and sediment transport regime in two Mediterranean intermittent rivers[END_REF] and high rates of biogeochemical transformations and ecosystem respiration [START_REF] Acuña | Drought and postdrought recovery cycles in an intermittent Mediterranean stream: structural and functional aspects[END_REF][START_REF] Romaní | Microbial Availability and Size Fractionation of Dissolved Organic Carbon After Drought in an Intermittent Stream: Biogeochemical Link Across the Stream-Riparian Interface[END_REF][START_REF] Von Schiller | Chapter 3.2 -Nutrient and Organic Matter Dynamics in Intermittent Rivers and Ephemeral Streams[END_REF][START_REF] Hale | Dynamic stream network intermittence explains emergent dissolved organic carbon chemostasis in headwaters[END_REF]. The expansion of the flowing stream network during wet periods, furthermore, leads to a more direct connection between the hillslopes and the stream, resulting in shorter travel times [START_REF] Van Meerveld | Expansion and contraction of the flowing stream network alter hillslope flowpath lengths and the shape of the travel time distribution[END_REF] and the potential bypassing of riparian buffer strips [START_REF] Wigington | Stream network expansion: a riparian water quality factor[END_REF].

Even though it is the repeated presence and absence of flowing water that shapes temporary stream ecosystems and the onset and duration of flow in headwater temporary streams affect water quantity and quality in downstream perennial streams, there are very limited hydrological data for temporary streams. Temporary streams are generally not included in stream monitoring networks, and where present, they are sometimes only operated seasonally as the dry period is not considered interesting for water management [START_REF] Peters | ZeroFlow: A PUB (Prediction in Ungauged Basins) Workshop on Temporary Streams -Summary of Workshop Discussions and Future Directions[END_REF]. Even in experimental headwater catchments, gauging stations are usually placed at the point of perennial flow. As a result, temporary streams are largely underrepresented in hydrological studies and monitoring networks [START_REF] Benstead | An expanded role for river networks[END_REF][START_REF] Snelder | Regionalization of patterns of flow intermittence from gauging station records[END_REF][START_REF] Godsey | Dynamic, discontinuous stream networks: hydrologically driven variations in active drainage density, flowing channels and stream order[END_REF]. The lack of gauging of temporary streams has to be kept in mind when datasets are compiled to determine the abundance and variation in temporary stream dynamics. The catalogue of temporary streams in Europe collected as part of the SMIRES initiative [START_REF] Sauquet | A catalogue of European intermittent rivers and ephemeral streams[END_REF] highlights the high variation in their hydrological response.

Still, systematic analyses of the spatial patterns in the onset and cessation of flow or trends therein are difficult due to the lack of data (Tramblay et al., in press). For example, only 7% of the U.K. benchmark network of near-natural catchments that are considered suitable for the analysis of trends in streamflow are non-perennial [START_REF] Harrigan | Designation and trend analysis of the updated UK Benchmark Network of river flow stations: the UKBN2 dataset[END_REF].Only 10% of the more than 4000 gauging stations in France with daily discharge data available in the national HYDRO database (http://www.hydro.eaufrance.fr) are likely naturally intermittent (Figure 1b). The fraction was highest (22%) for gauging stations with a catchment area ≤10 km 2 . The ONDE (Observatoire National des Etiages) network was designed by the French Biodiversity Agency, https://ofb.gouv.fr/) to complement the hydrometric network and reports the hydrological state (flowing water, standing water in isolated pools, dry streambed) for 3350 tributary streams at least five times per year (once per month between May and October). Most sites (85%) are located on streams with a catchment area ≤100 km 2 and 20% of the sites have a catchment area ≤10 km 2 (Figure 1a). For almost half (49%) of the sites, there was at least one observation of no flow prior to January 2020 (Figure 1a). Hydrometric challenges in measuring very low flows may mean that the stream is considered to be flowing, even though it is dry (or vice versa). Furthermore, the reported zero flows often reflect a data issue rather than a real measurement of zero flow [START_REF] Zimmer | Zero or not? Causes and consequences of zero-flow stream gage readings[END_REF]. This means that data on zero flows cannot be used without looking at their quality codes. For example, analysis of the data in the aforementioned French HYDRO database showed that a notable proportion of the 730 gauging stations for which an occurrence of zero flow was reported, are unlikely to be temporary streams (Figure 1b). This (admittedly fairly subjective) assessment involved data screening of no-flow occurrence (seasonality of null values, consistency with historical droughts, etc.), recession curves (changes in discharge before and after the sequence of null values, etc.), additional information from nearby ONDE sites, and information from the hydrometric services (personal communication). For publicly available datasets knowledge about specific practices is needed before the discharge data can be used to compute statistics on no-flow events. For example, Environment Canada used to denote flows smaller than 1 l s -1 as zero [START_REF] Peters | ZeroFlow: A PUB (Prediction in Ungauged Basins) Workshop on Temporary Streams -Summary of Workshop Discussions and Future Directions[END_REF]. We, therefore, recommend that especially for stations in temperate climates, where intermittence is not frequently observed, a label "temporary stream" should be added to the metadata of gauging stations to avoid any ambiguity.

Nevertheless, this status needs to be updated regularly, particularly after droughts.

For many applications (e.g., understanding ecological processes and biogeochemical cycling) it is crucial to know whether there is flowing water, or pools with standing water, or if the streambed is dry [START_REF] Gallart | A novel approach to analysing the regimes of temporary streams in relation to their controls on the composition and structure of aquatic biota[END_REF][START_REF] Bonada | Conservation and Management of Isolated Pools in Temporary Rivers[END_REF]. Gauging stations are designed to measure flow, not to provide information on the presence of pools. Pools of standing water are common and provide important refugia during dry periods (e.g., [START_REF] Marshall | Go with the flow: the movement behaviour of fish from isolated waterhole refugia during connecting flow events in an intermittent dryland river[END_REF]. At 39% of the sites in the ONDE network pools were observed at least once. For comparison, dry streambeds were observed at least once for 34% of the sites.

In light of the difficulties in gauging temporary streams, the costs associated with establishing and maintaining gauging stations, and the extremely high spatial variation in the occurrence of flow along the channel, new approaches to obtain data on the state of temporary streams are being tested.

While field mapping provides the most detailed spatial data [START_REF] Wigington | Stream network expansion: a riparian water quality factor[END_REF][START_REF] Malard | Flood-pulse and riverscape dynamics in a braided glacial river[END_REF][START_REF] Doering | Ecosystem expansion and contraction dynamics along a large Alpine alluvial corridor (Tagliamento River, Northeast Italy)[END_REF][START_REF] Godsey | Dynamic, discontinuous stream networks: hydrologically driven variations in active drainage density, flowing channels and stream order[END_REF][START_REF] Jensen | Headwater stream length dynamics across four physiographic provinces of the Appalachian Highlands[END_REF][START_REF] Sefton | Visualising and quantifying the variability of hydrological state in intermittent rivers[END_REF], it is difficult to do in headwater catchments during rainfall events because conditions can change quickly. Some studies have used drones with cameras for the mapping [START_REF] Spence | Deployment of an unmanned aerial system to assist in mapping an intermittent stream[END_REF][START_REF] Borg Galea | Mediterranean intermittent rivers and ephemeral streams: Challenges in monitoring complexity[END_REF][START_REF] Calsamiglia | Evaluating functional connectivity in a small agricultural catchment under contrasting flood events by using UAV[END_REF] but this is difficult for very small headwater streams, where vegetation is dense, or during intense rainfall events. Other studies have used low-cost electrical resistance [START_REF] Blasch | New Field Method to Determine Streamflow Timing Using Electrical Resistance Sensors[END_REF][START_REF] Goulsbra | A new approach to the application of electrical resistance sensors to measuring the onset of ephemeral streamflow in wetland environments[END_REF][START_REF] Bhamjee | Ephemeral stream sensor design using state loggers[END_REF][START_REF] Sherrod | A Low-Cost, In Situ Resistivity and Temperature Monitoring System[END_REF][START_REF] Chapin | Robust, low-cost data loggers for stream temperature, flow intermittency, and relative conductivity monitoring[END_REF][START_REF] Paillex | High stream intermittency in an alpine fluvial network: Val Roseg, Switzerland[END_REF] or temperature [START_REF] Ronan | Field study and simulation of diurnal temperature effects on infiltration and variably saturated flow beneath an ephemeral stream[END_REF][START_REF] Constantz | Heat as a tracer to determine streambed water exchanges[END_REF] sensors to determine the onset and cessation of flow. The sensor networks developed by [START_REF] Bhamjee | Monitoring ephemeral headwater streams: a paired-sensor approach[END_REF] and [START_REF] Assendelft | A Low-Cost, Multi-Sensor System to Monitor Temporary Stream Dynamics in Mountainous Headwater Catchments[END_REF] even allow differentiation of standing water (pools) and flowing water. Even though the initial tests of these sensors are promising, their use has yet to become commonplace, likely due to the need to invest in sensor development and maintenance. Aerial photographs, images from Google Street View and interviews with inhabitants have been used determine the medium-term state of temporary streams in populated areas [START_REF] Gallart | TREHS: An open-access software tool for investigating and evaluating temporary river regimes as a first step for their ecological status assessment[END_REF]. Physical and biological indicators can also be used to determine the duration of the flowing state for temporary streams [START_REF] Fritz | Classifying Streamflow Duration: The Scientific Basis and an Operational Framework for Method Development[END_REF].

Crowdsourcing or citizen science is an alternative approach to obtaining data on the state of temporary streams [START_REF] Kampf | Testing the waters: Mobile apps for crowdsourced streamflow data[END_REF]. Visual observations at a range of locations can lead to data with a relatively high temporal resolution (Figure 2) or to obtain detailed maps of the presence of flow along rivers [START_REF] Turner | Wet/Dry Mapping: Using Citizen Scientists to Monitor the Extent of Perennial Surface Flow in Dryland Regions[END_REF][START_REF] Allen | Citizen scientists document long-term streamflow declines in intermittent rivers of the desert southwest, USA[END_REF]. Although initial analyses suggest high interrater agreement [START_REF] Seibert | Wasserdaten sammeln mit dem Smartphone -Wie können Menschen messen, was hydrologische Modelle brauchen[END_REF], the accuracy and usefulness of these data still need to be determined. Furthermore, the involvement of the public is a challenge, particularly for national or international projects for which it is more difficult to organize local outreach events to raise awareness of the project and the importance of temporary stream observations. The involvement of the public can be a challenge as dry streams are valued less than flowing streams [START_REF] Armstrong | What's a Stream Without Water? Disproportionality in Headwater Regions Impacting Water Quality[END_REF] and because small streams are often overlooked. Citizen science helps to increase public awareness on environmental issues, and short-term, large-scale projects are particularly well-suited for this [START_REF] Pocock | Developing a Strategic Framework to Support Citizen Science Implementation in SEPA[END_REF]. To obtain repeated data from many sites, it is useful to engage environmental management agencies. When they include the quick citizen science-based approaches in their regular monitoring, a large number of additional data points can be collected. For example, the French authority SR3A in charge of water management for tributaries to the Rhone River used the CrowdWater approach (www.crowdwater.ch) to map the presence of flow in temporary streams (Figure 3). Together with the data from the ONDE network, these observations contributed to realtime monitoring of the state of the rivers and supported water restriction measures. Similarly, in the UK, 1050 observations were submitted for 145 spots between January 2019 and July 2020 using the CrowdWater app, mostly on chalk streams in the south-east of England. These observations complement surveys conducted by the Environment Agency. Knowledge of the patterns of intermittence for these groundwater-fed streams helps the agency to identify the impacts of abstractions and other stressors, track droughts, and inform ecological flow requirements. Hydrological models [START_REF] Williamson | Classification of Ephemeral, Intermittent, and Perennial Stream Reaches Using a TOPMODEL-Based Approach[END_REF][START_REF] Ward | Simulation of dynamic expansion, contraction, and connectivity in a mountain stream network[END_REF][START_REF] Yu | Quantifying spatial and temporal patterns of flow intermittency using spatially contiguous runoff data[END_REF][START_REF] Partington | What Triggers Streamflow for Intermittent Rivers and Ephemeral Streams in Low-Gradient Catchments in Mediterranean Climates[END_REF], topographic data [START_REF] Prancevic | Topographic Controls on the Extension and Retraction of Flowing Streams[END_REF] and statistical approaches [START_REF] Snelder | Regionalization of patterns of flow intermittence from gauging station records[END_REF][START_REF] Russell | A Spatially Explicit Model for Mapping Headwater Streams[END_REF][START_REF] González-Ferreras | Mapping the temporary and perennial character of whole river networks[END_REF][START_REF] Beaufort | A classification approach to reconstruct local daily drying dynamics at headwater streams[END_REF][START_REF] Jaeger | Probability of Streamflow Permanence Model (PROSPER): A spatially continuous model of annual streamflow permanence throughout the Pacific Northwest[END_REF][START_REF] Konrad | Estimating minimum streamflow from measurements at ungauged sites in regions with streamflow-gauging networks[END_REF][START_REF] Durighetto | Intraseasonal Drainage Network Dynamics in a Headwater Catchment of the Italian Alps[END_REF] have been used to predict where streams are temporary and can be used to determine where additional data on the state of temporary streams may be most useful. However, to train and validate these models, more observations of the state of temporary streams and stream network dynamics are needed. Data on the presence or absence of flowing water in different tributaries or the total flowing stream length can be used to calibrate hydrological models [START_REF] Stoll | Explicit simulations of stream networks to guide hydrological modelling in ungauged basins[END_REF] or to validate the simulations of the stream network from physically-based coupled surface-subsurface flow models. The comparison of observations and simulations is less direct for conceptual (i.e., bucket-type) models, but observations can be used indirectly in model calibration or validation because they provide information on storage dynamics.

The basic approach, in this case, is to compare the average or typical stream conditions to the dynamics of the simulated (groundwater) storage.

Previous model studies have focused on climate change and other human impacts on flow intermittence, particularly for Mediterranean catchments [START_REF] Jaeger | Climate change poised to threaten hydrologic connectivity and endemic fishes in dryland streams[END_REF][START_REF] Reynolds | Modeled intermittency risk for small streams in the Upper Colorado River Basin under climate change[END_REF][START_REF] Pumo | Climate change effects on the hydrological regime of small non-perennial river basins[END_REF][START_REF] Querner | Simulating streamflow variability and aquatic states in temporary streams using a coupled groundwater-surface water model[END_REF][START_REF] Tzoraki | Assessing the flow alteration of temporary streams under current conditions and changing climate by Soil and Water Assessment Tool model[END_REF][START_REF] Girolamo | Hydrology under climate change in a temporary river system: Potential impact on water balance and flow regime[END_REF] or globally (e.g., [START_REF] Döll | How is the impact of climate change on river flow regimes related to the impact on mean annual runoff? A global-scale analysis[END_REF]. They predict a shift from perennial to intermittent flow regimes and an increased duration of the dry state, which will impact freshwater ecosystems (e.g., [START_REF] Cipriani | Impact of climate change on aquatic ecosystems along the Asse river network[END_REF][START_REF] Jaeger | Climate change poised to threaten hydrologic connectivity and endemic fishes in dryland streams[END_REF]. Observations of trends in flow persistence in headwater streams can provide important information and an early warning of how the dynamics of larger streams may change due to climate or land-use change because even small changes can cause them to switch from being perennial to temporary.

We add our call for more studies on temporary streams to those of similar commentaries [START_REF] Larned | Emerging concepts in temporary-river ecology[END_REF][START_REF] Datry | Recent perspectives on temporary river ecology[END_REF][START_REF] Kampf | Testing the waters: Mobile apps for crowdsourced streamflow data[END_REF][START_REF] Shanafield | Science gets up to speed on dry rivers[END_REF]. In particular, we call on hydrologists and citizens to observe, sense and report the hydrological state of the aqua temporaria incognita. These data will improve our understanding of these unique streams and the impacts of climate and land use change and water management on them, both directly and through the testing and refinement of hydrological models. Without these data, it is as if we are trying to complete a puzzle on how headwater catchments function and how water affects ecological processes, while the majority of the puzzle pieces are hidden under the carpet. Recent studies provide some information on how many of the pieces are hidden, but our knowledge is so limited that we do not even know what is printed on them. This makes it impossible to complete the puzzle of our landscape and how it functions. Developing that understanding will not only expand our knowledge about temporary streams but will also entail a fundamental rethinking of how water is connected to landscapes. That is because the current understanding of high flows, when much of the water leaves the landscape, has not included a large component of the land-water interface during and after these high flowsaqua temporaria incognita.

Figure 1 .

 1 Figure 1. Frequency distribution of the number of sites in the ONDE network (a) and the number of French gauging stations available in the HYDRO database (b) as a function of catchment area. The sites for which at least once no flowing water (i.e., dry streambed or standing water in isolated pools) was observed and stations on potential intermittent streams are indicated by the filled area. For the gauging stations (b), all stations that include an occurrence of zero flow are shown with a dashed line (all data until 01.01.2020).

Figure 2 .

 2 Figure 2. Examples of eight-month time series of observations of the hydrological state of temporary streams made with the CrowdWater app for a site in Portugal (a) and a site in Switzerland (b). Note that the two Figures show a different period. Source: https://www.spotteron.com/crowdwater/spots/89106 (a) and https://www.spotteron.com/crowdwater/spots/245853 (b).

Figure 3 .

 3 Figure 3. The number of observations made with the CrowdWater app (Seibert et al., 2019) for temporary streams (colored circles,) by August 9, 2020, as well as the locations of the gauging stations (black triangles) and the ONDE sites (grey squares). The inset shows the location of the area in France and the Rhone river.Background elevation data from https://land.copernicus.eu/imagery-in-situ/eu-dem/eu-dem-v1-0-and-derivedproducts. The stream network was obtained from http://www.sandre.eaufrance.fr/.

Acknowledgements

We thank the SMIRES (Science and Management of Intermittent Rivers & Ephemeral Streams) COST Action CA15113 for facilitating our many discussions on temporary streams, Willy Bertin and Gaela Le Bechec from SR3A and Judy England from the Environment Agency for information on how and why they use the CrowdWater app, and Helena Ramos Ribeiro and Auria Buchs for collecting the data shown in Figure 2.


Fahed Zaïri
email: fahed.zairi@gmail.com

Mélissa Moulart

Christian Fontaine

Fahmi Zaïri

Vincent Tiffreau

Régis Logier

Relevance of a novel external dynamic distraction device for treating back pain

Keywords: Back pain, Actuators, Traction, Disc pressure, Anybody Modeling System

Introduction: Low back pain is a common, expensive and disabling condition in industrialized countries. There is still no consensus for its ideal management. Believing in the beneficial effect of traction, we developed a novel external dynamic distraction device. The purpose of this work is to demonstrate that external distraction allows limiting the pressure exerted on orthostatism on the lower intervertebral discs.

Materials and methods:

We firstly used the Anybody Modeling System, which is a validated musculoskeletal software. The device has been implemented in the model and we calculated the lower disc pressure while traction forces were applied by the four actuators. Secondly, we performed an anatomical study using a Biomet cadaver placed in a sitting position. The two belts of the device were tightened at the base of the thorax and on the waist, respectively. A pressure sensor was placed in the lower discs under fluoroscopic control through a Jamshidi needle. The intradiscal pressure was then measured continuously at rest while applying a traction force of 20 Kg.

Results: Both numerical and cadaveric studies have demonstrated a decrease in intradiscal pressures after application of a traction force by the external device. Using the numerical model, we have shown that tensile forces below 50 Kg in total were sufficient. The application of superior forces is useless and potentially deleterious.

Conclusion:

External dynamic distraction device is able to significantly decrease the intradiscal pressure in a sitting or standing position. The effects are obtained using "reasonable" traction forces, inferior to those applied by most conventional traction tables.

However, the therapeutic effects need to be proven using clinical studies.

Introduction

Back pain is becoming an increasing concern in most industrialized countries. Indeed, it is estimated that 80% of the adult population suffers or will suffer from low back pain with an annual prevalence of 30% [START_REF] Andersson | Epidomiological features of chronic low-back pain[END_REF][START_REF] Jayson | ABC of work related disorders: Back pain[END_REF] . This frequency more than tripled in Europe and in the US between 1980 and 2000 [START_REF] Andersson | Epidomiological features of chronic low-back pain[END_REF][START_REF] Jayson | ABC of work related disorders: Back pain[END_REF] . Chronic low back pain, whose evolution is by definition longer than 3 months, concerns only 20% of these patients, that is to say 10 to 20% of the general population with a preferred age ranged between 35 and 45 years. In industrialized countries, low back pain is a major public health problem because it generates significant costs in healthcare and is one of the leading causes of disability and sick leave [START_REF] Maher | Non-specific low back pain[END_REF][START_REF] Hartvigsen | What low back pain is and why we need to pay attention[END_REF] . From a physiopathological point of view, chronic low back pain is due to degenerative changes in the intervertebral disc and articular joints [START_REF] Borenstein | Mechanical low back pain: A rheumatologist's view[END_REF][START_REF] Bogduk | On the definitions and physiology of back pain, referred pain, and radicular pain[END_REF] . The lower segments are the more frequently and severely affected in clinical practice (L4L5 and L5S1 segments). Many risk factors have been highlighted such as smoking, overweight, sedentary lifestyle, and some occupations that require heavy loads [START_REF] Taylor | Incidence and risk factors for first-time incident low back pain: A systematic review and meta-analysis[END_REF][START_REF] Wippert | Risk and protective factors in the clinical rehabilitation of chronic back pain[END_REF] . The exaggerated mechanical stresses are clearly identified to accelerate the degeneration and thus aggravate the symptoms.

In spite of its frequency and its social impact, there is no consensus for the treatment of chronic low back pain [START_REF] Foster | Prevention and treatment of low back pain: Evidence, challenges, and promising directions[END_REF][START_REF] Buchbinder | Low back pain: A call for action[END_REF] . The surgery failed to show its effectiveness in the absence of radicular pain, except in very particular situations [START_REF] Copeland | Surgical versus nonsurgical treatment for back pain[END_REF][START_REF] Gaudin | Considerations in spinal fusion surgery for chronic lumbar pain: Psychosocial factors, rating scales, and perioperative patient education -A review of the literature[END_REF] . In the same way, the majority of conservative treatments have not shown clear evidence of their effectiveness [START_REF] O'connell | Clinical guidelines for low back pain: A critical review of consensus and inconsistencies across three major guidelines[END_REF] . Among them, lumbar traction is a commonly used method to treat patients with low back pain with or without sciatica. In the UK and the US, lumbar traction is used by 41 and 77% of outpatient rehabilitation providers, respectively [START_REF] Thackeray | The effectiveness of mechanical traction among subgroups of patients with low back pain and leg pain: A randomized trial[END_REF][START_REF] Oh | The impact of manual spinal traction therapy on the pain and Oswestry disability index of patients with chronic back pain[END_REF] . Thus, there is a discordance between the lack of evidence-based recommendation and how lumbar traction is regarded in current clinical practice, which is explained by the great heterogeneity of practices and the methodological problems found in most clinical studies [START_REF] Wegner | Traction for low-back pain with or without sciatica[END_REF] . Believing in the beneficial effect of traction, we developed an external distraction device. The purpose of this work is to demonstrate that external distraction allows to limit the pressure exerted on orthostatism on the lumbar spine and mainly on the intervertebral disc.

Materials and Methods

Numerical assessment a. "Anybody" model

Anybody Modeling System (AMS) is a validated musculoskeletal software revolving around inverse dynamic simulations able to calculate individual muscle forces, joint contact forces and pressures [START_REF] Bassani | Validation of the AnyBody full body musculoskeletal model in computing lumbar spine loads at L4L5 level[END_REF] . Each body part is implemented using validated cadaveric or anatomical data ensuring high accuracy and anatomical fidelity of the model. Every bones, joints, muscles, ligaments and tendons are represented (Figure 1). Finally, AMS allows the importation of computer assisted design (CAD) components in order to study their effects in interaction with the body. Note that anthropometric measures can be modified.

b. Exoskeleton model

The aim of the exoskeleton is to apply vertical traction forces to reduce pressure on the lumbar spine in the upright position (Figure 2A). In order to preserve both proper spinal alignment and freedom of motion, the device is composed of two sets of actuators positioned on both side of the body. The traction is produced by these four actuators, each of them represented by a pivot joint between the upper part and the lower part of the device (Figure 2B). These joints are activated by electric motors thanks to worm screws. At the end of each of them, there is a ball joint to preserve the motion of the trunk. Each actuator can generate a force of 8 Kg, ensuring a maximum distraction force of 32 Kg.

c. Implementation of the exoskeleton in Anybody

As the simulation does not include the skin, we used the torso on which the exoskeleton was designed to position the device in the simulation. The torso used for the design of the prototype is a standard morphology transmitted by a local orthoprosthesist. Then, a skeleton has been fitted inside the CAD of the body with size adaptation relative to the body dimension. The skeleton was obtained from a free CAD on Grabcad (Figure 3). The skeleton inside the CAD was then fitted to the Anybody Model to position the torso, therefore fitting the position of the exoskeleton in the simulation. Finally, the connection points were visually fitted to the solid part of the pelvis and the thorax.

Therefore, it was necessary to simplify the model (Figure 2C). First, a perfect model was applied for the actuators (no friction) to simplify the calculations. Secondly, as the AMS software performs the simulations by inverse dynamics operations, we had to constrain the DOF. The top ball joint has been replaced by a universal joint to avoid rotation in the actuators along the z axis. Similarly, the pivot joint between the two parts of each actuator has been replaced by a slide joint to limit rotation. Finally, the connection between the body and the belts has been simplified as a housing to prevent movement

d. Study protocol

For this study, we have chosen morphometric data to get closer to the European average for a man. We parameterized a height of 1.76 meters, a weight of 74 kilograms and a lumbar disc area of 19.8 cm 2 . The intradiscal pressure (P) was estimated by dividing the force (F) with the corresponding disc area applying a correction factor. The correction factor considered the non-uniform load distribution in the disc. In accordance with the Nachemson study [START_REF] Nachemson | Lumbar intradiscal pressure. Experimental studies on post-mortem material[END_REF][START_REF] Nachemson | Mechanical stresses on lumbar disks[END_REF] as confirmed by Brinckmann et al. [START_REF] Brinckmann | Change of disc height, radial disc bulge, and intradiscal pressure from discectomy. An in vitro investigation on human lumbar discs[END_REF] and Cripton et al. [START_REF] Cripton | In vitro axial preload application during spine flexibility testing: Towards reduced apparatus-related artefacts[END_REF] , we applied a correction factor of 0.66. Thus, the final formula to calculate the intradiscal pressure was P=F/(S×0.66). Even if the actuators of our device can not develop more than 8kg of traction force, we decided to realize extreme simulations of 0 to 100 Kg (25 Kg per actuator) in order to define the optimal traction force and to analyze changes of disc pressure during extreme traction.The measurements were performed for the two lower discs L4L5 L5S1, which are most frequently affected in clinical practice. We also decided to analyze the activity of the lumbar muscles to detect and measure any possible reaction contractile activity.

Anatomical study a. Cadaver characteristics

The full body of a 62-year-old woman (height 1.62 cm, weight 70 Kg) was used for the study.

The cadaver has been prepared using the Biomet liquid. She had no history of spinal surgery.

Radiological evaluation confirmed a disc high superior to 10 mm from L3 to S1. Similarly, there was nosclerosis of the endplates or voluminous osteophytes, overall confirming the absence of severe disc degeneration at those levels.

b. Measurement technique

A table of maintenance has been made for this study. The corpse is kept in a sitting position, in rectitude avoiding any support under the arms, likely to reduce the weight of the body. Due to the elasticity of tissues, we waited thirty minutes before starting any measurement to reach the plateau phase. A Jamshidi needle has been inserted percutaneously through a strict midline posterior approach. The placement has been performed under a strict anteroposterior and lateral fluoroscopic control (Figure 4). On the lateral radiograph, the disc was divided into three zones: posterior, medial and anterior (Figure 4C). The needle was placed respectively in the three zones beginning with the posterior, then median and finally anterior. Due to the bony overlapping of the iliac wing, it was difficult to certify the proper positioning of the needle at L5S1 level and measurements were thus made on the L3L4 and L4L5 discs. The device consists of two belts that were tightened at the base of the thorax and on the waist (over the iliac crests), respectively; no rigid attachment to the cadaver was performed, to maintain the potential effects of slips that could be encountered in therapeutic condition.

c. Measurement protocol

After having placed the needle, a pressure sensor was inserted through. We used a Gaeltec sensor coupled to the Picolog analysis software. Once the sensor introduced, awaiting time of 5 minutes was respected in order to reach an equilibrium, because of the minimal tissue lesions induced by the puncture. We then performed five measurements for each zone and for each disc. For each measurement the same protocol has been respected. After recording the base pressure, we applied a force of 5 kg per actuator (total of 20 kg) for 2 minutes. We then recorded the baseline pressure for 5 minutes to regain balance before starting a new measurement. We then performed a prolonged recording to measure the evolution of the disc pressure when a prolonged distraction is applied. We placed the sensor in the middle of the disc and after observing a latency of 5 minutes, we measured the disc pressure while a force of 5 Kg per actuator for 30 minutes was applied.

Results

Numerical assessment: AMS a. Development of the intradiscal pressure

For the L4L5 disc, it was found that the disc pressure gradually decreased from 0.41 to 0.2 MPa for a total traction force close to 50 Kg. By increasing the tensile force, paradoxically, a rise in the disc pressure was observed (Figure 5A).

For the L5S1 disc, the intradiscal pressure decreased from 0.42 to 0.17 MPa for a total traction force close to 60 Kg. Similarly, by increasing the traction, we observed a rise in the discal pressure that reached 0.2 MPa for 100 Kg of traction (Figure 5B).

Under the conditions of this simulation, the optimum total traction force seems to be close to 50 kg. In addition to safety concerns, a higher traction force seems useless or even deleterious.

b. Muscles activity

At the lumbar level, we have three powerful and stabilizing muscles that are inside and outside, the multifidus, the longissimus and the iliocostal, acting in compression [START_REF] Park | Paraspinal muscles of patients with lumbar diseases[END_REF] . The activity of these muscles was recorded in parallel with the disc pressure (Figure 6). We found that there was little muscle activity for traction forces below 40 Kg. Beyond 50 Kg, we observed a significant and progressive increase in the activity of these three muscles, exceeding twice the basis activity for the longissimus and three times the basis activity for both multifidus and iliocostal.

Anatomical study

As stated previously, the measurements were performed on L3L4 and L4L5 discs, as L5S1 disc was not easily identifiable on the fluoroscopic control, because of bony overlapping of the iliac wing, which did not allow confirming the proper positioning of the sensor with certainty.

For the L3L4 disc, we measured a significant decrease in intradiscal pressure during the distraction phase which remained stable. This decrease was reproducible in the five completed registrations. Standardized results are shown in Figure 7. We found that the decrease in pressure was greater in the middle and at the back of the disc, whereas it was less significant at the front of the disc. Indeed, the pressure drop reaches 43.96% in the middle, 27.82% in the back and only 17.90% in the front of the disc. For the L4L5 disc, we also measured a significant drop in pressure after activation of the actuators. This decrease was as stable and reproducible in the five recordings made. We found that the pressure drop was more significant at the back and the middle of the disc while it was minimal at the front of the disc. The normalized averages obtained for each zone of each disc are shown in Figure 8.

We also found that the pressure drop obtained under the effect of traction was durable over time. Indeed, the decrease in pressure recorded in the middle of the L3L4 disc remained significant (up to 40%) beyond thirty minutes of continuous traction (Figure 9).

Discussion

Low back pain is a common, expensive and disabling condition in industrialized countries [START_REF] Andersson | Epidomiological features of chronic low-back pain[END_REF][START_REF] Jayson | ABC of work related disorders: Back pain[END_REF][START_REF] Maher | Non-specific low back pain[END_REF][START_REF] Hartvigsen | What low back pain is and why we need to pay attention[END_REF] .

The pathophysiology is complex but the exaggerated mechanical stresses were clearly identified as a main pejorative factor. The existing therapeutic solutions are multiple but to date the prognosis is still often unfavourable, reflecting the need for new therapeutic tools [START_REF] Moseley | Innovative treatments for back pain[END_REF] .

Among the existing solutions, traction is very popular but no study has been able to demonstrate its clinical effectiveness in the medium or long term [START_REF] Wegner | Traction for low-back pain with or without sciatica[END_REF] . However, experimental studies have shown that traction tables are likely to increase the height of the intervertebral disc and even reduce the conflicts between the disc and the nerve roots in case of associated sciatica [START_REF] Chow | Mechanical effects of traction on lumbar intervertebral discs: A magnetic resonance imaging study[END_REF][START_REF] Kuo | Spinal traction promotes molecular transportation in a simulated degenerative intervertebral disc model[END_REF][START_REF] Chung | Herniated lumbar disks: Real-time MR imaging evaluation during continuous traction[END_REF][START_REF] Tadano | Lumbar mechanical traction: A biomechanical assessment of change at the lumbar spine[END_REF] . To date, the lack of evidence of effectiveness is likely to be related to a lack of technical solution, rather than a lack of concept. In this perspective, we developed an exoskeleton to obtain a distraction in orthostatism. Thus, the traction can be applied more prolonged, on a subject in a position of function and potentially in motion and in activity insofar as the actuators allow the maintenance of the amplitudes of movement. The objective would be to reduce the mechanical stress exerted on the lumbar disc (L4L5 and L5S1 being the most frequently affected), while maintaining the activity of the patient (recreation or professional). It would also aim to limit the muscular deconditioning caused by the inactivity or rigid contention belts sometimes proposed.

In our study we demonstrated that an external distraction was able to induce a significant decrease in intradiscal pressure. Using the numerical model, we have shown that tensile forces below 50 Kg in total were sufficient. The application of superior forces is useless because it is accompanied by a deleterious increase of the disc pressure. This effect can be explained in large part by the reflex muscular activity (recorded in our study) but also by the elastic properties of the ligament and tendinous structures not represented in the anybody model. Indeed, the developers stated: "when building the lumbar spine model the original idea was to include ligaments as well. But at some point, we decided not to include the ligaments in the model, because of lack of readily available information about the mechanical properties and slack lengths. We were in fear that ligaments with wrong properties might give worse results than excluding them". Note that most traction tables currently apply forces up to 100 Kg in clinical practice [START_REF] Tadano | Lumbar mechanical traction: A biomechanical assessment of change at the lumbar spine[END_REF][START_REF] Farajpour | Effects of different angles of the traction table on lumbar spine ligaments: A finite element study[END_REF] . The absence of representation of the ligaments is not the only limit of this model. The skin is also not represented in AMS, and the device had to be attached to the skeleton. As a result, the sliding of the device on the skin and the soft tissues is not taken into account, which may increase the physiological effects observed. In addition, the disc pressure is not directly measured but calculated according to the defined surface and by means of a chosen correction factor, which can be a source of approximation. The biomechanics of the soft tissue should be considered in further more accurate models [START_REF] Derrouiche | Osmo-inelastic response of the intervertebral disc[END_REF][START_REF] Derrouiche | A chemo-mechanical model for osmo-inelastic effects in the annulus fibrosus[END_REF] . In the current study, it was therefore necessary to carry out direct measurements "in vivo" in order to limit these biases.

Unlike cadavers prepared with formaldehyde, we used a BIOMET cadaver, in order to preserve much of the elasticity of the tissues. This cadaveric study has demonstrated that the application of external traction significantly reduces disc pressure. This pressure drop is significant, reproducible and durable over time, as demonstrated during the prolonged recording of thirty minutes, allowing to appreciate the potential therapeutic effects. Note that we measured larger effects in the front and middle of the disc while the pressure drop was lower at the front of the disc at the L3L4 and L4L5 levels. This is probably related to a postural effect. Indeed, the corpse was sitting slightly leaning forward, which can induce a slight inversion of curvature (lower lordosis) and increase the stresses exerted on the front of the disc [START_REF] Tadano | Lumbar mechanical traction: A biomechanical assessment of change at the lumbar spine[END_REF][START_REF] Farajpour | Effects of different angles of the traction table on lumbar spine ligaments: A finite element study[END_REF] .

This cadaveric study has certain limitations. First, muscle activity is non-existent and its effect cannot be measured. In addition, it would have been useful to perform measurements on several subjects to confirm reproducibility. However, the combination of a simulation on a validated model and a cadaveric study with direct measurements, makes it possible to validate the effect of external distraction on the decrease in intradiscal pressure. The therapeutic effect of this device deserves to be carefully studied, and for this an observational clinical study is currently ongoing.

Conclusion

In summary, external dynamic distraction device is able to significantly decrease the intradiscal pressure in a sitting or standing position. The effects are obtained using "reasonable" traction forces, inferior to those applied by most conventional traction tables.

However, the therapeutic effects need to be proven using clinical studies.

Legends

Figure 1 :

 1 Figure 1: Anybody Modeling System default model.

Figure 2 :

 2 Figure 2: Conceptual design of the device (A). The exoskeleton is composed of two belts tightened at the base of the thorax and on the waist (over the iliac crests), respectively, and linked by four actuators able to apply a traction force of 8 Kg each. The kinematic of the exoskeleton is schematized (B). A simplified kinematic of the exoskeleton has been implemented (C).

Figure 3 :

 3 Figure 3: Implementation of the exoskeleton in Anybody.

Figure 4 :

 4 Figure 4: Anatomical protocol. The corpse is placed in a sitting position under strict fluoroscopic control (A). The Jamshidi needle is inserted percutaneously under anteroposterior and lateral fluoroscopy (B). On the lateral view (C) the disc is divided into three zone (posterior, middle, anterior). On the AP view, we ensure the strictly median placement (D).

Figure 5 :

 5 Figure 5: Evolution of the intradiscal pressure in L4L5 (A) and L5S1 (B) as a function of the total traction force applied.

Figure 6 :

 6 Figure 6: Evolution of the activity of the three main erector muscles of the back as a function of the total traction force applied.

Figure 7 :

 7 Figure 7: Evolution of the intradiscal pressure in each zone of the L3L4 disc. Standardized results.

Figure 8 :

 8 Figure 8: Evolution of the intradiscal pressure (average of the five measurements) for each zone of the L3L4 (A) and L4L5 discs (B).

Figure 9 :

 9 Figure 9: Prolonged recording in the middle of the L3L4 disc. There is a prolonged and significant decrease in pressure when traction is maintained.

Figure 1 :Figure 2 :

 12 Figure 1: Anybody Modeling System default model

Figure 3 :

 3 Figure 3: Implementation of the exoskeleton in Anybody.Different steps: As the simulation does not include the skin, we used the torso on which the exoskeleton was designed to position the device in the simulation. Then, a skeleton has been fitted inside the CAD of the body with size adaptation relative to the body dimension. The skeleton inside the CAD was then fitted to the Anybody Model to position the torso, therefore

Figure 4 :

 4 Figure 4: Anatomical protocol. The corpse is placed in a sitting position under strict fluoroscopic control (A). The Jamshidi needle is inserted percutaneously under anteroposterior and lateral fluoroscopy (B). On the lateral view (C) the disc is divided into three zone (posterior, middle, anterior). On the AP view, we ensure the strictly median placement (D).

Figure 5 :

 5 Figure 5: Evolution of the intradiscal pressure in L4L5 (A) and L5S1 (B) as a function of the total traction force applied.

Figure 6 :

 6 Figure 6: Evolution of the activity of the three main erector muscles of the back as a function of the total traction force applied.

Figure 7 :Figure 8 :

 78 Figure 7: Evolution of the intradiscal pressure in each zone of the L3L4 disc. Standardized results.

Figure 9 :

 9 Figure9: Prolonged recording in the middle of the L3L4 disc. There is a prolonged and significant decrease in pressure when traction is maintained.

Acknowledgments

The authors would like to thank the Clinical Anatomy Laboratory for providing the cadaver.

The authors are grateful to Maurice Demeulaere, Fabien Descamps, Antoine Vanholber and Franck Stevendart for their assistance during the experiments.

Conflict of interests

None declared.