Stéphane Abel
email: stephane.abel@cea.fr

Massimo Marchi

Deciphering the Structure of the Gramicidin A Channel in Presence of AOT Reverse Micelles in Pentane using Molecular Dynamics Simulations

Keywords: Gramicidin A, SDS micelle, DMPC membrane, AOT reverse micelle, Molecular Dynamics Simulation

Structural studies of proteins and in particular integral membrane proteins (IMPs) using solution NMR spectroscopy approaches is challenging due not only to their inherent structural complexities, but also to the fact that they need to be solubilized in biomimetic environments (such as micelles), which enhances the slow molecular reorientation. To deal with these difficulties and increase the effective rate of molecular reorientation, encapsulation of the IMPs in the aqueous core of the reverse micelle (RM) dissolved in a low viscosity solvent has been proved to be a viable approach. However, the effect of the reverse micelle (RM) environment on the IMPs structure and function is little known. To gain insights into these aspects, this paper presents a series of atomistic unconstrained molecular dynamics (MD) of a model ion channel (gramicidin A, gA) with RMs formed with the anionic surfactant di-acyl chain bis(2-ethylhexyl) sodium succinate (AOT) in pentane at a water-to-surfactant molar ratio (W0) of 6. The simulations were carried out with different protocols and starting conditions and for a total of 2.4 µs and were compared with other MDs performed with the gA channel inserted in models of SDS micelle or DMPC membrane. We show here that in presence of AOT RMs, the gA dimer did not look like the "dumbbell-like" model anticipated by experiments, where the C-term parts of the gA are capped with two RMs and the rest of the dimer protected from the oil solvent by the AOT acyl chains. In contrast, the MD simulations reveal that the AOT, Na + and water formed two well defined and elongated RMs attached to the C-term ends of the gA dimer, while the rest is in direct contact with the pentane.

The initial β 6.3 secondary structure of the gA is well conserved and filled with 6 -9 waters, as in SDS micelle or DMPC membrane. Finally, the water movements inside the gA is strongly affected by the presence of RMs at each extremity and no passage of water molecules through the gA channel is observed even after a long simulation period, whereas the opposite was found for gA in SDS and DMPC.

I. Introduction

Integral Membrane proteins (IMPs), which can only be released from the membrane by disruption of the membrane, represent roughly 30 % of the proteins total coded from the genome of the eubacterial, archaean, and eukaryotic organisms [START_REF] Wallin | Genome-Wide Analysis of Integral Membrane Proteins from Eubacterial, Archaean, and Eukaryotic Organisms[END_REF]. To maintain them stable and functional for biophysical studies outside their native lipidic environments, a common strategy is to solubilize these proteins in membrane-mimetic environments (e.g. detergent micelles, amphipathic polymers, nanodiscs) (2)(3)(4), which makes their investigation even more complex. Due to their crucial role in many biological functions, diseases or in drugs target, structural and functional characterization of IMPs is a major research topic and still today remains a challenge.

Among the numerous biophysical approaches available to study soluble and IMPs, nuclear magnetic resonance (NMR) spectroscopy of protein solutions is a powerful technique and can yield not only information on the structure and dynamics of proteins [START_REF] Danmaliki | Solution NMR Spectroscopy of Membrane Proteins[END_REF], but also details on their interactions with small ligands [START_REF] Furukawa | Quantitative Analysis of Protein-Ligand Interactions by NMR[END_REF]. However, a fundamental limitation of solution NMR spectroscopy is its size limitation caused by the slow molecular reorientation time (τm) and shorter NMR signal relaxation times (T2), which limits its use to proteins of molecular weights around 30 -35 kDa [START_REF] Yu | Extending the Size Limit of Protein Nuclear Magnetic Resonance[END_REF]. To escape these fundamental barriers, advances in both hardware and experimental designs have been made possible in the past by the use of transverse relaxation optimized spectroscopy (TROSY) [START_REF] Pervushin | Attenuated T2 Relaxation by Mutual Cancellation of Dipole-Dipole Coupling and Chemical Shift Anisotropy Indicates an Avenue to NMR Structures of Very Large Biological Macromolecules in Solution[END_REF] [START_REF] Liu | Simultaneous NMR Assignment of Backbone and Side Chain Amides in Large Proteins with IS-TROSY[END_REF] and perdeuteration [START_REF] Gardner | The Use of 2H, 13C, 15N Multidimensional NMR to Study The Stucture and Dynamics of Proteins[END_REF] [START_REF] Del Amo | Quantification of Protein Backbone Hydrogen-Deuterium Exchange Rates by Solid State NMR Spectroscopy[END_REF] that allows structural studies of proteins with molecular mass of above 100 kDa (see for instance, refs. [START_REF] Riek | Solution NMR Techniques for Large Molecular and Supramolecular Structures[END_REF] [START_REF] Fernandez | TROSY in NMR Studies of the Structure and Function of Large Biological Macromolecules[END_REF]). However, these approaches can be sometimes limited by the increased cost and decreased yields from protein perdeuteration and also by the loss of the richness of information provided from the abundance of 1 H nuclei.

A solution proposed by the Wand's group (see [START_REF] Kielec | Reverse Micelles in Integral Membrane Protein Structural Biology by Solution NMR Spectroscopy[END_REF] [START_REF] Babu | Characterization, and NMR Spectroscopy of Encapsulated Proteins Dissolved in Low Viscosity Fluids[END_REF] [START_REF] Fuglestad | Reverse Micelle Encapsulation of Proteins for NMR Spectroscopy[END_REF] [START_REF] Dodevski | Optimized Reverse Micelle Surfactant System for High-Resolution NMR Spectroscopy of Encapsulated Proteins and Nucleic Acids Dissolved in Low Viscosity Fluids[END_REF] [START_REF] Pometun | Cold Denaturation of Encapsulated Ubiquitin[END_REF] [START_REF] Dodevski | Optimized Reverse Micelle Surfactant System for High-Resolution[END_REF] [START_REF] Peterson | High-Resolution NMR Studies of Encapsulated Proteins in Liquid Ethane[END_REF] [START_REF] Jorge | Characterizing Protein Hydration Dynamics Using Solution NMR Spectroscopy[END_REF] [START_REF] Kielec | A Method for Solution NMR Structural Studies of Large Integral Membrane Proteins: Reverse Micelle Encapsulation[END_REF] [START_REF] Peterson | Novel Surfactant Mixtures for NMR Spectroscopy of Encapsulated Proteins Dissolved in Low-Viscosity Fluids[END_REF] [START_REF] Lefebvre | NMR Spectroscopy of Proteins Encapsulated in a Positively Charged Surfactant[END_REF] [START_REF] Nucci | High-Resolution NMR Spectroscopy of Encapsulated Proteins Dissolved in Low-Viscosity Fluids[END_REF] and references cited therein) in order to decrease significantly the rotational correlation time is to encapsulate proteins within the water cavity of reverse micelles dispersed in low viscosity apolar solvents. In sub-critical regime, the viscosity of some apolar solvents is several orders lower than that of water. The low-viscosity fluids can be short chain alkanes (e.g. pentane [START_REF] Van Horn | Use of Reverse Micelles in Membrane Protein Structural Biology[END_REF] [START_REF] Nucci | High-Resolution NMR Spectroscopy of Encapsulated Proteins Dissolved in Low-Viscosity Fluids[END_REF], ethane [START_REF] Dodevski | Optimized Reverse Micelle Surfactant System for High-Resolution[END_REF] [START_REF] Peterson | High-Resolution NMR Studies of Encapsulated Proteins in Liquid Ethane[END_REF]), noble gases (Xe) [START_REF] Meier | Reverse Micelles Dissolved in Supercritical Xenon: An NMR Spectroscopic Study[END_REF], or liquid CO2 [START_REF] Gaemers | NMR of Biomolecules in Low Viscosity, Liquid CO2[END_REF].

RM systems have been originally employed decades ago with Terahertz (29) [START_REF] Levinger | Ultrafast Energy Migration Pathways in Self-Assembled Phospholipids Interacting with Confined Water[END_REF] or NMR spectroscopy techniques to study the structure and dynamics of nucleic acids [START_REF] Workman | Stabilization of RNA Oligomers through Reverse Micelle Encapsulation[END_REF] or the proteins [START_REF] Kielec | Reverse Micelles in Integral Membrane Protein Structural Biology by Solution NMR Spectroscopy[END_REF] [START_REF] Babu | Characterization, and NMR Spectroscopy of Encapsulated Proteins Dissolved in Low Viscosity Fluids[END_REF] [START_REF] Fuglestad | Reverse Micelle Encapsulation of Proteins for NMR Spectroscopy[END_REF] [START_REF] Dodevski | Optimized Reverse Micelle Surfactant System for High-Resolution NMR Spectroscopy of Encapsulated Proteins and Nucleic Acids Dissolved in Low Viscosity Fluids[END_REF] [START_REF] Pometun | Cold Denaturation of Encapsulated Ubiquitin[END_REF] [START_REF] Van Horn | Use of Reverse Micelles in Membrane Protein Structural Biology[END_REF] [START_REF] Dodevski | Optimized Reverse Micelle Surfactant System for High-Resolution[END_REF] [START_REF] Peterson | High-Resolution NMR Studies of Encapsulated Proteins in Liquid Ethane[END_REF] [START_REF] Jorge | Characterizing Protein Hydration Dynamics Using Solution NMR Spectroscopy[END_REF] [START_REF] Kielec | A Method for Solution NMR Structural Studies of Large Integral Membrane Proteins: Reverse Micelle Encapsulation[END_REF] [START_REF] Nucci | High-Resolution NMR Spectroscopy of Encapsulated Proteins Dissolved in Low-Viscosity Fluids[END_REF] [START_REF] Meier | Reverse Micelles Dissolved in Supercritical Xenon: An NMR Spectroscopic Study[END_REF] [START_REF] Gaemers | NMR of Biomolecules in Low Viscosity, Liquid CO2[END_REF] and their solvation water [START_REF] Jorge | Characterizing Protein Hydration Dynamics Using Solution NMR Spectroscopy[END_REF] in a well-defined confined environment. The structural properties of the RMs (such as their size and shape) or their capacity to solubilize biomolecules depend significantly on a multitude of parameters such as the chemical nature of the surfactant [START_REF] Fanun | Microemulsions: Properties and Applications[END_REF] [START_REF] Walde | Phospholipid-Based Reverse Micelles[END_REF], the co-surfactant [START_REF] Palazzo | The Role of the Cosurfactant in the CTAB/Water/n-Pentanol/n-Hexane System: Pentanol Effect on the Phase Equilibria and Mesophase Structure[END_REF] or the water content inside the RM water pool (i.e. the water-to-surfactant ratio, W0) [START_REF] Fayer | Analysis of Water in Confined Geometries and at Interfaces[END_REF] [START_REF] Moilanen | Water Dynamics in Large and Small Reverse Micelles: From Two Ensembles to Collective Behavior[END_REF]. Consequently, the choice of a type of the RMs surfactant in NMR experiments to stabilize the protein in its native state is crucial.

For instance, the widely used anionic di-acyl chain bis(2-ethylhexyl) sodium succinate (AOT, Figure 1a) is found not a suitable (i.e.. denaturing) environments for many globular proteins [START_REF] Nicot | Proteins as Invited Guests of Reverse Micelles: Conformational Effects, Significance, Applications[END_REF] [START_REF] Larsson | Interactions of Native and Modified Cytochrome C with a Negatively Charged Reverse Micellar Liquid Interface[END_REF] [START_REF] Naoe | Higher Order Structure of Proteins Solubilized in AOT Reverse Micelles[END_REF] [START_REF] Ternström | Unfolding and Inactivation of Cutinases by AOT and Guanidine Hydrochloride[END_REF]. In that context, systematic studies were carried out to provide a library of surfactants and optimized RM systems able to successfully solubilize a large variety of soluble [START_REF] Dodevski | Optimized Reverse Micelle Surfactant System for High-Resolution NMR Spectroscopy of Encapsulated Proteins and Nucleic Acids Dissolved in Low Viscosity Fluids[END_REF] [START_REF] Pometun | Cold Denaturation of Encapsulated Ubiquitin[END_REF] or membrane [START_REF] Van Horn | Use of Reverse Micelles in Membrane Protein Structural Biology[END_REF] [START_REF] Nicot | Proteins as Invited Guests of Reverse Micelles: Conformational Effects, Significance, Applications[END_REF] [START_REF] Nicot | Limited Proteolysis of Myelin Basic Protein in a System Mimetic of the Myelin Interlamellar Aqueous Space[END_REF] [START_REF] Kielec | Reverse Micelles in Integral Membrane Protein Structural Biology by Solution NMR Spectroscopy[END_REF] [START_REF] Binks | Structural Parameters of the Myelin Transmembrane Proteolipid in Reverse Micelles[END_REF] proteins in their native states. These include, for instance, the cationic cetyl trimethyl ammonium bromide (CTAB) [START_REF] Lefebvre | NMR Spectroscopy of Proteins Encapsulated in a Positively Charged Surfactant[END_REF], the nonionic n-dodecyl tetra ethylene glycol (C12E4) [START_REF] Peterson | Novel Surfactant Mixtures for NMR Spectroscopy of Encapsulated Proteins Dissolved in Low-Viscosity Fluids[END_REF], the zwitterionic mono-acyl chain lauryl dimethylamine oxide (LDAO) [START_REF] Dodevski | Optimized Reverse Micelle Surfactant System for High-Resolution NMR Spectroscopy of Encapsulated Proteins and Nucleic Acids Dissolved in Low Viscosity Fluids[END_REF] or the decyl-1-rac-glycerol (10MAG) [START_REF] Fuglestad | Reverse Micelle Encapsulation of Proteins for NMR Spectroscopy[END_REF] [START_REF] Dodevski | Optimized Reverse Micelle Surfactant System for High-Resolution NMR Spectroscopy of Encapsulated Proteins and Nucleic Acids Dissolved in Low Viscosity Fluids[END_REF].

In case of IMPs, studies [START_REF] Kielec | Reverse Micelles in Integral Membrane Protein Structural Biology by Solution NMR Spectroscopy[END_REF] [START_REF] Van Horn | Use of Reverse Micelles in Membrane Protein Structural Biology[END_REF] showed that the surfactant could form a "dumbbell-like" or "shower cap" structures where the IMP, in its native fold, forms a "bridge" between two RMs where the protein hydrophilic domains are in direct contacts with the aqueous core of the RMs and its transmembrane domains are in direct contact the surfactant alkyl chain. These results were obtained for IMP of simple or complex structure such as gramicidin A (gA) in presence of AOT [START_REF] Van Horn | Use of Reverse Micelles in Membrane Protein Structural Biology[END_REF] [START_REF] Salom | Characterization of Gramicidin A in an Inverted Micellar Environment. A Combined High-Performance Liquid Chromatographic and Spectroscopic Study[END_REF] or the KcsA potassium channel with CTAB/hexanol surfactants [START_REF] Kielec | Reverse Micelles in Integral Membrane Protein Structural Biology by Solution NMR Spectroscopy[END_REF].

MD studies of soluble proteins with RMs are very rare (see for instance refs [START_REF] Senske | Local Chemistry of the Surfactant's Head Groups Determines Protein Stability in Reverse Micelles[END_REF] [START_REF] Abel | Molecular Dynamics Simulations of Cytochrome c Unfolding in AOT Reverse Micelles: The First Steps[END_REF] [START_REF] Tian | Simulations of the Confinement of Ubiquitin in Self-Assembled Reverse Micelles[END_REF]) and despite their great interests, in our best knowledge, no MD has been performed with IMPs. Consequently, the influence of the reverse micellar environment on the IMPs structures and functions is not well described at the atomic level. 2a). In biological membranes, it could dimerize in two major conformations: head-to-head left-handed dimer (or h-h dimer) [START_REF] Urry | The Gramicidin A Transmembrane Channel: A Proposed (L,D) Helix[END_REF] and in double-stranded helical dimers (or ds-dimer) (75) depending on the local environments such as the membrane thickness [START_REF] Mobashery | The Conformational Preference of Gramicidin Channels Is a Function of Lipid Bilayer Thickness 1[END_REF]. The h-h dimer, also termed the "pore" conformation adopting a β 6.3 sheet secondary structure allows the uncontrolled passage of monovalent ions (such as Na + ) across the membrane, thus compensating the ion gradient with the cytoplasm (77) [START_REF] Hladky | Ion Transfer across Lipid Membranes in the Presence of Gramicidin A. I. Studies of the Unit Conductance Channel[END_REF]. Due to this effect, gA is considered a natural bacterial antibiotic, especially against gram-positive bacteria (79) [START_REF] Koeppe | Engineering the Gramicidin Channel[END_REF]. High performance liquid chromatographic, circular dichroism [START_REF] Salom | Characterization of Gramicidin A in an Inverted Micellar Environment. A Combined High-Performance Liquid Chromatographic and Spectroscopic Study[END_REF] and high resolution NMR spectroscopy [START_REF] Van Horn | Use of Reverse Micelles in Membrane Protein Structural Biology[END_REF] experiments suggested that in presence of AOT, the two C-term portions of gA h-h dimer point toward the AOT RM aqueous core as it would in membrane, whereas the N-term is directed toward either the solvent or to form an N-N gA-dimer, thus bringing two RM entities together.

In this work, MD simulations of gA were carried out in AOT and pentane at W0 = 6 to mimic the experimental condition of Van Horn et al. [START_REF] Van Horn | Use of Reverse Micelles in Membrane Protein Structural Biology[END_REF] and by using different strategies to build the gA_AOT complex. These MD simulations were further compared with other MDs carried out with a gA dimer inserted in models of sodium dodecyl sulfate (SDS, Figure 1b) micelle and 1,2dimyristoyl-sn-glycero-3-phosphocholine (DMPC, C14:0, Figure 1c) membrane. This phospholipid was chosen since it forms a membrane with a lipid hydrophobic thickness around 25.0 Å, close to the length of the gA dimer, while retaining the structure and the activity of the gA channel [START_REF] Basu | Ion Channel Stability of Gramicidin A in Lipid Bilayers: Effect of Hydrophobic Mismatch[END_REF].

The outline of the remainder of this paper is as follows: in the section 2, we describe the construction of the different systems and the simulation methodology. In section 3, we discuss the MD results and, in particular, the aggregation process of the AOT in the presence of a gA dimer and the influence of the RM environment on the secondary structure and activity gA channel. These results were also compared with those obtained for gA micelle and membrane environments. Finally, some concluding remarks are presented in section 4.

Material and Methods

Construction of the models

As an initial structure for the gA dimer, we used in this work its 3D structure in its h-h conformation obtained from solution state 2D-NMR spectroscopy in SDS micelle (PDB entry: IJNO) [START_REF] Townsley | Structures of Gramicidins A, B, and C Incorporated into Sodium Dodecyl Sulfate Micelles[END_REF]. Each gA monomer was modeled by the CHARMM36m force field [START_REF] Huang | CHARMM36m: An Improved Force Field for Folded and Intrinsically Disordered Proteins[END_REF] with the amino-acid side chains in their typical protonation states at pH 7 and CMAP corrections for L-and D-amino acids [START_REF] Ingólfsson | Gramicidin A Backbone and Side Chain Dynamics Evaluated by Molecular Dynamics Simulations and Nuclear Magnetic Resonance Experiments. I: Molecular Dynamics Simulations[END_REF]. The AOT, lipids and pentane, were modeled with the parameters available in CHARMM27/36 force fields (82) [START_REF] Abel | Molecular Modeling and Simulations of AOT-Water Reverse Micelles in Isooctane: Structural and Dynamic Properties[END_REF]. The initial configurations of the micelle and membrane systems were built with a pre-assembled approach with the Micelle and the Membrane Builder modules in CHARMM-GUI (www.charmm-gui.org), respectively (84) [START_REF] Jo | CHARMM-GUI 10 Years for Biomolecular Modeling and Simulation[END_REF]. In case of the gA-SDS system (named gA_SDS, hereafter), we approximated the solution-state NMR experiments (26)(64) by building a system with 1 gA dimer in the h-h conformation, 65 SDS and Na + counter ions and 15000 TIP3P water [START_REF] Jorgensen | Comparison of simple potential functions for simulating liquid water[END_REF] leading to a concentration in SDS around 0.24 M. In case of the DMPC membrane system (named gA_DMPC, hereafter), we inserted a gA h-h dimer in a DMPC membrane composed of 90 lipids in each leaflet [START_REF] Kim | Influence of Hydrophobic Mismatch on Structures and Dynamics of Gramicidin A and Lipid Bilayers[END_REF]. The membrane was then solvated with 10000 TIP3P water (i.e. 55.5 waters per lipid) and 27 Na + and Cl -ions to have a NaCl concentration of 0.15 M. We point out that to improve the sampling and examine the influence of the simulation starting conditions on these MD results, we replicated these simulations twice (named hereafter, replica) with different random seeds.

For gA in presence of AOT RMs we used two protocols for building the initial structure of the systems, named hereafter self-assembled" (SA) and "pre-assembled" (PA). These two protocols were intended to evaluate the influence of the starting conditions and the validity of the "bridging model" by Van Horn et al. [START_REF] Van Horn | Use of Reverse Micelles in Membrane Protein Structural Biology[END_REF]. In the SA approach (named gA_AOT_SA, hereafter), the initial configuration of the system was built by placing 1 gA h-h dimer in the center of a simulation box with a side length of 78 Å immersed in 75 AOT molecules and Na + counterions, 450 waters and 8500 randomly oriented pentanes. To evaluate the influence of the starting conditions, the coordinates of the initial systems were randomly generated twice. In the PA protocol, instead, we constructed a system (named gA_AOT_PA, hereafter) with two identical spherical RM with each 40

AOT and its sodium counterions and with 225 waters inside, corresponding to a radius of 17 Å from the center of the sphere. This system was created with the Packmol program (88)(89) and using the same approach described in ref. [START_REF] Abel | Molecular Modeling and Simulations of AOT-Water Reverse Micelles in Isooctane: Structural and Dynamic Properties[END_REF]. These two RMs were then carefully minimized in gas phase with the CHARMM default simulation parameters. After these steps, we placed the gA dimer in its h-h conformation between the two RMs, forming a bridge. To avoid steric clashes between the AOT and gA dimer atoms, two AOT molecules in each RM were removed leading to a final complex with 1 gA and 76 AOT (and Na + ) and 450 waters, corresponding to a system with W0 = 5.92 (Figure S1). The overall complex was then fully minimized again, in gas phase and then inserted in cubic box with 8500 pentane with a dimension of 125 Å for each side. We observed that the resulting complex tends to be unstable during the course of the MD production and that the gA dimer could unfold very quickly during the simulation production (see below) if we use a too fast equilibration protocol (i.e. similar to that of for gA_AOT_SA system, see section 2.2). To avoid this problem, we have tested different protocols and propose a "slow" equilibration protocol with twenty steps (listed in Table S1) for a total a total of 18 ns where the simulation time step are slowly increase from 0.05 to 2 fs and the position restraints applied to each components of the gA_RM complex (i.e.; gA, AOT, Na + and water) are gradually decreased from 5000, 1000 and 0 kJ.mol -1 .nm -1 . In advance to the results presented below, we show that with this protocol, the overall secondary structure of the gA dimer is maintained during the course of the MD run as in the SDS and DMPC environments. The properties of the systems equilibrated using a "fast" and "slow" protocols (named hereafter "gA_AOT_PA_fast" and "gA_AOT_PA_slow", respectively) will be discussed together throughout this paper with the other systems. Finally, the characteristics of all the simulated systems are summarized in Table 1.

MD simulation details

All the simulations discussed in this paper were carried out with the CPU and GPU versions of GROMACS MD package (v2018.2) (90) in the NPT ensemble at a pressure adjusted to 1.015 bar and at 293.15 K for the direct and reverse micelles and at 303.15 K for the membrane [START_REF] Kim | Influence of Hydrophobic Mismatch on Structures and Dynamics of Gramicidin A and Lipid Bilayers[END_REF]. The former and the latter temperatures were chosen to mimic the NMR experimental conditions of Van Horn et al. [START_REF] Van Horn | Use of Reverse Micelles in Membrane Protein Structural Biology[END_REF] and to be above the DMPC membrane gel (Lβ)/liquid (Lα) transition (297 K (91)), respectively. To prepare the systems prior to the simulation productions, each system was minimized and equilibrated by using the CHARMM-GUI default simulation parameters. In brief, we used the steepest descent algorithm and an energy tolerance lower than 1000 kJ mol -1 nm -1 . In the second stage, the minimized systems were equilibrated at constant volume and temperature (NVT)

with T = 293.15 K (inverted and direct micelles) and T = 303.15 K (membrane) for 400 ps using the Berendsen thermostat with a coupling constant τT = 0.1 ps. This stage was followed by another equilibration stage carried out in the NPT ensemble using an isotropic (inverted and direct micelles) or semi-isotropic (membrane) pressure coupling schemes for 1 ns with the temperature and pressure controlled with the Nosé-Hoover thermostat (92,93) (τT = 0.1 ps) and Parrinello-Rahman barostat (τp = 5.0 ps and T = 45.10 -6 bar -1 ) [START_REF] Hoover | Canonical Dynamics: Equilibrium Phase-Space Distributions[END_REF][START_REF] Nosé | A Molecular Dynamics Method for Simulations in the Canonical Ensemble[END_REF]. During these stages the restraint forces on the gA_AOT_SA, gA_AOT_PA_fast, SDS and DMPC complexes were reduced 1000 to 0 kJ mol -1 nm -1 . In case of the gA_AOT_PA_slow system, as mentioned in the previous section, we used the equilibration stages summarized in Table S1. The electrostatic interactions were computed with the particle mesh Ewald (PME) method ( 96) with a cutoff of 12 Å and the van der Waals interactions smoothly were switched off at 10 -12 Å by a force-switching function [START_REF] Steinbach | New Spherical-Cutoff Methods for Long-Range Forces in Macromolecular Simulation[END_REF]. The P-LINCS algorithm [START_REF] Hess | LINCS: A Parallel Linear Constraint Solver for Molecular Simulation[END_REF] was applied to restrain bond lengths to their equilibration values. Periodic boundary conditions and a time step of 2 fs for integration of the equations of motion were used with the neighbor list updated every 100 fs. Finally, we performed production runs for 230/260 and 230/250 ns for the gA_SDS and gA_DMPC systems and 695/750 and 500/600 ns foe the gA_AOT_SA and gA_AOT_PA systems, respectively. We point out that all the production simulations were carried out without any harmonic restraints applied to the gA dimer and therefore its conformation during the MD was only driven by the force field and the environment. Ultimately, the atomic configurations of each system were collected every 10 ps and the structural and dynamical properties analyzed with different GROMACS tools and programs developed by us. Some representative snapshots of the systems extracted from the simulations are shown in Figure 3.

Results and discussion

Aggregation process of the AOT reverse micelle with the gA dimer

We begin by examining how the gA dimer interacts with the AOT RMs depending on how the gA_AOT complex were constructed. In particular, we investigate the aggregation dynamics of the two gA_AOT_SA systems constructed with the self-assembled approach. We show in Figure 4 where 𝑁 𝑐 (𝑡) is the number of clusters at the time t. The number of clusters at the beginning and

𝐺(𝑡) = 𝑁 𝑐 (𝑡) -𝑁 𝑐 (∞) 𝑁 𝑐 (0) -𝑁 𝑐 (∞) (1)
end of the simulation are then 𝑁 𝑐 (0) and 𝑁 𝑐 (∞), respectively. These values are obtained with the gmx clustsize tool of GROMACS (v2018.2) and the results are shown in Figure 5. We consider that a molecular contact between two AOT molecules exists if any of the distances between atoms on different molecules is less 2.8 Å. This cutoff distance was chosen by comparing the results (i.e. number of clusters) given by gmx clustsize with the number of clusters manually counted from selected PDB snapshots extracted at different simulation times. As shown in Figure 5, the details of the AOT aggregation for the two systems is somewhat different, but in both cases the cluster formation occurs with two different time scales: One fast (time tfast < 70 ns) and one slow ( 70< tslow < 700 ns). Previously, we found that the formation of RMs of AOT in isooctane can be fitted by a stretched exponential defined as 𝐺(t) = e (-x/t) α (100). For the new system here, we found that the best fit is obtained for the two MDs by a double exponential function:

𝐺(𝑡) = 𝑐 1 𝑒 (-𝑥/𝜏 1 ) + (1 -𝑐 1 )𝑒 (-𝑥/𝜏 2 ) (2)
Here, 𝜏 1 and 𝜏 2 are the characteristic relaxation times for the fast and slow processes, respectively (101)(102) and c1 parameter represents the waning number of AOT clusters. The fitting parameters of the double-exponential function are given in the Table in Figure 5. We first notice that the faster kinetics has a timescale similar for the two simulations with a relaxation time 𝜏 1 = 2 -2.7 ns. This faster aggregation process is followed by a relaxation time of 𝜏 2 = 16.9 -35.4 ns. It worth mentioning that the small variations of the numbers of AOT cluster in the two simulations after 300 ns are artefacts due to the choice of the cutoff used to estimate a contact between two AOT molecules, not due to the escape of AOT molecules in the oil phase.

The calculation of several structural parameters of the RMs obtained at the end of all the gA_AOT simulations are computed with our code OpenTRJ (103) and reported in Tables 2 and3. This includes the AOT aggregation numbers (NAOT), the RM water-surfactant ratio values (W0), the overall RM radius of gyration (𝑅 𝑔 𝑅𝑀 ) and finally, the eccentricities of the clusters, 𝑒 𝑅𝑀 = √(1 -𝑐 𝑅𝑀 2 /𝑎 𝑅𝑀 2 ). Here, 𝑒 𝑅𝑀 ranges from 0, for a sphere, to 1, for a rod or a disk). To obtain this latter parameter, the lengths of the major (𝑎 𝑅𝑀 ) and minor (𝑐 𝑅𝑀 ) semi-axis were computed from the whole RM moments of inertia, including the AOT, water and Na + ion (104). We notice that at the end of the gA_AOT_SA MD simulations, the two complexes are formed with two elongated RMs with different sizes and shapes and one free RM for gA_AOT_SA_rep1. The free RM (not attached to the gA dimer) in the gA_AOT_SA_rep1 have NAOT and W0 values of 37 and 6.7, respectively, and is more spherical (𝑒 𝑅𝑀 = 0.61) than the two RMs attached to the gA channel. As shown by the Figure 3d and the values reported in the Tables 2 and3 In the case of the two gA_AOT_PA systems constructed with the "pre-assembled" the fast and slow protocols, Figures S3 and S4 clearly show different behaviors at the beginning of the simulation. In the case of the gA_AOT_PA_fast complex, the micelle conformation significantly changes from its initial state during the course of the MD and deforms in less of 10 ns of simulation.

We observe that during this period 6 AOT molecules are transferred from one end of the gA dimer to another (not shown). The central region of the gA dimer is already in contact with the pentane after less than 20 ns of simulation. Thus, after only 50 ns, the gA dimer forms a bridge that connects the two extremities of the single elongated RM. Finally, after 300 ns of the MD production, the structure of the gA is no longer similar to the initial "bridging model" of Van Horn et al. [START_REF] Van Horn | Use of Reverse Micelles in Membrane Protein Structural Biology[END_REF] but looks like that of the final structure obtained for the gA_AOT_SA_rep2 complex constructed with the SA protocol. In the case of gA_AOT_PA_slow system (Figure S4), the conformation of the complex at the beginning of the production run is already similar to the gA_AOT_SA systems due to the longer equilibration period of 13 ns without constraints on the RMs. Also, during the course of the equilibration and production periods we did not observe transfer of AOT molecules between the two extremities of the grA. The final configurations of these two systems are shown in Figures 3e andf and the structural parameters of the two RMs attached to gA dimer at the end of these simulations are also reported in Tables 2 and3. They show that the bonded RMs in the gA_AOT_SA_fast system have similar W0, 5.7 and 5. and discussions in the next sections).

Interaction between gA dimer with surfactants and phospholipids.

As mentioned in the introduction, attempts to address the structural characteristics of the gA_AOT complex in pentane using pulsed-field-gradient NMR diffusion and NOESY experiments have been made in the past [START_REF] Van Horn | Use of Reverse Micelles in Membrane Protein Structural Biology[END_REF]. From our atomistic simulations, it is possible to examine the interactions between the AOT and the gA dimer during the course of the simulations and evaluate how the RM mimics the hydrophobic environments found in the SDS micelle or in the DMPC membrane. As mentioned in the previous section, our MD simulations did not reproduce the "bridging model" proposed in ref. [START_REF] Van Horn | Use of Reverse Micelles in Membrane Protein Structural Biology[END_REF]. Indeed, we find that the two C-term portions of the gA dimer interact with the AOT, whereas the rest of the gA peptide (including the N-term) interact with the pentane. To describe this in details, we computed the numbers of average contacts between each residue of the gA two dimers and the different parts of the AOT (i.e. sulfonate headgroup and the two ethyl hexyl chains), water and pentane (Figures 6c andS7-8). These results were compared with those obtained for headgroups and tails of SDS and DMPC obtained in other simulations (Figures 6ab andS5-6). In all these cases, to define a molecular contact between a protein residue and a surfactant group we use a distance criterion with a cutoff of 4 Å (105,106) between atoms on different groups. In this context, we consider the formyl (CHO) and the ethanolamine (EAM) groups at the N and C extremities of each gA monomer as two residues (named hereafter CHO0 and EAM16, respectively). From these figures, we notice that the patterns of the contacts in all simulations present similar features and do not change significantly between the two simulated replicas. In particular, we find that in the first half of the histograms delimited by residues (CHO0-DVal8), each gA monomer is in contact with the SDS and DMPC alkyl chains and, of course, with the water that fills the gramicidin channel. This is not unexpected as these residues are mostly buried in the micelle or in the membrane hydrophobic cores. In the case of the AOT RMs, these 9 residues are only in direct contact with pentane and the channel water, not with the AOT acyl-chain as suggested in ref. [START_REF] Van Horn | Use of Reverse Micelles in Membrane Protein Structural Biology[END_REF]. We also notice that the residues of the second half of the histograms (i.e. those between residues Trp9-EAM16) share some contacts not only with the SDS sulfonate and with the headgroup and tail of the DMPC, but also with the water inside the channel and at the micelle or membrane surfaces. The hydration of these residues increases depending on their localizations on the micelle and membrane surfaces. In the AOT RM systems, the numbers of contacts of these residues with the AOT headgroup, alkyl chain, the water inside the RM and the pentane depends on the size (and the water content) of the RM but also on the folding of the gA monomer as it is shown by Figures S7-8. When the RM is large enough the contacts between the gA C-term residues are very small.

Stability of the gA channel in the different environments.

We will now turn our attention to the structural variations of the gA dimer in presence of the SDS, DMPC and AOT RMs. We first examine the stability of the gA conformation by computing the root mean square deviation of the backbone atoms (or RMSDBC) of the gA dimer during the course of the eight simulations using as a reference the minimized structure of the gA dimer in each environment (Figure 7). Consistent with previous MD results ( 63 S9g).

To investigate the degree of similarity between the gA dimer in the four simulations and the initial β 6.3 helix secondary structure, we have also computed additional structural parameters such as endto-end distance, dEndToEnd, (64), the rise per turn, drise, and the average distances between the amide protons of selected residues of the gA dimer. As seen on Figures 8 and9 8d and9d). This behavior is particularly noticeable for the gA_AOT_SA_rep2 simulation where we observe a long increase (100 ns) of the gA dimer end-to-end distance, which reverse later in the simulation. The corresponding average drise and dEndToEnd values computed from the last 200 ns in the two MD are in range of 24.8 Å and 4.8-5.0 Å/turn, respectively. These are similar that the values found for the gA in the micelle and the membrane. In the case of dEndToEnd and drise computed for the overall gA dimer in the gA_AOT_PA simulations, we observed, again, some differences between the two systems caused by the partial unfolding of the gA dimer in one of the simulations. For the gA_AOT_PA_fast system, the dEndToEnd and drise are larger (32.0 Å and 5.8 Å/turn, respectively) than those found for the gA in the micelle, membrane or in two gA_AOT_SA replicas whereas the dEndToEnd and drise values are found similar (26.0 Å and 4.8 Å/turn) in the gA_AOT_PA_fast system.

Finally, we have computed the average distances between the intra chain NH distances between the residues X at the positions i and i-6 on the gA peptidic chain and the inter chain NH distances between the 6 first N-term residues (i.e. Val1, Gly2, Ala3, DLeu4, Ala5 and DVal6) of each gA monomer. These intra and inter chain NH distances in the canonical β 6.3 helical structure are six residues apart, lie on the same face of the β 6.3 helix, less than 5 Å of one another (Figure S10), and can be monitored by NMR spectroscopy [START_REF] Van Horn | Use of Reverse Micelles in Membrane Protein Structural Biology[END_REF]. The results of these calculations with those obtained from the 1JNO PDB structure are reported in Tables S2 and3 and confirm once again, that in presence AOT surfactant (except for gA_AOT_PA_fast) the β 6.3 helical structure is well conserved as in the SDS micelle or in the membrane, consistent with the NMR experiments (26).

The gA channel activity in the different environments.

The structural analysis of the gA dimer discussed in the previous sections clearly indicates that the initial β 6.3 -helix conformations is mostly retained in all environments (with the notable exception of gA_AOT_PA_fast) allowing the possible passage of water along the channel in the single-file configuration as it is shown in the representative snapshots provided in Figure 3. To examine further this behavior, for the eight systems investigated here, we computed different structural properties for the pore and its water. These included the pore average radius and profile, the water average pore occupancy, its conformation and also the movement of the water inside the pore. For these calculations, the center of mass of the group formed by the formyl and Val1-Gly2-Ala3 residues of the gA dimer was centered at each time frame at z = 0 Å and the gA channel aligned along the z axis of the simulation box.

Dimension of the channel pore. The gA channel radius (Rz) along the pore z axis for each simulations was calculated with the program HOLE (108) interfaced with GROMACS by using the gmx_clusterByFeatrues code (109). For these calculations, we included all the gA hydrogen atoms and used the van-der-Waals atomic radii of Bondi (110). The average pore radius profiles as well as the pore lining surfaces along z axis for the five simulated systems and the 1JNO PDB structure are shown in Figure 10. In this Figure, the gA pore radius profiles of all systems are quite similar and show only constrictions (a smaller radius) at the entrance to the channel (minimum radius of 1.0 -1.4 Å) in particular for the gA channels bonded to the AOT RM. In the central region of the pore (defined with respect to z as -10 Å < z < -2 Å and 2 Å < z < 10 Å) the radius profile is mostly flat with an average radius around 1.7 Å, whereas the inter-monomer junction region (defined as -2 Å < z < 2 Å and including the formyl-Val1-Gly2-Ala3 residues) has a slightly wider radius of 1.8 Å.

Water occupancy and configuration inside the gA channel. As shown in the previous section, the radius of each gA pore is broad enough to allow the water to permeate through the channel in a single file configuration (73)(111). To examine the water occupancy and configuration inside the pore in more details, we computed the time evolution of the number of penetrating water molecules in the channel pore by considering that one water is in the pore if the distance between the water oxygen and the peptide backbone (i.e. N, C, O, and CA atoms that lie the pore lumen) are less than 3.5 Å. This cutoff value was chosen by comparing the number of water oxygen given by the gmx select tool of GROMACS (v2018.2) with those counted manually from PyMOL GUI from selected PDB at identical simulation times. These calculations (not shown) reveal that in the β 6.3 helix conformation the average numbers of water molecules that reside near the gA dimer "entire" residues backbone atoms depend on the channel environment and configuration. Averaged over the last 200, for gA_SDS, gA_DMPC, gA_AOT_SA_rep1_2 and gA_AOT_SA_slow systems, we find a number of waters of 18.9 ± 1.0, 15.9 ± 0.2, 16.2 ± 0.3 and 13.3 ± 1.2, respectively. For sake of comparison, the value obtained in the DMPC simulation is close to the value obtained by Basu et al. [START_REF] Basu | Ion Channel Stability of Gramicidin A in Lipid Bilayers: Effect of Hydrophobic Mismatch[END_REF] for gA in the same environment (around 13). By limiting the calculations to the central region of the pore (i.e. where the gA residues are in the β 6.3 helical conformation, i.e., -10 Å < z < 10 Å), we obtained that the gA channel in the SDS/DMPC, gA_AOT_SA and gA_AOT_SA simulations are continually filled during the course of the MD simulations with on average a single file of 6 -9 waters and sometime with one or two Na + ions located at the two channel entrances (Figure 11). For the gA channel in the gA_AOT_PA_fast simulation due to the shorter length of the pore folded region, we counted a water file with, on average, only 6 -7 waters (Figure 11e).

We also notice that when the Na + ions reside near the channel entrances and aligned along the channel axis, they are coordinated with the carbonyl moiety of the Trp11, 13 and 15 or Dleu10 residues of the gA channel, which are considered preferential Na + binding sites of the gA (49) [START_REF] Woolf | The Binding Site of Sodium in the Gramicidin A Channel: Comparison of Molecular Dynamics with Solid-State NMR Data[END_REF].

It is worth to mention that during the course of all the MD simulations, passage of the Na + ion through the gA channel is very rare (less than two/three events in case of the two gA_SDS/DMPC simulations and one/two events in the case gA_AOT_SA and PA simulations). This result may be the consequence of the use, here, of a force field with fixed charges that does tend to overestimate the Na + permeation free energy barrier (64) [START_REF] Peng | Accurate Evaluation of Ion Conductivity of the Gramicidin A Channel Using a Polarizable Force Field without Any Corrections[END_REF].

The distribution of distances between neighboring water oxygen atoms along the pore z axis for the eight channels were also computed (not shown). They show that the first neighbors oxygen distances are almost constant and, on average, equal to 2.9 ± 0.2 Å, which is similar to that found in the bulk (2.8 Å) (112). Moreover, as shown by the snapshots in Figure 11, the water dipole O-H orientation depends on the presence of the Na + in the channel and thus favors formation of different types of hydrogen-bonded (HB) between neighboring water molecules and the porelining carbonyl groups [START_REF] Pomès | Structure and Dynamics of a Proton Wire: A Theoretical Study of H+ Translocation along the Single-File Water Chain in the Gramicidin A Channel[END_REF]. Most of the time, when the water molecules form a single file of water in the pore lumen (Figure 11b), the water molecules are involved in three HBs and adopt the 20 S11). As it can be seen, there are striking differences in water mobility between the SDS/DMPC and AOT simulations. In the RM simulations (Figures 12c andS11), the water inside the central region of the pore remains in the pore and does not leave it even after a long period (> 300 ns) of simulation time and therefore the total number of waters crossing events (i.e. the water flux) is zero in these systems. In that case, the water molecules fluctuate around a well-defined position spaced of 3.0 Å inside the gA channel and may be dictated by the small deformation of the channel of the hydrogen bonds between the lining residues during the course of the simulations, in agreement with other simulations (e.g. ( 57)(113)). This zero-flow can be explained by the absence of motive force in a preferred direction, regardless of the fact that there are RMs of different sizes at each side of the gA channel. Another reason for this lack of water flux are that the two channel entries are sporadically blocked by Na + ions, or the AOT headgroup or deformed by the RMs. In case of the water movement inside the gA channel in the SDS/DMPC systems Figures 12a-b show larger concerted motion during the simulation time period where the water reside in the pore before to exit it. We point out that in the SDS and DMPC simulations, the water channel entry/exit is also reduced due not only to the absence of chemical/electrical gradients in the simulation, but also to the possible presence of DMPC/SDS headgroups (since the gA is slightly smaller than the membrane/micelle thickness (see Figures 3a-b)) or to the presence of ions at each side of the channel entry during some periods of the simulation.

Conclusions

To the best of our knowledge, this work presents the first detailed MD simulation study of the structural properties of a simple, but relevant, model of integral membrane protein, the gA channel, simulated in presence of AOT reverse micelles in pentane with a W0 = 6, this mimicking the experimental conditions of Van Horn et al [START_REF] Van Horn | Use of Reverse Micelles in Membrane Protein Structural Biology[END_REF]. The results of these simulations carried out with multiple replicas obtained from different starting conditions for a total of 2.4 µs were compared with other MDs, where the gA channel was inserted in models of SDS micelle and DMPC membrane. The main findings of this paper are given in the Figures 3c-f. They show that, in agreements with 1 H NMR experiments [START_REF] Van Horn | Use of Reverse Micelles in Membrane Protein Structural Biology[END_REF], the gA homodimer forms a stable water channel in presence of RMs of AOT in pentane as in micelle or membrane. Remarkably, our simulations, performed with different starting conditions and protocols, show that the gA dimer in AOT RMs did not look like to the "bridging model" where the gA dimer interacts only with AOT acyl chains and not with the solvent, as anticipated by experiments [START_REF] Van Horn | Use of Reverse Micelles in Membrane Protein Structural Biology[END_REF]. Indeed, in our simulations the AOT, Na + and water form two well defined elongated RMs bonded to the gA channel with different sizes, shapes and similar W0 values in the range of 9.5 to 19.3 Å, 0.66 to 0.92 and 4.7 to 7.0, respectively. If the RMs are only bonded to the C terminal parts of the gA dimer, the central region of the gA dimer is found in direct contact with the pentane and does not interact with the AOT acyl chains.

Concerning the secondary structure of the gA dimer, our simulations clearly show that the β 6.3 helix conformation of the dimer in RM is stable and well conserved during the course of the simulations as in SDS or DMPC, independently of how the model is initially constructed. Thus, the gA channel forms a water channel with constricted entries and an average central pore radius around 1.7 Å. This pore is filled by a single line of 6 -9 waters, where each water forms an extended HB network with the neighboring waters and the pore-lining carbonyl groups. Finally, the analysis of the water movement through the gA channel shows differences depending on the three environments. For instance, we find that in the RM environment the water molecules do not leave the gA channel even after a long simulation period (t > 300 ns) and that these molecules fluctuate around a well-defined position. This contrasts with our findings for the water inside the pore in the DMPC and SDS simulations, which show highly correlated water movements prior to exit the pore. Notwithstanding, also in these systems, the water flux is reduced probably due to the DMPC/SDS headgroups or ions that during the simulation times can sporadically block the channel entries. Author Information:

List of Figures

List of Tables

Corresponding Author:

  Figure 3d). In case of the gA_AOT_SA_rep1 simulation, the aggregation process is quite similar,

  , the structures of the RMs attached to the gA dimer have different characteristics between the two runs with NAOT, (W0) and 𝑅 𝑔 𝑅𝑀 values in range of 9 -47, (4.7 -6.0) and 9.5 -19.3 Å respectively. The computed semi-axis and eccentricity values confirm that they are mostly ellipsoidal without a clear symmetry of revolution and therefore cannot be characterized as either prolate (aRM > bRM = cRM) or oblate (aRM = bRM > cRM) ellipsoids, with 𝑒 𝑅𝑀 values between 0.63 -0.92.

  )[START_REF] Ingólfsson | Gramicidin A Backbone and Side Chain Dynamics Evaluated by Molecular Dynamics Simulations and Nuclear Magnetic Resonance Experiments. I: Molecular Dynamics Simulations[END_REF], we notice that the gA RMSDBC values are stable, low and, on average, around 0.75 Å for the two replicas during the course of the SDS and DMPC simulations, indicating that the initial dimer conformation of the gA dimer is well conserved (see also discussion below) in the direct micelle and the membrane (Figures 7ab). These results are also consistent with the conservation of the initial β-helix secondary structure (FigureS9a-d). In the case of the gA dimer in AOT RMs, the RMSDBC depicted in the Figure7c-dshow different behaviors as a function of time depending on how the gA_AOT complexes were initially built. For the gA_AOT_SA_rep1 and rep2 systems, the RMSDBC values are found slightly higher (< 1.5 Å) and fluctuate more than the corresponding RMSDBC computed from the SDS and DMPC simulations, reflecting notable short-lived variations of the β-sheet secondary structure, in particular in the middle and C-term of one monomer (see FigureS9e-f). In case of the gA in the gA_AOT_PA_fast simulation (Figure7din indigo color), the change of the gA secondary structure leads to a rapid increase of the RMSDBC value of up to 3.0 Å during the first 30 ns, which remains stable at this value until the end of the simulation. This occurs only in the Cterm region of one dimer, delimited by the residues TRP13-DLeu14-Trp15-Eth and located at the surface of the RM_2 aggregate, whereas the rest of the gA dimer keeps its initial secondary structure until the end of the run (FigureS9h). In case of the gA_AOT_PA_slow the RMSDBC value is found close to the other systems  1.0 Å, thus confirming the stability of the gA secondary structure during the course of the simulation, (Figures7d and

  following organization: Donation of one HB to a channel backbone carbonyl, donation of the other HB to a neighboring water, and acceptance of one HB from the other neighboring water. Such a coordination results in a continuous chain of water-water hydrogen bonds directed along the z axis of the channel. These overall results are consistent, for instance, with the MD results of those, for instance, of Pomes and Roux[START_REF] Pomès | Structure and Dynamics of a Proton Wire: A Theoretical Study of H+ Translocation along the Single-File Water Chain in the Gramicidin A Channel[END_REF] or Kalathingal et al.(113) of a cytotoxic 48-residues channel (polytheonamide B of the marine sponge Theonella swinhoei) that shares same structural features (i.e. presence of D-and L-amino acids and a native β 6.3 -helical conformation) with the gA channel (114).Water movements inside the pore. To complete our analysis of the pore activity, we have also examined the effect of the presence of two RMs at each of the gA dimer extremity on the water movement through the channel and compared with the results obtained for the SDS/DMPC simulations. To do this, we defined a region of L  20 Å and 14 Å for the gA_AOT_PA_fast indicated by two bars in Figures12 and S11, corresponding approximatively to the length of the β folded region of the pore, and track only the oxygen water that traverse this region at different simulation time in the two directions. Representative examples of theses water movements are shown in Figure 12 and were computed from different simulation time windows for the gA_SDS, gA_DMPC and gA_AOT_SA_rep1 simulations (other examples for gA_AOT_PA_rep2 and gA_AOT_PA_fast and slow simulations are shown in Figure

Figure 1 :

 1 Figure 1: Chemical structures of the surfactant and lipid used in this study. (a) bis(2-ethylhexyl) sulfosuccinate sodium salt (AOT), (b) Sodium dodecyl sulfate (SDS), (c) 1,2-dimyristoyl-sn-glycero-3-phosphocholine (DMPC).

Figure 2 :

 2 Figure 2: (a) Amino acid sequence of the three variants of gramicidin (gA) synthetized by Bacillus brevis (a). All the amino acid side chains are either hydrophobic (Ala, Leu and Val) or amphipathic

Figure 3 :

 3 Figure 3: Representative snapshots extracted from the MDs with gA in the (a) SDS micelle (tsim = 230 ns), (b) DMPC membrane (tsim = 230 ns), (c) gA_AOT_rep1 (tsim = 690 ns), (d) gA_AOT_rep2 (tsim = 750 ns), (e) gA_AOT_PA_slow (tsim = 600 ns) and (f) gA_AOT_PA_fast (tsim = 500 ns). The gA dimer is colored according to the nomenclature shown in Figure 1. The AOT and SDS alkyl chain are shown in grey, whereas the oxygen, sulfur, phosphorus and hydrogen atoms are in red, yellow, green and orange spheres, respectively. The sodium and chloride ions are in blue and red colors. Only the water within 5.0 Å of the gA channel are shown in the SDS, DMPC snapshots. Each RM in the Figures c and d are labelled "RM_X" according to their structural characteristics given in Tables2 and 3. The figures were rendered with PyMOL (99).

Figure 4 :

 4 Figure 4: Representative snapshots vs. time of the aggregation process of AOT monomers into reverse micelles in presence of the gA channel extracted from the gA_AOT_rep2 run. gA dimer, AOT are shown in sticks and blue, red and grey colors. Each gA monomer are represented as βsheet and in blue and red colors as in Figure 3. Water and Na + ions are represented as van der Waals spheres. Pentane is not shown for visual clarity. The black line shows the box limits. The figures were rendered with PyMOL (99).

Figure 5 :

 5 Figure 5: Function G(t) for the two gA_AOT trajectories showing the formation of the AOT RMs in presence of the gA dimer. The double-exponential fitting values for each gA_AOT simulations are reported in the inset and drawn with cyan, brown colors, see main text for details. To distinguish from the two results, in the figure, an offset of 0.2 is added to G(t). In the inset, factors of 2.9 multiply the function, respectively. c1, 1 and 2 are the vanning number and time constants (in ns), respectively.

Figure 6 :

 6 Figure 6: Average molecular contacts between the gA residues and the different parts of the (a) SDS, (b) DMPC and (c) the AOT RM systems and the water and pentane solvents. Note that the molecular contacts were averaged out from the contributions of the two monomers and the two replicas of SDS/DMPC and four AOT simulations. The same figures for each gA monomer and simulations are given in Figures S5-8 in the SI. The standard errors in less than 5 % of each type of molecular contacts.

Figure 7 :

 7 Figure 7: Time series of the root mean square deviation for the gA dimer backbone atoms (RMSDBC) for each replica (a) SDS micelle (black/grey), (b) DMPC membrane (blue/cyan), (c) gA_AOT_SA_rep1_2 (red/orange) and (d) gA_AOT_PA_fast/slow simulations (indigo/magenta).

Figure 8 :

 8 Figure 8: Time series of the end-to-end distance, dEndToEnd for the gA dimer and each replica (a) SDS micelle (black/grey), (b) DMPC membrane (blue/cyan), (c) gA_AOT_SA (red/orange) and (d) gA_AOT_PA_fast/slow simulations (indigo/magenta). The dEndToEnd values were computed from the distance between the first carbon atoms of the Ethaloamine group bonded to the last Trp15 residue of each gA monomer.

Figure 9 :

 9 Figure 9: Time series of the rise per turn; drise for the gA dimer and each replica (a) SDS micelle (black/grey), (b) DMPC membrane (blue/cyan), (c) gA_AOT_SA (red/orange) and (d) gA_AOT_PA_fast/slow simulations (indigo/magenta).

Figure 10 :

 10 Figure 10: Average pore radius profiles (a) and the pore lining surfaces (right panel) of the gA channel in the three environments computed from the last 150 and 300 ns of the gA_SDS/DMPC and gA_AOT simulations, respectively. The pore size was evaluated along the pore axis z and are colored in green, black/maroon, blue/cyan, red/orange and indigo/magenta for the gA 1JNO structure, gA_SDS, gA_DMPC, gA_AOT_SA and gA_ AOT_fast and slow simulations. Each pore lining surfaces are colored according to their minimal (in blue) and maximal (in red) values of the pore radius and were rendered with PyMOL (99). The gA dimer are represented with the same color code than in Figure 2.

Figure 11 :

 11 Figure 11: Representative snapshots obtained at the end of each MD of the positions and orientations of the water molecules inside the gA channel in (a) SDS micelle (rep1), (b) DMPC membrane (rep1), (c) gA_AOT_SA_rep1, (d) gA_AOT_SA_rep2, (e) gA_AOT_PA_fast and (f) gA_AOT_PA_slow simulations. The pore backbone atoms are in sticks with the Trp9, 11 (hidden), 13 and 15 highlighted in yellow. Na + ions are in van-der-Waals blue sphere. The black dots indicate hydrogen bonds between neighboring waters or the carbonyl oxygen atoms of the pore lining residues. The Figures were rendered with PyMOL (99).

Figure 12 :

 12 Figure 12: Representative picture of the spontaneous travels with time of the water oxygens within along the z-axis of the channel for (a) gA_SDS, (b) gA_DMPC and (c) gA_AOT_SA_rep1 simulations. Note that for the gA_SDS, we only show two water molecules and consider only the first 50 ns of the MD simulation for visual clarity. In case of the gA_DMPC we show only the first 4ns whereas of the gA_AOT RM_rep1, the last 300 ns of the MD when the gA_AOT form a stable complex (see main text). The calculations were done every 100 ps. The calculations for the two others gA_SDS and gA_DMPC replicas were similar and shown in the Figures. The water channel colored according to their relative positions in the pore are shown in the left of each figure. The two bars indicate the 20 Å long region in which water movement is analyzed and corresponds approximatively to the length of the folded region of the gA channel in these simulations. The same plots for the other gA_AOT systems are given in Figure S11 in the SI. The Figures were rendered with PyMOL (99).

Stéphane Abel -Figure S1 :

 AbelS1 Figure S1: Molecular snapshot of the gA bridging two RM particles construct with the PA protocol. Each RM contains 38 AOT (and Na + ions) and 225 waters (see main text and Figure 3) for details.

Figure S2 :

 S2 Figure S2: Representative snapshots vs. time of the aggregation process of AOT monomers into a reverse micelle in presence of the gA channel extracted from the gA_AOT_SA rep2 run. gA dimer, AOT are shown in sticks and blue, red and grey colors. Each gA monomer are represented as β-helix and in blue and red colors as in Figure 3 in the main text. Water and Na + ions are in van-der Waals spheres. Pentane is not shown for visual clarity. The black line shows the box limits. The figures were drawn with PyMOL (1).

Figure S3 :

 S3 Figure S3: Same legend as Figure S2 for the gA_AOT_PA_fast system constructed with the preassembled protocol. The figures were drawn with PyMOL (1).

Figure S4 :

 S4 Figure S4: Same legend as Figure S3 for the gA_AOT_PA_slow system. The figures were drawn with PyMOL (1).

Figure S5 :

 S5 Figure S5: Average molecular contacts between the gA residues, the SDS parts and water for the two gA_SDS systems.

Figure S6 :

 S6 Figure S6: Average molecular contacts between the gA residues and the DMPC parts and water for the two gA_DMPC systems.

Figure S7 :

 S7 Figure S7: Average molecular contacts between the gA residues, the AOT parts headgroup, water and pentane for the two gA_AOT_SA systems.

Figure S8 :

 S8 Figure S8: Average molecular contacts between the gA residues, the AOT parts headgroup, water and pentane for the gA_AOT_PA_fast and slow systems.

Figure S9 :

 S9 Figure S9: Time evolution of the gA dimer secondary structure obtained in different simulations. (a-b) In two gA_SDS, (c-d) gA_DMPC, gA_AOT_SA_rep1_2 (e-f) and gA_AOT_PA_fast and slow systems (g-h), respectively. The final structure of gA dimer are depicted according to the color convention given in Figure 1 of the article are shown in right side of each figure. The black box highlights the region permanently unfolded region of the gA monomer in the gA_AOT_PA_fast simulation. See main text for details.

Figure S10 :

 S10 Figure S10: Snapshots of the gA dimer with colored beads showing the distance between the different amide pairs (i.e. DVal1A↔DVal6B (red), Gly2A↔Ala5B (blue), DAla3A↔DLeu4B (green), DLeu4A↔Ala3B (yellow), Ala5A↔Gly2B (magenta) and DVal6A↔Val1B (cyan) in black dashed lines. The figures were drawn with PyMOL (1).

Figure S11 :

 S11 Figure S11: Representative picture of the spontaneous travels with time of the water oxygens within along the z-axis of the channel for (a) gA_AOT_SA_rep2 (b) gA_AOT_PA_fast and (c) gA_AOT_PA_slow simulations. We only considered in the calculations and the plotting the time where the gA_AOT RMs forms a stable complex (see main text). The representation of gA channel dimer with the channel water colored according to their relative positions in the pore are shown in the left of each figure. In panels (a-c) two bars indicate the 20 Å long region in which water movement is analyzed and corresponds approximatively to the length of the folded region of the gA channel in these simulations. This distance is reduced to 14 Å for gA_AOT_PA_slow system (panel b) due to the partial folding of the gA dimer.

  Herein, we carried out unconstrained explicit MD simulations of RMs with the gA channel, which is a simple and at the same time a relevant model of IMP and is also well characterized experi-

	mentally	(43)(47)(48)(49)(50)(51)(52)(53)(54)(55)(56)	and	theoretically

(48)(57)(58)(59)(60)(61)(62)(63)(64)(65)(66)(67)(68)(69)(70)(71)(72)(73) in various lipidic environments. gA is formed with two identical hydrophobic peptides of 15 amino acids each with alternative -D or -L conformations (Figure

Table 1 :

 1 Summary of the simulated system. The system size is given by the lengths of the simulation box vectors LX × LY × LZ.1 These simulations were performed twice with different random seeds.2 In each membrane leaflet.3 These RM systems were built with the SA and PA protocols and simulated twice with different random seeds and different equilibration protocols.

	System	Lipid	NL	NH2O Na + /Cl -natoms	Systems size (in Å)	tsim (ns)
	gA_SDS 1	SDS	65	15000	65/0	48347	77.9 x 77.9 x 77.9	230/260
	gA_DMPC 1	DMPC 90 2	10000	27/27	51846	72.2 x 72.2 x 95.0	230/250
	gA_AOT_SA 3 AOT	75	450	75/0	151352 119.7 x 119.7 x 119.7 695/750
	gA_AOT_PA 3 AOT	76	450	76/0	151418 119.7 x 119.7 x 119.7 500/600

Table 2 :

 2 List of the structural parameters of the simulated RMs at the beginning (tsim = 0 ns) and the end of the two gA_AOT simulations. NAOT, W0 and 𝑅 𝑔 𝑅𝑀 (in Å) are the numbers of AOT per RM, the water-to-surfactant ratio and the radius of gyration, respectively. The values separated by "/" are for different replicas. The values in bold are for the RMs bonded to the gA dimer. These latter values were computed with the OpenTRJ code (103) at the end of each simulation. See alsoTable 3 for the computed shape parameters of each RM.

	System	gA_AOT_SA	gA_AOT_PA_fast/slow
	Nc(0)	43/43	38/38
	Nc (∞)	3/2	2/2
	NAOT (∞)	37; 29; 9/47; 28	44; 32/38;38
	W0(∞)	6.7; 5.1; 4.7/6.0; 5.6	6.0; 5.7/5.9;5.7
	𝑅 𝑔 𝑅𝑀 (∞)	14.8;

13.7; 9.5/19.3; 13.4 18.1; 14.3/16.6;14.9Table 3 :

 3 Shape parameters of the whole RM (including AOT, water and ions). aRM, bRM, cRM and e are the lengths of semi-axis (in Å) computed from the whole RM moment of inertia and the eccentricity, respectively. As in Table 2, the values in bold are for the RMs bonded to each gA monomer. Clark, G. N. I.; Cappa, C. D.; Smith, J. D.; Saykally, R. J.; Head-Gordon, T. The Structure of Ambient Water. Mol. Phys. 2010, 108, 1415-1433. (113) Kalathingal, M.; Sumikama, T.; Mori, T.; Oiki, S.; Saito, S. Structure and Dynamics of Solvent Molecules inside the Polytheonamide B Channel in Different Environments: A Molecular Dynamics Study. Phys. Chem. Chem. Phys. 2018, 20, 3334-3348. (114) Hamada, T.; Matsunaga, S.; Fujiwara, M.; Fujita, K.; Hirota, H.; Schmucki, R.; Güntert, P.; Fusetani, N. Solution Structure of Polytheonamide B, a Highly Cytotoxic Nonribosomal Polypeptide from Marine Sponge. J. Am. Chem. Soc. 2010, 132, 12941-12945. The Supporting Information is available free of charge on the ACS Publications website at DOI  Additional Figures and Tables (PDF)

		System Biol. 1981, 59, 155-171.		gA_AOT_SA_rep1
	(112)	RM	aRM	bRM	cRM	𝑒 𝑅𝑀
		RM_1	16.6	18.7	22.0	0.66
		RM_2	13.4	18.3	20.5	0.76
		RM_3	10.4	12.7	13.4	0.63
		System		gA_AOT_SA_rep2
		RM	aRM	bRM	cRM	𝑒 𝑅𝑀
		RM_1	14.3	17.0	37.1	0.92
		RM_2	15.1	16.4	20.0	0.66
		System		gA_AOT_PA_fast
		RM	aRM	bRM	cRM	𝑒 𝑅𝑀
		RM_1	14.9	16.2	34.0	0.90
		RM_2	14.4	17.3	22.7	0.77
		System		gA_AOT_PA_slow
		RM	aRM	bRM	cRM	𝑒 𝑅𝑀
		RM_1	14.4	16.0	30.2	0.88
		RM_2	15.7	18.2	23.0	0.73

Acknowledgment

SA acknowledges Dr. Rajendra Kumar for his kind help for computing the channel radius with its gmx_clusterByFeatures program. This work was granted for access to the HPC resources of [CINES/IDRIS/TGCC] under the allocation 2018-A0040710413 made by GENCI. We also be

10000

Table S1: Proposed twenty steps (or "slow") equilibration protocol for the gA_AOT system constructed with PA approach and depicted in the Figure S1 after the minimization stage.

X_PR, t, tequ are the position restraints applied to the X group, the time step and the duration of the equilibration of period, respectively. (++), (+) and (-) are for the position restraint forces of 5000, 1000, and 0 kj.mol -1 , respectively.

III. References


Peter Reinhardt
email: peter.reinhardt@sorbonne-universite.fr

Ilya Popov

Tchougréeff Minimum

Andrei L Tchougrée

Minimum Atomic Parameter basis sets for elements 1 to 54 in a Hartree-Fock setting

come

Introduction. Why a new set of Slater-like functions to describe atoms?

Quantum chemical modeling starts from nite sets of basis functions so that quality and eciency of electronic-structure calculations depend on the properties of the latter. This is emphasized in almost every handbook of quantum chemistry (see e.g. Ref. 1) and manifests itself in a wide variety of basis sets used for dierent objects and their properties. 24 Cur- rently, the dominating motivation of quantum-chemical modeling is the accurate numerical reproduction of dierent properties by using heavy computations. This approach answers questions like What numerical value has a quantity of interest? However, a purely numerical approach cannot provide any explanation. A rarely recognized deciency of this approach is that the experimental data to be numerically reproduced are not available with sucient accuracy or even denition.

5 An alternative to this numerical approach which almost disappeared from the literature aims at a conceptual explanation of experimental trends by using semi-quantitative models. Although this second way sacrices the numerical accuracy to some extent, these approaches are to be considered as complementary to the former. They, however, need as well basis sets as a starting point, although the requirements dier, as formulated for instance in Ref. 6. In fact, only the following seems to be mandatory for the qualitative or at best semi-quantitative approach: Each basis function should reect the nature of the problem, such that a good representation of the orbitals can be achieved by a limited (small) number of functions.

6 If this requirement is satised in a consistent way, such that the small number of functions and, thus, of the parameters dening them, is really physically substantiated, then one can hope to be able to fulll another desideratum, underlined in Ref. 6, namely to provide basis sets for a good fraction of the Periodic Table ,   and additionally have a regular and sensible dependence of the parameters and basis orbitals on the only true parameter the atomic number, i.e. the nuclear charge Z. [START_REF] Van Den Broek | Das α-Teilchen und das periodische System der Elemente[END_REF][START_REF] Moseley | The High-Frequency Spectra of the Elements. Part II[END_REF] This is what we would like to achieve in the present contribution.

The paper is organized as follows. In Section I we provide the theoretical scheme of constructing minimally parameterized sets of atomic orbitals with one single exponent for each n subshell while keeping the correct nodal structure. Next, in Section II, we describe the results of solving the energy-optimization problem for atoms up to Z = 54 (Xe) in a restricted-open-shell Hartree-Fock setting, and discuss the obtained results in Section III, in order to draw conclusions (Section IV. A compilation of the data can be found in the accompanying Supplementary Material, and all purely technical details are collected in the Appendix.

I. THEORY

A. The construction procedure Minimal Slater-based basis sets of the form

ϕ n m ( r) = N n P n (2ξ n r)e -ξ n r × Y m (θ, φ) (1)
are dened so that the correct nodal structure of the orbitals (orthogonality between dierent functions of same ) is maintained by the radial polynomial P n (2ξ n r) the form originally used by Fock and coworkers 9 for the Be atom and for closed 2s

2 2p
6 shell ions (like Na + )

and then independently rediscovered and employed for the entire 2 nd period 10 leaving the exponents ξ n to be the only free parameters of the basis, as many as subshells. Writing the normalized radial functions as:

R 1s (r) = N 1s e -ξ 1s r R 2s (r) = N 2s (1 + a 2s (2ξ 2s r)) e -ξ 2s r R 3s (r) = N 3s 1 + a 3s (2ξ 3s r) + b 3s (2ξ 3s r) 2 e -ξ 3s r R 4s (r) = N 4s 1 + a 4s (2ξ 4s r) + b 4s (2ξ 4s r) 2 + +c 4s (2ξ 4s r) 3 e -ξ 4s r . . . (2)
(and similarly for higher values of the angular momentum) and evaluating the general overlap integral

ˆ∞ 0 r n e -ξ 1 r e -ξ 2 r r 2 dr = (n + 2)! (ξ 1 + ξ 2 ) n+3 (3)
we arrive at an iterative scheme as a 2s orbital should be orthogonal to 1s, 3s orthogonal to both 2s and 1s etc, leading together with the normalization to nconditions for the ncoecients of the polynomial P n (2ξ n r). The coecients depend parametrically on the ensemble of ξ's of the functions of lower n. For instance for the 4s orbital we obtain 3 equations for the 3 coecients a 4s , b 4s and c 4s , in a systematic and generalizable manner using reduced exponents k ij = ξ i ξ i +ξ j with one of the indices i and j being s4 and the other running from 1s to 3s.

mize the total energy with respect to them. Technically, the total energy for an (open-shell) atom is written, using a non-relativistic Hamiltonian for the atomic system with one or several open shells.

1113 As no orbital expansion coecients have to be determined by mini- mization (they come immediately from the orthonormality condition), one single evaluation of E Ψ = Ψ| Ĥ|Ψ is needed for each set of exponents. The determination of the exponents through the minimization of the total energy is then left to a simplex procedure 14 since the gradients used originally 10 become too cumbersome for the elements beyond the 2 nd period. Note that the present procedure is indeed a Hartree-Fock scheme, optimizing the total energy with respect to the parameters of a single-determinant wavefunction, however without the construction and diagonalization of a Fock matrix of the commonly employed Roothaan-Hall procedure.

II. RESULTS

A. Energy minimization and orbital exponents

For the elements Z = 1 -54 (H Xe) we optimized the exponents for atomic ground states as given by the Aufbau principle, without paying attention to exceptions like Cr, Cu, and most of the 4d transition elements. For the same elements we have from the literature the highly optimized multi-exponent Bunge basis sets, 15 which, in terms of total energy, are very close to numerical Hartree-Fock energies. 16 Our basis sets yield total energies slightly higher than the Hartree-Fock limit, about 0.30.7 %, which is an amazingly small dierence, given the number of free parameters in our calculations. Set aside the proposition of Fock et al. 9 , closest to our approach and covering a large part of the periodic table, are the single- zeta basis sets of Clementi and Raimondi, 17 with the same number of exponents as ours, but constructed as optimized linear combinations of nodeless Slater functions. The number of free parameters is thus about twice as large as ours (exponents and expansion coecients with orthonormality constraints), and the dierence to total energies obtained with Bunge's basis sets about a factor two reduced with respect to ours (see Figure 1).

Overlap integrals between core orbitals in the Bunge sets and respective MAP orbitals are found to be larger than 0.99 (> 0.999 for the 1s shell), and still larger than 0.90 for valence orbitals (for more details see below). One exception is the Pd atom for which the basis set of Bunge et al had been determined using the physical 4d 10 5s 0 ground state without the need to describe a 5s orbital.

Remarkably, in most cases the exponents of valence s-, pand d-orbitals, determined from the energy minimization procedure, follow, as functions of the atomic number, roughly linear trends as shown in Figure 2, corroborating the semi-empirical argument given in the Introduction.

A closer look reveals that despite almost perfect linear trends particularly for smaller n the deviations from the trends of Figure 2 are not statistical, but systematic. Furthermore, the dependency of orbital exponents with given (n ) on the nuclear charge is piecewise linear, clearly visible for instance for the 4s exponents in the three groups 1920 (K and Ca, 4s as valence), 2130 (ScZn, lling of the 3d inner shell) and 3154 (4s being outside of completely lled shells). In Table I we t the dierent regions identied for all orbitals to the linear expression ξ(Z) = a Z + b. The corresponding linear trends are shown in Fig. 3 for 4s and 5s-exponents.

These linear ts, obtained for small Z, may then serve as starting points for the optimization procedure for higher Z, being as reduced as Slater's screening model from 1930, 18 I.

setting. It also comes up in the studies on ab initio Gaussian atomic basis sets [START_REF] Jensen | Atomic orbital basis sets[END_REF][START_REF] Nagy | Basis sets in quantum chemistry[END_REF][START_REF] Jensen | Unifying General and Segmented Contracted Basis Sets. Segmented Polarization Consistent Basis Sets[END_REF] and may be a general property of all basis-set-expanded atomic orbitals representing, more or less realistic, the atomic one-electron states.

Additionally, we see that the positions of the nodes with the same number, as counted from the nucleus as functions of Z, follow fairly smooth curves, at least for 2s to 4s. They appear in groups, almost coinciding for the respective nodes, starting at the Z values at which further nodes for same and dierent n appear.

We observe an evident cusp for the 3rd root of the 4s polynomial between Z = 30 and Z = 31, at the end of the rst transition series. The 4s-exponent is almost constant in this range, but orthogonality has to be achieved with evolving 1s, 2s, and 3s exponents.

Beyond Z = 31 the 4s exponent varies as well, which adds a dierent component to the orthogonality constraint, without successive llings of inner shells. The same is observed for the 5s at Z = 39 (Y) and Z = 48 (Cd) for the lling of the 4d shell.

III. DISCUSSION

When developing the atomic Hartree-Fock theory in the 1930-1960's (Slater, Roothaan, Clementi) the focus was then on molecular applications, where a purely spherical atom is 23 Our minimally parametrized basis sets are not intended to replace the modern ones designed for molecular applications with numerical precision. They, however, may serve as a set of physical parameters for more qualitative or semi-quantitative studies aimed at establishing and explaining trends in the chemical behavior.

A. Linear trends for the exponents and deviations from them

We may start the analysis of the linear relations described in section II in terms of hydrogen-like orbitals (ξ = Z/n) or Slater's parameter set [START_REF] Slater | Atomic Shielding Constants[END_REF] ξ

(n ) S = Z * n * = (Z -σ n ) n * , ( 4
)
where σ is composed of the screening by the equivalent and inner electrons, and n * is an eective principal number. As result one obtains as well piecewise linear relations of the As one can see from Table I, the slopes in the core ranges deviate very little from the 1 /n rule where n is the true principal quantum number the one determining the number of radial nodes of the function. On the other hand, in the valence ranges the slopes for nontransition elements are systematically smaller than 1 /n most signicantly for n = 2, less pronounced for n = 3; in either case, however, not precisely following the Slater prescription of 0.65/n (which almost exactly shows up for 2p valence orbitals only). At n = 4 the situation inverts and the slopes of the linear ts of the MAP exponents as functions of Z remarkably exceed 1 /n in the valence range (named non-transition in Table I). This feature is partially absorbed by the eective principal quantum numbers, n * < n, introduced by Slater. However, using n * instead of n still produces very small slopes in the valence range:

the estimate of 0.65/n * based on Slater rules is smaller than the MAP ts in Table I, for which the ratio of slopes between core and valence trends is more like 1 than 0.65.

The linear relations for valence ns exponents of transition elements show very small slopes (0.03 and 0.02, respectively for n = 3 and 4), ultimately meaning that in those series the eld felt by s-electrons does not change signicantly with the change of the nuclear charge as expected from Slater's rules ascribing to each electron in the 3d-shell the screening power of 0.85 relative to the 4s electrons. This results in the slopes of 0.15/n * (0.05 and 0.04, respectively) for the latter exponents as functions of Z. This weak increase correlates to the experimental ionization potentials of the 3d transition-metal series, spanning a narrow interval of 1.2 eV.

B. Preliminary physics and chemistry implications

As we mentioned in the Introduction, the proposed MAP setting 10 for the atomic orbital basis in its primitive form dating back to 1930-ies 9 has the advantage that its parameters orbital exponents have a direct physical signicance. Thus, one can hope to relate them more or less directly to some observable characteristics of atoms and eventually those of molecules and crystals/materials. For instance atomic radii may be extracted. Experimentally, these are generally derived by tting some data on interatomic separations observed in various contexts.

The MAP exponents provide already the size of an atom due to the asymptotically exponential decay of the electron density.

24 A simple measure for an atomic radius is the average radial position r of the most diuse orbital. This quantity can be evaluated in any basis set, and we can compare directly orbitals of the Hartree-Fock limit to the MAP orbitals characterized by one single exponent. Figure 6 gives these data for all elements up to Xenon.

The coincidence is again quite striking, atomic radii obtained with the two orbital sets For the MAP orbitals we plot as well estimates of the atomic radii as given by eq. 6 (red bullets), see text; the respective interpolating lines are eyeguides. are nearly superposable, but for the 3d and 4d transition metals.

Within the MAP philosophy we may develop another way to estimate the atomic size via the the position of the outermost density maximum of its valence shell. For either value of the highest power of the P n polynomial is n -1 so that the dominant contribution to the outermost wave is a simple Slater monomial of the same degree. We can check by tting a Slater monomial to our MAP orbitals, by maximizing the overlap integral between the two radial functions, dependent on the Slater exponent. We see (Figure 7 the outermost maximum of the radial density r 2 R 2 n (r) is given by the expression

r max (Z) = n ξ n = n 2 Z -σ n (Z) (5)
where n can be safely understood as the element's period (row) in the Periodic Table .   We add to this value a spherical layer of the thickness of ξ -1 n so that ultimately the intrinsic atomic radius reads

r at (Z) = n + 1 ξ n (6)
where n refers to the outermost subshell, in particular 4s or 5s for the transition elements.

The resultant data are added to Figure 6. As one can see, this simple estimation generally follows the same trend as the estimates based on r .

C. Nodes of the radial MAP functions

Let us come back to the result of section II B, by looking more in detail at the polynomials P n (2ξ n r) which maintain the orthogonality between dierent functions of the same .

From the hydrogen atom we have the orthogonality relations involving Laguerre polynomials as

N n N n ˆ∞ 0 r 2 +2 e -2Zr( 1 n + 1 n ) × × L 2 +1 n--1 2Zr n L 2 +1 n --1 2Zr n dr = δ n,n (7)
as the hydrogenic orbital is written as

R n (r) = N n r L 2 +1 n--1 2Zr n e -Zr /n . ( 8
)
Setting ξ n, = Z/n (or a screened nuclear charge Z * ) we arrive at a generalized orthogonality constraint within the MAP series, provided that the exponents for a given angular quantum number derive from one single Z (or Z * ):

N n N n ˆ∞ 0 r 2 +2 e -(ξ n, +ξ n , ) r × × P n, (2ξ n r)P n , (2ξ n r) dr = δ n,n (9)
This is only approximately the case for the Hartree-Fock derived MAP radial parts we nd ξ n, (Z) ≈ Z /n + b n with a b n nearly independent of Z, see Table I. The higher Z, the smaller the eect of b n , which we can formulate as

Z → ∞ : P n (2ξ n r) ∝ L 2 +1 n--1 2 Z n r .
However, ignoring the screening (i.e. b n ) comes back to studying a bare, hydrogenoid, atom only, ignoring the poly-electronic situation to describe. Better seems to introduce an average screening for all shells through an eective Z * , averaged over all n, leading to the average setting

P n (2ξ n r) -→ L 2 +1 n--1 2 Z * n r .
As a consequence, the nodes of the MAP polynomials coincide on average with the corresponding nodes of the associated Laguerre polynomials. The Laguerre polynomials in their turn can be approximated through Bessel functions of the rst kind:

25,26

P n (2ξ n r) ∝ L 2 +1 n--1 2Z * r n ∝ J 2 +1 2 2Z * r 2Z * r 2 +1
.

From this we can immediately read the position of the roots of the MAP polynomials:

r n (k) ≈ j 2 2 +1,k 8 Z * , (11)
where j 2 +1,k is root number k of the Bessel function J 2 +1 of the rst kind, independent of n and inversely proportional to the nuclear charge, corroborating the graphs of Figure 4.

The Bohr from the origin, respectively. We nd 0.0431757, 0.168321 0.429621, 0.989825 Bohr as nodes of the 5s orbital which is close for the rst node, but more and more o for the subsequent ones. The expected positions are reported in Figure 5.

Of course, the reason for this result is that for the inner shells the screening is fairly correct, and as well the approximation of the Laguerre polynomial via a Bessel function.

For the outer nodes, screening via the chosen Z * becomes quite wrong and the approximation through the Bessel function less reliable.

IV. CONCLUSIONS

Minimally parameterized exponential basis sets of the form eq. ( 1) originating yet from 1930-ies 9 are dened so that the correct nodal structure (the orthonormality) of the orbitals is maintained through a polynomial depending only on the exponential parameters. Within this setting each (sub)shell can be described by only one single parameter -its orbital exponent. Assuming this we recently explored 10 atoms of the 2 nd row of the Periodic Table ,   and extended the approach up to (Z = 54), i.e. including the 5-th period.

The (sub)shell-specic orbital exponents ξ

(n ) MAP are derived through optimizing the total energy written in the Hartree-Fock approximation. Overlaps of the obtained orbitals with the corresponding ones from purely numerical approaches or much longer expansions 15 are quite large (> 90 %) and may serve as basis for atomic projection and subsequent analysis of the results of PAW/DFT calculations.

28 The details of the spatial structure of atomic states, like position of the nodes close to the nucleus or the asymptotic behavior for large r are fairly well reproduced.

Despite the restricted number of free parameters the proposed basis sets result in good total energies, and exponents follow simple piecewise linear relations as functions of the nuclear charge Z. The switches between dierent linear trends occur either at the end of rows of the Periodic System or when within a row a new (sub)shell starts being lled.

FIG. 1 :

 1 FIG. 1: Total energies (left scale) and dierence to Bunge's large basis sets in green and Clementi and Raimondi's 17 single-ζ basis sets in light blue (right scale).

FIG. 2 :

 2 FIG.2:The optimal exponents as functions of the atomic number Z. The visually almost strait lines coming up to surface up to 3d are in fact misleading: see the text for details.

FIG. 3 :

 3 FIG.3: Piecewise linear trends for 4s (left) and 5s (right) orbital exponents. Yttrium (Z=39) seems an outlier, probably not arrived at the optimum set of exponents. Parameters of all linear ts are assembled in TableI.

FIG. 4 :

 4 FIG.4: Position of the nodes of the radial functions 2s to 5p. Multiple lines of the same color refer to the subsequent nodes counted from the origin, r = 0, of the same orbital.

11 FIG. 5 :

 115 FIG.5: Plot of the dierent radial functions of s symmetry for a Iodine atom (Z = 53). We clearly see the close coincidence of the rst and the second orbital node. Black stars correspond to the estimates described by eq. (11), see Section III C.

FIG. 6 :

 6 FIG.6: Atomic radii obtained from Bunge's (violet upright crosses) and the MAP orbitals (green skew crosses), dened as the expectation value r of the outermost (partially) lled orbital or subshell. For the MAP orbitals we plot as well estimates of the atomic radii as given by eq. 6 (red

FIG. 7 :

 7 FIG. 7: Dierence between MAP exponents and Slater exponents (upper panel), and radial expectation values (lower panel), as function of Z.

TABLE I :

 I Piecewise linear ts ξ(Z) = a Z + b for s, p, and d exponents together with R 2 values.For the 4s and 5s orbitals and the rst two elements the character is in quotes as the same orbital is still a valence orbital in the corresponding transition-element series.

	Orbital	character	Z-range a	b	R 2
	1s	core	354	0.9999 -0.28 1.0000
	2s	valence	34	0.4017 -0.44
	2s	core	554	0.5046 -1.41 0.9999
	2p	valence	510	0.3319 -0.45 0.9984
	2p	core	1154 0.5046 -1.41 0.9999
	3s	valence	1112 0.3282 -2.63
	3s	core	1354 0.3354 -2.92 0.9985
	3p	valence	1318 0.2878 -2.52 0.9986
	3p	core	1954 0.3326 -3.47 0.9980
	4s	valence	1920 0.3135 -5.01
	3d	valence	2130 0.2616 -3.23 0.9890
	3d	core	3154 0.3669 -6.42 0.9999
	4s	transition	2130 0.0311 +0.69 0.9959
	4s non-transition 3136 0.3074 -7.53 0.9937
	4s	core	3754 0.2571 -5.50 0.9996
	4p	valence	3136 0.2778 -7.09 0.9947
	4p	core	3754 0.2580 -6.17 0.9996
	4d	valence	3948 0.2756 -8.54 0.9900
	5s	valence	3738 0.2812 -9.37
	5s	transition	3948 0.0207 +0.68 0.9990
	5s non-transition 4954 0.2686 -11.05 0.9973

  above estimate basically repeats the old observation based on the WKB treatment of the hydrogen atom (see Ref. 27). For an Iodine atom (Z = 53) we have the s exponents ξ 1s to ξ 5s as 52.6862, 25.2486, 14.8828, 8.1798, and 3.1897, respectively, instead of the hydrogenic Z/n series 53, 26.5, 17.666, 13.25, 10.6. Minimizing the quadratic deviation from the MAP exponents for a common Z * leads to a value of Z * = 49.845, or an average screening of a

	little more than 3 electrons.
	Via the roots of the Bessel functions j 2 1,k /8 amounts to 1.8352, 6.1523, 12.9375, 22.1901
	and thus expected positions of the nodes are 0.036818, 0.123428, 0.259552, and 0.445178

Acknowledgments

This work has been conceived during the stay of P.R. and A.L.T. in Aachen at R. Dronskowski's group whose kind hospitality is gratefully acknowledged. Calculations have been performed mainly in Paris, with the support of the research federation IP2CT to which P.R.

Implementation

From a given set of exponents, the linear equations to obtain the coecients of the polynomials are solved, and a corresponding input for Clementi's atomic Hartree-Fock program 29 is produced. The procedure is easily embedded into a simplex algorithm, as described in standard textbooks of numerical optimization problems.