Patrick Trouillas
email: patrick.trouillas@unilim.fr

Gabin Fabre

Srinivas Banala

Roderich Süssmuth

Keywords: Force field, lipid bilayer, membrane, partitioning, DMPC, molecular dynamics, COSMOmic Antioxidant, Lipid bilayer membrane, Molecular dynamics, Lipid peroxidation polyphenols, antioxidants, membrane, molecular dynamics Germany

Lipid peroxidation is a major deleterious effect caused by oxidative stress. It is involved in various diseases such as atherosclerosis, rheumatoid arthritis and neurodegenerative diseases. In order to inhibit lipid peroxidation, antioxidants must efficiently scavenge free radicals and penetrate inside biological membranes. Lipocarbazole has recently been shown to be a powerful antioxidant in solution. Here, we show its powerful capacity as lipid peroxidation inhibitor. Its mechanism of action is rationalized based on molecular dynamics simulations on a biomembrane model, quantum calculations and experimental evaluation. The role of the lipocarbazole side chain is particularly highlighted as a critical chemical feature responsible for its antioxidant activity.

My first and biggest thanks are addressed to Patrick Trouillas. Patrick, I do not think I can express how lucky I was to meet you and to have you as my advisor and mentor. During all these years, you arranged the perfect stimulating conditions for this thesis. You taught me everything you know, and it has been a real pleasure to work by your side every day. Thank you for your academic and scientific qualities, for your help, your patience, and your generosity. I will always be grateful for your support and for the sacrifices you made to help me.

I would like to express my sincere gratitude to Michal Otyepka for cosupervising this thesis. I truly appreciated your hearty welcome in your lab, your city and your country. Thank you for all the stimulating and fruitful discussions that pushed me further, and for giving me the opportunity to obtain this double degree. I sincerely wish that we continue collaborating in the future.

My next thanks go to all members of the jury, starting with Martin Hof. I am truly honored that you accepted to be part of this jury. I really appreciate your kindness, your expertise and your insightful discussions.

I am much obliged to Claire Rossi, twice. First, thank you for accepting to review this work. Second, thank you for this very efficient collaboration that really meant a lot to me. It has been a real pleasure to meet you, and I hope we will again have the opportunity to work together.

Frédéric Bonté, it is an honor to have you as a member of my jury. It is always stimulating to open new perspectives with industrial groups and to witness the applications of our fundamental research.

Bertrand Liagre, I am grateful that you accepted to be in my jury. Thank you for supporting me all along my studies and for teaching me academic and scientific rigor.

I would like to thank Pavel Banáš and Luca Muccioli for evaluating this work. Your expertise in MD simulations always brought interesting discussions. Thank you both for your constant help during my studies.

My next thanks are addressed to Karel Berka for your huge help on all aspects of my research. I really appreciate your availability to answer my questions, the new ispiring ideas and projects that you brought that contributed to this work. Also, on the personnal side it is a pleasure to know you as a friend. Thank you for your warm welcome inside and outside the lab. I wish you all the best for the future.

I thank all the students, post-docs and senior researchers of the department of physical chemistry in Olomouc: Veronika, Martin and Martin, Vojtěch, František, Mikulaš, Vašek, Petr and Petr, Marie, Petra, Žofie, Sylva and all the others. Special thanks to Marketa for the stimulating discussions and collaborations. It is a real pleasure to know you and to work with you. Děkuji mockrát.

Next I thank my fellow theoretical chemists in Limoges for your help and the friendly environment during my Ph.D studies: Benjamin, Tahani, Imene, Thomas, Nicolas, Michal and Michal, and people from all around the world who I had the chance to meet in the lab. In particular, I would like to express my sincere gratitude to my friend Florent. Thank you for enlightening me to theoretical chemistry and for all your support these last years. I wish you the best for the bright scientific carrier that you deserve.

My next thanks go to the members of my research team. Vincent Sol, I am grateful that you encouraged me to work on this research topic and that you allowed me to teach. Stéphanie Lhez, thank you for making my teaching experience a pleasure. I also thank all the students for the nice moments we shared: Olivier, Mark Arthur, Chloé, Amandine, Benjamin, Shihong, Jean-Pierre, Florent, Manu, and all the others. I would also like to thank the members of the U850 team: Benjamin, Bastien, Sébastien, Danko, Khadija and all the others. Thank you for the friendly atmosphere and the fun we had. Thanks to you, it is always a pleasure to go to work. Of course, special thanks to my bro Pierre-André for these nine years together. Here is not enough space to thank you for all these nice moments, but you know what I mean! I wish you the best for your personnal and professional lives.

My sincere thanks to all collaborators from other cities and countries that helped me build this research topic: Stefan Knippenberg, Burkhard Bechinger, Roderich Süssmuth and the people I already mentioned above. Thanks for all the insightful and motivating discussions.

Next, I would like to thank my family for supporting me. This is it, after all these years I am finally leaving school! Last but not least, a huge thank to my wife and the love of my life. Gaelle, thank you for being immensely kind and generous and for supporting me. I am really lucky to have you. Figure 19: Experimental partition coefficients plotted against the respective calculated values (upper panel) with parameters of the linear fit, i.e. coefficient of determination, R 2 , slope (a) (standard deviation in bracket) and intercept (b). Slopes significantly differing from 1.0 and intercepts from 0.0 significantly on the probability level of 0.975 are highlighted in red. Each data point is labeled by a number, which corresponds to the number of the molecule in Table 3. The fitting parameters for GROMOS 43A1-S3 (G43A1-S3) and GAFFlipids re-calculated by omitting outliers (acetone and 2-nitrotoluene, in GAFFlipids and GROMOS 43A1-S3, respectively) are shown in blue. The bar charts (lower panel) depict the mean differences and the mean absolute differences Table 1: 4: Simulation parameters. R coulomb is a short-range electrostatic cut-off, long-range electrostatics are evaluated by PME, R vdw is Lennard-Jones cut-off, in case of switching off the Lennard-Jones interactions, the switching begins at R vdw-switch . In case of CHARMM36 and Slipids, we tested the structural parameters also using different cut-off lengths (in brackets, not affecting the total CPU time in this table). CPU hours/project display the total CPU hours for the calculations -for obtaining the topologies and 30 ns z-constraint simulations for MD simulations and for DFT calculations and final free energy profile calculation in case of COSMOmic. The detailed CPU times are in Table S2. .......................60 Table 5: Mean differences and mean absolute differences of water/lipids ΔG wat and penetration ΔG pen barriers with respect to data obtained from Slipids FF. The values in brackets show the differences with excluded outlier (2- 8: Characterization of the interaction between a series of polyphenols and DOPC bilayer, as given by their positions, orientations, dipole moments, and free enthalpies of partition (ΔG part ) and penetration (ΔG pen ) are reported. .........................................................83 Table 9: Association energies and enthalpies (kcal.mol -1 ) calculated as the difference in energy (enthalpy) between the most stable complex and the isolated fragments, in the gas phase and in PCM-type benzene and water solvents. Negative values indicate that the association is thermodynamically favored compared to the pair of isolated fragments quercetin and vitE. .

List of Figures

List of Tables

Introduction

During their journey through the organism, the vast majority of drugs interact with membranes; they can even interact with membrane at the end of their road i.e. their action site. Some drugs never cross membranes (e.g., in the gastrointestinal tract) while others can easily diffuse through lipid bilayers, or even accumulate inside them. This has a strong impact of their mechanism of action, and so their activity and toxicity. Thus, drug capacity to incorporate/cross lipid bilayer membranes is of prime importance. How is it possible to evaluate drug-membrane interactions? Is there a universal rule that relates a drug chemical structure to its capacity to cross or to affect membranes? Can we predict the mechanism of action of a drug in which interaction with membrane is involved? Is it possible to reveal the intimate details of these interactions? If yes, can we predict and design the perfect drug? In this manuscript, we tackle some of these fascinating questions.

Because experimental techniques sometimes require expensive and long investments to get partial information, in silico molecular modelling has become a powerful alternative to tackle these issues. In the past few years, molecular dynamics (MD) has opened many perspectives, providing an atomistic description of the related intermolecular interactions, efficiently supporting experimental data. Using MD simulations, we have explored the capacity of several compounds (polyphenols, vitamins E and C, plantazolicin, carprofens) to incorporate lipid bilayer membranes. The different compounds were studied according to their different biological functions, namely (i) antioxidant activity against lipid peroxidation, (ii) antimicrobial activity with the possibility of trans-membrane pore formation, or (iii) inhibition of enzymes involved in Alzheimer's disease. In order to rationalize their mechanism of action, their location and orientation in membranes were assessed; accumulation and permeation capacities were also evaluated lipid bilayers were assessed.

Having in mind a predictive purpose in drug design, the accuracy of MD simulations relies in particular on quality of the in silico membrane models. By ensuring correlation between experimental and theoretical data, methodological improvements have been implemented on membrane. In particular, force field selection, xenobiotic parameterization and bilayer constitution have emerged as crucial factors for a correct prediction of drug-membrane interactions. The latter issue (composition) strongly impacts membrane penetration and lipid mixtures have recently been built in silico and properly parameterized; the role of cholesterol has deserved a particular attention. This manuscript is divided into four chapters that provide a gradual description of molecular interactions between drugs and membranes. At this point, it should be noted that henceforth the generic term 'drug' refers to natural but also synthetic or semisynthetic compounds exhibiting therapeutic actions or protective effects on human health.

To tackle drug-membrane interactions, a comprehensive description of membranes is required. This is the subject of Chapter I. Following a general introduction on importance of membranes in biology and therapeutics, their complex composition is detailed in section I.2. Then, specific the biophysical properties related to lipid bilayers are depicted in section I.3. Along this chapter, the focus is made on some notions that have been developed in our publications, such as structure of phospholipids, membrane phases and domains.

Chapter II focuses on MD simulation methodologies. The aim of this chapter is dual: to describe the underlying parameters and equations of MD simulations; and to serve as a manual on how to set up and analyze them. It begins with a short introduction to molecular modelling and the particular advantages of MD simulations. Then, it focuses on molecular mechanics and the notion of force field, which holds every parameter of interatomic interactions. Practical considerations guiding the choice of force fields and drug parameterization are given. In section II.3, the algorithms operating time evolution in MD are presented, such as temperature and pressure regulation. Another interesting ability of MD simulations is the calculation of free energy profiles, also denoted as potentials of mean force. They are related to partition coefficients, permeability and ultimately bioavailability, and thus are a powerful tool in pharmacology. The different methodologies and their respective advantages are reported in section II.4. The last section of this chapter refers to the analysis of MD simulations and to the different quantities that can be measured from the coordinate trajectories obtained by MD simulations. Thanks to the atomic resolution of MD simulations, various micro-and macroscopic data can be statistically analyzed e.g., molecular location, orientation, diffusion coefficients, and order parameters. The relationship between theoretical and experimental data is emphasized.

Chapter III reviews the literature on drug-membrane interactions evaluated by MD simulations. The advantages of in silico models and a brief history of the topic are stated in section III.1. Then (section III.2), the current knowledge on mechanisms of action of several therapeutic classes is highlighted. These classes were chosen because biological activities of the related drugs involve their interaction with lipid bilayers. Some drugs can affect the biophysical properties of membranes such as thickness, fluidity or order. As exemplified in section III.2, these structural modifications may be directly or indirectly correlated to their biological and pharmaceutical activities. Therefore, the evaluation of membrane perturbation by drugs at the molecular level may help rationalizing their mechanism of action. Namely, are reviewed (i) non-specific membrane effects of anesthetics and β-blockers; (ii) perturbation of lipid bilayers by NSAIDs, and (iii) inhibition of lipid peroxidation by antioxidants. Eventually, the limitations of MD simulations and perspectives for methodological improvements are discussed in section III.3.

Finally, Chapter IV lists some of the publications I co-authored during my Ph.D. studies. In a methodological study reported in section IV.1, we benchmarked five force fields and an alternative method (COSMOmic) for their accuracy in prediction of partition coefficients of small organic molecules. Next (sections IV.2-4) three publications of antioxidants are collected. Antioxidants are essential to inhibit deleterious effects of an overproduction of reactive oxygen species, which in membranes induces lipid peroxidation. To efficiently inhibit lipid peroxidation, antioxidants must penetrate lipid bilayers. MD simulations have appeared a valuable tool to determine antioxidant position in membranes. First we assessed antioxidant activity of lipocarbazole (section IV.2), a natural compound of bacterial origin. The combination of its free radical scavenging capacity, its affinity to membranes and its location rather deep in lipid bilayers rationalized its powerful antioxidant activity measure experimentally. Section IV.3 is the evaluation of twelve polyphenols well-known for their powerful antioxidant capacity. The relationship between their structures and their interaction with membranes is established. Section IV.4 better rationalizes collaborative effects between three antioxidants namely vitamin E, vitamin C and quercetin by that occur in non-covalent associates formed in membrane. These supra-molecular assemblies may have applications in the research of cocktails of antioxidants. The existence of such non-covalent assemblies in lipid bilayers opens many perspectives for other drugs. As biological membranes are complex mixtures of various lipids, we also investigated the influence of membrane composition. Section 0 is the evaluation of positions and orientations of two carprofen derivatives with lipid bilayers being in two phases and at two temperatures. These compounds are known to inhibit a transmembrane protein involved in Alzheimer disease, so their location in lipid bilayers is matter of crucial importance. Interestingly, lipid bilayer composition has a strong influence on position and orientation of these compounds. This work confirms that (i) membrane composition should be identical when comparing experimental and theoretical data and should reflect similar physical characteristics, and (ii) simple bilayer models -whether in vitro or in silico -may not be sufficient to account for the complexity of real biological membranes. Finally, as a perspective we assessed the ability of an antimicrobial compound (namely plantazolicin) to form pores in bilayers.

Chapter I. Molecular dynamics simulations of membranes

I.1. Introduction

Lipid bilayer membranes are essential to life, as they constitute boundaries of all biological cells. They accomplish numerous vital functions. Because living organisms are made of cells and that membranes encompass cells, the vast majority of drugs have to interact with lipid bilayer membranes, either during their action or on the way to their site of action. Understanding molecular mechanisms of interaction of drugs with membranes may help improvement of their activity and decrease of their toxicity. In this perspective, molecular dynamics (MD) simulations of drugs with membranes has become an attractive tool to address issues that are not readily accessible by experimental techniques. Both atomic and fs time resolution are accessible with MD simulations.

To study processes occurring in lipid bilayer membranes, this chapter will first detail their constitution in section I.2 and then their physical properties in section I.3. The relation to MD simulations is emphasized.

I.2. Membrane composition

Among all living organisms, a wide range of membrane sizes, shapes, functions and composition exist. In biology, the word 'membrane' refers to the envelope that delimits intraand extra-cellular compartments. Here, it is worth noting to distinguish envelopes from membranes, namely other layers than a bilayer constitute an envelope. For instance, Grampositive bacterial cell wall is constituted of a lipid bilayer and peptidoglycan and anionic polymers. This manuscript essentially focuses on lipid bilayer membranes. Thus throughout text, as a misusage, the word 'membrane' may sometimes refer to lipid bilayer membrane. It should be noted that organelles such as Golgi apparatus, endoplasmic reticulum, mitochondria or nucleus are also delimited by lipid bilayer membranes.

Three main components compose membranes: lipids, proteins and carbohydrates (Figure 1). However, the membrane composition may dramatically differ in nature and respective concentrations, depending on species and types of cells (Figure 2). The current computer power does not allow simulating membrane models sufficiently large to fully describe a biological membrane, at least at a relevant sampling time. It is however possible to model one or a few proteins embedded in a lipid bilayer. In this work, we mainly focus on interactions of drugs and natural compounds with lipid bilayers.

I.2.1. Lipid composition

In most of membranes, lipids constitute half the membrane (Figure 2), being major components of the membrane that drive physical properties and strongly influence protein functions. The lipid bilayer is a selective barrier allowing certain xenobiotics to diffuse from extra to intra-cellular compartments (passive diffusion).

In mammals, the main lipids constituting membranes are phospholipids, sterols and glycolipids. Phospholipids and glycolipids consist of a polar hydrophilic head group and lipophilic lipid chains (tails). Thus, they are amphiphilic. This is also true to a lesser extent for sterols that have a polar OH group. This property drives the spontaneous association of lipids either into micelles or lipid bilayers. Sterols are of conical shape that favors micelles, whereas phospholipids and glycolipids are mostly cylindrical, promoting bilayers. To minimize edge effects, bilayers can adopt spherical shapes. For instance, liposomes and cell membranes are lipid bilayers adopting a spherical-like shape. For the MD simulations, pure phospholipid bilayers have been mainly used, as they provide a relevant representation of lipid bilayers constituting biological membranes. In the past few years, models of lipid mixtures have also been developed to mimic some selective interactions with different lipids. Although they can be more realistic, they usually require longer simulation time to reach a correct sampling of molecular motions.

I.2.1.1 Phospholipids

Phospholipids are the principal lipid components of membranes. They all are constituted similarly i.e. one to four lipid tails (generally two, mostly fatty acids); a central platform (glycerol or sphingosine); and a polar head group consisting of a phosphate moiety and a polar alcohol (Scheme 1 and Figure 4). Phospholipids with a glycerol backbone are called phosphoglycerides, whereas sphingolipids are those having a sphingosine backbone.

Ceramides are special kinds of sphingolipids that do not bear a polar head group other than OH groups of sphingosine. Fatty acids are important moieties of phospholipids. They vary in length (i.e. number of carbon atoms quoted C, mainly an even number) and in the number of unsaturation, quoted D. The most common fatty acids have 16 or 18 carbon atoms. In their natural form, double bonds are cis, but trans forms can also be synthetized by dehydrogenation. Unsaturated fatty acids have one saturated methylene group in between two double bonds. The first double bond is often either 3 or 6 atoms from the end of the lipid chain (quoted n-3 and n-6 fatty acids, respectively, also called ω-3 and ω-6). The fatty acid composition of a lipid bilayer has a dramatic impact on its physical properties. The length of lipid chains is obviously directly correlated to membrane thickness. The number of cis unsaturations influences chain packing, and thus melting temperature T m (i.e., transition temperature between gel and fluid phases, see section I. 3.4). For instance, oleate chains (18:1) spread on larger lateral space and are more disordered than stearate chains (18:0) at similar temperature. Therefore, T m is lower in phospholipids containing oleate than stearate chains.

The other part defining a phospholipid is its polar head groups, consisting of a phosphate moiety linked to an alcohol. Alcohols found in natural phospholipids are represented in Figure 5.

Head groups can provide different charges to phospholipids. For instance, phosphocholines (PC) are globally neutral (one negative charge on phosphate moieties and one positive on choline moieties). Similarly, phosphor ethanolamines (PE) are neutral, whereas phosphoserines (PS), phosphoglycerol (PG) and phosphoinositol (PI) are negatively charged. The proportion of charged lipids and their distribution in both leaflets can induce a surface potential that may dramatically affect biological processes [START_REF] Goormaghtigh | Evidence of a specific complex between adriamycin and negatively-charged phospholipids[END_REF][START_REF] Fong | Stabilization of acetylcholine receptor secondary structure by cholesterol and negatively charged phospholipids in membranes[END_REF][START_REF] Van Dijck | Negatively charged phospholipids and their position in the cholesterol affinity sequence[END_REF].

For convenience reasons, the nomenclature of phospholipid heads and tails can be reduced to acronyms. For phosphoglycerides, the first two letters stand for fatty acids and the last two for the head group, see Table 2 for examples.

Abbreviation Chemical name DMPC 1,2-dimyristoyl-sn-glycero-3-phosphocholine DPPC 1,2-dipalmitoyl-sn-glycero-3-phosphocholine DSPC 1,2-distearoyl-sn-glycero-3-phosphocholine POPC 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine DOPC 1,2-dioleoyl-sn-glycero-3-phosphocholine DLPC 1,2-dilinoleoyl-sn-glycero-3-phosphocholine DOPE 1,2-dioleoyl-sn-glycero-3-phosphoethanolamine POPS 1-palmitoyl-2-oleoyl-sn-glycero-3-phospho-L-serine DPPG 1,2-dipalmitoyl-sn-glycero-3-phospho-(1'-rac-glycerol) 18:0 SM N-stearoyl-D-erythro-sphingosylphosphorylcholine
Membrane composition in lipids is highly variable according to species, individuals, and types of cell, organelles nutritional intakes and many other factors. Some lipids are only present in bacteria (e.g., lipid A [START_REF] Raetz | Discovery of new biosynthetic pathways: the lipid A story[END_REF]), others are almost exclusively found in mitochondrial membranes (e.g., cardiolipin [START_REF] Schlame | Mitochondrial cardiolipin in diverse eukaryotes[END_REF][START_REF] Nowicki | Cardiolipin synthase of Arabidopsis thaliana[END_REF]). Additionally, all types of head group and fatty acids can recombine to yield more than 1,000 phospholipids [START_REF] Van Meer | Membrane lipids: where they are and how they behave[END_REF]. Therefore, the analysis of phospholipids often distinguishes head group from fatty acid composition. In the excellent review of van Meer et al., the synthesis, transport and distribution of lipids is reported in eukaryotic and yeast membranes [START_REF] Van Meer | Membrane lipids: where they are and how they behave[END_REF].

Figure 7: Phospholipid distribution in mammals and yeast plasma membrane. The sterol / phospholipid ratio is mentionned for mammals (cholesterol, CHOL) and for yeast (ergosterol, ERG). Adapted from ref. [START_REF] Van Meer | Membrane lipids: where they are and how they behave[END_REF].

As typical examples, red blood cell membranes have also been extensively studied [START_REF] Dodge | Composition of phospholipids and of phospholipid fatty acids and aldehydes in human red cells[END_REF][START_REF] Verkleij | The asymmetric distribution of phospholipids in the human red cell membrane. A combined study using phospholipases and freeze-etch electron microscopy[END_REF]. In this case, the proportion of PC and SM is higher in the outer leaflet whereas PE and PS are preferentially found in the inner leaflet. In the brain, both length of fatty acids and type of phospholipids varies when considering neurons, astrocytes or gangliocytes, with a higher proportion of SM in the myelin sheath surrounding axons [START_REF] Norton | The lipid composition of isolated brain cells and axons[END_REF].

In another example, the most external layer of skin -stratum corneum -acts like a barrier to external molecules and microorganisms. It consists of a stack of dead cells embedded in a lipid matrix constituted of a strong proportion of ceramides with long fatty acid chains [START_REF] Menon | The structure and function of the stratum corneum[END_REF][START_REF] Paloncýová | Rationalization of Reduced Penetration of Drugs through Ceramide Gel Phase Membrane[END_REF].

Ceramides are not phospholipids stricto sensu as they lack the phosphate group; they are classified as sphingomyelins.

I.2.1.2 Sterols

Sterols constitute an important part of membranes as they represent up to 50 % of membrane lipids. They are steroid lipids with bulky, conical shape. The OH group at position 3 is oriented towards the water phase in bilayers, embedded in the polar head groups. A particular case of sterols is cholesterol (Figure 8) that is found in all eukaryotic cells but is absent in prokaryotic cells. It is the predominant sterol in animals. Plants contain large amount of cholesterol but also significant proportion of sitosterol [START_REF] Behrman | Cholesterol and Plants[END_REF]. In yeasts, ergosterol (Figure 8) is the dominant sterol. Sterols are known to influence the order of membranes. When the bilayer is in its fluid phase, sterol addition increases the order. Conversely, it will slightly decrease the order of a 'gel' phase membrane while also increasing the fluidity. Sterols thus favor the liquid ordered phase (L o ) defined as a highly fluidity and low-ordered phase [START_REF] Petruzielo | Phase behavior and domain size in sphingomyelin-containing lipid bilayers[END_REF][START_REF] Berkowitz | Detailed molecular dynamics simulations of model biological membranes containing cholesterol[END_REF].

Regarding the abundance of cholesterol in biological membranes, the incorporation of this constituent in membrane model has appeared mandatory to predict membrane behaviors.

The effects of cholesterol on lipid bilayers have been studied by MD simulations in the past years. Since the first investigations in the 1990's, atomistic [START_REF] Berkowitz | Detailed molecular dynamics simulations of model biological membranes containing cholesterol[END_REF] and coarse grained [START_REF] Müller | Biological and synthetic membranes: What can be learned from a coarse-grained description?[END_REF][START_REF] Venturoli | Mesoscopic models of biological membranes[END_REF][START_REF] Waheed | Phase Transitions in Coarse-Grained Lipid Bilayers Containing Cholesterol by Molecular Dynamics Simulations[END_REF] models of cholesterol-containing lipid bilayer models help rationalize the influence of cholesterol on membrane physical properties and raft formation. We have used membrane mixtures containing cholesterol in section IV.5.

I.2.1.3 Glycolipids

The third family of lipids constituting membranes is glycolipids. They are based on sphingosin linked with an ether bond to one or several sugar moieties e.g., cerebroside (Scheme 1 and Figure 4). Complex residues, including gangliosides, bear longer sugar chains with up to seven sugar units. They are always distributed asymmetrically in all eukaryotic membranes, as carbohydrate moieties are only present in the outer layer. A typical function of glycolipids is cell signaling. They are for instance responsible for the blood type ABO antigens.

Another example related to glycolipids is lipopolysaccharide, a constituent of the outer membrane of Gram-negative bacteria (Scheme 2). It consists of a complex lipid part (lipid A) and a long polysaccharide chain. As it is an endotoxin, it induces a stimulation of the immune system in hosts infected with Gram-negative bacteria.

Scheme 2: Gram negative bacteria cell wall MD simulations of glycolipid-containing membranes have been performed [START_REF] Demarco | Molecular Dynamics Simulations of Membrane-and Protein-Bound Glycolipids Using GLYCAM[END_REF][START_REF] Demarco | Presentation of Membrane-Anchored Glycosphingolipids Determined from Molecular Dynamics Simulations and NMR Paramagnetic Relaxation Rate Enhancement[END_REF][START_REF] Kapla | Molecular Dynamics Simulations of Membranes Composed of Glycolipids and Phospholipids[END_REF][START_REF] Patel | Structure and Dynamics of Glycosphingolipids in Lipid Bilayers: Insights from Molecular Dynamics Simulations[END_REF], although glycolipids are less extensively studied than phospholipids or sterols. The principal limitations are the limited availability of force field parameters and the variability in glycolipid structures.

I.2.2. Proteins

Proteins are essential components of biological membranes, as they represent from 18 to 75 % of membrane mass (Figure 2). Myelin sheath membranes have the least amount of protein, while specialized organelles such as mitochondria or chloroplast have the highest amount. Proteins are responsible for most of biological processes occurring in or close to membranes. These functions are highly cell type dependent; for instance protein receptors transmit chemical information through membrane, pumps and transporters carry molecules and ions in and out; enzymes transform molecules; pores are channel allowing molecules to pass through; and peptides can signal immune system or cell attachment.

Protein attachment to lipid bilayers is mainly related to their secondary structures (Figure 1). In general, α-helices are constituted of hydrophobic aminoacids. This provides such part of proteins a lipophilic character particularly adapted to anchors proteins rich in α-helices into transmembrane domains (TMD). Some TMDs are constituted of β-sheets forming a β-barrel, for instance in porins. Proteins can also be attached in membranes via modified amino acids that act as lipid anchors, such as palmitylcystein.

Although transporters and pores span the whole membrane thickness, some smaller proteins are just adsorbed on the surface or in one layer of the membrane. They are often associated to transporters or pores and act as regulators or signal transductors.

I.3. Membrane physical properties

Understanding the physical properties of membranes first requires a thorough knowledge of their composition. Among other descriptors, size, thickness, curvature, area per lipid, lipid phases, order, domains and rafts characterize membranes. All of the descriptors are directly affected by any biological processes and drug interaction with membranes. In this section, some usual values of these descriptors are given, showing relationships between them if any. The methodological aspects on how to evaluate them by MD simulations are given in the next chapter, section II.5.4.

I.3.1. Thickness

Thickness is the most obvious membrane property. Usual thickness for a biological membrane is around 5 nm. Naturally, thickness varies according lipid chain length but also lipid phase which is related to chain order; namely if order increases, chains are straighter and bilayers are thicker. For instance, with the same composition, a 'gel' phase membrane is thicker than a fluid phase. Other factors such as membrane composition or temperature influence lipid phase and thus thickness indirectly. Definitions and methods to evaluate membrane thickness are given in section II.5.4.1.

I.3.2. Curvature

For small liposomes, membrane curvature directly depends on liposome diameter. However, for larger membranes as plasma membranes, local curvatures can be induced by an asymmetric lipid composition in both leaflets. Local curvatures are dynamic processes which are inherent of lipid structure and flexibility [START_REF] Kamal | Measurement of the membrane curvature preference of phospholipids reveals only weak coupling between lipid shape and leaflet curvature[END_REF]. Indeed, lipids with a cylindrical shape spontaneously form flat bilayers whereas phospholipids with short lipid tails and large head group exhibit conical shape and induce (positive) curvature. Similarly, lipids with long, unsaturated chains and small head groups induce (negative) curvature [START_REF] Kamal | Measurement of the membrane curvature preference of phospholipids reveals only weak coupling between lipid shape and leaflet curvature[END_REF]. A simple approach to predict membrane curvature was proposed by Israelachvili et al. [25]: the lipid packing parameter 𝑃:

𝑃 = 𝑣 𝑎𝑙
where 𝑣 is the molecular volume, 𝑎 the head group cross section area and 𝑙 the lipid length.

Figure 9: Influence of the lipid packing parameter P on the lipid bilayer curvature Reproduced from ref. [START_REF] Yandrapalli | Lipid domains in HIV-1 assembly[END_REF].

Although spontaneous curvature does play a role in liposomes, in biological membranes local curvatures do not drive association of lipids with similar shapes [START_REF] Kamal | Measurement of the membrane curvature preference of phospholipids reveals only weak coupling between lipid shape and leaflet curvature[END_REF]. Therefore, other processes are most likely responsible for membrane curvature in biological membranes [START_REF] Cooke | Coupling between Lipid Shape and Membrane Curvature[END_REF][START_REF] Sorre | Curvature-Driven Lipid Sorting in Biomembranes[END_REF].

Membrane fusion can occur between biological membranes due to strong local curvatures. They are driven by membrane proteins (e.g., SNAREs) [START_REF] Martens | Mechanisms of membrane fusion: disparate players and common principles[END_REF][START_REF] Mcmahon | Membrane curvature and mechanisms of dynamic cell membrane remodelling[END_REF].

Membrane curvature decreases or increases lipid head group packing, thus facilitating or hindering drug penetration, respectively. It was shown that addition of free fatty acids and lysolipids induces membrane curvature and lowers permeability barriers [START_REF] Jespersen | Lipids, curvature stress, and the action of lipid prodrugs: Free fatty acids and lysolipid enhancement of drug transport across liposomal membranes[END_REF].

I.3.3. Area per lipid

Area per lipid is, as a self-explanatory definition, the average area occupied by one lipid in a bilayer. Although membrane curvature relates to head group packing, area per lipid is associated with lipid packing. As for bilayer thickness and curvature, it depends on lipid composition, temperature, and subsequently fluidity and lipid order. For instance, addition of cholesterol in PC bilayers reduces the area per lipid [START_REF] Edholm | Areas of Molecules in Membranes Consisting of Mixtures[END_REF][START_REF] Lindblom | Lipid lateral diffusion and membrane heterogeneity[END_REF][START_REF] Kheyfets | Area per Lipid in DPPC-Cholesterol Bilayers:Analytical Approach[END_REF]. Besides, area per lipid also influences drug partitioning [START_REF] Wennberg | Large Influence of Cholesterol on Solute Partitioning into Lipid Membranes[END_REF].

I.3.4. Fluidity, order and lipid phase

Membrane fluidity can be defined by the lateral diffusion coefficient of lipids, 𝐷 𝐿 . The lateral diffusion of lipids was first observed in 1970 by observing the movement of antigens in cells resulting from the fusion of mouse and human cells [START_REF] Frye | The Rapid Intermixing of Cell Surface Antigens After Formation of Mouse-Human Heterokaryons[END_REF]. Since then, the fluid mosaic model of membrane has been widely accepted, that is that membrane proteins are embedded in fluid matrix of phospholipids [START_REF] Singer | The fluid mosaic model of the structure of cell membranes[END_REF]. It can now be efficiently evaluated by fluorescence techniques such as fluorescence recovery after photobleaching (FRAP). We must note that the higher 𝐷 𝐿 , the more fluid the membrane. Fluidity is increased vs temperature, but it also depends on membrane composition. Saturated lipids facilitate lipid packing compared to unsaturated ones, and thus lower fluidity. Typical values of 𝐷 𝐿 are ranging from 10 -7 to 10 -8 cm 2 s -1 for fluid phase membranes and from 10 -8 to 10 -9 cm 2 s -1 for more ordered membranes [START_REF] Lindblom | Lipid lateral diffusion and membrane heterogeneity[END_REF]. At a macroscopic scale, this means that lipids move rapidly towards cell surface. For instance, a lipid can diffuse to the opposite side of a bacterial cell within one second. In the time scale of MD simulations (ca. a few µs), lipid diffusion is clearly observed.

Another property is the diffusive capacity of lipids to cross the membrane and change leaflet; this process is usually called 'flip-flop'. However, because polar head groups have to cross the hydrophobic core, this process takes place much more slowly 10 -15 s -1 [START_REF] Contreras | Transbilayer (flip-flop) lipid motion and lipid scrambling in membranes[END_REF].

Fluidity can also be described by the order parameter. It is a dimensionless quantity including both the preferential orientation defined by an angle 𝜃 and the related deviations from the average value:

𝑆 = 3 2 〈cos 2 𝜃〉 - 1 2
where brackets mean time average. To simplify, 𝜃 can be considered as the angle between the C-H bonds of the lipid tails and membrane normal vector. 𝑆 can take values in the [-0.5 ; 1] interval. When lipid tails are completely disordered, 𝑆 equals 0; whereas a perfectly ordered membrane (i.e. parallel lipid chains) bears a value of -0.5 [START_REF] Berkowitz | Detailed molecular dynamics simulations of model biological membranes containing cholesterol[END_REF]. A value of 1 means a perfectly ordered membrane with the orientation of C-H bonds parallel to the reference axis, that is with the lipid chains perpendicular to the reference axis. A more detailed explanation of order parameters is given in section II.5.4.3.

Both fluidity and order characterize lipid phases. There exist three phases: 1) The liquid disordered phase (L d , also called fluid phase, liquid crystalline phase or L α ) is characterized by a high fluidity and low order. This is the state of pure phospholipid membranes above their gel-to-liquid transition temperature, Tm. It is favored by lipids with short and/or unsaturated lipid chains that possess cis double bonds as for all unsaturated natural lipids.

2) The liquid ordered phase (L o ) features a high fluidity and a high order. This state is typically present in binary or ternary mixtures containing cholesterol.

3) The gel phase also referred as solid crystalline phase, S o or L β , is characterized by low fluidity and high order. It is favored by saturated and long lipid chains that can all align to form an ordered crystal and having high Tm values, thus favoring the existence of the S o phase at room temperature.

Experimental and theoretical studies have repeatedly reported influence of temperature and membrane composition on partition coefficients, thus suggesting influence of lipid phase [START_REF] Wennberg | Large Influence of Cholesterol on Solute Partitioning into Lipid Membranes[END_REF][START_REF] Luxnat | Partition of chlorpromazine into lipid bilayer membranes: the effect of membrane structure and composition[END_REF]. It should be noted that there is no trivial correlation between lipid phase and drug permeability. Numerous factors such as membrane composition, as well as drug polarity or size participate in membrane permeability [START_REF] Wennberg | Large Influence of Cholesterol on Solute Partitioning into Lipid Membranes[END_REF].

I.3.5. Rafts and domains

As mentioned above, lipid composition of both leaflets of a cell membrane is different, thus transmembrane or 'vertical' segregation occurs. Likewise, lateral segregation can also occur in lipid mixtures. In other words, the lateral distribution of a given lipid type may be inhomogeneous. Hence, 'patches' of membrane with a specific composition float in a mixture of a different lipid composition. This is the raft hypothesis, first suggested in the 1970's [START_REF] Stier | Spin labels as enzyme substrates Heterogeneous lipid distribution in liver microsomal membranes[END_REF][START_REF] Israelachvili | Physical principles of membrane organization[END_REF], refined in 1997 [START_REF] Simons | Functional rafts in cell membranes[END_REF] and now widely evidenced and accepted [START_REF] Lingwood | Lipid rafts as a membrane-organizing principle[END_REF][START_REF] Simons | Model Systems, Lipid Rafts, and Cell Membranes1[END_REF]. Rafts are defined as 'small (10-200 nm) heterogeneous, highly dynamic, sterol-and sphingolipidenriched domains that compartmentalize cellular processes. Small rafts can sometimes be stabilized to form larger platforms through protein-protein and protein-lipid interactions' [START_REF] Pike | Rafts defined: a report on the Keystone symposium on lipid rafts and cell function[END_REF]. Depending on their size, they can also be termed micro-or nanodomains [START_REF] Petruzielo | Phase behavior and domain size in sphingomyelin-containing lipid bilayers[END_REF].

Chapter II. Theory and Methods

II.1. Introduction -Born-Oppenheimer approximation

In the end of the 19 th century, Max Planck observed phenomena that were not explained by classical physics. He postulated that the emission of light was not described by a continuum energetic spectrum but as a series of discrete energy quanta. This was further generalized by Einstein, de Broglie and others, light and other particles being considered both as waves and corpuscles (wave-particle duality). So quantum mechanics was born. In 1926, Erwin Schrödinger formulated an equation to describe the motion of electrons: [START_REF] Schrödinger | An Undulatory Theory of the Mechanics of Atoms and Molecules[END_REF]

𝑯𝛹 = 𝐸𝛹

where H is the Hamiltonian operator, E the energy of the system and Ψ the wave function containing all information on the studied particles. This time-independent Schrödinger equation is the basis of most of theoretical chemistry and molecular modeling methods of calculation. It describes any system of particles, including molecules made of electrons and nuclei. As the mass and the momentum of electrons are very different from those of nuclei, the Born and Oppenheimer approximation [START_REF] Born | Zur Quantentheorie der Molekeln[END_REF] proposes to decouple and treat separately both motions of electrons and nuclei. This gives rise to two main families of methods in theoretical chemistry. When treating electrons (quantum chemistry, QM), they are considered in a field of fixed nuclei; when treating nuclei (molecular mechanics, MM), they are considered in a field of forces representing all interactions including those of electrons.

Treating electrons allows accessing valuable chemical information. An accurate evaluation of ground and transition state energies allows tackling precisely the reactivity. Among other applications, this can drive the choice of catalyzers or activating groups; predict light absorption and emission properties in the UV-visible range thus rationalizing colors of molecules and sunscreen capacity; rationalize NMR data, polarizability, ionization potential and electron affinity.

The exact solution of the Schrödinger equation is only accessible for systems containing one electron. It means that for any polyelectronic system, approximations are required. Their study has really started to be developed in the 70's with the developments of computational facilities. From this period, according to Moore's law, computer power has doubled every 18 months. This has allowed tackling bigger and bigger molecular systems together with better and better precision. Over the past decades the number of methods of calculations has dramatically increased, which have been classified as semi-empirical, Hartree-Fock (HF), post-HF and density functional theory (DFT) methods. The accuracy of quantum chemistry calculations depends on the method used. In general, the less approximation, the better accuracy; however this is generally correlated with dramatic increase of computational time.

Therefore the challenge of a theoretical chemist is always to choose the best compromise between accuracy and computational time.

Whatever the method of calculations, QM can only treat relatively small molecular systems (up to 100, 1000 and 10000 atoms with post-HF, DFT and semi-empirical methods, respectively). In order to treat bigger systems, MM appears as an adequate alternative. In this thesis, we are mainly dealing with modelling lipid bilayers, which contain several tens of thousands of atoms. Therefore, the emphasis is made on MM and molecular dynamics (MD) methods. Specific QM-related topics are only developed in the respective methods' section of each project in Chapter III.

This chapter first describes principles and equations behind MM, MD and force fields. The analysis of MD simulations is then detailed in order to get valuable data to study the interaction of drugs with lipid bilayers.

II.2. Molecular mechanics

MM treats motion of particles within the classical (Newton) physics instead of quantum physics. Atoms are pictured as balls and bonds as springs. MM is tremendously less computational demanding than QM and thus allows treating molecular systems containing tens of thousands of atoms. Each interaction between all atoms is parameterized empirically.

There are parameters for the interactions between atoms separated by one or several bonds (i.e., bonded interactions) and between atoms through space independently of the existence of bonds between them (i.e., non-bonded interactions). The collection of all these parameters is one part of the so-called force field. The second part of the force field is the set of equations that are used to calculate the potential energies and ultimately the forces by derivation. The function assembling this set of equations to calculate the total energy in the case of a force field can be given as: Interactions between atoms separated by bonds consist of three types. Bond stretching describes energetic dependence vs length 𝑟 𝑖𝑗 of a bond between atoms 𝑖 and 𝑗 by a harmonic potential 𝑉 𝑏𝑜𝑛𝑑 (Figure 10):

𝑉 𝑡𝑜𝑡𝑎𝑙 =
𝑉 𝑏𝑜𝑛𝑑 (𝑟 𝑖𝑗 ) = 1 2 𝑘 𝑖𝑗 𝑏 (𝑟 𝑖𝑗 -𝑏 𝑖𝑗 ) 2
where 𝑘 𝑖𝑗 𝑏 is the force constant of the "spring" and 𝑏 𝑖𝑗 the equilibrium bond length.

The energetic variation vs the angle 𝜃 𝑖𝑗𝑘 formed by three atoms are also usually described by a harmonic potential 𝑉 𝑎𝑛𝑔𝑙𝑒 :

𝑉 𝑎𝑛𝑔𝑙𝑒 (𝜃 𝑖𝑗𝑘 ) = 1 2 𝑘 𝑖𝑗𝑘 𝜃 (𝜃 𝑖𝑗𝑘 -𝜃 𝑖𝑗𝑘 0 ) 2
where 𝑘 𝑖𝑗𝑘 𝜃 is the force constant and 𝜃 𝑖𝑗𝑘 0 the equilibrium angle.

Finally, dihedral angles are defined as the angle 𝜑 𝑖𝑗𝑘𝑙 formed by four atoms around a bond. One must distinguish (i) proper dihedrals, usually used to describe the periodic rotation around a single bond, and (ii) improper dihedrals, useful in out-of-plane motion in rings.

𝑉 𝑑𝑖ℎ𝑒𝑑𝑟𝑎𝑙 = 𝑉 𝑝𝑟𝑜𝑝𝑒𝑟 + 𝑉 𝑖𝑚𝑝𝑟𝑜𝑝𝑒𝑟

The most common form of proper dihedral potential 𝑉 𝑝𝑟𝑜𝑝𝑒𝑟 is periodic (𝑉 𝑝𝑟𝑜𝑝𝑒𝑟-𝑝𝑒𝑟𝑖𝑜𝑑𝑖𝑐 ):

𝑉 𝑝𝑟𝑜𝑝𝑒𝑟-𝑝𝑒𝑟𝑖𝑜𝑑𝑖𝑐 (𝜑 𝑖𝑗𝑘𝑙 ) = ∑ 𝑘 𝜑,𝑛 (1 + cos(𝑛𝜑 -𝜑 𝑠,𝑛 )) 𝑛
where 𝑘 𝜑 is the force constant, 𝑛 the periodicity and 𝜑 𝑠 the equilibrium dihedral angle. Depending on substituents of the two central atoms of the dihedral, potential can be represented as a weighted sum of cosines functions, for instance for alkanes. In this case, the Ryckaert-Bellemans expression of the potential can be more appropriate:

𝑉 𝑝𝑟𝑜𝑝𝑒𝑟-𝑅𝐵 (𝜑 𝑖𝑗𝑘𝑙 ) = ∑ 𝐶 𝑛 (cos(𝜓)) 𝑛 5 𝑛=0
Improper dihedrals are non-periodic and often described by a harmonic potential:

𝑉 𝑖𝑚𝑝𝑟𝑜𝑝𝑒𝑟 (𝜉 𝑖𝑗𝑘𝑙 ) = 1 2 𝑘 𝜉 (𝜉 𝑖𝑗𝑘𝑙 -𝜉 0 ) 2
Dihedral angles are important parameters in the conformation of a molecule as the rotation of a single bond can dramatically impact the tridimensional structure and thus on molecule properties. In some case, it has to be carefully reparameterized. For instance, the dihedral angle between π-conjugated rings of quercetin influences its antioxidant properties, and has been reparameterized within the Gromos force fields (see IV. 4.3.1).

II.2.2. Non-bonded interactions

Whether two atoms are within a molecule or not, non-bonded interactions always affect them. They are weak and rather long distance (several Ångströms) interactions being responsible for numerous physical, chemical or biological processes. For instance, ebullition temperature, chemical reactivity and the existence of lipid bilayer membranes are directly driven by non-bonded interactions.

There are three main types of non-bonding interactions: electrostatic (or Coulombic), Van der Waals and hydrogen bond interactions. Electrostatic interactions (between charges) are described by the classical Coulombic interaction given by:

𝑉 𝑒𝑙𝑒𝑐𝑡𝑟𝑜𝑠𝑡𝑎𝑡𝑖𝑐 (𝑟 𝑖𝑗 ) = 1 4𝜋𝜀 0 𝑞 𝑖 𝑞 𝑗 𝜀 𝑟 𝑟 𝑖𝑗
where 𝑟 𝑖𝑗 is the distance between atoms 𝑖 and 𝑗, 𝑞 𝑖 is the charge of atom 𝑖, and 𝜀 𝑟 is the relative dielectric constant.

The van der Waals interaction corresponds to the dispersion interactions (fluctuating dipolefluctuating dipole and higher order interactions). It is usually well described by the Lennard-Jones potential 𝑉 𝐿𝐽 : 𝑉 𝐿𝐽 (𝑟 𝑖𝑗 ) = 𝐶 𝑖𝑗 [START_REF] Menon | The structure and function of the stratum corneum[END_REF] 𝑟 𝑖𝑗 12 -𝐶 𝑖𝑗 [START_REF] Schlame | Mitochondrial cardiolipin in diverse eukaryotes[END_REF] 𝑟 𝑖𝑗 6

including the attractive dispersion term (-𝑟 -6 ) and the repulsive term (𝑟 -12 ) at short distances. The parameters 𝐶 𝑖𝑗 (12) and 𝐶 𝑖𝑗 (6) are defined per type of atom. For instance, an aromatic carbon atom does not exhibit the same van der Waals interaction as an aliphatic one. Therefore, each atom type (as defined depending on the chemical function to which they belong) bears atomic coefficients 𝐶 𝑖 (12) and 𝐶 𝑖 (6) , which allow calculating 𝐶 𝑖𝑗 (12) and 𝐶 𝑖𝑗 [START_REF] Schlame | Mitochondrial cardiolipin in diverse eukaryotes[END_REF] for a given 𝑖𝑗 bond. Depending on the force field, these combinations are obtained by either geometric or arithmetic averages.

Hydrogen bond interactions have a strong electrostatic character, therefore being often sufficiently well described by the electrostatic potential. However, some force fields of the AMBER or OPLS-AA families also include a specific description of hydrogen or halogen bonding [START_REF] Paton | Hydrogen Bonding and π-Stacking: How Reliable are Force Fields? A Critical Evaluation of Force Field Descriptions of Nonbonded Interactions[END_REF][START_REF] Jorgensen | Treatment of Halogen Bonding in the OPLS-AA Force Field: Application to Potent Anti-HIV Agents[END_REF]. This is performed either by fine-tuning electrostatic and van der Waals parameters for halogens and polar hydrogens or by introducing special potentials for these atom types.

Whereas bonded interactions are roughly proportional to the number of molecules, nonbonded interactions scale to 𝑁(𝑁 -1)/2, 𝑁 being the number of atoms. Thus, as the system size increases, the number of non-bonded interactions to compute becomes very large. However, the non-bonded interactions between two atoms further apart than 1.4 nm is small. Therefore, in order to reduce computational time, the non-bonded interactions between two atoms separated by more than a certain distance (i.e., a cutoff) are not directly evaluated. Depending on the force field, this cutoff can be set from 0.9 to 1.4 nm for Coulombic and van der Waals interactions. In order to account for small contributions to the potential energy of atoms further than the cutoff, long range corrections have to be applied. They consist of an additional potential energy term calculated in 3D space. The most popular implementation of long-range correction to electrostatic potential is the Particle-mesh Ewald (PME) [START_REF] Darden | Particle mesh Ewald: An N⋅log(N) method for Ewald sums in large systems[END_REF]. In this method, point charges are transformed into Gaussian charge potentials and reported on a 3D grid. This grid or mesh is then Fourier transformed, so that the Poisson equation calculating the potential from the charges can be solved much more easily in reciprocal space. PME has the advantage of scaling to 𝑁 • log 𝑁 thanks to efficient Fourier transform algorithms.

II.2.3. Force field resolution

Some phenomena (e.g., membrane crossing, protein folding) occur over long time periods (i.e., from several µs to several hours). Current computer power may prevent describing long biological processes with regular MD simulations if one wants avoid years of calculation. Two options are however available to tackle such processes: (i) lower the spatial resolution of the simulation by decreasing the number of atoms or (ii) lower the time resolution by increasing the time step.

The former possibility (decreasing the number of atoms) is not about reducing dimensions of the system but rather about omitting certain atoms however, keeping accuracy. The general idea is to regroup several atoms into one pseudo-atom instead of deleting some. For instance, aliphatic hydrogens are not as crucial as polar hydrogens involved in H-bonds. Thus, in an aliphatic carbon chain, two hydrogen atoms bonded to one carbon atom are regrouped into one "CH 2 " pseudo-atom. The force fields using this trick are called united atom force fields (UA-FF). Along this line, the number of atoms in a lipid chain is divided by 3, greatly improving the speed of membrane simulations. The parameterization of the pseudo-atoms is achieved by defining new atom types which have refined bonded and nonbonded parameters. Parameterization is usually achieved by fitting these parameters to experimental values.

The latter possibility to speed up the simulation is to increase time step (see section II.3.1). However, time step is fixed to a minimum so that the fastest vibration can be properly simulated. In UA-FF, polar hydrogens are still present and have the fastest vibration, thus the time step cannot be larger than for all-atom force fields (AA-FF). However, by regrouping not only hydrogens with heavy atoms (e.g., C, N, O…) but also heavy atoms together, hydrogens can be completely omitted and time step can be increased 10 to 20 times. This also considerably lowers the number of atoms and additionally increases the speed of simulations. The clustering of heavy atoms into large beads is the principle of coarse-grained force fields (CG-FF). The number of atom types and Hamiltonian parameters are also reduced, e.g. there is no explicit description of dihedral angles. Accounting for all these simplifications, one would expect that sampling time is increased several tens of times. However, as many degrees of freedom are neglected, friction between atoms is reduced and events may artificially occur much faster than they should do. For instance, the MARTINI 2.0 CG-FF simulates events approximatively 4 times too fast, thus the effective speed up is actually 2 to 10 times, depending on the system [START_REF] Marrink | The MARTINI Force Field: Coarse Grained Model for Biomolecular Simulations[END_REF]. The effective speed-up can be measured by comparing diffusion coefficients to experimental values. In summary, the lower the resolution of a given force field, the faster the simulation but the less detailed the structural information. Also, result reliability directly depends on the quality of force field parameterization, regardless resolution.

II.2.4. Force field versions and specificities

A number of force fields exist and are divided in families and versions. They differ by resolution (i.e., AA-FF, UA-FF and CG-FF), form of the Hamiltonian, bonded and non-bonded parameters, cutoffs, long-range corrections or time step. Their parameterization is either based on high quality QM calculations or done by fitting all parameters to reproduce experimental values such as densities, heats of vaporization, free energy of hydration, membrane structural parameters or Ramachandran plots of protein dihedral angles [START_REF] Jorgensen | Potential energy functions for atomic-level simulations of water and organic and biomolecular systems[END_REF]. Consequently, each force field is more or less accurate in reproducing given experimental values. The OPLS/AA force field was shown to be well-adapted for organic liquids [START_REF] Jorgensen | Potential energy functions for atomic-level simulations of water and organic and biomolecular systems[END_REF][START_REF] Caleman | Force Field Benchmark of Organic Liquids: Density, Enthalpy of Vaporization, Heat Capacities, Surface Tension, Isothermal Compressibility, Volumetric Expansion Coefficient, and Dielectric Constant[END_REF][START_REF] Kaminski | Evaluation and Reparametrization of the OPLS-AA Force Field for Proteins via Comparison with Accurate Quantum Chemical Calculations on Peptides †[END_REF].

For proteins, the recent revisions of AMBER, CHARMM and OPLS/AA force field yield much acceptable results [START_REF] Ponder | Force Fields for Protein Simulations[END_REF][START_REF] Beauchamp | Are Protein Force Fields Getting Better? A Systematic Benchmark on 524 Diverse NMR Measurements[END_REF]. The accurate description of DNA is still a tricky issue and it is even more delicate for RNA [START_REF] Cheatham | Twenty-five years of nucleic acid simulations[END_REF]. The force fields for nucleic acids are still currently under improvements [START_REF] Šponer | Molecular dynamics simulations of G-DNA and perspectives on the simulation of nucleic acid structures[END_REF][START_REF] Condon | Optimization of an AMBER Force Field for the Artificial Nucleic Acid, LNA, and Benchmarking with NMR of L(CAAU)[END_REF][START_REF] Krepl | Reference Simulations of Noncanonical Nucleic Acids with Different χ Variants of the AMBER Force Field: Quadruplex DNA, Quadruplex RNA, and Z-DNA[END_REF][START_REF] Krepl | Can We Execute Stable Microsecond-Scale Atomistic Simulations of Protein-RNA Complexes?[END_REF], even though recent modifications of AMBER force fields provide the most accurate results at the moment. Numerous force fields also exist to simulate lipids and lipid bilayers. Among them, Berger lipids [START_REF] Berger | Molecular dynamics simulations of a fluid bilayer of dipalmitoylphosphatidylcholine at full hydration, constant pressure, and constant temperature[END_REF], the GROMOS family (e.g., 43A1-S3 [START_REF] Chiu | An improved united atom force field for simulation of mixed lipid bilayers[END_REF] and 53A6 [START_REF] Oostenbrink | Validation of the 53A6 GROMOS force field[END_REF]), CHARMM36 [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], GAFFlipids [START_REF] Dickson | GAFFlipid: a General Amber Force Field for the accurate molecular dynamics simulation of phospholipid[END_REF], LIPID11 [START_REF] Skjevik | LIPID11: A Modular Framework for Lipid Simulations Using Amber[END_REF], LIPID14 [START_REF] Dickson | Lipid14: The Amber Lipid Force Field[END_REF] and Slipids [START_REF] Jämbeck | Derivation and Systematic Validation of a Refined All-Atom Force Field for Phosphatidylcholine Lipids[END_REF][START_REF] Jämbeck | An Extension and Further Validation of an All-Atomistic Force Field for Biological Membranes[END_REF][START_REF] Jämbeck | Another Piece of the Membrane Puzzle: Extending Slipids Further[END_REF] can be cited. Some studies benchmarked different properties of lipid force fields, e.g., membrane structural parameters [START_REF] Piggot | Molecular Dynamics Simulations of Phosphatidylcholine Membranes: A Comparative Force Field Study[END_REF][START_REF] Siu | Biomolecular simulations of membranes: Physical properties from different force fields[END_REF]. We evaluated the capacity of force fields to reproduce membrane partitioning and crossing, and found that Slipids is the most accurate force field (see reference [START_REF] Paloncýová | Benchmarking of force fields for molecule-membrane interactions[END_REF] or section IV.1).

Accuracy of a force field is important, but when it comes by introducing new terms to be computed, computational time is increased. Thus, there is a delicate balance between quality of results (e.g., free energy differences, structural parameters) and computational cost.

Polarizable force fields are a successful example of improvement. They allow modifying partial charges along MD simulations as a function of atom environment. When the dynamic description of polar species is crucial (e.g., solvation energy, evaluation of pK a values, MD of ion channels) the polarization is essential [START_REF] Warshel | Polarizable Force Fields: History, Test Cases, and Prospects[END_REF]. However, the computational cost is twice larger than for fixed-charge force fields. This is why non-polarizable force fields are still more commonly employed.

To simulate complex systems including together a lipid bilayer and protein or a protein and DNA, a straightforward solution would be to use the most accurate force field for each part of the system. However, mixing force fields is sometimes simply infeasible, as it involves mixing parameters that are not consistently developed with each other and that are not always designed to use similar potentials. Some force fields are however compatible e.g., Slipids [START_REF] Jämbeck | Derivation and Systematic Validation of a Refined All-Atom Force Field for Phosphatidylcholine Lipids[END_REF][START_REF] Jämbeck | An Extension and Further Validation of an All-Atomistic Force Field for Biological Membranes[END_REF][START_REF] Jämbeck | Another Piece of the Membrane Puzzle: Extending Slipids Further[END_REF], LIPID11 [START_REF] Skjevik | LIPID11: A Modular Framework for Lipid Simulations Using Amber[END_REF], LIPID14 [START_REF] Dickson | Lipid14: The Amber Lipid Force Field[END_REF] and GAFF are generally compatible with Amber force fields, whereas Berger lipids [START_REF] Berger | Molecular dynamics simulations of a fluid bilayer of dipalmitoylphosphatidylcholine at full hydration, constant pressure, and constant temperature[END_REF] is compatible with GROMOS force fields.

II.2.5. Water models

In MD simulations, the solvent is often treated explicitly, i.e., each solvent molecule is described. As water is the solvent of choice for biological studies, several force fields have been developed for water. They differ by the number of sites, i.e., the number of points interacting with each other. For simple models like the Simple Point Charge model (SPC) [START_REF] Berendsen | Interaction Models for Water in Relation to Protein Hydration[END_REF], the Extended SPC (SPC/E) [START_REF] Berendsen | The missing term in effective pair potentials[END_REF] or TIP3P [START_REF] Jorgensen | Comparison of simple potential functions for simulating liquid water[END_REF], three points are defined, each on the center of each atom. For more sophisticated models, additional sites are added. In TIP4P [START_REF] Jorgensen | Comparison of simple potential functions for simulating liquid water[END_REF], the fourth site represents a charge delocalized from the oxygen atom, while in TIP5P [START_REF] Mahoney | A five-site model for liquid water and the reproduction of the density anomaly by rigid, nonpolarizable potential functions[END_REF] they stand for both oxygen lone pairs. Force field compatibility must also be insured with water models. SPC and SPC/E are compatible with GROMOS, whereas TIP3P and TIP4P are compatible with AMBER and OPLS, respectively.

II.2.6. Creation of topologies for small molecules

To describe molecules other than lipids, aminoacids, nucleic acids, carbohydrates and water, by force fields, specific parameters are sometimes required. This can be done 'by hand' using high-level QM calculations or automatized or semi-automatized by producing topologies (i.e., force field parameters) for any type of organic molecules. There are programs producing such topologies for a given family or force field, i.e., PRODRG [START_REF] Schüttelkopf | PRODRG: a tool for high-throughput crystallography of protein-ligand complexes[END_REF] for GROMOS [START_REF] Lemkul | Practical Considerations for Building GROMOS-Compatible Small-Molecule Topologies[END_REF]; Antechamber and GAFF [START_REF] Wang | Development and testing of a general amber force field[END_REF] for Amber; and CGenFF [START_REF] Vanommeslaeghe | CHARMM general force field: A force field for drug-like molecules compatible with the CHARMM allatom additive biological force fields[END_REF][START_REF] Vanommeslaeghe | Automation of the CHARMM General Force Field (CGenFF) I: Bond Perception and Atom Typing[END_REF] or SwissParam [START_REF] Zoete | SwissParam: A fast force field generation tool for small organic molecules[END_REF] for CHARMM. Recently, a similar program was developed for CG models using the MARTINI force field [START_REF] Bereau | Automated Parametrization of the Coarse-Grained Martini Force Field for Small Organic Molecules[END_REF].

Here we focus on producing topologies for the GROMOS force fields using PRODRG. While the parameters provided by this program are rather suitable, it is known that the corresponding partial charges are erroneous [START_REF] Lemkul | Practical Considerations for Building GROMOS-Compatible Small-Molecule Topologies[END_REF]. Thus, partial charges have to be recalculated using QM methods. DFT calculations of the electrostatic potential (ESP) followed by a restrained fit of electrostatic potential (RESP [START_REF] Bayly | A well-behaved electrostatic potential based method using charge restraints for deriving atomic charges: the RESP model[END_REF]) is the method of choice to generate partial charges in AMBER force fields and has also appeared suitable for GROMOS. The Duan (DFT-based) method [START_REF] Duan | A point-charge force field for molecular mechanics simulations of proteins based on condensed-phase quantum mechanical calculations[END_REF] using B3LYP/cc-pVTZ and implicit solvation with ε = 4 should be preferred over the one proposed by Cornell et al. [90] (HF/6-31G(d) in gas). In order to obtain accurate and reproducible partial charges with RESP, the R.E.D III program [START_REF] Dupradeau | tools: advances in RESP and ESP charge derivation and force field library building[END_REF] was used. It allows to average charges on multiple conformations and orientations for a given molecule. Generation of conformers was performed either by a short MD simulation in vacuum followed by conformation clustering, or by using Confab [START_REF] O'boyle | Confab -Systematic generation of diverse low-energy conformers[END_REF], a program generating conformers systematically.

Another point which is underestimated by PRODRG is the influence of aromaticity on dihedral parameters. Thus for some molecules, dihedral angles between conjugated moieties have to be reparameterized (see IV.4.3.1 for details).

II.2.7. Periodic boundary conditions

When simulating the motions of a molecule surrounded by solvent molecules, the size of the system is limited in space, usually in a cubic box. However, a problem arises at the boundaries of this box for the treatment of interactions. Considering vacuum outside the box leads to severe artifacts. One elegant way to solve this problem is to consider a periodic representation, mimicking the same box at each of face of the primary box. Following this procedure, called periodic boundary conditions, a molecule can escape the box by one face but reenter by the opposite face. If the box is large enough for molecules not to interact with their own copies, periodic boundary conditions mimic an infinite system.

II.3. Molecular dynamics

II.3.1. Principle -Integrators

MM allows computing potential energies (𝑉) and forces acting on all particles of a system in a given state. The movement of molecules vs. time (trajectory) can also be obtained by solving Newton's second law of motion (𝐹 = 𝑚𝑎), which can be expressed in its differential form as:

- 𝑑𝑉 𝑑𝒙 𝒊 = 𝑚 𝑖 𝑑 2 𝑥 𝑖 𝑑𝑡 2
where -

𝑑𝑉 𝑑𝒙 𝒊
is the force acting on an atom of mass 𝑚 𝑖 along a coordinate 𝑥 𝑖 .

From a given structure it is possible to calculate potential and force acting on each atom and thus to deduce their speed and positions after a given time step. This constitutes a loop that is repeated at each new position. To accurately describe molecular motion, a time step smaller than the fastest atomic movements must be used. Since the fastest vibration of a hydrogen atom is approximately 13 fs, a 2 fs time step is generally used. Additionally, since biologically relevant processes (e.g. protein rearrangement, DNA folding or permeation of molecules through membranes) occur at least in the range of hundreds of nanoseconds, hundreds of millions of steps are required. This only became possible in the last few years thanks to the new developments of computers, hardware, algorithms and software [START_REF] Klepeis | Long-timescale molecular dynamics simulations of protein structure and function[END_REF].

The algorithm propagating velocities 𝑣 and coordinates 𝑥 along time is the integrator. It solves Newton's equation in its simple form:

𝐹(𝑡) = 𝑚 ∆𝑣 ∆𝑡
The leap-frog algorithm is one of the most common integrators. It computes new coordinates every ∆𝑡 step and the velocities at 𝑡 + 1 2

∆𝑡:

𝑣 (𝑡 + 1 2 ∆𝑡) = 𝑣 (𝑡 - 1 2 ∆𝑡) + ∆𝑡 𝑚 𝐹(𝑡) 𝑥(𝑡 + ∆𝑡) = 𝑥(𝑡) + ∆𝑡 𝑣(𝑡 + 1 2 ∆𝑡)

II.3.2. Temperature regulation

Temperature 𝑇 is the macroscopic result of molecular agitation, and as such, atom velocities 𝑣 𝑖 and kinetic energy 𝐸 𝑘𝑖𝑛 :

𝐸 𝑘𝑖𝑛 = 1 2 ∑ 𝑚 𝑖 𝑣 𝑖 2 𝑁 𝑖=1 = 1 2 𝑁 𝑑𝑓 𝑘𝑇
where 𝑘 is Boltzmann's constant and 𝑁 𝑑𝑓 is the number of degrees of freedom, which usually equals 3×(number of atoms) in MD simulations.

Direct use of MD simulations leads to use the NVE ensemble (constant number of particles, volume and total energy). However, in NVE ensemble, temperature is allowed to fluctuate. In order to reproduce experimental or biological conditions, the NVT ensemble (constant temperature instead of constant total energy, also called canonical ensemble) appears more adapted. In NVT, temperature is thus regulated by scaling velocities. Several algorithms exist among with two are detailed now.

A simple possibility would be to immediately scale the velocities of all atoms as soon as a temperature deviation is measured. However, such rough solution induces non-negligible artifacts. Berendsen's thermostat [START_REF] Berendsen | Transport Properties Computed by Linear Response through Weak Coupling to a Bath[END_REF] corrects temperature divergence smoothly and slowly according to an exponential decay:

𝑑𝑇 𝑑𝑡 = 𝑇 0 -𝑇 𝜏
where 𝑇 0 is the reference temperature and 𝜏 is the decay time constant. It mimics the presence of one external temperature bath, and it is very efficient at relaxing the system to the desired temperature. However, it quenches the variations of the kinetic energy and thus does not sample a proper canonical ensemble.

The Nosé-Hoover's thermostat [START_REF] Nosé | A molecular dynamics method for simulations in the canonical ensemble[END_REF][START_REF] Hoover | Canonical dynamics: Equilibrium phase-space distributions[END_REF] enables proper canonical ensemble simulations by introducing an external "heat bath" variable 𝑝 𝜉 directly into the equation of motion:

𝑑 2 𝑥 𝑖 𝑑𝑡 2 = 𝐹 𝑖 𝑚 𝑖 - 𝑝 𝜉 𝑄 𝑑𝑥 𝑖 𝑑𝑡
The variable 𝑝 𝜉 varies as according to:

𝑑𝑝 𝜉 𝑑𝑡 = (𝑇 -𝑇 0 )
and:

𝑄 = 𝜏 𝑇 2 𝑇 0 4𝜋 2
As opposed to Berendsen's thermostat which allows quick convergence to the reference temperature 𝑇 0 , Nosé-Hoover's thermostat induces temperature oscillations that slowly converge to 𝑇 0 . Therefore, the latter is recommended only for pre-equilibrated simulations. Additionally, Nosé-Hoover's thermostat can be coupled to its own external heat bath controlled by another Nosé-Hoover's thermostat. These so-called Nosé-Hoover chains allow a better exploration of the phase-space.

An interesting possibility offered by MD simulations is the ability to define various thermostats and reference temperatures for different parts of a given system. Such a procedure can be used either (i) to ensure that each group is simulated at the proper temperature, or (ii) to heat or cool down a part of the system e.g. to heat a solute while preventing solvent boiling.

Another problem is how to define a temperature while starting with a set of fixed atomic coordinates. In the leap-frog integrator, the velocities of atoms at time 𝑡 0 -1 2 ∆𝑡 can be generated according to a Maxwell-Boltzmann distribution:

𝑝(𝑣 𝑖 ) = √ 𝑚 𝑖 2𝜋𝑘𝑇 exp (- 𝑚 𝑖 𝑣 𝑖 2 2𝑘𝑇
)

II.3.3. Pressure regulation

Pressure coupling is very similar to temperature coupling. The relationship between temperature and velocities is replaced by the relationship between pressure and box size scaling. Moreover temperature bath is replaced by a "pressure bath". The Berendsen's barostat scales the coordinates according to an exponential decay relationship. The counterpart of Noosé-Hoover's thermostat is Parrinello-Rahman's barostat [START_REF] Parrinello | Polymorphic transitions in single crystals: A new molecular dynamics method[END_REF][START_REF] Nosé | Constant pressure molecular dynamics for molecular systems[END_REF] that can be used to generate a proper NPT ensemble.

II.4. Potential of mean force

When a system is at thermodynamic equilibrium between two states A and B, the equilibrium is described by the constant 𝐾 𝑒𝑞 as the ratio between the populations of states A and B:

𝐾 𝑒𝑞 = [𝐴] [𝐵]
The free energy difference ∆𝐺 between the two states is related to 𝐾 𝑒𝑞 by:

∆𝐺 = -𝑘 𝐵 𝑇 ln 𝐾 𝑒𝑞
where 𝑘 𝐵 is Boltzmann's constant and 𝑇 the absolute temperature. Therefore, the sign of ∆𝐺 determines which state is favored, classically describing thermodynamic balances.

There are two methods to theoretically evaluate free energy differences. The first one is to obtain the relative populations of the states A and B (i.e.,

𝑃 𝐴 𝑃 𝐵

) by statistical analysis of a free unbiased MD simulation. However, obtaining relevant statistical data on the relative populations suppose that the system is ergodic, i.e. time average is equivalent to ensemble average. In other words, this method could be applied if every state of the system is sufficiently sampled during the time scale of the free unbiased simulation. Although it can be possible for fast events such as bond and angle vibrations, the time scales to simulate a protein folding or membrane crossing by a drug are several orders of magnitude longer. Therefore, with the current computer resources it is impossible to sample all states of a membrane crossing within affordable computer time.

As a consequence, rare events occurrences must be accelerated to evaluate free energy differences. To do so, resolution can be lowered (e.g., CGFF) or temperature increased. However, this could lead to lose precision or displace equilibrium, respectively. The general approach to tackle this issue is to calculate the work along the reversible path connecting the two states. Then, integrating the mean forces 〈𝐹(𝑥)〉 𝑡 along this path allows reconstruction of the potential of mean force (PMF, i.e., free energy profile):

∆𝐺(𝑥) = -∫〈𝐹(𝑥)〉 𝑡 𝑑𝑥
A profile not only provides ∆𝐺 between two states but also the description of the transformation between these two states.

To calculate the work along the path, one has to i) know this path in the multi-dimensional space and ii) force the system to sample the states along this path. Various methods were developed in this purpose. In the following sections, the methods used in this work are detailed.

II.4.1. Pulling

When drugs cross lipid bilayers, the free energy is obtained vs the distance between drug and bilayer center along membrane normal (𝑧 axis). In order to sample the states along this path in a reasonable computational time, the drug has to be forced to fully explore this path.

One simple way is to create starting structures where the drug is manually placed at different positions along the 𝑧 axis by translating coordinates. However, steric clashes can arise from overlapping atoms of drug, membrane or solvent. Another way is to use a pulling MD simulation, in which the drug is slowly forced by a bias potential to move along the 𝑧 axis. To allow some flexibility of drug movements, the bias potential is often described by a harmonic potential:

𝑉 𝑝𝑢𝑙𝑙𝑖𝑛𝑔 = 𝑘(𝑧 𝑖 -𝑧 0 ) 2
The drug velocity along the 𝑧 axis should be as small as possible to avoid deformations at the lipid bilayer surface and to allow sufficient relaxation.

After pulling which forces membrane crossing, a set of positions is obtained at regular intervals. Independent MD simulations are then assessed for each position allowing sampling the close conformational space within a window. At each step, a bias potential is applied to maintain the drug in the window. The potential can be either a harmonic restraint in umbrella sampling or a constraint in z-constraint simulations.

II.4.2. Umbrella sampling

A harmonic restraint -also called umbrella potential, 𝑈 ′ (𝑧) -is applied to maintain drug in windows, close to its initial value. After sufficient sampling for each window (usually several tens of nanoseconds), the biased probabilities 𝑃 ′ (𝑧) to find the system in a given state along the path are computed. The unbiased free energy is then reconstructed:

𝐺(𝑧) = -𝑘 𝐵 𝑇 ln 𝑃 ′ (𝑧) -𝑈 ′ (𝑧) + 𝐹
where 𝐹 is the free energy shift to unbias the free energy. It should be noted that it is undetermined, it depends directly on 𝑈 ′ (𝑧) and thus it is different for each window.

The weighted histogram analysis method (WHAM) [START_REF] Kumar | THE weighted histogram analysis method for free-energy calculations on biomolecules. I. The method[END_REF] determines the values of 𝐹 by combining the biased probabilities 𝑃 ′ (𝑧) (also called histograms) of all windows. This algorithm therefore reconstructs the unbiased PMF. From a practical viewpoint, window spacing, sampling time, width and strength of the bias affect quality of the final PMF and should be determined beforehand.

II.4.3. Z-constraint method

The z-constraint method is similar to umbrella sampling. It is however simpler as the restraint is replaced by a constraint, removing the degrees of freedom for the drug to move in the 𝑧 dimension, while it is still free to move in the 𝑥𝑦 plane. The force to constrain the drug at its original position is monitored in each window. After sufficient sampling, the forces are averaged for each window and then integrated along the 𝑧 axis, from outside to the center of the membrane [100]:

∆𝐺(𝑧) = -∫ 〈𝐹(𝑧)〉 𝑡 𝑧 𝑜𝑢𝑡𝑠𝑖𝑑𝑒 𝑑𝑧
Window spacing and sampling time are also of importance, although larger spacing and smaller sampling times yield comparable results to umbrella sampling [100]. It is particularly adapted to evaluate free energy profile of membrane crossing, see Chapter III.

II.4.4. Metadynamics

Metadynamics uses a different approach to bias MD simulations. Instead of having several simulations (windows) along the 𝑧 axis, the bias is constructed dynamically along one MD simulation [101,102]. At regular intervals, a repulsive Gaussian potential is added at the reaction coordinates 𝑠 𝑖 (𝒓), so that the already visited coordinates are discouraged. Several reaction coordinates can be considered simultaneously, they are called collective variables (CV). The total repulsive potential 𝑉(𝒔, 𝑡) at time 𝑡 is then equal to the sum of all repulsive Gaussian functions:

𝑉(𝒔, 𝑡) = 𝜔 ∑ exp (-∑ (𝑠 𝑖 (𝒓) -𝑠 𝑖 (𝒓(𝑡 ′ ))) 2 2𝜎 𝑖 2 𝑑 𝑖=1 ) 𝑡 ′ =𝜏 𝐺 , 2𝜏 𝐺 ,… 𝑡 ′ < 𝑡
where 𝜎 𝑖 is the Gaussian width of the 𝑖th CV, 𝜔 is the Gaussian height, 𝑑 is the number of CV, and 𝜏 𝐺 is the rate at which the Gaussian functions are added.

When the system initially lies in an energy minimum on the potential energy surface (PES), the bias potential gradually "fills" this minimum until the nearest energy barrier is passed. Once all minima of the explored space are filled with repulsive bias, free diffusion in the explored space is observed. At the end of the procedure, the bias potential represents the negative of the estimated free energy profile:

𝐺(𝒔, 𝑡) ≈ -𝑉(𝒔, 𝑡)
Despite the straightforward theoretical background of this methodology, practical use of metadynamics is delicate [101]. The first reason is that 𝜎 𝑖 , 𝜔 and 𝜏 𝐺 must be chosen prior to the MD simulation while these parameters drastically influence outcomes. Too high 𝜎 𝑖 and 𝜔 values imply loss of precision, whereas too low 𝜎 𝑖 and high 𝜏 𝐺 valueslead to very slow convergence. Moreover, low 𝜏 𝐺 values (i.e. frequent deposition) introduce artifacts because the system does not have sufficient time to relax. A correct estimate of 𝜎 𝑖 is generally half the standard deviation during an unbiased MD simulation. As for the other parameters, a determination by a trial-and-error procedure is often necessary. The second reason of delicate usage of metadynamics is the difficulty to choose one or more CV to get reliable free energy estimates. Although being a common problem with other bias methods such as umbrella sampling, it is particularly true in metadynamics, i.e. in the case of free diffusion is never reached even with optimized parameters or when hysteresis are observed in the free energy profiles. Bad choices of CV are often responsible for these drawbacks. In these cases, the optimal energy path is on a multi-dimensional PES, and one must include more CV to follow it. For instance, in the case of a drug crossing the membrane, the drugmembrane distance is one evident CV, but drug rotation may play an important role e.g. in the surrounding of phospholipid head group region. CV should allow connecting initial and final states with intermediates, and they should describe slow events [101].

It is clear that the time required sampling the PES increases exponentially with the number of its dimensions, i.e., the number of CV in metadynamics. This method performs well up to 2 to 3 CV in regard to other methods [101]. Depending on the system, or if more CV are required, the sampling may be inefficient in a reasonable computational time. In these cases, enhanced sampling methods can be coupled to metadynamics. For instance, multiple walkers method runs several metadynamic simulations in parallel that share a common bias potential; parallel tempering method is similar to replica exchange, where parallel simulations exchange coordinates at various temperatures; and bias exchange metadynamics also runs several replicas that explore each CV and that exchange biases on CV at regular intervals.

II.4.5. COSMOmic

The COSMOmic divides membrane and surrounding solvent in layers (usually 1 Å thick) along the 𝑧 axis (perpendicular to membrane surface). For each 𝑛 𝑡ℎ layer, 𝑚 molecule orientations are evaluated to properly sample the PES. The partition function of a given drug, say 𝑋, in a membrane M is:

𝑍 𝑀 𝑋 = ∑ ∑ exp (- 𝜇 𝑀 𝑋 (𝑟 𝑖 , 𝒅 𝑗 ) 𝑘 𝐵 𝑇 ) 𝑚 𝑗=1 𝑛 𝑖=1
where 𝑟 𝑖 is the drug position and 𝒅 𝑗 its orientation.

The probability to find the drug in the layer 𝑖 is then:

𝑝 𝑀 𝑋 (𝑟 𝑖 ) = 𝑍 𝑀 𝑋 (𝑟 𝑖 )
𝑍 𝑀 𝑋 and the free energy profile:

𝐺 𝑀 𝑋 (𝑟 𝑖 ) = -𝑘 𝐵 𝑇 ln 𝑝 𝑀 𝑋 (𝑟 𝑖 )
In summary, COSMOmic presents an elegant alternative to the computationally expensive methods based on MD simulations. Only one DFT calculation per molecule specie (e.g., drug, phospholipid, solvent) is required; the subsequent evaluation of the free energy profiles only requires minutes of computational time as opposed to weeks or months with MD-based methods. We have showed that it is as precise as the best force fields to evaluate water/membrane partition coefficients (see section IV.1). Therefore, COSMOmic is suitable for screening drug partition coefficient in fluid lipid bilayer membranes. Although it is possible to evaluate membranes that are constituted of lipid mixtures, one should keep in mind that COSMOmic is designed under the assumption of fluid phase, and that entropic contributions of liquid ordered or gel phase can be poorly evaluated as the environment is described implicitly.

II.5. Analysis of molecular dynamics simulations

II.5.1. Sampling times

Coordinates and velocities of tens of thousands of atoms for millions of steps can generate an incredible amount of data. These data have to be carefully analyzed to make sure that no unrealistic behavior happen and to extract valuable information. Some results of the simulation can be related to macroscopic or microscopic experimental data. However, to obtain relevant averaged values of physical parameters, a huge amount of data may not be sufficient and sometimes simulations have to be prolonged.

The time required to properly sample a rare event depends on many parameters, in particular structural re-arrangements requiring free energy of activation. For instance, i) large molecules usually reach their equilibrium position in the bilayer much slower than small compounds; ii) crossing lipid bilayer in gel or liquid ordered phase is longer than in fluid phase.

It is possible to estimate the permeability 𝑃 (in cm/s) of molecules through membranes thanks to the following equation:

𝑃 = 𝐾𝐷 ∆ 𝑥
where 𝐾 is the partition coefficient, 𝐷 the diffusion coefficient (in cm 2 /s) and ∆ 𝑥 the membrane thickness [105].

II.5.2. Distances -positions

One of the most valuable results yielded by MD simulations of drugs interacting with lipid bilayers is the drug equilibrium position. Namely, we usually estimate distance along the 𝑧 axis between drug and center of mass of the membrane * . An interesting possibility is to evaluate the position of either the center of mass or of specific moieties of the drug. Then, the difference of position of two moieties is related to the general orientation of the drug along the 𝑧 axis.

When drug reach a proper equilibration position, its averaged value must be weighted by its standard deviation, which is often for small drugs about 2-3 Å. Naturally, it is also possible to follow other positions or distances of interest, whether they are projected in one, two or three dimensions.

Positions obtained from MD simulations can be related to various experimental techniques [106]. Among them are (i) small angle X-ray and neutron scattering; (ii) electron spin * Even though center of mass motion of the system is removed periodically, the COM of the membrane still fluctuates in the 𝑧 direction. Thus, in order to properly evaluate the distance distribution, one must take care not to calculate a 𝑧 distance by the simple difference of time averaged positions but to evaluate the average of 𝑧 distances along time. Evaluation of distances can also be indirectly related to other observable quantities such as radial distribution functions, quantification of hydrogen bonding or root mean square deviation (RMSD) from a reference structure.

II.5.3. Angles -orientations

The evaluation of the distribution of angles between two vectors in MD simulations can be performed for a variety of cases, e.g., (i) for angles between 3 consecutive atoms, (ii) for dihedral angles, (iii) for orientation of bonds or molecules in regard to given axes or (iv) for plane normal vectors to evaluate coplanarity. Only few experimental methods can assess orientation. Among them, 2 H-NMR can obtain order parameters and from that estimate orientation (see section II.5.4.3).

II.5.4. Membrane parameters

To ensure that a lipid bilayer model reflects experiment requires comparison of a maximum of membrane structural parameters. Bilayer thickness, area per lipid, and order parameter are among the most that should be correctly predicted by a model.

II.5.4.1 Bilayer thickness

Bilayer thickness is one of the characters that characterize lipid phase. Although its principle is simple, its precise definition requires a particular attention from a practical point of view, because there is no obvious boundary between lipids and solvent. Hence, several definitions of lipid bilayer thickness have been proposed [109].

Head-head thickness (𝐷 𝐻𝐻 ) is defined as the distance between peaks of electron density, that can be obtained by MD simulations or X-ray experiments [109][110][111] and should correspond to the distance between head groups of both layers.

Similarly, phosphate-phosphate thickness can be defined as the distance between the COM of phosphate moieties. Hydrocarbon chain thickness 𝐷 𝐶 is also valuable, however experimentally it is only accessible indirectly by evaluation of partial head group thickness 𝐷 𝐻1 [109,111].

𝐷 𝐶 = 𝑉 𝐶 𝐴 = 𝐷 𝐻1 - 𝐷 𝐻𝐻 2
where 𝐴 is the area per lipid. 𝐷 𝐻1 is the measure of the distance between phosphate and the average lipid chain boundary, which is estimated using molecular modelling [109].

Likewise, the steric bilayer thickness 𝐷 𝐵 ′ is defined as a function of steric head group thickness 𝐷 𝐻 ′ .

𝐷 𝐵 ′ = 2(𝐷 𝐶 + 𝐷 𝐻 ′ )
On the other hand, Luzzati thickness 𝐷 𝐵 is based on volume per lipid 𝑉 𝐿 and area per lipid 𝐴:

𝐷 𝐵 = 2 𝑉 𝐿 𝐴
In MD simulations, 𝐷 𝐵 can be evaluated as a function of box volume 𝑉 𝐵 and water phase volume:

𝐷 𝐵 = 2 (𝑉 𝐵 -𝑛 𝑊 𝑉 1𝑊 ) 𝐴
where 𝑛 𝑊 is the number of water molecules and 𝑉 1𝑊 is the volume of one water molecule. Luzzati thickness can be evaluated in neutron diffraction experiments by the high contrast between protonated lipids and deuterated water [111].

All these thicknesses are all valid and different definitions. Luzzati thickness is convenient to calculate in MD, however experimental methods can measure of thickness definitions. Thus, it is important to correlate thicknesses matching the same definition.

II.5.4.2 Area per lipid

In MD simulations, area per lipid 𝐴 𝐿 is simply calculated by dividing the 𝑥𝑦 area of the box by the number of lipids in one layer. It can be derived experimentally by measurement of Luzzati thickness 𝐷 𝐵 and precise evaluation of 𝑉 𝐿 .

It should be noted that aside from bilayer composition and temperature differences, there are large discrepancies between experiments (Figure 12). Disparities between force fields have also been observed, mostly because some were fitted to reproduce non-consistent experimental data. Some force fields include constant surface tension, introducing an additional potential to fix 𝐴 to a desired value. However, it should be noted that the discrepancy between all measurements, experimental and theoretical, has been reduced or elucidated in the past years [112]. Recent force fields such as Slipids have been capable to reproduce 𝐴 very accurately for a wide range of lipid bilayers and temperatures [START_REF] Jämbeck | Derivation and Systematic Validation of a Refined All-Atom Force Field for Phosphatidylcholine Lipids[END_REF][START_REF] Jämbeck | An Extension and Further Validation of an All-Atomistic Force Field for Biological Membranes[END_REF][START_REF] Jämbeck | Another Piece of the Membrane Puzzle: Extending Slipids Further[END_REF].

II.5.4.3 Order parameters

Order parameter is a dimensionless quantity related to both a preferential orientation of an angle 𝜃 and deviations around this orientation. It is defined as:

𝑆 = 3 2 〈cos 2 𝜃〉 - 1 2
where brackets mean time average. When a molecule has no preferential orientation, as for a molecule tumbling in homogeneous solvent, the order parameter S value equals 0.

However, it important to stress that a 0 value also corresponds to the perfect order and 𝜃 = 54.7356° (i.e., the so-called magic angle). This angle is obtained by solving the above equation with 𝑆 = 0 and considering no average (i.e., perfect order). Therefore, an orientation can only be translated in terms of order parameter from dynamical collection of data obtained by a sufficient sampling. Conversely, an orientation can only be obtained from an order parameter in the case of a perfect ordering, with a sign uncertainty arising from the square root of cos 2 𝜃.

Order parameters are a valuable tool in the study of lipid bilayer membranes as they can be precisely evaluated by solid state NMR (ssNMR) of deuterated compounds (lipids or drugs).

The order parameters obtained are related to the angle 𝜃 between carbon-deuterium bonds and the axis of the magnetic field. Quadrupolar splitting ∆ 𝜈𝑞 is measured and it is related to 𝑆 by the following relationship:

∆ 𝜈𝑞 = 3 2 𝑆 𝑒 2 𝑞𝑄 ℎ
where

𝑒 2 𝑞𝑄
ℎ is the quadrupolar coupling constant depending on the atom type (e.g., aliphatic, aromatic or polar C-D bond).

In this work we have used order parameter in two applications: (i) to confirm the order in lipid chains, which is related to phase; and (ii) to evaluate the orientation of deuterated drugs in lipid bilayer (see section IV.5).

Chapter III. Interaction of drugs with membranes

III.1. Introduction

A comprehensive understanding of mechanisms of action of drugs is crucial to control their efficacy, side effects, and toxicity. This requires a complete description, at an atomic scale, of all intermolecular interactions between drugs and biological targets. Such atomistic description has become a crucial step in drug design aiming at increasing activity while reducing toxicity. Drug targets are as diverse as (i) the gastro-intestinal tract; (ii) other organs; (iii) extra-or intracellular compartments; (iv) proteins, nucleic acids, carbohydrates or lipids. Interactions with biological membranes are a key step in pharmacology, as they directly affect drug delivery and accumulation. Models mimicking cell membranes are thus of great interest to predict both pharmacodynamics and pharmacokinetics behaviors. Depending on drugs and many environmental conditions, membrane crossing can occur by passive diffusion (simple diffusion through lipid bilayers following stochastic events to equilibrate concentration gradients); facilitated diffusion (e.g. through pores made of transmembrane proteins); active transport (membrane proteins that influx or efflux compounds often via large conformational changes which can be activated e.g. by ATP binding and hydrolysis). It should be noted that this concerns not only drugs but also any xenobiotic including natural products from diet (polyphenols, vitamins) or compounds synthetized by cells (hormones…).

There is a gamut of experimental methods that tackle xenobiotic-membrane interactions including advanced fluorescence spectroscopy and microscopy, solid-state nuclear magnetic resonance (SSNMR), surface plasmon resonance (SPR), and neutron scattering. Most of the time these methods are employed with liposomes of various sizes for in vitro assays. Applying these biophysical techniques often provide fragmented information on drug penetration, partioning and orientation. For instance, average location in the inner or outer parts of lipid bilayers, diffusion coefficients, partition coefficients, membrane permeability, or partial and indirect information on orientation. Large screening is often virtually impossible because the experiments take long and are rather expensive. Alternatively, in silico molecular modelling has gained much interest over the past decade. Although limited to simple models since their conceptualization in the XX th century, relevant biological models have been developed and successfully used in the last years [113][114][115]. The exponential growth of computing power is also driving the development towards accuracy at reasonable computational time. Molecular dynamics (MD) simulations of lipid bilayer membranes evaluate drug-membrane interaction at both atomic and femtosecond resolutions, which is unreached by experimental methods. The atomic resolution allows a precise description of drug position and orientation in membrane, highlighting intermolecular interactions driving penetration; the femtosecond resolution allows correct evaluation of thermodynamic quantities as far as the system is sampled over a sufficiently long time scale. The pioneer studies that investigated interactions of small solutes with lipid bilayer membranes were published in the mid 1990's on benzene [116,117] and water [118] permeation. Since then, the number of publications has rocketed and we are convinced that the constant improvement of the computational facilities as well as the methods of calculation is paving the way towards a systematic usage of MD simulations in drug discovery [119].

For drug-protein interactions in pharmacology one could refer to the review written by Salsbury [120]; whereas influence of drugs on chemical or physical properties of lipid bilayers have been thoroughly reviewed by Kopeć et al. [121]. In this review we would like to stress how the recent advances in MD simulations of drug interaction with membranes allow considering MD as a new pharmacological tool. One of the main advantages of MD simulations is the ability to describe the preferred location and orientation of drugs in membranes at an atomic resolution. Positions and orientation have been successfully predicted for a wide variety of drugs being often compared and in agreement with experimental data obtained on biomimetic membrane models [121]. More than providing locations and orientations, MD simulations enable evaluation of Gibbs energy profiles along a chosen direction, the z-axis perpendicular to membrane surface allows in particular to follow membrane crossing. Membrane affinity and partition coefficient can be extrapolated from these profiles, as corresponding to the Gibbs energy difference (∆𝐺 𝑝𝑎𝑟𝑡 ‡ ) between the bulk water phase and the minimum Gibbs energy in the lipid bilayer. Membrane permeability that is related to passive diffusion can also be evaluated from these profiles as being the highest Gibbs energy barrier to cross membrane (∆𝐺 𝑝𝑒𝑛 ‡ ), which is most of the time the energy required to cross membrane center (see section II.5.1 for more details).

III.2. Drug penetration in lipid bilayer

III.2.1. Anesthetics

Anesthetics are known to act mainly by disrupting the physical properties of lipid bilayers. A distinction has to be made between local anesthetics (LA) and general anesthetics (GA) although similarities in the mechanism of action can be observed in some cases.

LAs, including articaine, lidocaine, prilocaine, tetracaine or phenytoin, bear ionizable amino groups and may exist in both charged or uncharged states at physiological pH. Location of LA in lipid bilayer membrane was extensively studied with free MD simulations and Gibbs energy profiles [122][123][124][125][126][127][128][129][130][131][132]. All charged forms partition in the polar head group region, in contact with water molecules. Conversely, the uncharged forms penetrate bilayers, preferentially locate below this region head groups, and cross membrane by passive diffusion [122][123][124][125][126][127][128][129][130][131][132]. For instance, the (∆𝐺 𝑝𝑒𝑛 ‡ ) values of benzocaine and phenytoin were ca. 5 and 10 kcal mol -1 , respectively, allowing a relatively fast passive diffusion (Figure 14A). Even though the location of LA in membrane is well rationalized, their mechanism of action is still under debate. The first hypothesis was a direct inhibition of voltage-sensitive sodium and potassium channels [133][134][135], however it has clearly appeared that their interaction with lipid bilayer participate to the mechanism of action. As shown both experimentally and passive diffusion facilitated diffusion active transport endocytosis theoretically by MD simulations, LA increase lipid bilayer fluidity and decrease lipid order (Figure 14) [122,127,130,132]. MD simulations succeeded at rationalizing fluidity increase, as attributed to increase of the dipole electrostatic potential and decrease of order parameters in the lipid core [128,132]. Membrane modifications associated to LA lipid bilayer penetration most likely occur in the surrounding of ion channels, thus affecting ion exchanges. Another interesting mechanism suggested that LAs affect nerve pulse propagation, related to a thermodynamic soliton (i.e., a single impulsion) through the neuron membrane [136,137]. Decrease in membrane fluidity due to LA lowers compressibility and dissipate solitons, stopping transmission of neuronal impulse. GAs constitute a family that gather a wide range of chemical structures having similar effects. It includes small gases as neon or NO; fluorinated molecules (e.g., halothane, desflurane); amphiphilic compounds (e.g., benzodiazepines, propofol, ketamine); or 1-alkanols. As for LAs, the mechanism of action of GAs is controversial and should be elucidated at a molecular level [138]. Although it is accepted that GA action involve post-synaptic ligandgated channels [139], the exact mechanism by which GA affect these channels is a matter of discussion. Two main hypotheses have been investigated, namely the receptor and the membrane hypotheses.

The receptor hypothesis proposes that some GAs can bind to ion channels and receptors in the neuronal membrane. It is supported by the fact that one isoflurane isomer is more active than its enantiomer, even though their membrane affinities are identical [138], thus refuting a membrane-only hypothesis. Therefore, MD simulation studies have been performed to identify binding sites of GA. For instance, in agreement with experiments, MD simulations showed that n-alcohols and inhaled anesthetics modulate the potassium channel Shaw-2 by allosteric effects [128]. Similarly, the binding sites of isoflurane were identified inside voltagegated sodium channels [140]; the authors hypothesized that the binding sites were identical to that suggested for LAs. As another example based on QM/MM calculations and MD simulations, xenon was predicted to inhibit competitively NMDA receptors [141]. Other receptors were experimentally identified as related to GA action, however the lack of crystalline structure prevent proper docking and MD studies [138].

The membrane-mediated mechanism hypothesis has also been extensively explored. As for LA, GA activity is correlated to their greater affinity to oil than to water, as stated by the Meyer-Overton rule. Although this rule still stands for most of GAs, some exceptions have been reported (e.g., 1-alkanols are more potent than predicted by this rule) [138]. A possible mechanism establishing the relationship between membrane affinity and ion channel function is the modification of lateral pressure profile of lipid bilayer [139]. Indeed GA insert in the polar head group region of lipid bilayer. Doing so it increases lateral pressure at the membrane-water interface. Such surface tension modifications can induce significant to dramatic conformational changes of ion channels, possibly closing extremities of the pore [139]. For instance, MD simulations confirmed X-ray diffraction studies showing that, at therapeutic concentration, ketamine inserts in membrane at the lipid/water interface without affecting neither membrane thickness nor area per lipid but inducing significant changes in lateral pressure profile that could affect ion channels [142].

Interestingly, there exists a clinical phenomenon observed in anesthesia that is called pressure reversal. It is defined as the cessation of anesthesia by hyperbaric pressure, typically between 80 and 200 atm. The activity of GA being related to two biophysical properties, namely their affinity for membrane and the ambient pressure, it is tempting to relate these pressure effects to lateral membrane pressure. Despite the fact that pressure reversal was also observed for non-GA compounds, it has been consistently described for most of GAs. Several MD studies investigated the relationship between lateral pressure and pressure reversal to link up both effects at the molecular level. For instance, xenon was shown to disorder lipid bilayers and to increase area per lipid at atmospheric pressure [143]. At high pressures xenon location is restricted to membrane center and lipids are packed and ordered as in the absence of xenon [143,144]. Thus, the disordering effect of xenon disappears with pressure increase, in correlation with pressure reversal. Similar mechanisms were recently simulated for chloroform, halothane, diethyl ether and enflurane [145]. By comparing MD simulations at low and high pressure in the presence or absence of GA molecules, Fábián et al. ruled out some molecular descriptors as factors rationalizing pressure reversal such as location in membranes of some GA moieties or orientation of lipid head groups and tails [145]. Again it appeared that the GA global location in membrane was influenced by pressure. As opposed to xenon, these compounds were found to locate in the middle of membrane at atmospheric pressure and below polar head groups at high pressure [145]. Alternatively, GA aggregation in membrane has been proposed to participate in pressure reversal. Indeed, MD simulations revealed that under the conditions of pressure reversal, halothane molecules could aggregate [146,147]. The authors hypothesized that aggregation may reduce proportion of free halothane molecules for receptor binding, thus explaining the lower activity of GA under these conditions.

III.2.2. β-blockers

Beside their β-adrenergic blocking acitivities, β-blockers can also affect lipid bilayer properties, especially the non-selective β-blockers such as alprenolol, oxprenolol and propranolol. This mechanism of action is known as non-specific membrane effect, and it results in anesthetic [148] . This result agrees with the modification of lateral pressure observed for LAs, and could explain the anesthetic effects of β-blockers. However, this mechanism still requires confirmation, in particular paying much attention to i) differences that could be observed on membrane modification from one β-blocker to another, and ii) membrane composition.

III.2.3. Non-steroidal anti-inflammatory drugs

Non-steroidal anti-inflammatory drugs (NSAIDs) are widely used drugs acting by inhibiting cyclooxygenases (COX). The most common drugs in this family are ibuprofen, aspirin and naproxen. Their main side effect is gastro-intestinal (GI) toxicity, in particular GI ulceration.

Several studies highlighted that this toxicity is not COX-related, but that it could result from perturbations of the phospholipid barrier of the mucosa [151]. Interaction between NSAIDs and lipid bilayers is mainly driven by the fact that they all bear a carboxylic acid moiety that can be deprotonated according to pH (e.g. pK a values of ibuprofen, aspirin and naproxen are 4.5-5.2, 3.5 and 4.2, respectively). The position of NSAIDs in lipid bilayers was evaluated experimentally. Considering physiological pH, the negatively charged forms of NSAIDs are predominant and have been seen to locate just below the polar head groups of lipid bilayers as seen theoretically and experimentally [152,153]. MD simulations were also carried out on the neutral form. In this case, NSAIDs were predicted much deeper in bilayer, lying between lipid chains close to the bilayer center [152,[154][155][156]. Although the deprotonated form predominates in water, reprotonation events are likely to occur in the polar head group region, which would allow the subsequent neutral form to relocate deep in the bilayer. When being inserted in between lipid chains, deprotonation is unlikely due to absence of water molecules. Therefore, crossing of the bilayer core most likely occurs for NSAIDs in their neutral form.

Although no clear mechanism of GI toxicity of NSAIDs has been identified yet, several studies have suggested relationship with their ability to intercalate in between phospholipids.

MD simulations have highlighted perturbation of bilayer structure by various NSAIDs, namely they were shown to induce membrane thinning and fluidizing effects [152,156,157]. Such alteration may increase membrane permeability to H + , which may participate in rationalizing GI toxicity (Figure 15) [157].

III.2.4. Antioxidants

Antioxidants have been extensively studied for their beneficial effects on human health.

Although the exact role of dietary antioxidants and even more antioxidant supplementation is still under debate, their application in cosmetics and food preservation is clear. From a medical viewpoint, only organ conservation is requiring antioxidant usage. Indeed in transplantation, the organ is subjected to severe damages induced by ischemia-reperfusion that should be limited as much as possible. Among other processes, antioxidants are capable of inhibiting lipid peroxidation (LPO). LPO is a chain reaction that degrades cell membrane bilayer structure, endangering cell survival [158]. This reaction can be inhibited by antioxidants thanks to their capacity to (i) scavenge free radicals, and (ii) to locate deep inside lipid bilayer membrane where LPO occurs. Whereas theoretical prediction of the former property requires QM calculations [159-162], the latter can be theoretically tackled by MD simulations. Polyphenols, as the prototypical quercetin antioxidant, were shown to locate below phospholipids' head groups, where they can inhibit the initiation stage of LPO. These theoretical results agreed with experimental studies [163][164][165][166].

Methylation or hydrophobic moieties drive antioxidants deeper penetration in membrane, allowing a better antioxidant activity during the propagation stage of LPO. This is particularly exemplified by α-tocopherol having one phenolic OH group, adjacent methyl groups and a long lipid tail that penetrates deeper lipid bilayers than most of polyphenols. Its active OH group locates underneath polar head groups but the rest of the structure lies in between lipid chains [167]. Additionally, it is able to 'flip-flop' from one to the other leaflet, increasing contact with lipid chains [167]. Another example is lipocarbazole, a bacterial compound, bearing a similar chemical structure than α-tocopherol i.e. a polar phenolic OH group and a lipid tail [168]. Argenteane, a natural antioxidant from nutmeg, also exemplifies the role of the methyl group, in the guaiacyl moiety, combined with an apolar linkage [169]. For these three compounds, the combination of their efficient free radical scavenging capacity and their ideal location in lipid bilayers (i.e., in contact with lipid chains) allow an efficient LPO inhibition (see section IV.2, IV.3, and IV.4 for details).

Conversely, polar groups drive antioxidants in polar head groups of lipid bilayers. For instance, catechins bear several OH groups and were clearly shown to partition in between phospholipid head groups; this interaction being driven by H-bonding [170,171]. Ascorbic acid is also a very common antioxidant that partitions outside lipid bilayer, in contact with the water phase [167].

MD simulations have appeared powerful to predict positioning of antioxidants, which is directly related to their biological activities. However, membrane complexity should be systematically investigated. For instance α-tocopherol's depth of penetration in bilayer strongly dependent of lipid composition [172][173][174]. Also collaborative effects between antioxidants are likely to occur inside lipid bilayer and may increase their total antioxidant activity [175,176]. MD simulations have allowed better understanding of synergism between vitamin E, vitamin C and polyphenols, which may occur within non-covalent complexes that are formed in between the lipid chains just below the polar head group region [167]. The existence of such non-covalent antioxidant association was predicted theoretically and further confirmed by fluorescence quenching of vitamin E in the presence of polyphenols [167].

III.3. Limitations and perspectives

MD simulations have definitely become a pharmacological tool, supporting the fragmented knowledge on drug-membrane interactions therefore making possible rationalization of action, bioavailability and toxicity at an atomistic level. Thanks to a dramatic increase of computer power over the past decades, MD simulations have become efficient at predicting partition coefficients of small drugs in simple bilayers; passive diffusion coefficients; drug mechanisms of action in membranes. Although MD simulations cannot evaluate all possible mechanisms for a single drug, the simulations with lipid bilayer models have been successfully used to address some underlying molecular mechanisms, so paving the way to a global understanding of their actions.

Although the MD-based predictions agree most of the time with experimental data performed on similar membranes, there are still a series of drawbacks that must be carefully considered if one aim at predicting behaviors under actual biological conditions. First, a major limitation of MD simulations is the time of sampling, inherent to time resolution (usually 2 fs for all-atom simulations), system size and affordable computing time. Two decades ago the Gibbs energy profiles were limited to very small molecules within a few ns time-scale for sampling [118]; this has considerably be improve and the time scale available nowadays is ca. 10 3 time longer. Even though, when molecular flexibility is important, the time required for a proper sampling of the entire conformational space is hardly reachable. Other techniques have been developed to overcome this sampling issue.

First, coarse-grain force fields can use longer time step by lowering resolution; for instance they rationalize membrane crossing of antimicrobial peptides [177][178][179][180][181][182] or larger drugs such as paclitaxel [183]. Second biased MD simulations can be used to enhance and fasten sampling including metadynamics. Metadynamics was used in the case of ibuprofen and showed a conformation transition of ibuprofen from trans to cis in membrane center [184].

Although metadynamics is an elegant solution to comprehensively explore a given conformational space, the inherent too many parameters are a severe limitation that must be well chosen before the simulation.

Second, the quality of force fields can sometimes be a limiting factor to reach accuracy. We have showed that the choice of force field can significantly influence agreement between theoretical and experimental partition coefficients of small drugs (see ref [START_REF] Paloncýová | Benchmarking of force fields for molecule-membrane interactions[END_REF] and section IV.1). Additionally, while lipid force fields are generally available, drug parameterization has to be parameterized properly and systematically. A particular attention has to be paid to the description of partial charges in drugs, which have been shown of crucial importance to describe interactions with the different regions of membranes [START_REF] Lemkul | Practical Considerations for Building GROMOS-Compatible Small-Molecule Topologies[END_REF]126].

Third concerns another inherent limit to classical MD simulation, namely absence of any chemical reactivity as electron motion is not treated explicitly. MD simulations are thus appropriate for drugs that do not chemically react with membranes. However, drugs often bear ionizable moieties and are subject to chemical variations according to pH. The influence of protonation state is delicate to evaluate experimentally, in particular because pH in the surrounding of the polar head group of membranes is a complex issue. Classical MD simulations thus have the inconvenient of not allowing a drug to freely change its protonation state; only evaluation of the different charge states separately is possible. Indeed, using different simulations with fixed charge forms, MD can provide an atomic rationalization of the drug protonation/deprotonation events required for membrane crossing. This was for instance performed for ibuprofene [155] or vitamin C [167]. It should be noted that attempts of constant-pH simulations exist, but are currently marginal [185]. The ultimate methodology may reside in quantum dynamics, where electrons are explicitly defined while the system dynamically evolves over time. These methods may provide accurate results and dynamically describe electronic processes. However, they are incredibly expensive in terms of computational resources and they are currently limited to few atoms and short time scale [186]. Nevertheless, they might open considerable possibilities for pharmacology in a "distant" future.

Forth, that is last drawback but not least, most of theoretical studies performed so far have considered bilayer models made of a single lipid type, often pure phosphatidylcholine. However, lipid bilayer membranes are mixtures of lipids e.g. cholesterol is an essential component of biological membranes; and MD simulations evidenced its dramatic influence on drug partitioning [START_REF] Wennberg | Large Influence of Cholesterol on Solute Partitioning into Lipid Membranes[END_REF]187]. Adding cholesterol increases lipid chain ordering and may induce phase transition from L d to L o phase under certain conditions of concentration and temperature [START_REF] Petruzielo | Phase behavior and domain size in sphingomyelin-containing lipid bilayers[END_REF]. Sphingolipids are also important components of biological membranes, as they can favor L d /L o phase coexistence and thus formation of domains [START_REF] Petruzielo | Phase behavior and domain size in sphingomyelin-containing lipid bilayers[END_REF]. We showed in section IV.5 that temperature and membrane composition can dramatically influence drug position and orientation, its partition coefficient and membrane permeability. Therefore, it is crucial to use in silico models as close as possible to biological environment in its extreme complexity and diversity.

Chapter IV. Publications

IV.1. Benchmarking of Force Fields for Molecule-Membrane Interactions

Foreword

The exploration of drug-membrane interactions relies on a robust methodology. In particular, the affinity of drugs for membranes (i.e., their partition coefficient) is a critical descriptor of their bioavailability. Here, the capacity of different force fields to predict partition coefficients was benchmarked for eleven molecules against experimental data.

Paloncýová M, Fabre G, et al. J Chem Theory Comput. 2014;10: 4143-4151 [75].

Authors

Markéta Paloncýová † , Gabin Fabre †, ‡ , Russell H. DeVane ║ , Patrick Trouillas †, §, ⊥ , Karel Berka* † ,and Michal Otyepka* † Studies of drug-membrane interactions witness an ever-growing interest, as penetration, accumulation and positioning of drugs play a crucial role in drug delivery and metabolism in human body. Molecular dynamics simulations complement nicely experimental measurements and provide us with new insight into drug-membrane interactions, however, the quality of the theoretical data dramatically depends on the quality of the force field used.

We calculated the free energy profiles of eleven molecules through a model dimyristoylphosphatidylcholine (DMPC) membrane bilayer using five force fields, namely Berger, Slipids, CHARMM36, GAFFlipids and GROMOS 43A1-S3. For the sake of comparison, we also employed the semi-continuous tool COSMOmic. High correlation was observed between theoretical and experimental partition coefficients (log K). Partition coefficients calculated by all-atomic force fields (Slipids, CHARMM36, and GAFFlipids) and COSMOmic differed by less than 0.75 log units from the experiment and Slipids emerged as the best performing force field. This work provides the following recommendations i) for a global, systematic and high throughput thermodynamic evaluations (e.g. log K) of drugs COSMOmic is a tool of choice due to low computational costs; ii) for studies of the hydrophilic molecules CHARMM36 should be considered; and iii) for studies of more complex systems, taking into account all pros and cons, Slipids is the force field of choice.

IV.1.1. Introduction

In nature, biomembranes make selectively permeable walls separating inner and outer cell environments, or inner organelles and cytosol [188]. They play a key role in the control of active transport and passive permeation of endogenous or exogenous compounds [189][190][191]. Hence, the molecular interaction of xenobiotics (e.g. drugs and pollutants) with biomembranes is of major importance for understanding their flux through tissue and targeting in the human body [192][193][194]. Biomembranes are complex supramolecular systems, which mostly consist of lipids arranged as bilayers. They also contain proteins attached or embedded in the membrane bilayer [195]. The xenobiotics may interact with all these constituents during their passage through the membrane. Interactions of xenobiotics with the membrane-anchored cytochrome P450 represents a typical example of the complexity of membrane trafficking [196][197][198].

Basic features of the interaction of xenobiotics with biomembranes are known from experimental observations [199]. However the understanding is fragmented and the molecular picture is often missing. Molecular dynamics (MD) simulations have appeared as an alternative way to gain insight into structural features [200] and thermodynamics of interaction of guest molecules with biomembranes [201][202][203][204][205][206][207][208][209]. MD follows motions of all atoms of molecular system and generates a wealth of information having extremely fine resolutions both in time (subpicosecond) and space (atomic). This provides MD a major advantage with respect to all other techniques to tackle the interaction of xenobiotics with biomembranes, which nicely complements observations from the experimental techniques.

On the other hand, the quality of MD simulations is heavily limited by the underlying empirical potential, also termed force field (FF), and affordable sampling, i.e. duration of MD simulation [203,210,211]. In other words, inaccurate FF parameters may lead to biased structural or thermodynamic membrane parameters, hence, developed FFs are tested to determine the level of agreement with experimental observations. ). However, the accurate description of not only membrane structural parameters, but also molecular interactions between guest molecules and biomembranes requires highly advanced FFs. For even more complicated goals like membrane protein studies they should also achieve a properly balanced description of structural and dynamical features of proteins. To this end, advanced FFs compatible with advanced protein FFs would be a promising tool to describe the behavior of guest molecules within realistic complex biomembranes.

To simulate thermodynamics of the interaction between a guest molecule and membrane with MD is computationally demanding as they require robust sampling and in turn accumulation of long simulation times [203,207] This study aims at a critical analysis of molecule-membrane interaction, as evaluated by free energy profiles, which were derived from z-constraint MD simulations. In the test set, eleven organic compounds were included, having a broad range of affinities for dimyristoylphosphatidylcholine (DMPC) bilayers and also bearing common organic functional groups. Five advanced FFs dedicated to biomembrane simulations have been evaluated, including Berger, Slipids, CHARMM36, GROMOS43A1-S3 and GAFFlipids; for the sake of comparison COSMOmic has been also employed. Based on free energies, the partition coefficients were calculated for each molecule and each FF, and were compared to the available experimental data in order to investigate the performance of individual FFs for drugmembrane interactions.

IV.1.2. Methods

IV.1.2.1 Small molecule parameterization

A set of eleven molecules was selected for which experimental partition coefficients to DMPC membrane were available (Table 3) [227]

[242]

IV.1.2.2 MD simulation parameters

Fully hydrated membrane patches -bilayers, were prepared with 36 DMPC lipids in each monolayer surrounded by 0.15 M NaCl solution to mimic the physiological conditions (Figure 16). The bilayers were then equilibrated and the simulation setup was tested against the experimental structural membrane properties [243]. The simulation setup was then used for the z-constraint simulation (see all specific simulation parameters for all FFs in Table 4). The bilayer normal was oriented parallel to the z-axis and the origin of the axis was set in the middle of the bilayer. All MD simulations were performed by the GROMACS 4.

IV.1.2.3 Z-constraint simulation

Two drug molecules were initially placed in the simulation box: one in the middle of the membrane and another on the top of the simulation box i.e. into the water phase. The system was left for 500 ps to equilibrate and then both molecules were pulled in the same direction along the z-axis with a pulling rate of 0.05 nm•ns -1 and a harmonic force constant of 500 kJ•mol -1 •nm -2 . The initial structures for z-constraint simulations were separated from this pulling simulation. In each simulation box two drug molecules were placed, one in each monolayer. The windows for z-constraint simulations were chosen with separating distance of 0.3 nm, whenever possible.

Z-constraint simulations constrain a distance between different groups and monitors the required force applied on the molecule to keep this distance. The averaged force is then used to calculate the free energy profile also called potential of mean force (Eq. 1):

∆𝐺(𝑧) = -∫ 〈𝐹 ⃗ (𝑧)〉 𝑡 𝑑𝑧 𝑧 ′ 𝑜𝑢𝑡𝑠𝑖𝑑𝑒 , (Eq. 1)
where 〈𝐹 ⃗ (𝑧)〉 𝑡 is the force applied on the molecule in order to keep it at a given depth z. We constrained the two molecules in a box and monitored the applied force separately. Over the last years, we have systematically optimized the simulation protocol for free energy profile calculation in order to minimize the computer time cost [204]. Several authors have identified that the selection of an initial structure can slow the convergence of free energy profiles, especially in area of head groups [203,207,253]. The z-constraint simulation converges quicker compared to umbrella simulation, even when the initial structure is unequilibrated [203]. As it was also successfully used earlier [253], the amount of simulation windows was halved by adding two solute molecules in one simulation box. Free energy profile by zconstraint simulation allows a window of 0.3 nm. This significantly reduces the computingtime cost. It should be noted that cut-off lengths and water models dramatically influence computational time (Table 5). The z-constraint simulations were run for 30 ns per simulation window and the convergence of free energy profiles was monitored. The initial 15 ns of constraint simulation were left for equilibration and the free energy profiles were calculated from the last 15 ns. In the case of too slow convergence, the window lengths were extended to 50 ns (see the Supporting Information Table S1).

IV.1.2.4 COSMOmic free energy profile calculation

To increase the precision of COSMOmic calculations, 30 DMPC bilayer structures obtained from S-lipids simulation were used; this approach was successfully applied in earlier works [204,228]. The geometries and σ-profiles of DMPC, water and guest molecules were obtained by DFT/COSMO calculations at the BP/TZVP level of theory [254,255]. A single conformation as a result of geometry optimization was used. Free energy profiles were calculated at 310 K. Using the COSMOmic software [224] from the COSMOtherm 13 package, the bilayers were separated into 50 layers [256]. A total of 162 orientations of the solute molecules were used for each membrane to produce individual free energy profiles.

The final free energy profile was averaged over the individual free energy profiles of all the DMPC bilayer structures.

IV.1.2.5 Log K calculation

The free energy profiles obtained with MD (all FFs) and COSMOmic were analyzed and the partition coefficients were calculated according to an implemented method of COSMOmic [224,228] that removes the need for setting a membrane border and which is independent on the system size (Eq. 2):

𝐾 = ∫ (𝑒 -𝛥𝐺(𝑧) 𝑅𝑇 - 𝜌 (𝑧) 𝑤𝑎𝑡𝑒𝑟 𝜌 (𝑛) 𝑤𝑎𝑡𝑒𝑟 )𝑑𝑧 × 𝐴𝑃𝐿 𝑀 𝑙𝑖𝑝𝑖𝑑𝑠 𝑚 𝑢 𝑛 0 , (Eq. 2)
where 𝛥𝐺(𝑧) stands for a free energy at depth 𝑧, 𝜌 (𝑧) 𝑤𝑎𝑡𝑒𝑟 stands for water density at depth 𝑧 and 𝜌 (𝑛) 𝑤𝑎𝑡𝑒𝑟 stands for density of bulk water. The multiplying factor converts the partition coefficient into units used in experimental works kg (lipid) /L (water) . APL is the area per lipid, 𝑀 𝑙𝑖𝑝𝑖𝑑𝑠 is the molecular weight of lipids and 𝑚 𝑢 is the atomic mass constant.

IV. ), and in terms of the parameters of the linear log K exp vs. log K calc fit (Eq. 3):

log 𝐾 𝑒𝑥𝑝 = 𝑎 • log 𝐾 𝑐𝑎𝑙𝑐 + 𝑏, (Eq. 3)

which was constructed by the least-square method. The significance of the slope differing from 1 and intercept differing from 0 were evaluated at the probability level of 0.975. We analyzed the outliers of log 𝐾 predictions based on a Williams plot [257] and identified acetone in GAFFlipids, 2-nitrotoluene in GROMOS 43A1-S3 and 2-nitrotoluene and hexachlorobenzene in Berger. Due to the limited number of molecules investigated, we included the outliers in our analysis. However for analysis in a given FF, we excluded the outliers. We also analyzed the predictability of proper ordering of molecules according to their lipophilicity based on Spearman's rank correlation coefficient. Further we analyzed the heights of free energy barriers -the water/lipids barrier ΔG wat , the membrane center penetration barrier ΔG pen and the free energy at various membrane depths -and compared them to the values from Slipids that provided log K calc in the best agreement with experimental data. Additional relevant structural characteristics of fluid membranes are the order parameters of lipid tails [258]. The average order parameters were monitored (i.e. both sn1 and sn2 chains were averaged, Figure 18) during both unbiased and z-constraint simulations. Slipids, Berger and GROMOS43A1-S3 FFs reproduced the order parameters as best (mean absolute differences equal to 0.012, 0.013, and 0.015, respectively). On the other hand, the order parameters calculated by GAFFlipids and CHARMM36 were slightly overestimated (mean absolute differences 0.031 and 0.035). These findings agree with a recent work by Piggot et al. [200], comparing structural parameters of DPPC and POPC; the calculated order parameters of lipid tails in the plateau region below the head groups were the lowest with Berger, followed by GROMOS 43A1-S3 and CHARMM36. It should be noted that in the original publication of GAFFlipids [221] the order parameters were also slightly overestimated. However, DMPC membranes were in fluid phase with all FFs, for the full simulation time.

IV

The structural features of the DMPC membrane are sensitive to the simulation setup, especially cut-offs and water models. So, we used the setup suggested by the developers of each FF and when necessary we optimized the setup to acquire structural parameters best agreeing with the experimental data (see Table 4). As expected, from the point of view of computational time, the united atom FFs (i.e. Berger and GROMOS 43A1-S3) were the most efficient (Table 3 and Table S2). There were also differences among the all-atom FFs, the most effective being GAFFlipids due to a very short cut-off (0.8 nm). Slipids take advantage of uncharged carbons and hydrogens in the middle of aliphatic tails, while CHARMM36 was the slowest among all tested FFs, because of the long cut-off used, and the CHARMM modified TIP3P water model. In order to use parameters compatible with AMBER ff99SB FF for proteins, we also carried out Slipids simulations with 1.0 nm cut-off and tested CHARMM36 simulations with a 1.2 nm cut-off. In this case, the DMPC bilayer structural parameters stayed reasonably close to the experimental values (data not shown).

Decreasing the cut-off is an attractive way to increase performance for future simulations on larger membrane systems.

IV.1.3.2 The calculated partition coefficients agreed with the experimental values

Membrane/water partition coefficients were calculated by eq. 2 and compared with the experimental values (Table 3 and Table S2). The relative ranking of the molecules according to their partition coefficients, which was evaluated by the Spearman's rank order correlation coefficient, was reproduced best by Slipids and CHARMM36 (Table S2). The differences in ranking appeared for the medium lipophilic molecules for both CHARMM36 and Slipids, while both FFs ranked all lipophilic molecules properly. CHARMM36 ranked adequately even the most hydrophilic molecules (log K < 0.5) while Slipids ranked well all molecules with log K higher than 1.7. The COSMOmic approach also ranked properly the lipophilic molecules and performed just a little worse than Slipids and CHARMM36. GAFFlipids, Berger, and GROMOS 43A1-S3 showed worse ranking performance over the whole lipophilicity scale (Table S2). It should be stressed that all FFs and COSMOmic reproduce the right ranking of affinities to DMPC membrane, with α = 0.05 statistical significance.

The absolute predicted values of the partition coefficients log K calc also agreed with the corresponding experimental values log K exp (Figure 19, Table S3). The mean absolute difference with respect to log K exp of log K calc obtained with Slipids was 0.42 log unit, which is comparable with the experimental uncertainty for determination of log K exp . With this FF, the linear fit between calculated and experimental partition coefficients (cf. Eq. 3) led to a slope of effectively 1 and a y-intercept of effectively 0 (0.97(0.09) and -0.12(0.26), respectively, see Figure 19). CHARMM36 and COSMOmic exhibited similar performance (MAD 0.65 and 0.62, respectively), but CHARMM36 showed systematic shifts towards hydrophilic results (b = 0.59(0.22)), whereas COSMOmic towards hydrophobic results (b = -0.70(0.31)). GAFFlipids (MAD 0.74) gave one outlier (acetone) and GROMOS 43A1-S3 (MAD 1.08) gave 2nitrotoluene as outlier. When omitting the outliers the mean absolute differences dropped to more reasonable values, namely 0.50 and 0.78 for GAFFlipids and GROMOS 43A1-S3 FFs, respectively. The reason for the existence of these respective outliers has not been rationalized. Berger FF is known to overestimate lipophilicity of guest molecules and showed the largest deviation from experimental values [204]. In summary, taking the mean absolute differences and the linear fit of log K into consideration, the best performing FF among those tested here appears to be Slipids. However, the other FFs appear predictive enough, with the significant exception of Berger FF. Taking the predictive power (see also ref.

[204]) into consideration and regarding low computer cost, COSMOmic can be recommended for high throughput screening of interaction of small molecules, e.g. drugs, cosmetics, antioxidants, pollutants, pesticides and warfare agents with lipid bilayers.

IV.1.3.3 Properties of the free energy profiles

From the previous section, Slipids was taken as a reference, and the performance of the other FFs was tested in terms of water/lipids barrier ΔG wat and penetration barrier ΔG pen with respect to the corresponding values obtained with Slipids. The water/lipid barriers ΔG wat (that strongly correlates with log K calc , r 2 = 0.96) predicted by CHARMM36, GAFFlipids and COSMOmic were similar to those obtained with Slipids (Table 5, Table S4, Table S5, and Figure 20). GAFFlipids exhibited the lowest mean difference (MD 0.02 kcal/mol) and both GAFFlipids and CHARMM36 yielded the best mean absolute difference (0.72 kcal/mol, or even better -0.68 kcal/mol -when excluding the acetone outlier from GAFFlipids data set). ΔG wat values calculated by GROMOS 43A1-S3 exhibited a mean absolute difference of 1.65 kcal/mol; when removing the 2-nitrotoluene outlier from the dataset, the mean absolute difference dropped to 1.12 kcal/mol. Berger as expected predicted higher values of ΔG wat with a 2.09 kcal/mol mean absolute difference due to its over attractive Lennard-Jones interactions as we suggested earlier [204]. Concerning the mean difference of the penetration barrier ΔG pen , the best agreement with Slipids was achieved with CHARMM36 having a -0.15 kcal/mol mean difference and a 0.89 kcal/mol mean absolute difference. COSMOmic predicted ΔG pen values lower than Slipids with a -0.73 kcal/mol mean difference and a 0.91 kcal/mol mean absolute difference. The mean absolute difference calculated from GAFFlipids data was 1.33 kcal/mol (and 0.65 if acetone was excluded). The mean absolute differences calculated from GROMOS 43A1-S3 and Berger data were 1.28 and 1.06 kcal/mol, respectively. Though ΔG pen range is lower the range of ΔG wat with Slipids (5.8 and 8.7 kcal/mol, respectively), the relative mean absolute difference (with respect to Slipids) of ΔG wat of CHARMM36, COSMOmic and GAFFlipids is less than or equal to the mean absolute difference of ΔG pen . Therefore, CHARMM36, COSMOmic and GAFFlipids agreed with Slipids better for ΔG wat than ΔG pen . However, it must be stressed that in the case of GAFFlipids, the ΔG pen description was affected by the presence of one outlier (Figure 20). On the other hand, the mean absolute difference of both free energy barriers of CHARMM36 and COSMOmic compared to Slipids was lower than 1.0 kcal/mol. This confirms the ability of Slipids, CHARMM36 and COSMOmic to provide comparable and rather accurate predictions of the free energy barriers. The free energy profiles were also compared at different membrane depths calculated by all methods vs. the free energy profile from Slipids (Figure 21). The reference free energy value (ΔG = 0 kcal/mol) was set to water and the largest increase in the differences occurred at the water/membrane interface (2.5-1.5 nm from the membrane center). For COSMOmic, the maximum mean absolute difference (1.3 kcal/mol) was reached at 2.0 nm, dropped back to 0.7 kcal/mol at 1.75 nm, and slowly increased again to 1.2 kcal/mol in the middle of the membrane. With CHARMM36, it increased gradually up to 1.2 kcal/mol at the membrane center and the bump at the interface is less pronounced. GAFFlipids exhibited a slightly similar behavior with a mean absolute difference below 1.0 kcal/mol except at the center of the membrane. Berger and GROMOS 43A1-S3 failed in the description of the free energy profiles with respect to Slipids. Berger produced an excessively lipophilic description (i.e. too deep, Fig. S1) with a mean absolute difference reaching 2.9 kcal/mol in the center of the membrane. Concerning the united atom FFs, GROMOS 43A1-S3 is a better choice than Berger and all-atomic FFs and COSMOmic performed better than any of the united atoms FFs.

IV.1.4. Conclusion

This work compared the performance of five (two united atom and three all atom) FFs and the implicit COSMOmic method to reproduce the experimentally observed partition coefficients of eleven molecules into the DMPC membrane. Slipids appeared to be the most precise method, followed by COSMOmic, CHARMM36, GAFFlipids, GROMOS 43A1-S3 and Berger. COSMOmic and the all-atomic FFs performed well and reproduced the log K with a mean absolute difference lower than 0.8 log units. Perhaps a more relevant result is that Slipids, CHARMM36 and COSMOmic performed well in the prediction of free energy barriers; GAFFlipids predicted ΔG wat very well. In terms of computational time, COSMOmic is by far the best choice at predicting log K for fluid membranes. To study hydrophilic molecules, CHARMM36 is the only FF able to predict a correct ranking of lipophilicity. However, in the GROMACS software due the specific TIP3P water model required, CHARMM36 is the slowest, which might be limiting for larger systems, such as proteins and lipids. Taking all pros and contras into account, we recommend Slipids as the versatile FF for simulations of complex molecular systems containing lipid bilayers.

IV.2. Lipocarbazole, an efficient lipid peroxidation inhibitor anchored in the membrane

Foreword

Antioxidants play a major role in the prevention of lipid peroxidation. By using MD simulations, the structural properties responsible for efficient antioxidant activities can be elucidated. In this study, we rationalized the strong antioxidant activity of lipocarbazole as well as the function of its lipid side chain. Fabre G, et al. Bioorg Med Chem. 2015;23: 4866-4870 [168].

Authors

IV.2.1. Introduction

Oxidative stress is defined as an imbalance between production and regulation of reactive oxygen species (ROS), mainly free radicals. 22) is more active than ascorbic acid at scavenging DPPH (1,1-diphenyl-2-picrylhydrazyl) in methanol [263]. The DPPH assay is extensively used to provide a solid starting point to evaluate the capacity of a compound at scavenging free radicals by hydrogen and electron transfers [266]. An effective LPO inhibitor must i) scavenge efficiently free radicals from both thermodynamic and kinetic points-of-view, which is indeed well-related to free radical scavenging and ii) incorporate into lipid bilayer membranes. The combination of both features allows the compound to efficiently inhibit the LPO chain reaction [158]. Molecular modeling is a unique and powerful tool allowing the evaluation of these two properties at the atomic scale. In the present study, the thermodynamics of free radical scavenging reactions was obtained using quantum chemistry calculations for compound 1, whereas molecular dynamics (MD) simulations was performed to describe its incorporation, position and orientation into a lipid bilayer model. In order to evaluate the role of the lipid side chain, MD simulations (5µs total) were performed for 1 and its lipid-side-chain-less counterpart derivative, hereafter referred as "carbazole" (2).

IV.2.2. Material and methods

IV.2.2.1 Preparation of Large Unilamellar Vesicles and lipid peroxidation inhibition

A solution (100µL) of L-α-phosphatidylcholine from soybean (95%) at 40 mg/ml (Soy-PC, Avanti ® Polar Lipids inc.) was prepared in chloroform and was further evaporated under vacuum in a round bottom flask to produce thin Soy-PC film. Multilamellar vesicles (MLV) were produced by vortexing the thin film after hydration with distilled water. The MLVs were extruded through a 0.1 m double layer polycarbonate membrane using a Lipex extruder (Northern Lipids) to produce Large Unilamellar Vesicles (LUV). The preparation resulted in a aqueous solution of 2.50.10 -4 M Soy-PC LUV. The particle size was ranging from 90 to 110 nm, as determined using an N4plus submicron particle size analyzer (Beckman-Coulter).

The tested compounds of various concentrations (6.2 µL in methanol) were added to 500 µL of LUVs mixture prior to the lipid peroxidation initiation. Oxidative stress was generated by peroxyl radicals (R-OO  ) produced during AAPH (2,2'-Azobis (2-methylpropionamidine) dihydrochloride, Aldrich) degradation for 90 min at 37°C (125 µL, 2.5 mM in water). The oxidative stress effect was determined following the formation of conjugated dienes at  = 233 nm UV-visible absorption (Shimadzu UV-2401PC).

[267]

IV.2.2.2 Bond dissociation enthalpies

Over the past decade, free radical scavenging by polyphenols have been extensively investigated using quantum chemistry calculations [159,162,268,269]. The O-H bond dissociation enthalpy (BDE) was shown to be the major descriptor to predict free radical scavenging; the lower the BDE, the higher the capacity of H-atom transfer (HAT) from the antioxidant to the free radical, and the higher the antioxidant activity of the corresponding OH group. It perfectly and systematically correlates with DPPH scavenging. It is a thermodynamic intrinsic parameter calculated for all potentially labile chemical groups (mainly OH groups here) as the following difference in enthalpy (at 298 K):

BDE(Antiox-H) = H 298K (Antiox • ) + H 298K (H • ) -H 298K (Antiox-H) (1)
, H 298K (Antiox-H) being the enthalpy of the antioxidant and H 298K (Antiox • ) being the enthalpy of the radical formed after H atom abstraction.

Flavonoid derivatives and their corresponding aryloxyl radicals were found to be accurately described by density functional theory (DFT) calculations [162]. The B3P86 functional has been shown to be particularly well-adapted to evaluate the thermodynamics of the reaction between polyphenols and free radicals [159,162,161]. The 6-31+G(d,p) basis set is used since it provides very similar results compared to the larger and more computationally demanding 6-311+G(2d,3pd) basis set [161]; in particular the use of triple-ζ basis sets and the second diffuse function did not significantly enhance BDE predictions (difference lower than 1 kcal.mol -1 ). Geometries, energies including the zero-point correction (V) and enthalpies (H) at 298 K were determined at the (U)B3P86/6-31+G(d,p) level. Ground-state geometries were confirmed by a vibrational frequency analysis that indicated the absence of imaginary frequency.

The solvent effect was taken into account using the integral-equation-formalism polarizable continuum model (IEF-PCM) as implemented in Gaussian 09 [270]. Continuum models consider the molecular system embedded in a shape-adapted cavity surrounded by a dielectric continuum characterized by its permittivity (for water ε = 78.4). Calculations in water reproduce a polar physiological environment, while calculations in the gas phase and in benzene give a good approximation of non-polar conditions such as lipophilic membranes.

The implicit solvent model weakly influences the quantitative evaluation of phenolic BDE values but may slightly alter qualitative description, i.e., modifying the relative contribution of the different H atom donor groups [159]. All calculations were carried out by the Gaussian 09 software [270].

IV.2.2.3 Force field and membrane model

All MD simulations were carried out using the GROMACS package version 4.5.4 [271]. Two compatible united-atom force fields were used, namely GROMOS 53a6 [START_REF] Oostenbrink | Validation of the 53A6 GROMOS force field[END_REF] and Berger's [START_REF] Berger | Molecular dynamics simulations of a fluid bilayer of dipalmitoylphosphatidylcholine at full hydration, constant pressure, and constant temperature[END_REF] for water/hetero-molecules and phospholipids, respectively. The model of bilayer membrane consisted of 128 molecules of 1,2-dioleoyl-sn-glycero-3-phosphatidylcholine (DOPC) [272] surrounded by approximately 5400 water molecules (SPC/E model).

Phosphatidylcholines represent the main type of phospholipids in human membranes [START_REF] Van Meer | Membrane lipids: where they are and how they behave[END_REF].

The Na + and Cl -ions were added to the system by replacing water molecules using the Genion program, according to a regular physiological concentration C (0.9% = 0.154 mol L -1 ). Since Genion calculates the number of ions to be added according to the volume of solvent first obtained as the box volume. Since membrane is empty of water molecules, the corrected concentration C corr was re-calculated as:

𝑪 𝒄𝒐𝒓𝒓 = 𝑪. 𝒛 𝒃𝒐𝒙 -𝒛 𝒎𝒆𝒎𝒃𝒓𝒂𝒏𝒆 𝒛 𝒃𝒐𝒙 (2)

IV.2.2.4 Solute parameters

For the solutes (carbazole and lipocarbazole derivatives), the topologies were obtained from the PRODRG2 webserver [START_REF] Schüttelkopf | PRODRG: a tool for high-throughput crystallography of protein-ligand complexes[END_REF].

The partial charges defined by PRODRG2 webserver were significantly lower than those issued from the GROMOS force field and were shown to describe poorly the partitioning between aqueous and cyclohexane phases [START_REF] Lemkul | Practical Considerations for Building GROMOS-Compatible Small-Molecule Topologies[END_REF]. The restrained fit of electrostatic potential (RESP) [START_REF] Bayly | A well-behaved electrostatic potential based method using charge restraints for deriving atomic charges: the RESP model[END_REF] partial charges were alternatively used. RESP-type charges were successfully used in lipid bilayer simulations of several compounds [START_REF] Paloncýová | Benchmarking of force fields for molecule-membrane interactions[END_REF]100,167,169,273]. The ESP charges were obtained from B3LYP/aug-cc-pVTZ [START_REF] Duan | A point-charge force field for molecular mechanics simulations of proteins based on condensed-phase quantum mechanical calculations[END_REF] calculations obtained on geometries optimized at the same level, with Gaussian 09 software [270]. RESP fit was carried out with the Antechamber package of AMBER 11 [274].

IV.2.2.5 Free MD simulations

Several free simulations were carried out for every studied molecule, with different starting points (far from, close to and inside the lipid bilayer membrane). Energy minimization using the steepest-descent algorithm was performed before production simulations i.e., hundreds ns long MD simulations. The Leap-frog Verlet integrator was used with a 2 fs time step. The cut-off for electrostatic and vdW interactions were set to 14 Å. Long range electrostatic interactions were calculated with the particle mesh Ewald (PME) algorithm [START_REF] Darden | Particle mesh Ewald: An N⋅log(N) method for Ewald sums in large systems[END_REF]. Temperature and pressure were set to be as close as possible to physiological conditions, the temperature being maintained at 310 K by velocity rescaling [275] every 0.1 ps and the pressure being regulated anisotropically (each axis independently) at 1 atm by Parrinello-Rahman barostat [START_REF] Parrinello | Polymorphic transitions in single crystals: A new molecular dynamics method[END_REF] with a time constant of 5 ps and a 4.5 × 10 -5 bar -1 compressibility. All bonds were constrained by the LINCS algorithm (linear constraint solver) [276]. Periodic boundary conditions were used along the three axes.

IV.2.2.6 Free energy profiles

While free simulations provide an insight into the position and orientation of molecules, free energy profiles allow accessing additional information on i) the global energy minimum along the z-coordinate, defined orthogonally to the P-atom surface of membrane, ii) free energy barriers of membrane penetration and crossing, iii) free energy differences between inside and outside of the lipid bilayer that correspond to partitioning between lipid and water phases.

The free energy profiles were obtained with the z-constrained method [154,277-280], using the same parameters than for free simulations. This method defines a series of windows along the z-axis, in which a part of the molecule of interest is constrained at one specific z' coordinate. Several windows are defined along 40 Å, from the center of the water phase to the center of the membrane. The step between adjacent windows is related to the difference in the constrained position, which was 0.1 ± 0.02 nm in the present work. In this case, the carbazole moiety of both compounds 1 and 2 was constrained. In each window, the equilibrated starting structure was taken from a 20 ns long MD simulation, and was used to produce the 40 ns long MD simulation. The forces acting on a constrained molecule at a given 𝑧 ′ depth were averaged over time and integrated along the z-axis to build free energy profiles according to the following equation:

𝜟𝑮(𝒛) = -∫ 〈𝑭 ⃗ ⃗⃗ (𝒛 ′ )〉 𝒕 𝒅𝒛 ′ 𝒛 𝒐𝒖𝒕𝒔𝒊𝒅𝒆 (3)

IV.2.3. Results and Discussion

IV.2.3.1 Antioxidant capacity

The free radical scavenging capacity of a phenolic antioxidant (Antiox-H) occurs by hydrogen atom transfer (HAT) from one of the active group of this compound to the free radical (R • ), according to reaction [START_REF] Murray | Harpers Illustrated Biochemistry[END_REF]:

Antiox-H + R • ⟶ Antiox • + R-H (4)
Antiox-H is active when it possesses at least one sufficiently labile group to make reaction (1) thermodynamically favorable. Both OH and NH groups may provide providing efficient HAT capacity. The HAT capacity of the chemical groups is well evaluated by BDE. This parameter has been extensively evaluated and it systematically correlates with the DPPH scavenging activity [159,162,269]. Both O-H and N-H BDEs were calculated for both compounds 1 and 2 (Table 6). The N-H BDEs of compound 1 (85.9 kcal mol -1 in water) is much higher than that of the corresponding OH group (76.7 kcal mol -1 ). The absence of the lipid chain slightly increases the N-H BDE (87.1 kcal mol -1 ), but it still remains significantly higher than the O-H BDE (77.0 kcal mol -1 ). This clearly indicates for these compounds the major role of the OH group to scavenge free radicals with respect to the NH group; HAT from NH being unlikely. It must also be stressed that the O-H BDEs of both compounds 1 and 2 are similar, which unambiguously evidenced the same free radical scavenging capacity in solution for both compounds.

The crucial role of the OH group in free radical scavenging has been clearly shown for numerous polyphenols [159,162,167,269]. The O-H BDE obtained for compounds 1 and 2 is very similar to that obtained for the most active OH group of quercetin, a powerful free radical scavenger. A 76.7 kcal mol -1 BDE confirms the efficient free radical scavenging capacity, as already measured for compound 1 [263]. Such a low BDE allows to efficiently scavenging almost all free radicals R • of biological interest (e.g., • OH and LOO • ) or those used in antioxidant assays (e.g., DPPH). For these three free radicals, reaction ( 1) is predicted to be thermodynamically favorable with compounds 1 and 2 (H around -36, -8 and -3 kcal mol -1 , respectively † [161]).

The LPO inhibiting capacity was different for both compounds. Compound 1 is much more active than its lipidless counterpart ( 2), exhibiting a 6.25 times lower IC 50 (Table 7).

Compound 1 appeared even more active than reference antioxidants including α-tocopherol (vitamin E) and quercetin (Table 7). This high activity was obviously attributed to the high HAT capacity from the OH group of this compound, but this did not appear the only parameter as 2 was much less active. The role of the lipid chain clearly appears crucial to figure out the efficacy of this antioxidant as LPO inhibitor.

IV.2.3.2 Positioning and orientating compounds 1 and 2 in lipid bilayers

In order to investigate the position and orientation of both compounds 1 and 2, unbiased free simulations were carried out. From all MD simulations, starting either in the water phase or in the middle of the lipid bilayer, each compound reached a converged location within 40 ns.

The equilibrium distances of center of mass (COM) of the carbazole moiety with respect to † The H values were obtained from the difference between 77 kcal mol the center-of-membrane were 1.10 ± 0.05 nm and 1.53 ± 0.13 nm for 1 and 2, respectively (Figure 23). In both cases, the antioxidant-active moiety (i.e., the OH group) was mainly orientated towards the surface of the membrane but was inserted deeper inside the membrane for 1 compared to 2 (1.41 ± 0.07 nm and 1.66 ± 0.17 nm from the center of membrane, respectively). The OH group of 2 is embedded in the region of ester groups of DOPC, in close contact with water molecules and with a limited contact with lipid chains. In such a position and orientation, the inhibition of the propagation stage of LPO is unlikely, except if oxidized lipids adapt a snorkel-like shape [173], which rarely occurs for lipid peroxyl radicals [282]. The only way for compound 2 to inhibit LPO is to scavenge free radical initiators coming from the outer part of the membrane i.e., before the initiation stage. In this respect, compound 1 exhibits a totally different behavior; due to its deeper location, both inhibition of the initiation and the propagation stages are allowed. This makes compound 1 behaving as vitamin E i.e. as a powerful LPO inhibitor.

IV.2.3.3 Free energy profiles

To further rationalize the efficient antioxidant activity of compound 1, the free energy profiles of both compounds were calculated. A similar profile shape was obtained for both compounds (Figure 24), allowing to confirm that they are more stabilized inside the membrane rather than in the water phase. No energetic barrier was observed for membrane incorporation (Figure 24), confirming an easy penetration as expected for the fast localization observed during the free MD simulations i.e. within 40 ns. The convergence of free energy profiles was assessed (Figure S2). As for other small molecules like coumarin [100], time windows of 10 ns were sufficient for convergence. However, for bigger molecules longer windows are required [273]. The final profile of 1 was almost reached after 12 ns, but slowly fluctuated until 40 ns.

Three main differences appeared between both compounds. First, the free energy profile confirmed that 1 is inserted deeper inside the bilayer in agreement with the free MD simulations. Second, the difference in free energy between the optimized location and the water phase (partitioning) was dramatically lower for 1 (-25.0 kcal mol -1 ) compared to 2 (-8.93 kcal mol -1 ). Here we can conclude that the partition coefficient is higher for 1 than for 2, indicating that the concentration of the former compound in the membrane is much higher than that of the latter. Third, the barrier to cross the center-of-membrane is lower for 1 (2.82 kcal mol -1 ) than for 2 (5.84 kcal mol -1 ). Consequently, 1 has also a better capacity to flip-flop from one layer to the other. ‡ As LPO may also occur in the second layer, compound 1 can inhibit the propagation step more efficiently.

IV.2.4. Conclusion

In this study, the capacity of lipocarbazole, namely compound 1, to inhibit LPO and the role of its lipid chain in this process has been fully rationalized. Even if the lipid side chain does not improve the thermodynamic ability to scavenge free radicals, it plays a key role in ‡ This flip-flop was not observed during free simulations because it would have required MD simulations at much longer time scales. positioning 1 in membrane: i) it pulls the active OH group deeper in the lipid bilayer, increasing the contact with lipid chains affected by LPO, ii) it enhances the free energy difference between polar and non-polar phases, thus increasing the concentration of 1 inside the membrane, and iii) it lowers the central free energy barrier in the membrane, allowing a better contact between lipid chains of the membrane and the active OH group. These results are in very good agreement with experimental LPO inhibition.

This study also opens many perspectives for future works. The respective behavior of compounds 1 and 2 serve as a prototype in understanding the role of lipid side chains in compounds incorporating biological membranes. As can be seen, lipocarbazoles are more active than vitamin E as LPO inhibitors with potential use as natural antioxidants in food preservatives or as food supplementation. Furthermore, due to the structural analogy of compound 1 with carazostatin [265] ( 3), an effective in vivo antioxidant, compound 1 would probably be similarly bioavailable.

IV.3. Interaction of polyphenols with lipid bilayers membranes

Foreword

The successful rationalization of lipocarbazole antioxidant activity led us to go further and explore the interaction with membranes of a series of well-known polyphenols. The structural features influencing this interaction were analyzed. This work will be submitted for publication as soon as possible. Although their role as food nutrients is still under debate, it clearly appears that regular intake of exogenous antioxidant from fruit and vegetables is important to limit the deleterious effects of oxidative stress [284]. Moreover, the search for new efficient, safe and stable antioxidants has deserved much effort in cosmetics, food industry (food preservation) and even medicine (conditioning organ for transplantation). Even though antioxidants have been extensively investigated over the last decades, a thorough understanding of antioxidant mechanisms of action is still needed at an atomic level. Such understanding is required for an optimal and safe usage in all possible applications. Many of these applications are lipid media-based, requiring antioxidant actions on LPO. An efficient LPO inhibitor requires as least to be an efficient free radical scavenger. However, interaction and penetration within lipid assemblies, mainly lipid bilayers membranes, is of crucial importance. Over the last few years, molecular dynamics (MD) has appeared a promising tool complementary to experiments to tackle membrane penetration. MD simulations provide nothing less than an atomic picture of lipidantioxidant interaction, depth of penetration, orientation, energetic barrier of penetration and passive diffusion. The accuracy has been dramatically improved and screening of mediumsize database of potential antioxidants of compounds should soon be possible. Here, we provide an atomistic description of the capacity of twelve representative polyphenols to scavenge free radicals and to incorporate into lipid bilayer membranes. Structure property relationship is then proposed, highlighting a few key parameters affecting position, orientation or interaction free energy with lipid bilayer membranes.

IV.3.2. Materials and methods

IV.3.2.1 Molecular dynamic simulations

The membrane model is of a bilayer of 128 DOPC (1,2-dioleoyl-sn-glycero-3phosphocholine) molecules, solvated by water molecules in an 8 nm thick box. Na + and Cl - ions were included in the box at a 0.9% w/v physiological concentration. Lipids, polyphenols and ions were described by the Gromos43A1-S3 force field [285] whereas water was described by the SPC/E model [START_REF] Berendsen | The missing term in effective pair potentials[END_REF]. Polyphenol geometries were first optimized with density functional theory (DFT) at the B3PLYP/6-31+G(d,p) level. The corresponding topologies were then generated by the PRODRG webserver [START_REF] Schüttelkopf | PRODRG: a tool for high-throughput crystallography of protein-ligand complexes[END_REF]. Partial charges were recalculated with RESP (restraint fit electrostatic potential) from optimization made at the B3LYP/cc-pVDZ level in implicit solvent (ε = 4) [START_REF] Duan | A point-charge force field for molecular mechanics simulations of proteins based on condensed-phase quantum mechanical calculations[END_REF], using multiple conformations and orientations within the R.E.D (RESP and ESP charge Derive) software [START_REF] Dupradeau | tools: advances in RESP and ESP charge derivation and force field library building[END_REF]. All MD simulations and analyses were conducted using the GROMACS package version 4.5.4 [271]. MD calculations were integrated using a 2 fs time step and the leap-frog Verlet scheme. Electrostatic and van der Waals short-range interaction cutoffs were set to 1.4 nm. PME (particle mesh Ewald) was used for long-range electrostatic interactions. Temperature was kept constant at 310 K with the Nosé-Hoover thermostat (τ T = 0.5 ps). The Parrinello-Rahman barostat was used to maintain the pressure anisotropically at 1 bar (τ p = 5.0 ps, compressibility = 4.5.10 -5 bar -1 ). Periodic boundary conditions were used in every dimension. Bond constraints were handled by LINCS. Free MD simulations were 2-4 µs long so to ensure sufficient exploration of conformational space. The MD simulations were prolonged until equilibrium locations were reached and sampled during at least 1 µs. All values characteristic of penetration were averaged along the MD trajectory section for which equilibrated location was reached. The total sampling time for all calculations was equivalent to 36 µs.

IV.3.2.2 Free energy profiles

Free enthalpy profiles of membrane crossing was evaluated by COSMOmic [103,104]. This method based on an implicit description of the environment was shown particularly efficient at predicting partition coefficients as well as diffusion barriers [START_REF] Paloncýová | Benchmarking of force fields for molecule-membrane interactions[END_REF]. Briefly, COSMOmic calculates free energy profiles based on statistical thermodynamics and molecular polarity. This polarity is described by σ-profiles, which are histograms of partial charge as a function of molecular surface area. σ-profiles were obtained for polyphenols, DOPC and water molecules within the DFT formalism at the COSMO-BP/TZVP level with TURBOMOLE. Free enthalpy profiles were obtained at 310 K on one leaflet of a DOPC-like membrane, which was divided into 50 slices. Membrane-polyphenol interaction was properly sampled by using ca. 150 different orientations for each compound.

IV.3.3. Results and discussion

IV.3.3.1 General requirements for penetration depth of antioxidants

LPO is a chain radical reaction that is initiated by ROS generated either endogenously (mainly enzymatic processes) or exogenously (e.g., radiation, pollution, smoking) and reaching the lipid bilayer from the water phase. The initiator ROS (often hydroxyl radicals) diffuse in between lipid chains and abstract H-atom to form carbon-entered radicals. Fast O 2 addition yields peroxyl radicals ROO • , which can propagate from lipid chain to lipid chain by H-atom abstraction. Propagation is the limiting step and its inhibition is particularly efficient to block LPO. Therefore, the deeper antioxidant, the most efficient is role as LPO inhibitor.

Conversely, antioxidants located too close from membrane surface can only inhibit the initiation stage, therefore being less active. Here, 'deep' means that the antioxidant should ideally locate close to the production of ROO • i.e. close to lipid unsaturation.

Moreover, if antioxidants are sufficiently deep in the lipid bilayer they may act as vitamin E regenerators, as it is the case of ascorbic acid and certain polyphenols [167]. Lipid-soluble antioxidants such as vitamin E or quercetin, a prototypical polyphenol, efficiently inhibit LPO through synergetic effects [172,286]; we recently elucidated the molecular interaction favoring the cooperative effects between these two antioxidants [167,273].

The ultimate aim of such theoretical study is to predict the antioxidant activity against LPO. However, the direct correlation between the actual LPO inhibition activity by an antioxidant in the membrane and its location, orientation and partitioning is not straightforward as several other factors may influence this activity. For instance, the sheer radical scavenging activity, membrane composition, temperature, concentrations of pro-and antioxidants, or the presence of catalyzers may play a role. As all these parameters were not investigated here, such a correlation cannot be established instead trends are given. Moreover, to our knowledge no experimental study reported on the activity against LPO of all the polyphenols presented here within a single methodology. Several studies evaluating a few polyphenols are available, however as the experimental conditions differ from one study to another, the strict (quantitative) comparison of the reported activities cannot be made, again only trends can be provided.

Yet, only a few studies have reported on polyphenol location in membranes. Accurate experimental measurements often provide fragmented information and are rarely adapted for screening series of compounds, namely X-ray diffraction, neutron scattering, nuclear magnetic resonance (NMR), or electron paramagnetic resonance (EPR). Quenching fluorescence of probes characteristic of certain depth regions is probably the most widely used experimental technique, which provide general information of drug positioning but hardly distinguish the relative depth of penetration of the different moieties of a given compound. MD simulations have appeared adapted to support these experimental data, allowing nm and fs spatial and time resolutions, respectively. A few theoretical studies have successfully reported depth of penetration of polyphenols in agreement with experimental data [167][168][169]273,287].

IV.3.3.2 Polyphenols' membrane positioning

General location

All neutral polyphenols penetrate the lipid bilayer within a few hundreds of ns and lie at an equilibrated location just below the polar head groups, defining membrane surface. From free MD simulation the center of mass (COM) lies at 1.28 to 1.68 nm from membrane center (Table 8). The minimum obtained with the COSMOmic-based free enthalpy profiles agreed with the free MD simulations although compounds were systematically predicted less deep by a few Å (Table 8). Location of polyphenols below membrane is in rather good agreement with experimental studies [164,286,[288][289][290][291][292], which also confirms that polyphenols do not penetrate too deep in the membrane being often reported as non-efficient inhibitor of the LPO propagation stage [290,293].

Role of the number of OH groups

The number of substituted OH groups often distinguishes one polyphenol from another. It influences many of their chemical behaviors and biological activities. Galangin perfectly exemplifies the role of OH groups. Having no OH group on the B-ring, this compound penetrates deeper in the lipid bilayer with respect to quercetin, luteolin or catechin (Table 8). This comes from the importance of H-bonding interactions between phenolic OH-groups and the polar head group of the membrane, which drive location (Figure 26). This thus appears as a major descriptor to tackle polyphenols' depth of penetration.

Compound size

Molecular weight could have been a descriptor of membrane penetration. However, here we only studied rather low-molecular weight compounds and no significant penetration depth was observed between e.g. the flavonolignan silybin and the flavonoids (quercetin, luteolin or catechin), see Table 8. Although silybin was reported to locate rather in head groups than in lipid tails for DPPC membranes in gel phase, it is able to shorten the fluorescence lifetime of the deeply-located probe DPH in liquid-crystalline phase PC membranes [294].

Secoisolariciresinol was the deepest compound; such location is however most likely attributed to the presence of methoxy groups that decrease polarity, but also to its somewhat linear shape and flexibility facilitating insertion between lipid chains. The only point that must be properly addressed when working with bigger compounds is the sampling required for averaging all conformational re-arrangements. It must be long enough to allow penetration and stabilizing in a given energetic minimum. The case of high-molecular weight phenolic oligomers or polymers is beyond the scope of this article.

Role of the formal charge

The influence of formal positive or negative charge of polyphenol compounds was also evaluated. Polyphenols may indeed exist in different charge states including the +1 positive flavylium cation form of anthocyanin and anthocyanidins (e.g. cyanidin) or the negative charge of the deprotonated forms of polyphenols. Although the positive charge is often depicted on the O atom of C-ring, the charges obtained by quantum mechanics (QM) calculations and used in MD simulations is spread over the whole molecule, as expected to sufficiently stabilize flavylium cation at low pH. Surprisingly, the flavylium cation penetrates the membrane even deeper than the corresponding flavonoids (e.g. catechin or quercetin). This is rationalized by the charge delocalization that allows cyanidin bearing a dipole moment similar to that to be comparable to that of the neutral catechin (Table 8). A thorough analysis of MD simulations revealed that phospholipids were able to bend their head groups so to allow a direct contact between the flavylium cation and phosphate groups (Figure 27). This effect is driven by strong electrostatic interactions. Nevertheless, flavylium cation in membrane is mainly a figment of imagination; indeed flavylium cation only exist at very low pH, except when involved in copigmentation complexes [295], therefore anthocyanins can only approach membrane in their neutral or anionic forms [296,297]. Conversely, negative charge totally prevents penetration. For instance, the deprotonated form of rosmarinic acid (rosmarinate) is predominant at neutral pH. In contrast to the flavylium cation of cyanidin, the negative charge of rosmarinate is almost exclusively located on ester and carboxylate moieties, resulting in a much higher dipole moment than other polyphenols. As a result, rosmarinate locates in the polar head group region in interaction with the positive charge of the phosphatidylcholines. This location agrees with fluorescence experiments showing a change of generalized polarization of Prodan but not of Laurdan, the former being located in polar head groups while the latter is close to lipid tails [290].

Additionally, rosmarinate was also shown experimentally to spontaneously insert in lipid membranes [290].

Role of the sugar moiety

The sugar moiety pulls the compounds towards membrane surface, as the OH groups of the sugar moiety are prone to form H-bonding within the polar head group region with deep water molecules, carbonyl or phosphate moieties. This effect is clearly exemplified with the hypothetical case of the flavylium form of cyanidin-3-O-glucoside (COM is at 1.68 ± 0.34 nm with respect to 1.30 ± 0.20 nm in absence of sugar moiety). In the case of rutin, the sugar moiety is bulky, making the compound more embedded in the polar head group region but the flavonoid moiety is still at 1.4 ± 0.4 nm.

IV.3.3.3 Polyphenols' membrane orientation

Orientation of xenobiotics in membrane is a matter of utmost importance in understanding biological activities. It indeed provides location of active groups with respect to strategic locations in lipid bilayers including membrane protein active sites or specific moieties such as unsaturation of lipid chains. The latter sites are crucial to rationalize LPO inhibition as being the bilayer region where ROO • are produced and where the propagation stage occurs. Therefore more than a global location of the antioxidant polyphenols, one needs to know where do lie the most active OH groups as H-atom donors. It has indeed been shown that chemical structure affect polyphenol activity against LPO and that their orientation may contribute to their antioxidant activity [298]. Here, the average -angle between the longest axis of the polyphenol aglycone and the membrane surface normal (z-axis) was measured as being characteristic of orientation. It clearly appears that the compounds having many OH groups at different positions (catechin, cyanidin, luteolin, quercetin, resveratrol, secoisolariciresinol or taxifolin) display a preferred orientation that is somewhat perpendicular to membrane surface normal, in other words parallel to membrane surface (Figure 26). Such a location, in agreement with results reported so far [273,288,299], is driven by H-bonding interaction between the phenolic OH groups and the polar head group region. This means that most OH groups of the polyphenols are located at the same location than the COM (Table 8). In the case of galangin, having an inhomogeneous OH group distribution (i.e., no OH group on B-ring), the orientation is parallel to the membrane surface normal (-angle of 140 ± 20 degrees) somewhat parallel to the lipid chains. Such an orientation is clearly driven by the hydrophobic character of the unsubstituted B-ring, lying deeper in the bilayer. It should be noted that for glycosylated compounds the longest molecular axis is defined along the sugar and the aglycone moieties. However, the flavonoid moieties of rutin and cyanidin-3-Oglucoside are also oriented parallel to membrane surface. Thus, the slight twist of rutin reported on Table 8 stems from the location of the sugar moiety in polar head groups.

Interestingly, silybin exhibits a specific behavior; its size and conformational flexibility drives a more tilted orientation (i.e., -angle of 50 ± 24 degrees). The E-ring is oriented towards polar head groups thus locating the active C20-OH group [300] rather far from the lipid chains. This suggests that silybin is more prone to inhibit initiation than propagation of LPO.

IV.3.3.4 ∆G profile of polyhenols' membrane crossing

Computing free energy profile of drugs crossing lipid bilayers allows evaluation of their capacity to partition into membranes and to cross them by passive diffusion. From these profiles, free enthalpy of partition (ΔG part ) is calculated as the difference between free enthalpy in water and the minimum free enthalpy in membrane. Free enthalpy of penetration (ΔG pen ) is defined as the difference between minimum and maximum free enthalpies in membrane, the maximum being obtained at membrane center.

All polyphenols exhibited negative ΔG part values, ranging from -5.6 to -11.4 kcal mol -1 (Table 8), which indicate a major partitioning (ca. 99% according to COSMOmic calculations) inside membrane rather than in the water phase. Neutral polyphenols exhibited very similar ΔG part values (ΔG part values of -7.0 ± 0.9 kcal mol -1 ). Statistically significant higher membranes affinities were observed for charged compounds (ΔG part values of -10.5 ± 0.7 kcal mol -1 on average, p < 0.0001) due to strong electrostatic interactions between charged polyphenols and phospholipids head groups, namely i) choline moieties in the case of cationic polyphenols, and ii) sulfate or carbonyl groups in the case of anionic polyphenols.

ΔG pen values differ dramatically from one polyphenol to another, ranging from 4.3 to 36.3 kcal mol -1 , the two extreme values being obtained for galangin and rosmarinate, respectively. When considering a 1 µM polyphenol concentration, which could appear as the upper bioavailable limit expected in the organism after rich-in-polyphenol diets, such energetic barrier heights roughly correspond to one crossing event every 0.1 ms or billion years for galangin or rosmarinate, respectively. Thus, only galangin, secoisolariciresinol, silybin, quercetin, resveratrol, luteolin or catechin are candidates for passive diffusion across DOPC lipid bilayers. Conversely, rosmarinic acid, taxifolin, cyanidin, rutin, cyanidin-3-O-glucoside or rosmarinate cannot cross the membrane within a reasonable time. Therefore, it appears that the presence of OH or OCH 3 groups do not completely prevent membrane crossing, whereas both sugar, and charged moieties hinder passive diffusion, both moieties acting as solid anchor to the polar head group region. Interestingly, taxifolin bears a ΔG pen value of 17.0 kcal.mol -1 whereas it is twice as small for quercetin. Although the extension of the πconjugation has no real impact on the most stable location in lipid bilayer, it appears to significantly modify the capacity for passive diffusion. Indeed, extension of π-conjugation as in quercetin makes the compound more hydrophobic, as expected, thus reducing the energetic barrier at the membrane center, which is highly hydrophobic.

IV.3.4. Conclusion

The present work has highlighted that MD simulations is relevant to establish structure activity relationship of polyphenol penetration in membrane. The role of OH groups, sugar moiety, charge, size or π-conjugation has been thoroughly analyzed. The strength of intermolecular interaction (electrostatic and H-bonding) between polyphenols and polar head group region appear the driving force of location, partitioning and orientation. Although MD simulations have appeared mature to establish relevant structure activity relationship of membrane penetration, one must keep in mind that this study has been performed with a simple lipid bilayer model. Interaction of small molecules with membranes may be dramatically affected by biophysical properties such as lipid composition and temperature, both parameters affecting lipid phase. Therefore, attention must be paid when relating theoretical and experimental data and when aiming at using MD simulation as a predictive tool. Further systematic studies on membrane composition and phase should be carried out.

IV.4. Synergism of Antioxidant Action of Vitamins E, C and Quercetin Is Related to Formation of Molecular Associates in Biomembranes

Foreword

This joint experimental and theoretical study rationalizes the molecular mechanism of the cooperation between antioxidants in membranes by the existence of non-covalent associates. Fabre G, et al. Chem Commun. 2015;51: 7713-7716 [167]. vitamin C (ascorbic acid, vitC) and natural polyphenols are major antioxidants found in food. Depending on their bioavailability [301,[303][304][305][306][307], these antioxidants are known to be highly efficient ROS scavengers in different phases, namely vitE in membranes [172,173], vitC in plasma or cytosol [308] and flavonoids at the membrane/water interface [273,309]. When acting simultaneously, their overall antioxidant activity is synergistically enhanced [175,176,261,[310][311][312][313][314][315]. Free radical scavenging by vitE yields the corresponding αtocopheroxyl radical by hydrogen atom transfer (HAT), which in turn can be regenerated back to vitE by vitC [175,176,261,[310][311][312][313][314][315]. This synergistic effect has been shown enhanced by flavonoids [175,176,[310][311][312][313][314][315][316][317][318], which are efficient hydrogen atom donor antioxidants [162,319,320].

Authors

Here, we present a molecular description of the interaction between vitE, vitC and a representative flavonoid antioxidant, namely quercetin[321] (Figure 31), in lipid bilayer membranes. Using both in vitro and in silico models, the formation of mutual associations at the membrane/water interface is described for the first time. This description enables better rationalization of vitE regeneration by vitC, which is often enhanced in the presence of flavonoids. The penetration and positioning of vitC, vitE and quercetin in membrane was evaluated using a lipid bilayer model comprising DOPC molecules, as phosphatidylcholines are major components of biological membranes in plant and animal cells [START_REF] Van Meer | Membrane lipids: where they are and how they behave[END_REF]. Molecular dynamic (MD) simulations were used, which have been repeatedly shown to predict the positioning of small molecules in lipid bilayers in agreement with experimental data [START_REF] Paloncýová | Benchmarking of force fields for molecule-membrane interactions[END_REF]121,273,281,322]. The behavior of those three (non-interacting) antioxidants was evaluated by placing a single molecule in the lipid bilayer model during the MD simulations.

The simulations showed that vitE localizes below the membrane/water interface and can penetrate through the membrane center. The peak position of the C5-methyl group of vitE was found to be 1.5 ± 0.3 nm from the bilayer center (Figure 32A), which agrees with recent experimental data in DOPC bilayers (1.7 ± 0.4 nm) [173]. The OH group of vitE, which is responsible for free radical scavenging by HAT [323], was mainly located close to the lipid polar head groups, i.e., at the lipid/water interface suggesting inhibition of both the LPOinitiation (directly) and LPO-propagation (if the lipid chains adopt a transient snorkel-like shape [173,282]). Moreover, flip-flops may occur with an energetic barrier of 0.65 kcal.mol -1 as obtained by COSMOmic (Figure S1). This roughly corresponds to an occurrence of 1 flip-flop event every 1 µs at a 10 -6 µM concentration, in agreement with observations from our MD simulations. The flip-flop process is accompanied by the transient presence of an active OH group inside the lipid bilayer (Figure 32A) to scavenge the deeply buried peroxy radicals thus playing a role in direct inhibition of LPO-propagation. VitC is less buried in the lipid bilayer than vitE and resides in the outer layer close to the water phase (1.9 ± 0.3 nm) because of the lower lipophilicity of vitC with respect to vitE. Interestingly, the average location of quercetin and its aryloxyl radical formed under oxidative stress (1.7 ± 0.3 nm) was found to lie between that of vitC and vitE (Figure 32A). The flip-flop of quercetin is much less efficient than that of vitE, due to higher energetic barrier of 10.2 kcal.mol -1 (Figure S3), corresponding to a 1 s time-scale occurrence at 10 -6 µM.

Under physiological conditions (pH 7.4) and in an aqueous environment, vitC and quercetin are deprotonated (first pK a equal 4.2 and 5.7 in water for vitC and quercetin, respectively). As expected [324] the corresponding anions lies outside the membrane (Figure 32A) i.e. 2.5 ± 0.3 nm and 2.4 ± 0.2 nm for ascorbate and the phenolate form of quercetin (deprotonated at C-7), respectively. Acid-base equilibrium is likely to occur in the overlapping regions with the protonated forms (Figure 32A).

The lateral (x,y-plane) diffusion coefficients of vitC, quercetin and vitE were 17 ± 2, 17 ± 2 and 22 ± 5 x 10 -8 cm 2 .s -1 , respectively, as obtained from averaging MD trajectories (Table S7). These values are in agreement with the experimental self-diffusion coefficients of DOPC at 313 K (14.10 -8 cm 2 .s -1 ) [325], confirming that the MD simulation time was sufficient to allow correct sampling of all intermolecular motions. The diffusion coefficients along the z-axis were lower by one order-of-magnitude for the three antioxidants (Table S7), confirming rather extended residence time in the equilibrium locations.

According to the respective locations of the three studied antioxidants, quercetin may act i) by scavenging free radicals diffusing into the membrane like vitE, both quercetin and vitE being regenerated by vitC; and/or ii) as a vitE regenerator, thus enhancing the regeneration by acting in synergy alongside vitC. The active OH group of vitE overlapped that of the center of mass of vitC and quercetin in the head group region (Figure 32A) highlighting the proximity of the three antioxidants, so that the formation of mutual complexes seems likely, in the membrane layer close to the surface.

To confirm that such intermolecular association can be formed in the membrane, a series of 300 ns free MD simulations of the lipid bilayer containing several vitC, vitE and quercetin molecules was performed. This procedure allowed sufficient sampling of all possible noncovalent rearrangements and interactions (see Methodology section in Supplement Information). During the MD simulations, long-lasting (> 90% of the time) and close-contact pairs were observed, namely hetero-association complexes quercetin:vitE, quercetin:vitC and vitC:vitE, and self-association complexes quercetin:quercetin and vitE:vitE (Figure S4, Table S8). An extensive set of one hundred of 100-ns-long MD simulations quantified formation of self-and hetero-association, amounting to 27:45:28% for quercetin:quercetin, quercetin:vitE and vitE:vitE, respectively (Table S9). This does not significantly differ from a random distribution (25:50:25%); however, this should be interpreted with care, as the sampling is still quite limited despite all the effort. The driving force of such non-covalent association was thoroughly analyzed with quantum chemical calculations. Quercetin:quercetin, quercetin:vitE and vitE:vitE pairs were mainly held together by π-stacking interactions, whereas pairs involving vitC were stabilized only by intermolecular H-bonding. The stability of these non-covalent interactions was confirmed with density functional theory (DFT) augmented by an empirical dispersion term, namely B3P86-D2 recently re-parameterized to accurately evaluate stabilities of polyphenol non-covalent complexes [326]. Different intermolecular arrangements were predicted, namely head-tohead and head-to-tail, in which the importance of π-stacking (ring-to-ring distance of around 3.6 Å, as typical for π-stacking of aromatic rings [327]) and H-bonding was confirmed (see Figure 33 for the most stable geometries and Dataset S1 for all xyz geometries).

The in vacuo enthalpies of association ranged from -24.4 to -10.8 kcal.mol -1 (Table 1). The presence of aqueous environment lowered the absolute values of these association enthalpies by 10.0, 5.8, 8.0 and 14.2 kcal.mol -1 for quercetin:vitE, quercetin:vitC, vitC:vitE and vitE:vitE, respectively (Table 1). An entropy loss is expected accompanying formation of the non-covalent complexes, probably counterbalancing the strongly negative enthalpies of association. However, this entropy loss is most probably lower in the organized membrane phase with respect to vacuum [328] (see Methodology section in Supplement Information).

In any event, the quantum calculations confirmed that the associations are stabilized by a combination of intermolecular hydrogen bonding and π-stacking. According to this quantum evaluation, attractive forces definitely exist between the three antioxidants, favoring the formation of non-covalent (self-and hetero-) associations of antioxidants. ) calculated as the difference in energy (enthalpy) between the most stable complex and the isolated fragments, in the gas phase and in PCMtype benzene and water solvents. Negative values indicate that the association is thermodynamically favored compared to the pair of isolated fragments quercetin and vitE. An experimental confirmation was obtained from the fluorescence quenching of vitE embedded in DOPC liposomes in the presence of quercetin, added at increasing concentrations. VitE-containing liposomes were formed by addition of vitE to DOPC prior to liposome formation. These liposomes were then pelleted and re-suspended in buffer by a double ultra-centrifugation/re-suspension procedure so that non-inserted vitE molecules were discarded (see Materials and Methods section for details). Following this procedure, the measured vitE fluorescence (Figure 34A, condition: 0 µM of quercetin) was unambiguously assigned to vitE molecules embedded in the bilayer and not lying on the liposome surface.

System

With increasing quercetin concentration to the vitE-containing liposomes, a significant decrease in vitE fluorescence intensity was observed (Figure 34A). Quercetin did not exhibit any fluorescence when excited at 291 nm (excitation wavelength of vitE, Figure S5A) in both aqueous solutions and liposomes (Figure S5), therefore ruling out interference. The quercetin concentration-dependent fluorescence quenching thus suggests that i) quercetin molecules have the capacity to insert into the DOPC bilayer, and ii) quercetin:vitE complexes are formed.

The I 0 /I = f([quercetin]) Stern-Volmer plot is clearly non-linear and follows a quadratic function (Figure 34B). The quadratic evolution is unambiguously attributed to the presence of both static and dynamic quenching [329,330]. The linearity of [I 0 /I -1]/[quercetin] = f ([quercetin]) also confirms this concomitant quenching (Figure S6). The confirmed occurrence of static quenching supports the results of the MD simulations and indicates that quercetin penetrates the membrane and forms non-covalent complexes with vitE.

Our findings help to rationalize the results of previous experimental studies showing that addition of flavonoids synergistically increases the antioxidant activity of a vitE and vitC mixture in membranes [176,310,311]. The existence of non-covalent complexes between these antioxidants explains how pairs can dramatically improve LPO inhibition by increasing intermolecular contacts between antioxidants, enhancing recycling and subsequent synergic effects.

Indeed, from a thermodynamic point of view, the capacity of regeneration is confirmed by comparing the bond dissociation enthalpies (BDEs) of the most labile hydroxyl group of each antioxidant (Figure 31). The BDEs were calculated as 75.5, 78.7, and 78.7 kcal.mol -1 for vitE, quercetin (4'-OH group) and vitC, respectively § . These low values agree with previous experimental data that have been strongly supported theoretically [159,162,320], showing that these three compounds have a strong capacity to scavenge free radicals by HAT. The BDE values were similar for all three compounds, which indicates that HAT between the different antioxidants (native or oxidized) is thermodynamically allowable i.e., enabling the regeneration process. The only limitation to this process is thus the capacity of two antioxidants to come into contact. Here, we have shown that non-covalent interactions (mainly π-stacking and hydrogen bonding) drive this association process and put in close contact the active OH groups (see Figure 33 and Figure S4). This geometrically and statistically enables quercetin undergo HAT towards vitE to regenerate it. Because the BDEs of both compounds are rather close in energy, the reverse process (regeneration of quercetin by vitE) is likely as well, despite being less preferred. Due to π-stacking interactions between aromatic rings in a given complex, electron transfer between the two π-conjugated antioxidant partners is also likely to occur.

These effects would be even more enhanced in larger aggregates, e.g., in nanodomains (lipid rafts). VitE has already been experimentally shown to preferentially localize in lipid rafts [331]. Aggregation and formation of domains have also been evidenced at the membrane surface for catechin derivatives [332], but also inside the bilayer for quercetin [324] and curcumin [322]. § Interestingly, the ascorbate BDE was significantly lower (68.1 kcal.mol -1

), indicating that scavenging of the ascorbate radical by other antioxidants is unlikely. We hypothesize that the combination of the very low BDE of ascorbate and its location in the water phase could result in one-way extraction of free radicals from the membrane.

The average position of the non-covalent associations in the membrane was also evaluated. No significant location difference was detected between the antioxidants in the complexes and their respective individual partners, except for quercetin:vitE. Indeed, quercetin in these pairs exhibited a probability density with two peaks (Figure 32B). Although 50% of the quercetin molecules remained at a similar location to the individual molecules (1.7 ± 0.2 nm), 50% were pulled deeper into the membrane (1.3 ± 0.1 nm). This latter location allows the quercetin:vitE pair to span a larger part of membrane with respect to the non-interacting quercetin. This shift towards the center of the membrane may increase the capacity of quercetin to directly inhibit the propagation stage of LPO by scavenging lipid peroxy free radicals, which may also contribute to the synergetic effects.

We have presented a molecular insight into the synergism of vitE, vitC and polyphenols. Our results showed that vitE can reach vitC in the polar head group region of the membrane and form associations that favor its recycling. Quercetin can readily form non-covalent associations with vitE and vitC in membranes, therefore enabling regeneration of vitE and mediating vitE regeneration by vitC. Moreover, in these complexes, quercetin and vitE are better positioned in the membrane to inhibit LPO with respect to the separated fragments.

The occurrence of such associations should be systematically considered to support the research in new cocktails of collaborative antioxidants.

. Methodological comments

When placing several antioxidants in the bilayer (which is modeled in a periodically repeating box of 6×6×8 nm 3 volume) we definitely model quite high concentrations (one antioxidant molecule in the simulation box corresponds to a concentration of several mM). In our model, the lipid bilayer is too small and the time scale computationally available is limited to several µs. In a real situation, a few µM vitamin and polyphenol concentrations are expected in the organism. Here, we accelerate the diffusion and interacting processes by increasing the initial concentration, corresponding to antioxidant accumulation in the membrane, which has been indeed repeatedly suggested for vitamins and polyphenols.

Our theoretical data provide collective evidence that non-covalent aggregation is likely to occur in a biomembrane. All types of simulations (with DFT-D calculations, and free and constrained MD simulations) clearly show that with regard to the non-covalent association (e.g., quercetin:vitE), the inter-molecular interaction is sufficiently strong to give stable complexes inside lipid bilayers under physiological conditions. This is also supported by the (static) fluorescence quenching.

The occurrence of this association results from a complex phenomenon, in which enthalpy, entropy and solvation effects interplay. QM accurately evaluates the electronic energy of non-covalent association. However, the evaluation of the Gibbs energies of association at the QM level is a more delicate task. In such calculations, the entropy term consists of three major components arising from translational, rotational and vibrational modes. As these contributions are estimated under ideal gas, rigid rotor and harmonic approximations (which are valid only for molecules in gas phase), one might expect that the loss of translational entropy accompanying formation of the non-covalent complexes in the lipid bilayer is dramatically overestimated, as we indeed observed.

In any case, from the different data discussed here, we believe that the difference in Gibbs energy of association is negative but small, as already measured in hetero polyphenol association in water. Therefore the co-existence of non-covalent dimers and free antioxidants is much likely. The association will also be favored in the case of higher concentration, i.e., in the case of antioxidant accumulation in tissues of organs. Such accumulation (local concentration increase) has already been suggested for polyphenols. Here we show again, at a molecular level, that such accumulation may increase the global antioxidant status.

IV.4.3. Materials and Methods

IV.4.3.1 Molecular dynamics simulations

The membrane model consisted of a bilayer of 128 DOPC molecules solvated in ca. 5400 water molecules with 20 Na + and 20 Cl -ions (equivalent to 0.9% w/v physiological concentration). Antioxidant molecules and ions were described by the Gromos53a6 force field [START_REF] Oostenbrink | Validation of the 53A6 GROMOS force field[END_REF], water by the SPC/E model and DOPC molecules by the Berger force field [START_REF] Berger | Molecular dynamics simulations of a fluid bilayer of dipalmitoylphosphatidylcholine at full hydration, constant pressure, and constant temperature[END_REF], which is compatible with Gromos53a6. All MD simulations and analyses were conducted using the GROMACS package version 4.5.5 [271]. The MD calculations were integrated using a 2 fs time step and the leap-frog Verlet scheme. The electrostatic and Van der Waals short-range interaction cutoffs were set to 1.4 nm; PME (particle mesh Ewald) was used for long-range electrostatics. The temperature was kept constant at 310 K with the Nosé-Hoover thermostat (τ T = 0.5 ps), and the Parrinello-Rahman barostat was used to maintain the pressure anisotropically at 1 bar (τ p = 5.0 ps, compressibility = 4.5.10 -5 bar -1 ). Periodic boundary conditions were used in every dimension. Bond constraints were handled by LINCS.

Topologies for antioxidant compounds were obtained from the PRODRG webserver. Partial charges were recalculated using RESP and the model of Duan et al. [89]. The torsion angle χ between the B and C rings of quercetin was also carefully re-parameterized according to the recent accurate method developed for nucleic acids [333].

𝐸 dih,𝜒 solv = 𝐸 QM//QM,COSMO -𝐸 -𝜒 MM//MM,PB (1)
The following equation was used for describing the dihedral torsion potential V d :.

𝑉 d = 𝑘 𝜑 (1 + cos(𝑛𝜑 -𝜑 s )) (2)
The dihedral torsion constant k φ was determined to be 22.5 kJ.mol -1 (Figure S7).

Antioxidant molecules were initially placed in the water phase. The positions of individual compounds (ascorbic acid, ascorbate, quercetin, phenolate form of quercetin, α-tocopherol and α-tocopheroxyl radical) were obtained from simulations of single molecules interacting with the membrane.

The interaction between vitC, vitE and quercetin was assessed by high concentration simulations: 6 molecules of each antioxidant were placed in the water surrounding one bilayer, 3 of each compound close to either one or the other leaflet. The same starting structure was run 6 times for 300 ns to allow better statistical sampling. All molecules converged to their equilibrium positions before 150 ns. Therefore, only data for the last 150 ns were used for analysis. In the high concentration simulations, two molecules were considered to form pairs if the average distance between their centers of mass was less than 10 Å. This cutoff corresponded to the intermolecular distance for which pairs spend more than 90% of the time attached together (Table S8). It is important to note that the difference between this cutoff and distances typical for π-π interactions (3 to 4 Å) is because (i) the intermolecular distances were measured between the centers of mass and not as minimum distances between π -conjugated moieties, and (ii) thermal motion occurs. In order to further assess the relative populations of homo-and hetero-associates, one hundred 100 ns long simulations were conducted. Quercetin and vitE molecules were initially placed in both leaflets of the membrane at their equilibrium depth and at random x,y-positions. The concentration was varying from 2 to 4 molecules of each antioxidant per leaflet. The overall simulation time was 13.6 µs.

The free energy profiles of non-interacting antioxidants and their corresponding radicals were evaluated with COSMOmic [104] in a DOPC bilayer at 310 K (Figure S3).

IV.4.3.2 Quantum mechanics calculations

The potential energy surface of the various dimers was explored using our previously reported method [326]. Association energies (∆E) of the complexes were calculated as the difference in energy between the complex and the isolated fragments (Table 1). The association enthalpies (∆H) were calculated for all conformers of the complexes using a frequency analysis, a temperature of 298 K and a pressure of 1 atm. Negative values indicated that the association was thermodynamically favored compared to the pair of isolated compounds. The calculations were performed with B3P86-D2(s 6 =0.78)/def2-QZVP//B3P86-D2(s 6 =0.78)/def2-SVP. This re-parameterized DFT functional, which includes Grimme's dispersive term with the s 6 parameter adjusted to a value of 0. VitE and quercetin stock solutions were prepared in ethanol at 10 mM. Further dilutions were carried out directly in the lipid solution or in buffer to the desired concentrations. The samples containing vitE or quercetin were protected from light throughout their preparation.

Liposomes with embedded vitE were obtained by addition of vitE stock solution to a lipid chloroform solution with molar 4:1 lipid:vitE ratio (corresponding to a final vitE concentration of 50 µM in the liposome suspension). A dried DOPC/vitE-lipid film was formed from a chloroform solution by removing the organic solvent under a nitrogen stream, followed by a minimum of 2 h drying under vacuum. The dried lipid film was then hydrated in HBS buffer and extruded 19 times through 50 nm size calibrated polycarbonate membranes using a syringe-type extruder (Liposofast, Avestin Inc.) [335]. The liposomes were separated from free vitE molecules by ultra-centrifugation (200000 g for 2 h at 277 K; Beckman Coulter Optima L-100 XP ultracentrifuge) and re-suspended in buffer prior fluorescence measurement. This procedure was repeated twice in order to remove all non-inserted vitE molecules from the liposome pellet and final suspension. The hydrodynamic mean diameter of the liposomes after centrifugation was determined by quasi-elastic light scattering (Zetasizer, Nano-ZS, Malvern Instruments). Their size found to be homogenous with a diameter of 159 ± 2 nm.

The efficiency of the liposome separation from the non-inserted vitE was checked by adding vitE only after the liposome extrusion step. Pure DOPC liposomes were exposed to 50 µM vitE during 30 min. Again, to discard all non-inserted vitE molecules, the liposomes were submitted to the ultra-centrifugation/re-suspension steps prior to fluorescence measurement.

When vitE was incubated in the pre-formed pure DOPC liposomes (latter procedure), the measured vitE fluorescence was very weak (Figure 34A, control), whereas it was much higher when vitE was added prior to liposome formation (former procedure), see Figure 34A, condition: 0 µM quercetin. This shows that under the former procedure, the fluorescence was definitively assigned to vitE molecules embedded in the bilayer and not lying on the liposome surface, supporting the robustness of our cleaning step. Interestingly, this also shows that vitE does not insert easily into liposomes from a buffer solution under our experimental conditions.

To evaluate the effect of quercetin on the fluorescence signal of vitE, different concentrations of quercetin (from 25 µM to 100 µM) were added to the liposome suspension (lipid concentration of 200 µM) and incubated 30 min at 310 K. Emission fluorescence spectra of vitE or quercetin were recorded in buffer (Figure S5A&B and S4) and DOPC liposomes (Figure S5 and Figure 34A) on a Varian Cary Eclipse fluorescence spectrophotometer. Spectra were corrected by subtracting the baseline spectra of the corresponding blank vesicles. The excitation wavelengths of vitE and quercetin were 291 nm and 370 nm, respectively. However, the quercetin emission signal, either in buffer or incubated with liposomes, was too weak for further analysis (Figure S5B&C). Therefore, only data for vitE emission were employed.

IV.5. Position and orientation of carprofen derivatives in lipid-bilayer membranes: a joint theoretical and experimental study

Foreword

Carprofens are potential new drugs against Alzheimer's disease, as they can inhibit the synthesis of toxic peptides. This collaborative experimental and theoretical work (i) assesses carprofen position and orientation in lipid bilayers; and (ii) evidences membrane composition strong influence. This study will be submitted in the next few months.

Authors

Gabin Fabre Interestingly, non-steroidal anti-inflammatory drugs (NSAIDs) were shown to modulate γsecretase activity [339,340]. NSAIDs are able to inhibit production of Aβ peptides, but most of them require toxic concentrations to be active [340]. Carprofens, which belong to a NSAID family used in veterinary medicine, have been evaluated as an alternative and have appeared particularly efficient as selective and non-toxic inhibitors of γ-secretase [341]. They are capable of modulating production of amyloidogenic Aβ peptides (Aβ 40 and Aβ 42 ) to shorter or non-amyloidogenic peptides as Aβ 38 . Carprofen derivatives that were N-substituted by a lipophilic moiety have appeared 10 times more active when compared to the nonsubstituted ones [341]. Benzylcarprofen and sulfonylcarprofen (Figure 35) are chemically close to other active N-substituted carprofen derivatives [342] and might exhibit similar activities. The mechanism of γ-secretase activity modulation by carprofen derivatives has not been elucidated so far as neither intermolecular interactions nor tridimensional structures have been described yet. As speculated by Narlawar et al. [341] and suggested by their amphiphilic character, benzyl-and sulfonylcarprofen are most probably located in membrane, i.e., where APP cleavage into Aβ peptides proceeds. The first step towards understanding the mechanism of action of N-substituted carprofens is therefore to determine their interaction with lipid bilayer membranes. Their capacity to penetrate membrane as well their location and orientation in bilayer must be described with much accuracy. This would pave the way towards identification of targeted interaction sites. In this work, we provide a precise description of positions and orientations of carprofen molecules in membrane models as obtained by experimental and theoretical investigations.

IV.5.2. Results and discussion

Membrane composition and lipid phase are tightly interlinked, which influence thickness, area per lipid, diffusion constants and partition coefficients, among other physical-chemical parameters. Thus, the lipid bilayer model has to be carefully selected to be of biological relevance. To mimic neuronal membranes, a ternary mixture of 1-palmitoyl-2-oleoyl-snglycero-3-phosphocholine: sphingomyelin: cholesterol (POPC:bSM:Chol) 1:1:0.67 (mole/mole/mole) was chosen, as such a composition has appeared characteristic of neuronal cell membranes [START_REF] Norton | The lipid composition of isolated brain cells and axons[END_REF].

IV.5.2.1 Drug positioning

A thorough and accurate evaluation of depth of penetration, if any, requires using a set of collaborative techniques. Here we report the results obtained by neutron diffraction experiments and molecular dynamic (MD) simulations. Such joint experimental and theoretical approach may provide ca. 0.1 nm resolution.

It was shown that at 288 K the ternary lipid mixture  POPC:bSM:Chol 1:1:0.67  exhibits phase separation is structured as nanodomains (4-14 nm diameter) of both L d or L o phases [START_REF] Petruzielo | Phase behavior and domain size in sphingomyelin-containing lipid bilayers[END_REF], with a higher proportion of POPC in the L d phase and of bSM and Chol in the L o phase [343]. The situation is slightly different with pure synthetic SM (d18:1/18:0 N-stearoyl-Derythro-sphingosylphosphorylcholine) with respect to natural brain SM, as in this case the nanodomains are larger (75-100 nm diameter) [START_REF] Petruzielo | Phase behavior and domain size in sphingomyelin-containing lipid bilayers[END_REF]343]. In any event, the lipid bilayer models used for MD simulations bears a membrane surface of 5-7 nm in diameter and therefore nanodomains can hardly be observed. To account for the possible segregation in the two L d or L o domains, two models were built, namely the 1:1:0.67 POPC:bSM:Chol mixture at 288 K clearly mimic L o phase (i.e. with i) high order parameters [344,345]; ii) thicker membrane than L d phase [346]; iii) lower area per lipid than L d phase; and iv) low lateral diffusion coefficient (4-5 µm 2 .s -1 ) -whereas the pure DOPC bilayer mimics the L d phase domains. At 288 K, neutron diffraction measurements positioned the deuterated moieties of benzylcarprofen and sulfonylcarprofen at 2.64 nm and 1.99 nm from the membrane center, respectively. Such locations indicates that both drugs i) interact with the lipid bilayer and not only with water molecules, and ii) are lying just above the lipid headgroups, in contact with water molecules, sulfonylcarprofen being deeper in head groups. Free MD simulations performed at 288 K on the L 0 phase (POPC:bSM:Chol mixture) showed that the deuterated moieties of both compounds reached an averaged location at 2.5 ± 0.4 nm and 2.6 ± 0.4 nm from the membrane center, respectively. This means that when approaching from the water phase, both molecules interact with the low-headgroup density region in direct contact with the water phase (region 1 defined in the Marrink and Berendsen membrane model [118]) but do not penetrate membrane. Free enthalpy profiles of membrane crossing (Figure 37A) exhibited free enthalpy minima (-7.5 ± 1.4 and -8.6 ± 2.8 kcal mol -1 , respectively) at 2.6 and 2.4 nm from membrane center, respectively (Table 10). Both neutron diffraction experiments and MD simulations show preferential partition outside membrane in contact with polar head group. Such location is driven by the electrostatic interactions between the negative charge on the carboxylic acid moieties of both compounds and the positive charge on the choline moieties of lipids (Table 10 and Figure 36A &B). This force strongly constrains location of the whole compound, which fluctuates around the carboxylic acid anchor. In the case of benzylcarprofen, the compounds hardly penetrate deeper because the rest of the compound is hydrophilic. In the case of sulfonylcarprofen, the sulfonyl moiety spends some time below membrane surface, making the distribution profile broader (Table 10 and Figure 36). This capacity to twist around the carboxylic acid anchor is attributed to the polar character of the sulfur atom. This most likely rationalizes the slightly deeper location experimentally observed for the deuterated moiety of sulfonylcarprofen (vide infra). The free enthalpy profiles also suggest that membrane crossing by passive diffusion is unlikely due to the large free enthalpy barriers at the membrane center (46.7 ± 9.0 and 53.8 ± 6.4 kcal mol -1 , respectively).

Free MD simulations performed at 288K on the L d phase (DOPC) exhibited somewhat similar locations in the outer part of membrane in contact with polar head groups, i.e. at 2.2 ± 0.3 nm and 3.0 ± 0.4 nm from membrane center for both benzyl-and sulfonylcarprofen deuterated moieties, respectively (Table 10). However, a second minimum is observed suggesting the transient presence of the compounds deeper in the L d phase, in particular a clear minimum is observed at 1.2 ± 0.2 nm for sulfonylcarprofen. This confirms the role of the carboxylic acid moiety as an anchor in the polar head group and the capacity for twisting that is easier in L d phase, in particular for sulfonylcarprofen. Nonetheless, this inner location was not observed experimentally. Therefore it is reasonable to hypothesize that, at 288K, both molecules partition preferentially in the L o rather than in the L d phase.

At 310K, the ternary POPC:bSM:Chol mixture exhibits a higher miscibility but L o /L d phase separation is still present [343]. At this temperature, the sulfonylcarprofen CD 3 moiety revealed a peak at 1.58 nm with neutron diffraction, therefore suggesting deeper penetration than at 288 K, in region 3. For technical reasons, neutron diffraction experiment was not performed at 310 K for benzylcarprofen. The MD simulation study shows the existence of two different locations inner and outer (Table 10 and Figure 36 E&F). Here we confirm the role of the carboxylic acid as an anchor to the polar head groups. Thanks to the conformational flexibility allowed by higher temperature, both compounds can reach a deeper minimum. In the ternary mixture, although the deuterated moieties of both compounds were also present deeper in the bilayer, the COM of sulfonylcarprofen stayed anchored outside head groups (Table 10 and Figure 36F). Conversely, the COM of benzylcarprofen could penetrate deeper, at 1.7 ± 0.2 nm from membrane center. This difference was confirmed by Gibbs energy profiles (Figure 37B), the position of the minimum energy for benzylcarprofen being deeper than for sulfonylcarprofen.

In DOPC, due to the higher fluidity of the L d phase only the inner minimum is populated. The sulfonyl moiety even pulls the compound deeper in the bilayer, the deuterated moiety lying at 1.1 ± 0.2 nm. Gibbs energy differences between the water phase and the minimum in the membrane (ΔG wat ) is correlated with partition coefficients, namely the lower ΔG wat value, the higher affinity and partition coefficient. Therefore, by comparing ΔG wat values between L o and L d phases, the favored partitioning of carprofen derivatives in one or the other phase can be predicted. At 310 K, Gibbs energy profiles (Figure 37) suggest that benzylcarprofen partitions into L o phase. Conversely, sulfonylcarprofen has a higher affinity to L d phase, which is confirmed by the location in L d phase (i.e., below head groups) which agrees with experimental results, whereas it's the location obtained in L o phase does not agree.

Although the free enthalpy barrier of membrane crossing are lower at 310 K with respect to 288 K, the values are still rather high (ranging from 17.8 to 20.1 kcal mol -1 ), making passive diffusion slow or virtually impossible. It should be noted that these barriers were obtained for the deprotonated and negatively charged forms of carprofen derivatives. In the event of reprotonation in the head group region, neutral forms are most likely able to cross the lipid bilayer by passive diffusion, as this mechanism was repeatedly showed [167,347].

IV.5.2.2 Drug orientation

The thorough analysis of positioning revealed the crucial role of the carboxylic moiety as an anchor to the polar head group. This influence positioning but may also influence orientation.

To provide a clear picture of caprofen orientation in membrane is of crucial importance to tackle mechanism of interaction with γ-secretase. More than for location, orientation requires a set of converging data that we have assessed with deuterium solid-state nuclear magnetic resonance ( 2 H ssNMR) and MD simulations.

A clear difference in orientation was observed in MD simulations between molecules penetrating or not below the membrane head groups, regardless of the bilayer model (mixture or pure DOPC), temperature or carprofen derivative. To illustrate this difference, the vector 𝑉 1 ⃗⃗⃗ connecting the benzyl and the carbonyl moieties of each molecule was averaged along the simulations. Concerning the compounds at the outer location (i.e. center of mass above 2.0 nm from the center of the membrane) a unique orientation was not observed, except that 𝑉 1 ⃗⃗⃗ avoided parallel orientations with respect to membrane normal (Figure 38 A&B). A similar orientational distribution was also observed during the simulations at 288 K.

Concerning the caprofens at the inner location (i.e., distance of the center of mass lower than 2.0 nm from the center of the membrane) 𝑉 1 ⃗⃗⃗ was mainly parallel to membrane normal, that is aligned along with the fatty acyl chains of the lipids (Figure 38 C&D). In this case, the charged carbonyl moiety anchors the compound inside the polar head group region of the membrane whereas the lipophilic N-substituents are inserted deeper, between lipid chains.

Experiments were performed on the same ternary mixture than for neutron diffraction experiments (i.e., POPC:bSM:Chol 1:1:0.67 by mol) at 288 K i.e., the carprofen derivatives being mainly in region 1, on the membrane surface. The 2 H NMR spectra of deuterated compounds provided information on the order parameters of C-D bonds, which are related to the orientation of C-D bonds in regard to the magnetic field direction. Well-resolved peaks of sulfonylcarprofen suggested two defined orientations of the C-D 3 group with respect to the magnetic field direction, θ = ± 27.4°, or its complementary θ angle = ± 152.6°, which reflects the orientation (i.e the tilt angle) of the C-S bond relative to the magnetic field direction. The uncertainty in θ sign arises from the uncertainty in the sign of the quadrupole splitting.

From a specific narrow time scale (up to 400 ns) along the MD simulations, we have calculated the order parameters of the C-D bonds, S CD (Table 11). Large variations of the order parameters were observed along trajectories (data not shown). Again, regardless of membrane composition or temperature, compounds outside membrane (i.e. center of mass above 2.0 nm from the center of the membrane) exhibited low order with S CD parameters close to 0 on average. The low order arises from rotational averaging of the whole carprofen compounds as shown by the orientation of 𝑉 1 ⃗⃗⃗ and intramolecular rotational averaging.

However, experimental absolute values of order parameters were rather high for both compounds. This difference can be rationalized by the fact that in solid-state 2 H-NMR, membranes are not surrounded by bulk water but are stacked on top of each other. This suggests that there is no space for carprofen molecules to locate outside the membrane and that they are driven to their other minimum, i.e., below head groups. This hypothesis is supported by theoretical order parameters of compounds located below head groups that are in better agreement with experimental data, i.e. exhibiting higher values (Table 11).

IV.5.3. Conclusion

This study based on neutron diffraction, 2 H-ssNMR and molecular description from MD simulation rationalizes the position and orientation of two carprofen derivatives in different membrane phases and at two temperatures. Membrane composition and temperature influence phase properties of lipid bilayers and subsequently positioning and orientation of carprofen derivatives. At 288 K, both compounds were located outside the lipid bilayer head groups and probably partition into the L o phase. At 310 K, sulfonylcarprofen partitions deeper inside the bilayer in the L d phase, whereas benzylcarprofen rather stays outside head groups in the L o phase. Carprofen derivatives are thus present in two equilibrium positions, namely just above and below the surface of the bilayer. In the latter, interaction with both polar headgroups and lipid chains were observed. Although the exact orientation of carprofen derivatives lying on membrane surface was not precisely elucidated, some clear trends were shown i.e. either lying parallel-to-tilted with respect to the surface of the bilayer or parallel with respect to the lipid chains for outer-or inner-located compounds, respectively. Finally, the interaction of carprofen derivatives with neuron membrane is a complex matter depending on temperature and lipid phase and implying very different positions and orientations. Therefore, lipid bilayer phase should be carefully taken into account when evaluating the interactions of these derivatives with their protein targets.

IV.5.4. Methods

IV.5.4.1 Molecular dynamic simulations

Two membrane models were built. The first one is a ternary mixture of 60 POPC (1palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine), 40 sphingomyelin ((d18:1/18:0) N-stearoyl-D-erythro-sphingosylphosphorylcholine) and 10 cholesterol molecules. The second consisted of a bilayer of 128 DOPC (1,2-dioleoyl-sn-glycero-3-phosphocholine) molecules. Both lipid bilayer models were solvated using water molecules in an 8 nm thick box. Na + and Cl -ions were included in the box at a 0.9% w/v physiological concentration. Lipids, carprofen derivatives and ions were described by the Gromos43A1-S3 force field [285] and water by the SPC/E model. All MD simulations and analyses were conducted using the GROMACS package version 4.5.4 [271]. MD calculations were integrated using a 2 fs time step and the leap-frog Verlet scheme. Electrostatic and van der Waals short-range interaction cutoffs were set to 1.4 nm. PME (particle mesh Ewald) was used for long-range electrostatics. The temperature was kept constant at 288 K or 310 K with the Nosé-Hoover thermostat (τ T = 0.5 ps). The Parrinello-Rahman barostat was used to maintain the pressure anisotropically at 1 bar (τ p = 5.0 ps, compressibility = 4.5.10 -5 bar -1 ). Periodic boundary conditions were used in every dimension. Bond constraints were handled by LINCS. Every simulation was 300 ns long but 5 non-interacting carprofen molecules were placed in the bilayer models; this strategy allowed enhancing sampling equivalent to 1.5 s per simulation. The total sampling time for all calculations was equivalent to 15 µs.

Although carprofens usually exist as a racemic mixture, only the R-configured were modeled.

Even though S-carprofen exhibits higher COX-2 inhibition activity [348], the compound with the R-configuration exhibits higher bioavailability [349,350]. Carprofen derivative geometries were first optimized with B3PLYP/6-31+G(d,p), within the density functional theory (DTF) at the level. The corresponding topologies were then obtained from the PRODRG webserver. The partial charges were recalculated using RESP (restraint fit of electrostatic potential) and the model of Duan et al [START_REF] Duan | A point-charge force field for molecular mechanics simulations of proteins based on condensed-phase quantum mechanical calculations[END_REF].

The Luzzati thickness was calculated as 2 times the ratio between the volume per lipid (VPL) and the area per lipid (APL):

𝐷 𝐵 = 2 𝑉𝑃𝐿 𝐴𝑃𝐿 = 2 𝑉 𝐵 -𝑛 𝑊 . 𝑉 1𝑤 𝐴𝑃𝐿
VPL is obtained from the volume of the box (V B ), the number of water molecules (n w ) and the volume of one water molecule (V 1w ).

While free simulations provide an insight into the position and orientation of molecules, free energy profiles allow accessing additional information on i) the global energy minimum along the z-coordinate, defined orthogonally to the P-atom surface of membrane, ii) free energy barriers of membrane penetration and crossing, iii) free energy differences between inside and outside of the lipid bilayer that correspond to partitioning between lipid and water phases.

IV.5.4.2 Free energy profiles

The free energy profiles were obtained with the z-constrained method [154,277-280], using the same parameters than for free simulations. This method defines a series of windows along the z-axis, in which a part of the molecule of interest is constrained at one specific z' coordinate. Several windows are defined along 40 Å, from the center of the water phase to the center of the membrane. The step between adjacent windows is related to the difference in the constrained position, which was 0.2 ± 0.02 nm in the present work. In this case, the COM of carprofen derivatives was constrained. In each window, the equilibrated starting structure was taken from a 100 ns long pulling MD simulation at 0.05 nm ns -1 rate and 500 kJ mol -1 nm -2 harmonic restraint constant. Then, each window was run for 50 ns which was sufficient for all profiles to converge. The forces acting on a constrained molecule at a given 𝑧 ′ depth were averaged over time and integrated along the z-axis to build free energy profiles according to the following equation:

𝛥𝐺(𝑧) = -∫ 〈𝐹 ⃗ (𝑧 ′ )〉 𝑡 𝑑𝑧 ′ 𝑧 𝑜𝑢𝑡𝑠𝑖𝑑𝑒

IV.5.4.3 Calculation of order parameters

Order parameters S CD reflect order and orientation and can be straightforwardly obtained from measurement of the angle θ between the C-D vector and the bilayer normal:

𝑆 𝐶𝐷 = ⟨ 3𝑐𝑜𝑠 2 𝜃 -1 2 ⟩
where brackets mean time-averaged value of θ over the MD simulation.

S CD can also be obtained with a sign uncertainty from 2 H ssNMR quadrupolar splittings Δν Q :

|𝑆 𝐶𝐷 | = 2 3 1 𝑒 2 𝑞𝑄 ℎ 𝑆 𝑙𝑖𝑏 ∆𝜈 𝑄
where 𝑒 2 𝑞𝑄 ℎ ⁄ is the static quadrupolar splitting constant. To account for librational motions, the quadrupolar splitting constants were weighted by S lib = 0.88. These constants depend on the type of C-D bonds. The static quadrupolar splitting constant is 180 kHz for a CD 3 moiety or an aromatic C-D bond, and 167 kHz for methylene CD 2 moieties.

For benzylcarprofen, the order parameters of Cδ-D and Cε-D were considered equivalent due to rotational motions. In this case, S CD can be obtained from simulations using the following equation:

𝑆 𝐶𝐷 = ( 5 8 (1 -cos 2 𝜑 -cos 2 𝜓)) -( 1 8 cos 2 𝜑) -( 1 2 cos 2 𝜓)
where ψ is the angle between Cβ-Cγ and membrane normal, and φ the angle between the benzyl ring normal and membrane normal.

IV.6. A complete conformational analysis of plantazolicin

Foreword

This study exemplifies how MD simulations can describe more complex phenomena. Here, the conformation and the interaction with membrane of a large antibacterial compound were assessed. Pore formation was ruled out, opening the way to other assumptions for its mechanism of action. This work will also be submitted shortly.

IV.6.1. Introduction

Natural compounds are a perpetual source of potentially new drugs. Bacteria have been widely studied to provide wide variety chemical structures. Bacillus amyloliquefaciens FZB42 is a bacteria that produce many antimicrobial secondary metabolites: polyketides (bacillaene, difficidin, and macrolactin), lipopeptides (surfactin, fengycin, and bacillomycin D), and siderophores (bacillibactin) [351]. It was recently found that the FZB42 mutant RS6 is unable to produce the compounds mentioned above but also to produce plantazolicin (PZN), which is active on closely related gram-positive bacteria [352]. The chemical structure of PZN was recently fully elucidated (Figure 39) [353]. It belongs to the thiazole/oxazole-modified microcin (TOMM) family, which includes numerous compounds having various biological activities [354] e.g., thiostrepton (50S ribosome inhibitor), trunkamide (anti-cancer drug), microcin B17 (DNA gyrase inhibitor), goadsporin (secondary metabolism inducer), yersiniabactin (siderophore) or ritonavir (HIV-1 protease inhibitor). TOMM-type compounds are synthesized from peptides with high concentration in cysteine, serine and threonine, which undergo posttranslational cyclodehydration and dehydrogenation to form oxazole and thiazole rings.

The TOMM-type structural characteristics provide to PZN, two -conjugated, rigid and hydrophobic extended moieties separated by two central isoleucines, which connect both rigid moieties by flexible classical amino bonds (Figure 1). As opposed to the central part of the molecule, the two edges are polar, even charged at physiological pH. With regards to molecular size, the structural and conformational features of PZN may strongly impact its biological response. Due to the hydrophobic character of the main part of this polypeptide, interactions with lipid bilayer membrane have been hypothesized. However, the mechanism of action of PZN as active compound on closely related gram-positive bacteria has not been elucidated yet. A thorough analysis of it conformational feature appears the crucial initial stage to pave the way of a better understanding of its biological activities. For that purpose, quantum mechanics (QM) and molecular dynamics (MD) simulations, combined with NOESY-NMR data can provide an accurate atomistic description of its 3D structural feature [355], and subsequently shed light on its behavior in a biological environment e.g., in interaction with lipid bilayer membranes. PZN has too many degrees of freedom to allow a complete exploration of the conformational space at the quantum level. For such a molecular system, free MD simulations and simulated annealing calculations are much better adapted to evaluate all intramolecular rearrangements. The summary of the conformational study is reported on Figure 40. The simulations were carried out in water considering physiological conditions (pH 7.3), under which the C-and N-terminal amino acids of PZN are ionized i.e., with a total charge of +1 (Figure 39). Interestingly, when starting from a linear conformation (Figure 41A), PZN folded on itself during both free MD simulation and simulated annealing (SA) procedure in water (Figure 41B). The intramolecular folding appeared driven by hydrophobic effects, both conjugated moieties being parallel to each other along the MD simulations. In order to rationalize the importance of π-stacking, we computed the amount of interatomic distances accounting for π-stacking (see Figure 40 and Methodology section for more information); this folded monomer exhibited 718 and 646 contacts for free MD simulations and SA, respectively. Concomitantly, both terminal amino acids showed close contacts (Figure 41B). When starting from the folded conformer, no unfolding was observed in both free MD simulation and SA. MD simulations of the neutral PZN (Figure 39) were also performed, and folding was similarly observed, confirming that the driving force of the folding is mainly the stacking. To better rationalize this intramolecular interaction, the folded geometry was optimized at the quantum level. Using an appropriate methodology (i.e., properly taking noncovalent interactions into account, see Methods section), the folded conformer (Figure 41C) appeared 41.5 kcal.mol -1 more stable, in terms of electronic energy, than the unfolded conformer (Figure 41A). Such a huge intramolecular force was explained by i) π-stacking interaction between the two extended conjugated moieties, and in a minor extend by ii) electrostatic interactions between the polar terminal amino acids and iii) hydrogen bonding between the N-terminal arginine amines and the C-terminal carboxylate, forming a sixmembered pseudo-ring (Figure 41C). To be comparable with the conditions used for NMR measurements, MD simulations were also carried out in DMSO. In this solvent, the conformation of PZN was studied in both neutral and charged forms, starting either from the linear or the folded conformation. In all cases, no folding was observed, and the number of contacts involved in π-stacking was close to 0 (Figure 40). The molecule did not fold when starting from its linear conformation and did unfold when starting from a folded conformation. Therefore all intramolecular interactions appear strong enough to provide hydrophobic effects i.e. folding the molecule in water, but not enough to fold the structure in less polar environment as DMSO. This result is consistent with the absence of intramolecular NOESY-NMR contact between C-and Nterminal moieties in DMSO. However, NOESY-NMR contacts were observed between i) N,N-diMeArg 1 and Ile 8 and ii) 5-MeOxl 13 and Ile 7 . They imply head-to-center and tail-to-center contacts that can hardly occur within a monomer, as both π-conjugated moieties are planar and rigid. Therefore these contacts are likely to be intermolecular, namely occurring within dimers or oligomeric clusters.

IV.6.2.2 Conformation of PZN dimer

Different types of geometries for PZN dimers were guessed and their close conformational space was explored. The most obvious guess geometry was built from two linear PZN monomers aligned in a head-to-tail conformation favoring electrostatic interactions between terminal moieties (dimer A, see Figure 41D). The other two geometries were guessed from NOESY-NMR data, showing contacts between N-terminal moieties and central isoleucines. Dimer B was made of two linear PZN monomers aligned in a head-to-center way, and dimer C was made of folded PZN monomers with aligned in a head-to-center way as well (Figure 41E). Free MD simulations (100 ns) of these guessed representative geometries were achieved in an explicit DMSO environment. Both dimers A and B quickly dissociated to yield two independent linear monomers, while dimer C remained stable at 310 K all along the simulation. In the guess geometry of dimer C the number of interatomic distances involved in -stacking interactions (3150) was higher than in guess geometries of dimers A and B (1210 and 669, respectively). Moreover, the stabilizing non-bonding interactions during free MD simulations appear much more important in dimer C with respect to the folded monomer (2093 vs 646 contacts, respectively), thus explaining the stabilization of dimer C in DMSO while the folded monomer was not stable enough.

Therefore, dimer C appears a relevant candidate to agree with experimental NOESY-NMR data. Due to the high number of degrees of freedom (inducing many possible supramolecular re-arrangements), SA was also performed to better explore the potential energy surface of dimer C. The high temperatures in the SA procedure (up to 600 K) overcame the stabilization observed in free MD simulations and yielded two independent monomers. Therefore, all distances corresponding to NOESY contacts were then restrained between 4 and 5 Å with a 750 kJ.mol -1 .nm -2 force constant. More favorable conformations were observed increasing πstacking interactions (2131 contacts, Figure 42). This illustrates how non-bonding interactions (mainly -stacking) play a crucial role in the conformation of the PZN dimers.

These results perfectly agree with the experimental NOESY-NMR data obtained in DMSO and provide an elegant molecular picture of the intermolecular contacts. The existence of these dimers was also envisaged in water, as water can be a more realistic environment when dealing with the antimicrobial activity of PZN. Free MD simulations were performed from the same three starting dimer geometries, namely dimers A, B and C. In this case, the molecules were considered in their charged form. Contrary to what was observed when using DMSO, all dimers appeared stable along the simulations. Only dimer B was modified and adopted geometry similar to that of dimer A (linear head-to-tail). Again, to further investigate the conformational space around these dimers, SA simulations were performed. Interestingly, in water, all PZN dimers quickly converged to a confirmation similar to dimer C i.e., both monomers folded, in close contact to each other, and aligned head-to-center (the terminal aminoacids of one monomer -arginine or phenylalanine -being in contact with the central isoleucines of the other monomer). This further confirms that this conformation (Figure 41E) is the most stable for a dimer and that the SA procedure was able to explore extensively the conformational space.

IV.6.2.3 Interactions with membranes

As an efficient antibacterial compound, PZN could hypothetically interact with membranes. Indeed, one of the main targets of small cationic antimicrobial peptides is the bacterial membrane [356] which mainly consists of a phospholipid bilayer. According to the Marrink and Berendsen membrane model [118], four regions can be clearly defined in a bilayer namely 1) the low head group density in direct contact with the water phase, 2) the high head group density defining an intermediate region between the water and lipid phases, 3) the lipid chains, and 4) the center of the membrane with low lipid density. The former two are polar while the latter two (approximately 4 nm width for the whole bilayer depending on the lipid chain length) are apolar. Interestingly, PZN (5 nm long) fits in size and in polarity (polar edges and apolar central moiety) into the bilayer. That is, one can easily imagine that PZN is a transmembrane compound in its linear conformation. This hypothesis was fully confirmed by MD simulation i.e., when guessing a transmembrane compound, after 100 ns no significant change was observed confirming that the transmembrane arrangement was stable. The crossing of the charged terminal aminoacids through the apolar regions is unlikely, which prevents folding of PZN when being inserted into the bilayer. However, starting out of the membrane, the penetration has appeared unfeasible within the µs time scale. Even if no longer simulation times are accessible, a classical penetration process is unlikely, again due to high energetic cost required for the charged terminal aminoacids to cross the apolar regions. Therefore PZN can only be in a transmembrane position if it takes part of bilayer formation (or destruction) processes.

In order to further investigate the antimicrobial action mechanism of PZN, we also made the hypothesis of a multimeric membrane pore. Pore creation is a typical mode of action of antimicrobial peptides that form holes in the lipid bilayer, allowing free ion and molecule exchange between intra and extracellular compartments, leading to bacterial death. Here, we designed a hypothetical pore consisting of six transmembrane and linear PZN molecules (Figure 43A, see Methodologies of calculations for more details). After only 12 ns MD simulation, the pore completely collapsed and all water molecules present in the cavity defined between the six PZN molecules flew outside the membrane (Figure 43B). This collapsing is fully rationalized by the strong non-covalent intermolecular interactions previously described. Therefore, the capacity of PZN to form a classical pore in the bacterial membrane is highly unlikely, whatever the pore shape that can be guessed.

IV.6.3. Conclusion

The conformational space of PZN has been explored exhaustively; this paves the way to a better understanding of its antibacterial mechanisms of action. Our findings suggest that in water, PZN exhibits a strong capacity to fold on itself due to intramolecular but also intermolecular -stacking interactions. In DMSO (in which NMR-NOESY data were collected), PZN consists of independent monomers adopting a linear conformation. However favoring intermolecular interactions (e.g. increasing concentration), non-covalent dimers are formed thanks to intermolecular π-stacking. The most stable conformer of these dimers obtained from simulated annealing agrees with the observed NOESY-NMR data. Such ability of PZN for stacking by non-covalent interactions prevents the formation of classical pore in membrane. One can only imagine concentration and accumulation effects of PZN at the membrane surface and afterwards inside the membrane.

IV.6.4. Methods

IV.6.4.1 Force field and membrane model

All molecular dynamic simulations were carried out using the GROMACS package version 4.5.4 [357,358]. Two compatible united-atom force fields were used, namely GROMOS 53a6 [START_REF] Oostenbrink | Validation of the 53A6 GROMOS force field[END_REF] for DMSO and PZN, and the Berger force field [START_REF] Berger | Molecular dynamics simulations of a fluid bilayer of dipalmitoylphosphatidylcholine at full hydration, constant pressure, and constant temperature[END_REF] for phospholipids. The SPC/E model was used for water. Na + and Cl -ions were added to the water at a physiological concentration C (0.9% = 0.154 mol.L -1 ). MD was integrated using a 2 fs time step and the leap-frog Verlet scheme. Electrostatics and Van der Waals short-range interaction cutoff were set to 1.4 nm, Particle Mesh Ewald was used for long-range electrostatics. The temperature was kept constant at 310 K with Nosé-Hoover thermostat ( T = 0.5 ps), while the Parrinello-Rahman barostat maintained the pressure anisotropically at 1 bar ( p = 5.0 ps, compressibility = 4.5.10 -5 bar -1 ). Periodic boundary conditions were used in every dimension. Bond constraints were handled by LINCS. Energy minimization using steepest-descent algorithm was performed before production simulations. Then, 100 ns long MD simulations were performed (total simulation time 2.2 µs).

IV.6.4.2 New molecule parameters

Bonded and van der Waals (vdW) parameters for PZN compatible with Gromos 53a6 forcefield were obtained from the PRODRG webserver [START_REF] Schüttelkopf | PRODRG: a tool for high-throughput crystallography of protein-ligand complexes[END_REF]. The torsion angles between conjugated moieties were re-parameterized from quantum calculations [167]. All partial charges were assigned to atoms as obtained with the RESP method [START_REF] Bayly | A well-behaved electrostatic potential based method using charge restraints for deriving atomic charges: the RESP model[END_REF]. The ESP charges were obtained from B3LYP/aug-cc-pVTZ [359] calculations obtained on geometries optimized at the same level, with the Gaussian 09 software [270]. The RESP fit was carried out with the Antechamber package of AMBER 11 [274].

IV.6.4.3 Conformational analysis

To investigate the conformation of PZN, MD simulations and QM optimization were performed in both water and DMSO in order to match different experimental conditions. The MD simulations were carried out in a truncated octahedronal box filled with solvent molecules. The edges of the box were 7 Å from the solute. Although PZN possesses large πconjugated and rigid moieties, the C and N-terminal aminoacids as well as the central isoleucines allow numerous degrees of freedom. To completely investigate the complexity of the potential energy hypersurface of PZN within the simulations time scales, temperature may be artificially increased to cross the energetic barriers of all possible torsions. For that purpose a simulated annealing procedure was used. The temperature followed the hundred successive artificial heating/cooling (1ns) loops, a loop being in this case i) 100 ps heating from 310 K to 600 K, ii) 800 ps slow and linear cooling from 600 K to 310 K and iii) 100 ns equilibration at 310 K. The thermostat for PZN molecules was coupled independently from that of the solvent, and only PZN molecules were heated to avoid solvent evaporation.

Because of the large temperature fluctuations, the velocity rescaling algorithm was used as thermostat ( T = 0.1 ps). Unbiased free simulations in the same conditions but at constant temperature (310 K) were also conducted in order to be compared to simulated annealing simulations. The simulated annealing procedure appears sufficient and sometimes such a range of extreme temperatures is necessary to tackle the torsion angles responsible for the flexibility of PZN within simulations time scale.

In order to assess the quantity of π-stacking within PZN monomers or dimers, all distances between atoms involved in central π-conjugated moieties were computed. Then, the number of interatomic contacts was measured by counting distances smaller than a cutoff. The number of contacts involved in π-stacking was obtained by subtracting the number of contacts within a linear monomer showing no π-stacking to the total number of contacts. The cutoff (0.9 nm) was chosen as the number of contacts in the linear head-to-tail dimer is twice the one of the folded monomer.

QM calculations were performed at the density functional theory (DFT) level. Due to the importance of van der Waals interaction in the supramolecular re-arrangement of PZN, the ωB97-XD [360] functional was used with the 6-31+G(d,p) basis set to properly take dispersion effects into account. Solvent effects were treated with an implicit method, namely polarizable continuum model (PCM) in which the molecule is embedded in a shape-adapted cavity surrounded by a continuum dielectric characterized by its dielectric constant ( = 78.355 and 46.826 for water and DMSO, respectively).

IV.6.4.4 Membrane simulation

The model of membrane bilayer consisted of one hundred and twenty eight 1,2-dioleoyl-snglycero-3-phosphatidylcholine (DOPC) [361] surrounded by approximately five thousand four hundred water molecules. DOPC is a simple model for fluid membranes at body temperature. The transmembrane pore was created by placing 6 PZN molecules in a linear conformation and a head-to-tail disposition. The pore was merged with the phospholipids with the g_membed [362] program and then filled with water and ions.

Conclusion

Interactions between drugs or natural compounds and lipid bilayer membranes are essential to their activity, toxicity and bioavailability. Indeed, most drugs have to cross membranes before reaching their biological targets. The mechanism of action of some therapeutic classes is even directly or indirectly related to their behavior in lipid bilayers. Designing better, more available and less toxic drugs nowadays requires molecular rationalization of their mechanisms of action and establishments of solid structure-activity relationships. Because computational facilities have dramatically increased over the past decades as well as the quality of bilayer in silico models, MD simulations have appeared as an interesting alternative to experimental methods to tackle drug interaction with membrane. They have for instance been successfully applied to help rationalize the mechanism of action of several classes of compounds e.g., anesthetics, β-blockers, non-steroidal anti-inflammatory drugs, antioxidants and other classes targeting membrane proteins.

Supported by MD simulations, this work highlighted mechanisms of action of various active compounds, mainly lipid peroxidation inhibitors for which cooperative effects were shown to play crucial roles, and in a less extend antimicrobial compounds and drugs against Alzheimer disease. Although MD successfully revealed mechanisms that had not been previously elucidated, MD simulations have sometimes failed at describing some complex or subtle behaviors in membranes (e.g. orientation, diffusion). However, "shuttles" from theory to experiment and return is a virtuous circle allowing constant methodological improvements of in silico membrane models and MD methodologies. The next crucial step in methodological development is the use of lipid mixtures including cholesterol, as carprofens exemplified the strong influence of membrane composition and temperature (section IV.5). While simple bilayer models consisting of pure phosphatidylcholine may sometimes be appropriate to rationalize some mechanisms of action, more complex models are necessary to tackle the diversity of biological membranes. Indeed, the composition of biomembranes is highly variable from one cell line to another (e.g., skin or intestinal cells, healthy or cancer cells), or even within a cell between organelles or in lipid rafts. Therefore, simulations with complex lipid mixtures as close as possible to experimental or biological conditions are nowadays mandatory if one would use MD simulations as a predictive tool.

The ultimate goal of molecular modelling in pharmacology is prediction of biological activity for a wide range of compounds. This PhD work and the increasing literature in this field of research allow us to believe that in the next decade MD simulations will unanimously convince researchers and industrial partners to be a prime tool in chemistry, biochemistry, pharmacology and drug design.

Table S3: Logarithms of partition coefficients between DMPC membrane and water (Log K) measured experimentally (Exp.) and calculated by considered FFs and COSMOmic. The mean differences, mean absolute differences and Spearman`s rank order correlation coefficient are calculated with respect to experiment. Molecular interaction of natural compounds with lipid bilayer membranes:

Molecule

One of the key lockers to understand mechanisms of biological action of drugs and natural compounds is their capacity to incorporate/cross lipid bilayer membranes. In the light of demanding experimental techniques, in silico molecular modelling has become a powerful alternative to tackle these issues. In the past few years, molecular dynamics (MD) has opened many perspectives, providing an atomistic description of the related intermolecular interactions.

Using MD simulations, we have explored the capacity of several compounds (polyphenols, vitamins E and C, plantazolicin, carprofens) to incorporate lipid bilayer membranes. The different compounds were chosen according to their different biological functions, namely (i) antioxidant activity against lipid peroxidation, (ii) antimicrobial activity with the possibility of trans-membrane pore formation, and (iii) inhibition of enzymes involved in Alzheimer's disease. In order to rationalize their mechanisms of action, their position and orientation in membranes as well as their capacity to accumulate or permeate lipid bilayers were assessed.

Having in mind a predictive purpose in drug design for MD simulations, the accuracy of the results relies on the quality of the in silico membrane models. By ensuring relationships between experimental and theoretical data, methodological improvements have been proposed. In particular, force field selection, xenobiotic parameterization and bilayer constitution emerged as crucial factors to appropriately depict drug-membrane interactions. For the latter issue, lipid mixtures e.g., including cholesterol have been developed. Mots-clés : dynamique moléculaire, membranes lipidiques, composes naturels, antioxydants

  Appendices ........................................................................................................................ 148

Figure 1 :

 1 Figure 1: Schematic representation of a lipid bilayer cell membrane and its constituents. ....19

Figure 2 : 19 Figure 3 : 20 Figure 4 : 21 Figure 5 21 Figure 6 : 22 Figure 7 : 24 Figure 8 : 25 Figure 9 : 27 Figure 10 : 31 Figure 11 : 34 Figure 12 : 45 Figure 13 : 48 Figure 14 : 49 Figure 15 : 52 Figure 16 : 59 Figure 17 : 62 Figure 18 :

 2193204215216227248259271031113412451348144915521659176218 Figure 2: Protein to lipid ratio in various lipid bilayer membranes. Adapted from ref. [1]. ......19 Figure 3: Schematic representations of the different auto-associations of lipids. ..................20 Figure 4: Chemical structure of representative lipids ............................................................21 Figure 5 Phospholipid head group alcohols. OH group able to bond to phosphates are highlighted. ...........................................................................................................................21 Figure 6: Membrane models of 128 DPPC molecules, (A) below T m at 298 K and (B) above T m at 323 K. The difference in lipid chain ordering can be clearly distinguished. ...................22 Figure 7: Phospholipid distribution in mammals and yeast plasma membrane. The sterol / phospholipid ratio is mentionned for mammals (cholesterol, CHOL) and for yeast (ergosterol, ERG). Adapted from ref. [8]. .................................................................................................24 Figure 8: Chemical structures of cholesterol and ergosterol, two common sterol lipids. ........25 Figure 9: Influence of the lipid packing parameter P on the lipid bilayer curvature Reproduced from ref. [26]. ........................................................................................................................27 Figure 10: Interactions described in force fields. ...................................................................31 Figure 11: Different force field resolutions of a DOPC molecule: AA-FF, UA-FF and CG-FF. .............................................................................................................................................34 Figure 12: Review of experimental area per lipid for fluid DPPC at 50°C (black) and gelphase DPPC at 20°C (grey). Reproduced from [109]. ..........................................................45 Figure 13: Mechanisms of membrane crossing. ...................................................................48 Figure 14: Interaction of local anesthetics with membranes. A) Free energy profiles of the penetration of benzocaine and phenytoin through a POPC bilayer at 298 K (reproduced from ref. [126]). B) Slice of DPPC bilayer at 310 K, in presence of C) lidocaine, D) procaine and E) tetracaine showing disordering effect of LAs (adapted from ref [132]). .................................49 Figure 15: Hypothesis of the mechanism of gastric mucosa perturbation by NSAIDs. A) healthy mucosa, B) attack of H + ions (black) possible thanks to the membrane effects of NSAIDs (red). Reproduced from ref. [157]. ...........................................................................52 Figure 16: Structure of a dimyristoylphosphatidylcholine (DMPC) bilayer (background) with highlighted glycerol oxygens (red balls), choline nitrogens (blue) and phosphorus (dark yellow). The electron density profile (upper panel) contains labels for membrane thicknesses, i.e., head group to head group distance (D (HH) ), hydrocarbon core thickness (D (C) ) and Luzatti thickness (D (B) ) calculated as a ratio of volume per lipid (VPL) and area per lipid (APL). The free energy profile (lower panel) has highlighted water/lipids barrier ΔG wat , representing the affinity to the membrane, and penetration barrier ΔG pen . .......................................................59 Figure 17: Structural parameters of DMPC bilayer as predicted by MD simulations with various FFs compared to experimental values at 30°C shown as dotted lines [243]. APLarea per lipid, VPL -volume per lipid, D (HH) -electron -electron density peak distance, D (C)hydrocarbon core thickness, D (B) -Luzatti thickness. The error bars show the standard

64 Figure 20 : 66 Figure 21 : 67 Figure 22 : 70 Figure 23 : 75 Figure 24 : 76 Figure 25 : 80 Figure 26 : 83 Figure 27 : 85 Figure 28 : 85 Figure 29 : 87 Figure 30 : 87 Figure 31 : 90 Figure 32 : 91 Figure 33 : 92 Figure 34 : 93 Figure 35 :Figure 36 :Figure 37 :Figure 38 :Figure 39 :Figure 40 :Figure 41 :Figure 42 :Figure 43 :

 64206621672270237524762580268327852885298730873190329133923493353637383940414243 Figure 19: Experimental partition coefficients plotted against the respective calculated values (upper panel) with parameters of the linear fit, i.e. coefficient of determination, R 2 , slope (a) (standard deviation in bracket) and intercept (b). Slopes significantly differing from 1.0 and intercepts from 0.0 significantly on the probability level of 0.975 are highlighted in red. Each data point is labeled by a number, which corresponds to the number of the molecule in Table 3. The fitting parameters for GROMOS 43A1-S3 (G43A1-S3) and GAFFlipids re-calculated by omitting outliers (acetone and 2-nitrotoluene, in GAFFlipids and GROMOS 43A1-S3, respectively) are shown in blue. The bar charts (lower panel) depict the mean differences and the mean absolute differences. The patterned bars show values when excluding outliers. .............................................................................................................................................64 Figure 20: Water/lipid barriers ΔG wat and penetration barriers ΔG pen calculated by all FFs and COSMOmic vs. the values obtained with Slipids. .................................................................66 Figure 21: Mean absolute difference of free energy profile values with respect to Slipids as a function of distance from the middle of the membrane. .........................................................67 Figure 22: Chemical structures of lipocarbazole derivatives .................................................70 Figure 23: Partial densities of 1 (blue) and 2 (red). Densities of carbazole moieties are drawn in solid lines and those of active OH group in dotted lines. ...................................................75 Figure 24: Free energy profiles of 1 (blue) and 2 (red). .........................................................76 Figure 25: Chemical structures of studied polyphenols. ........................................................80 Figure 26: Representative snapshot of quercetin location below head groups, stabilized by Hbonds with carbonyl moieties of phospholipids and deep water molecules. ..........................83 Figure 27: Location of cyanidin in the membrane. The cationic form interacts with negatively charged phosphate groups, which P atoms are depicted as ochre spheres. .........................85 Figure 28: Representative snapshots of the interaction of the polyphenols studied with the lipid bilayer. ..........................................................................................................................85 Figure 29: Location and orientation of silybin. The active OH group on position 20 is close to polar head groups, on the left-hand side of the picture. ........................................................87 Figure 30: ∆G profile of polyhenols' DOPC membrane crossing as calculated with COSMOmic. .........................................................................................................................87 Figure 31: Antioxidant compounds evaluated in this study. The active antioxidant OH groups (prone to HAT) are shown in red. .........................................................................................90 Figure 32: Position of center of mass of vitC and quercetin, and the antioxidant OH group of vitE in DOPC. (A) individual molecules, (B) close contact pairs ............................................91 Figure 33: Geometries of the most stable associations as obtained from quantum DFT-D calculations. (A) quercetin:vitE, (B) vitC:vitE, (C) vitE:vitE, and (D) quercetin:vitC. ...............92 Figure 34: Fluorescence emission of vitE in liposomes with increasing concentrations of quercetin (0 to 100 µM). (A) Fluorescence spectra, (B) Stern-Volmer plot. VitE was excited at

22 Table 2 : 23 Table 3 :

 222233 Most common fatty acids found in membranes .......................................................Common phospholipids and their abbreviations ......................................................Molecules used in this study. The experimental partition coefficients (log K exp ) between water and DMPC are given from extensive dataset [227]. They are given as an average of experimental values in case of multiple source of individual partition coefficients (shown in brackets in the Method column). ...........................................................................58 Table

73 Table 7 :

 737 nitrotoluene in GROMOS 43A1-S3 and acetone in GAFFlipids). .................................................................................................65 Table 6: O-H and N-H BDE (kcal mol -1 ) for compounds 1 and 2 in the presence or absence of a PCM-type polar (water) or non-polar (benzene) solvent. ...................................................Lipid peroxidation inhibition for compounds 1 and 2, vitamin E and quercetin. IC 50 are given in µmol L -1 . ............................................................................................................74 Table

  ...............................................................................................................93 Table 10: Positions of COM and deuterated moieities of benzyl-and sulfonylcarprofen in POPC:SM:Chol and DOPC at 288 K and 310 K, from neutron diffraction and MD simulations. Free energy differences for membrane affinity (ΔG wat ) and crossing (ΔG pen ) under these various conditions are also reported. .................................................................................. 104 Table 11: Experimental and theoretical order parameters for the C-D bonds of both carprofen derivatives. ......................................................................................................................... 108

Figure 1 :

 1 Figure 1: Schematic representation of a lipid bilayer cell membrane and its constituents.

Figure 2 :

 2 Figure 2: Protein to lipid ratio in various lipid bilayer membranes. Adapted from ref. [1].

Figure 3 :

 3 Figure 3: Schematic representations of the different auto-associations of lipids.

Scheme 1 :Figure 4 :Figure 5

 145 Scheme 1: Representation of the different classes of lipids.

Figure 6 :

 6 Figure 6: Membrane models of 128 DPPC molecules, (A) below T m at 298 K and (B) above T m at 323 K. The difference in lipid chain ordering can be clearly distinguished.

Figure 8 :

 8 Figure 8: Chemical structures of cholesterol and ergosterol, two common sterol lipids.

  𝑉 𝑏𝑜𝑛𝑑𝑒𝑑 + 𝑉 𝑛𝑜𝑛-𝑏𝑜𝑛𝑑𝑒𝑑 𝑉 𝑏𝑜𝑛𝑑𝑒𝑑 = 𝑉 𝑏𝑜𝑛𝑑 + 𝑉 𝑎𝑛𝑔𝑙𝑒 + 𝑉 𝑑𝑖ℎ𝑒𝑑𝑟𝑎𝑙 𝑉 𝑛𝑜𝑛-𝑏𝑜𝑛𝑑𝑒𝑑 = 𝑉 𝑒𝑙𝑒𝑐𝑡𝑟𝑜𝑠𝑡𝑎𝑡𝑖𝑐 + 𝑉 𝑣𝑎𝑛 𝑑𝑒𝑟 𝑊𝑎𝑎𝑙𝑠 II.2.1. Bonded interactions

Figure 10 :

 10 Figure 10: Interactions described in force fields.

Figure 11 :

 11 Figure 11: Different force field resolutions of a DOPC molecule: AA-FF, UA-FF and CG-FF.

Figure 12 :

 12 Figure 12: Review of experimental area per lipid for fluid DPPC at 50°C (black) and gel-phase DPPC at 20°C (grey). Reproduced from [109].

Figure 13 :

 13 Figure 13: Mechanisms of membrane crossing.

Figure 14 :

 14 Figure 14: Interaction of local anesthetics with membranes. A) Free energy profiles of the penetration of benzocaine and phenytoin through a POPC bilayer at 298 K (reproduced from ref. [126]). B) Slice of DPPC bilayer at 310 K, in presence of C) lidocaine, D) procaine and E) tetracaine showing disordering effect of LAs (adapted from ref [132]).

Figure 15 :

 15 Figure 15: Hypothesis of the mechanism of gastric mucosa perturbation by NSAIDs. A) healthy mucosa, B) attack of H + ions (black) possible thanks to the membrane effects of NSAIDs (red). Reproduced from ref. [157].

  To date, numerous FFs have been developed for biomembranes, mostly focusing on structural and dynamical features of lipid bilayers. They were based on coarse-grained (e.g. MARTINI [212], SDK [213]), united-atom (e.g. Berger [214] and GROMOS 43A1-S3 [215]) and all-atom models (e.g. Slipids [216-218], CHARMM36 [219,220], GAFFlipids [221], LIPID11 [222], LIPID14 [223]

  [227]. The molecules were chosen to cover a wide range of partition coefficients (from -1.04 to 5.64 measured at temperatures from 20 to 40 °C) and to include common functional groups present in drugs such as hydroxyl, carbonyl, chloro, methyl, nitro and amino groups on aliphatic chains or aromatic benzene rings. The MD parameters of these molecules were prepared for individual FF, as recommended by their developers. Bonding and van der Waals parameters were taken from i) GAFF[229] for Slipids and GAFFlipids, ii) PRODRG[230] for Berger and GROMOS 43A1-S3 and iii) ParamChem[231,232] for CHARMM36. For CHARMM36, partial charges were also taken from ParamChem. Special attention was paid to the description of partial charges for Slipids, GAFFlipids, Berger and GROMOS 43A1-S3 FFs. For these FFs, the partial charges were derived using the restrained fit of electrostatic potential (RESP) procedure and the R.E.D. III software [233] using multiple conformations and multiple reorientations to ensure reproducibility of charge derivation, as ESP charges are sensitive to orientation[233,234]. Conformations were generated from 1 ns MD simulation in vacuum followed by clustering using the single linkage method. Only clusters representing more than 10% of the total number of conformations were taken into account. Then, energy minimization and electrostatic potential (ESP) charges were calculated for each conformation with Gaussian09 (rev. A02)[235] either according to theDuan model [236] (B3LYP/cc-pVTZ and PCM solvation in diethylether) for Slipids, Berger and GROMOS 43A1-S3 or according to the Cornell model [237] (HF/6-31G* in vacuum) for GAFFlipids. In recent work a polarization scheme was also applied [238], however we focused on a single set of partial charges for each system.

  5.1 software package with a 2 fs time step and periodic boundary conditions in all directions. Electrostatic interactions were treated by the particle-Mesh Ewald method [244] and bonds were constrained by the LINCS algorithm [245]. A Parrinello-Rahman barostat [246] was used for a semi-isotropic pressure coupling at 1 bar and compressibility of 4.5•10 -5 bar -1 and Nose-Hoover thermostat [247,248] at 310 K.

Figure 16 :

 16 Figure 16: Structure of a dimyristoylphosphatidylcholine (DMPC) bilayer (background) with highlighted glycerol oxygens (red balls), choline nitrogens (blue) and phosphorus (dark yellow). The electron density profile (upper panel) contains labels for membrane thicknesses, i.e., head group to head group distance (D (HH) ), hydrocarbon core thickness (D (C) ) and Luzatti thickness (D (B) ) calculated as a ratio of volume per lipid (VPL) and area per lipid (APL). The free energy profile (lower panel) has highlighted water/lipids barrier ΔG wat , representing the affinity to the membrane, and penetration barrier ΔG pen .

. 1 . 3 .

 13 Results and DiscussionIV.1.3.1 Structure of DMPC bilayer is well represented by all FFsDuring both unbiased and z-constraint simulations, most of the membrane structural parameters stayed reasonably close to experimental values[243], though most of the FFs produced a bilayer with thickness lower than that measured experimentally (Figure17). The values of area per lipid (APL) were reproduced reasonably well by all FFs considered here. The volume per lipid (VPL) predicted by GAFFlipids significantly differed from the other FFs. On the other hand, GAFFlipids showed head group distance (D (HH) ) and hydrocarbon thickness (D (C) ) in agreement with the experimental data. The Luzatti thickness (D (B) ), which depends on a ratio of VPL and APL (see Figure16for thickness explanation)[243], was again well reproduced by all other FFs but GAFFlipids (Figure17). In summary, all FFs tested in this study accurately reproduce the structural features of the DMPC bilayer reasonably well.

Figure 17 :

 17 Figure 17: Structural parameters of DMPC bilayer as predicted by MD simulations with various FFs compared to experimental values at 30°C shown as dotted lines [243]. APL -area per lipid, VPLvolume per lipid, D (HH) -electron -electron density peak distance, D (C) -hydrocarbon core thickness, D (B) -Luzatti thickness. The error bars show the standard deviation of data obtained from multiple simulations, all the graphs are scaled to show 20 % of deviation from experimental values.

Figure 18 :

 18 Figure 18: Order parameters experimentally measured (stars) and calculated by MD simulations with five FFs.

Figure 19 : 2 ,

 192 Figure 19: Experimental partition coefficients plotted against the respective calculated values (upper panel) with parameters of the linear fit, i.e. coefficient of determination, R 2 , slope (a) (standard deviation in bracket) and intercept (b). Slopes significantly differing from 1.0 and intercepts from 0.0 significantly on the probability level of 0.975 are highlighted in red. Each data point is labeled by a number, which corresponds to the number of the molecule in Table 3. The fitting parameters for GROMOS 43A1-S3 (G43A1-S3) and GAFFlipids re-calculated by omitting outliers (acetone and 2nitrotoluene, in GAFFlipids and GROMOS 43A1-S3, respectively) are shown in blue. The bar charts

Figure 20 :

 20 Figure 20: Water/lipid barriers ΔG wat and penetration barriers ΔG pen calculated by all FFs and COSMOmic vs. the values obtained with Slipids.

Figure 21 :

 21 Figure 21: Mean absolute difference of free energy profile values with respect to Slipids as a function of distance from the middle of the membrane.

Gabin

  Fabre, a,b Anne Hänchen c , Claude Calliste a , Karel Berka, b Srinivas Banala d , Michal Otyepka b , Roderich D. Süssmuth c and Patrick Trouillas d,e,*

Figure 22 :

 22 Figure 22: Chemical structures of lipocarbazole derivatives 1: R = (CH 2 ) 7 CH=CH(CH 2 ) 7 CH 3 2: "Carbazole": R = H 3: Carazostatin: R = n-C 7 H 15

  -1 (O-H BDE of 1 and 2) and the BDE of H-OH, LOO-H and DPPH-H (i.e., around 113, 85 and 80 kcal mol -1 , respectively).

Figure 23 :

 23 Figure 23: Partial densities of 1 (blue) and 2 (red). Densities of carbazole moieties are drawn in solid lines and those of active OH group in dotted lines.

Figure 24 :

 24 Figure 24: Free energy profiles of 1 (blue) and 2 (red).

  Advanced Technologies and Materials, Department of Physical Chemistry, Faculty of Science, Palacký University, tř. 17 listopadu 12, 771 46 Olomouc, Czech RepublicIV.3.1. IntroductionUnbalanced production of ROS (reactive oxygen species) is related to oxidative stress, for which long-term effects have been evidenced in aging and various diseases including cardiovascular, Alzheimer and liver diseases[260,283]. At the cellular level, ROS excess can oxidize lipids by lipid peroxidation (LPO), in particular degrading lipid bilayer membranes[158,261]. Antioxidants are known to regulate ROS excess; they are either endogenous enzymatic systems (e.g., catalase, superoxide dismutase and glutathione peroxidase) or exogenous small compounds (e.g. vitamin C, vitamin E, carotenoids and polyphenols).

Figure 25 :

 25 Figure 25: Chemical structures of studied polyphenols.

Figure 26 :

 26 Figure 26: Representative snapshot of quercetin location below head groups, stabilized by H-bonds with carbonyl moieties of phospholipids and deep water molecules.

Figure 27 :

 27 Figure 27: Location of cyanidin in the membrane. The cationic form interacts with negatively charged phosphate groups, which P atoms are depicted as ochre spheres.

Figure 28 :

 28 Figure 28: Representative snapshots of the interaction of the polyphenols studied with the lipid bilayer.

Figure 29 :

 29 Figure 29: Location and orientation of silybin. The active OH group on position 20 is close to polar head groups, on the left-hand side of the picture.

Figure 30 :

 30 Figure 30: ∆G profile of polyhenols' DOPC membrane crossing as calculated with COSMOmic.

Gabin

  Fabre a,b , Imene Bayach a , Karel Berka b , Marketa Paloncýová b , Marcelina Starok c , Claire Rossi c , Jean-Luc Duroux a , Michal Otyepka *,b and Patrick Trouillas *,b,d a LCSN EA1069, Faculté de Pharmacie, Université de Limoges, 2 rue de Dr. Marcland, 87025 Limoges, France b Regional Centre of Advanced Technologies and Materials, Department of Physical Chemistry, Faculty of Science, Palacký University, tř. 17 listopadu 12, 771 46 Olomouc, Czech Republic c FRE CNRS 3580, Génie Enzymatique et Cellulaire, Université de Technologie de Compiègne, CS 60319, 60203 Compiègne, France d UMR 850 INSERM, Univ. Limoges, Faculté de Pharmacie, 2 rue du Docteur Marcland, 87025 Limoges Cedex, FranceAbstractVitamins E, C and polyphenols (flavonoids and non-flavonoids) are major natural antioxidants capable of preventing damage generated by oxidative stress. Here we show the capacity of these antioxidants to form non-covalent association within lipid bilayers close to the membrane/cytosol interface. Antioxidant regeneration is significantly enhanced in these complexes.IV.4.1. CommunicationOver the last decades, natural antioxidants have attracted increasing interest, largely because they have been shown to exhibit preventive effects against various disorders caused by oxidative stress, including cardiovascular and neurodegenerative diseases, ageing and also certain cancers[260]. Despite recent progress in the field, there are still many open and fundamental questions concerning antioxidant mechanisms and biological targets, and their exact role in various pathologies is still under scrutiny[301]. A deep understanding of these mechanisms of action is mandatory for their safe and efficient usage in nutrition, health prevention, cosmetics and food preservation. Most of the known antioxidants are efficient scavengers of reactive oxygen species (ROS), which are overproduced during oxidative stress. Oxidation of lipids (namely lipid peroxidation, LPO) is a major process in oxidative stress, which is initiated by various endogenous (e.g., inflammation, enzymatic processes) or exogenous (e.g., radiation, smoking, pollution) effects. The propagation stage ofLPO [261] can be inhibited by lipophilic or amphiphilic antioxidants sufficiently incorporated in lipid bilayers[169,273]. In addition, hydrophilic and polar antioxidants are able to scavenge ROS that diffuse toward membranes, thus inhibiting the initiation stage of LPO. Vitamin E (α-tocopherol, henceforth referred to as vitE) [302],

Figure 31 :

 31 Figure 31: Antioxidant compounds evaluated in this study. The active antioxidant OH groups (prone to HAT) are shown in red.

Figure 32 :

 32 Figure 32: Position of center of mass of vitC and quercetin, and the antioxidant OH group of vitE in DOPC. (A) individual molecules, (B) close contact pairs

Figure 33 :

 33 Figure 33: Geometries of the most stable associations as obtained from quantum DFT-D calculations. (A) quercetin:vitE, (B) vitC:vitE, (C) vitE:vitE, and (D) quercetin:vitC.

Figure 34 :

 34 Figure 34: Fluorescence emission of vitE in liposomes with increasing concentrations of quercetin (0 to 100 µM). (A) Fluorescence spectra, (B) Stern-Volmer plot. VitE was excited at λ exc = 291 nm after incorporation into liposomes. The control condition was performed by incubation of vitE (50 µM) with vitE-free DOPC liposomes. Prior to quercetin addition, the non inserted VitE molecules were eliminated from the liposome suspension by double centrifugation and resuspension.



  PME: particle mesh Ewald  RESP: restrained fit of electrostatic potential  ROS: reactive oxygen species  SPC/E: extended single point charge model  vitC: vitamin C  vitE: vitamin E IV.4.2

Figure 35 :

 35 Figure 35: Chemical structures of carprofen derivatives.

Figure 36 :

 36 Figure 36: Position distributions of the center of mass (COM), the deuterated and carboxyl moieties of carprofen derivatives in POPC:SM:Chol (L o ) and DOPC (L d ) at 288 K and 310 K.

Figure 37 :

 37 Figure 37: Free energy profiles of benzyl-and sulfonylcarprofen membrane crossing in A) POPC:SM:Chol at 288 K, B) POPC:SM:Chol at 310 K, and C) DOPC at 310 K.

Figure 38 :

 38 Figure 38: Orientation distribution given as the vector 𝑽 𝟏 ⃗⃗⃗⃗ connecting benzyl to carbonyl moieties of each carprofen derivative for outer-located (A and B) and inner-located (C and D) compounds.Orientation analysis was averaged over the second half of all simulations, independently of membrane composition and temperature.

Figure 39 :

 39 Figure 39: Neutral and charged chemical structures of plantazolicin

Figure 40 :

 40 Figure 40: Summary of the conformational study showing schematic initial and final geometries. The number of contacts involved in π-stacking in the final geometry is quoted in brackets.

Figure 41 :

 41 Figure 41: Conformations of the final structures for simulations: A) linear monomer, B) folded monomer from free MD, C) folded monomer after QM optimization, D) head-to-tail dimer (dimer A), and E) head-to-center folded dimer (dimer C).

Figure 42 :

 42 Figure 42: Plantazolicin dimer as obtained from simulated annealing simulation with distance restraints.

Figure 43 :

 43 Figure 43: Snapshots of the pore after 1 ns equilibration (A) and after a 12ns MD simulation (B).

Figure S1 :

 S1 Figure S1: Mean difference of free energy values in various membrane depths in respect to free energy calculated in Slipids force field shows regions with increased (below zero) and reduced (above zero) affinity to that region.

Figure S5 :

 S5 Figure S5: Fluorescence emission spectra of quercetin. (A) 100 µM, λ exc = 291 nm, in HBS solution or incubated with vitE-free liposomes, (B) 100 µM, λ exc = 370 nm, in HBS solution or incubated with vitE-free liposomes, (C) 0 to 100 µM, λ exc = 370 nm, incubated with liposomes with embedded vitE.

Figure

  Figure S6: [I 0 /I -1]/[Q] = f([quercetin]) plot. The fluorescence intensity of vitE in liposomes was recorded at λ exc = 291 nm in the absence (I 0 ) and in the presence (I) of increasing concentrations of quercetin.

Figure S7 :

 S7 Figure S7: Parameterization of the dihedral angle χ of quercetin. (A) Dihedral angle χ in the quercetin structure. (B) Determination of dihedral angle potential according to Equation (2).

Figure S8 :

 S8 Figure S8: Fluorescence emission spectra of vitE exposed to increasing concentrations of quercetin (0 to 100 µM) in HBS solution. VitE was excited at λ exc = 291 nm.

  Keywords: molecular dynamics, lipid bilayer membranes, natural compounds, antioxidantsInteractions moléculaires des composés naturels avec les membranes lipidiquesUne des clés pour comprendre les mécanismes d'action biologiques des molécules naturelles et thérapeutiques est leur faculté à incorporer ou traverser les membranes lipidiques. Parce que les méthodes expérimentales sont parfois couteuses et réponde partiellement aux questions posés par les interactions composé-membrane, la modélisation moléculaire est devenue une sérieuse alternative. Les simulations de dynamique moléculaire ont ouvert de nombreuses perspectives ces dernières années en offrant la possibilité de décrire ces interactions intermoléculaires au niveau atomique. À l'aide de ces simulations, nous avons évalué la capacité de plusieurs composés (polyphénols, vitamines E et C, plantazolicine et carprofènes) à s'incorporer dans les membranes. Ces molécules ont été choisies pour leurs activités biologiques diverses, à savoir (i) activité antioxydante, précisément inhibition de la peroxydation lipidique, (ii) activité antibiotique et possibilité de former un pore transmembranaire, et (iii) inhibition d'enzymes impliquées dans la maladie d'Alzheimer. Leurs positions et orientations ainsi que leur capacité à s'accumuler ou à traverser les membranes ont été évaluées pour comprendre leurs mécanismes d'action. Dans le but d'utiliser les simulations de dynamique moléculaire en drug design, l'accent a été mis sur la précision des calculs, qui dépend de la qualité sous-jacente du modèle utilisé. En corrélant données expérimentales et théoriques, la méthodologie de nos modèles a été systématiquement revisitée. Le choix du champ de force, les paramètres des composés étudiés ainsi que la composition de la membrane sont en particulier apparus comme d'importants facteurs dans la description des interactions entre les molécules naturelles et thérapeutiques et les membranes. Des mélanges de lipides contenant du cholestérol ont notamment été utilisés et ont montré un impact significatif sur les résultats obtenus.

Table 1 :

 1 Most common fatty acids found in membranes

	C:D	Common name	n-x	Formula
	12:0	Lauric acid	-	CH 3 (CH 2 ) 10 COOH
	14:0	Myristic acid	-	CH 3 (CH 2 ) 12 COOH
	16:0	Palmitic acid	-	CH 3 (CH 2 ) 14 COOH
	18:0	Stearic acid	-	CH 3 (CH 2 ) 16 COOH
	14:1	Myristoleic acid	n-5 CH 3 (CH 2 ) 3 CH=CH(CH 2 ) 7 COOH
	16:1	Palmitoleic acid	n-7 CH

3 (CH 2 ) 5 CH=CH(CH 2 ) 7 COOH 18:1 Oleic acid n-9 CH 3 (CH 2 ) 7 CH=CH(CH 2 ) 7 COOH 18:2 Linoleic acid n-6 CH 3 (CH 2 ) 4 CH=CHCH 2 CH=CH(CH 2 ) 7 COOH 20:4 Arachidonic acid n-6 CH 3 (CH 2 ) 4 (CH=CHCH 2 ) 4 (CH 2 ) 2 COOH

Table 2 :

 2 Common phospholipids and their abbreviations

  COSMO-RS software (COnductor like Screening MOdel for Real Solvents)[103] calculates chemical potentials 𝜇 in liquids based on DFT calculations.COSMOmic [104] is a specialized version of COSMO-RS applied to micelles. It uses statistical thermodynamics instead of MD simulations to evaluate the free energy profile of a drug crossing a lipid bilayer.The underlying principle of COSMO-RS is the evaluation with DFT methods of charge density 𝜎 at the surface of a molecule. The distribution of the charge density, 𝑝(𝜎), is the socalled 𝜎-profile. It should be noted that the 𝜎-profile of a mixture of compounds is the weighted sum of the respective 𝜎-profiles. Using statistical thermodynamics, and knowing the 𝜎-profile of the solvent 𝑝 𝑠 (𝜎′), it is possible to calculate the 𝜎-potential 𝜇

𝑠 (𝜎) of a solute embedded in a solvent

[103]

. The 𝜎-potential corresponds to the affinity of the solvent 𝑠 to a solute surface having a polarity 𝜎. Then, the chemical potential for the whole solute 𝑋 in a solvent 𝑆 is obtained by integrating over the polarity 𝜎:

𝜇 𝑆 𝑋 (𝜎) = 𝜇 𝑐𝑜𝑚𝑏 𝑋 + 𝜇 𝑑𝑖𝑠𝑝 + ∫ 𝑝 𝑋 (𝜎) 𝜇 𝑠 𝑑𝜎

where 𝜇 𝑐𝑜𝑚𝑏 𝑋 and 𝜇 𝑑𝑖𝑠𝑝 are additional combinatorial and dispersive contributions

[104]

. For the special case of lipid bilayers, elastic deformation and 𝜁 potential contributions are also added

[104]

.

  resonance (ESR) for spin-labeled molecules [106]; (iii) steady state or time-dependent fluorescence techniques that measure quenching of fluorescent probes (e.g., laurdan, prodan) [106,107]; (iv) 1D or 2D nuclear magnetic resonance (NMR) in solution or in solid state with various nuclei ( 1 H, 2 H, 13 C, 19 F or 31 P); (v) atomic force microscopy (AFM) [108]; or (vi) Fourier transform infrared spectroscopy (FT-IR) [106].

  . The huge computer cost of MD simulations has motivated many researchers to develop less expensive approaches to estimate thermodynamic properties of molecule-membrane interaction. An example of such approaches is the COSMOmic [224] tool of the COSMOtherm program[225], which is based on the conductor-like screening model for realistic solvation (COSMO-RS) theory[226]. It was repeatedly shown that COSMOmic provides thermodynamics of molecule-membrane interactions in good agreement with experimental data[227,228]. On the other hand, this implicit approach loses the fine time insight into the interaction, which is provided by MD simulations.

Table 3 :

 3 Molecules used in this study. The experimental partition coefficients (log K exp ) between water and DMPC are given from extensive dataset[227]. They are given as an average of experimental values in case of multiple source of individual partition coefficients (shown in brackets in the Method column).

	Nr. Compound	log K exp	Method	Ref
	1	glycerol	-1.04	Ultracentrifugation	[239]
	2	methanol	-0.53	Ultracentrifugation	[239]
	3	acetone	0.06	Ultracentrifugation (0.02, 0.10)	[239]
	4	1-butanol	0.51	Ultracentrifugation (0.54)	[239]
				Non-depletion PA-SPME (0.45)	[240]
	5	benzylalcohol	1.14	Ultracentrifugation	[239]
	6	aniline	1.63	Non-depletion PA-SPME	[240]
	7	2-nitrotoluene	2.41	Non-depletion PA-SPME	[240]
	8	p-xylene	2.98	Non-depletion PA-SPME	[240]
	9	4-chloro-3-	3.34	Non-depletion PA-SPME	[240]
		methylphenol
	10 2,4,5-trichloroaniline	4.16	Non-depletion PA-SPME	[240]
	11 hexachlorobenzene	5.64	n-hexane passive dosing (5.43)	[241]
				PDMS sheet dosing (5.90)
				SPCE-PDMS passive sampling
				(5.59)

Table 4 :

 4 Simulation parameters. R coulomb is a short-range electrostatic cut-off, long-range electrostatics are evaluated by PME, R vdw is Lennard-Jones cut-off, in case of switching off the Lennard-Jones interactions, the switching begins at R vdw-switch . In case of CHARMM36 and Slipids, we tested the structural parameters also using different cut-off lengths (in brackets, not affecting the total CPU time in this table). CPU hours/project display the total CPU hours for the calculations -for obtaining the topologies and 30 ns z-constraint simulations for MD simulations and for DFT calculations and final free energy profile calculation in case of COSMOmic. The detailed CPU times are in TableS2.

	Force Field	R coulom b	R vdw	R vdw-switch	Bond constraints	Water model	RESP set method/basis	CPUh/ project
		nm
	Berger	1.0	1.0	-	All-bonds	SPC/E[249,	B3LYP/cc-	21,200
						250]	pVTZ
	GROMOS	1.0	1.6	-	All-bonds	SPC/E[249,	B3LYP/cc-	34,400
	43A1-S3					250]	pVTZ
	CHARMM36	1.4	1.4	0.8	H-bonds	CHARMM		145,200
		(1.2)	(1.2)			TIP3P[251]
	Slipids	1.0	1.5	1.4	All-bonds	TIP3P[252] B3LYP/cc-	71,300
		(1.0)	(0.9)	(0.8)			pVTZ
	GAFFlipids	0.8	0.8	-	H-bonds	TIP3P[252] HF/6-31G*	44,900
	COSMOmic							3

1.2.6 Statistical evaluation

  Predicted log 𝐾 𝑐𝑎𝑙𝑐 were compared to the log 𝐾 𝑒𝑥𝑝 experimental values in terms of mean

	difference	(1/𝑁 ∑ (log 𝐾 𝑐𝑎𝑙𝑐,𝑖 -log 𝐾 𝑒𝑥𝑝,𝑖 ) 𝑁 𝑖	)	and	mean	absolute	difference

(1/𝑁 ∑ |log 𝐾 𝑐𝑎𝑙𝑐,𝑖 -log 𝐾 𝑒𝑥𝑝,𝑖 | 𝑁 𝑖

Table 5 :

 5 Mean differences and mean absolute differences of water/lipids ΔG wat and penetration ΔG

	Force Field	ΔG wat		ΔG pen
		Mean difference	Mean absolute difference	Mean difference	Mean absolute difference
		kcal/mol	kcal/mol	kcal/mol
	Berger	1.94	2.09	0.14	1.06
	CHARMM36	-0.27	0.72	-0.15	0.89
	GAFFlipids	0.02 (0.14)	0.72 (0.68)	1.04 (0.33)	1.33 (0.65)
	GROMOS 43A1-S3	-0.34 (-1.07)	1.65 (1.12)	-0.35 (-0.29)	1.28 (1.31)
	COSMOmic	0.12	0.91	-0.73	0.91

pen barriers with respect to data obtained from Slipids FF. The values in brackets show the differences with excluded outlier (2-nitrotoluene in GROMOS 43A1-S3 and acetone in GAFFlipids).

  e.g., superoxide dismutase, glutathione peroxidase, and catalase) and small compounds (e.g., glutathione). Exogenous antioxidants, contained in food (e.g., vitamins C and E, polyphenols, carotenoids) or food supplementation, also contribute to the total antioxidant action. Long-term effects of oxidative stress have extensively been studied over the past years; they have been shown to be responsible for various diseases e.g.

	cardiovascular, Alzheimer and liver diseases [260]. In this context, lipid peroxidation (LPO) is
	one of the most important processes involving free radicals and it is directly implicated in
	various diseases such as atherosclerosis, rheumatoid arthritis and neurodegenerative
	diseases [158,261]. To discover new LPO inhibitors is of particular importance in order to
	prevent those diseases. From a clinical point of view extensive research deals with new
	antioxidants being able to decrease lesions induced by ischemia/reperfusion in organ
	transplantation [262]. It is also a challenge of major importance in cosmetics and food
	industries.
	A series of lipocarbazole derivatives was isolated from the bacterium Tsukamurella
	pseudospumae Acta 1857 [263]. These compounds were later synthetized by a series of
	metal-catalyzed reactions [264]. Due to their structural analogy with carazostatin [265], an
	effective in vivo antioxidant, the antioxidant capacity of lipocarbazole is under scrutiny in this
	article. It was found that lipocarbazole A3 (1) (Figure

The subsequent ROS overproduction can be induced by many factors including UV light, hypoxia, cytokines, chemotherapy and highenergy radiation

[259]

. Various endogenous antioxidant systems regulate ROS production, namely enzymes (

Table 6

 6

	: O-H and N-H BDE (kcal mol	-1 ) for compounds 1 and 2 in the presence or absence of a PCM-
	type polar (water) or non-polar (benzene) solvent.
	Compounds	1		2
	O-H BDE	N-H BDE	O-H BDE	N-H BDE
	In vacuo	79.5	86.5	80.0	88.3
	water	76.7	85.9	77.0	87.1
	benzene	78.9	88.4	79.3	88.4

Table 7 :

 7 Lipid peroxidation inhibition for compounds 1 and 2, vitamin E and quercetin. IC 50 are given in µmol L

	-1 .
	Compounds	IC 50
	1	0.16 ± 0.03
	2	1.00 ± 0.13
	vitamin E	0.80 ± 0.03
	quercetin	0.20 ± 0.03

Table 8 :

 8 Characterization of the interaction between a series of polyphenols and DOPC bilayer, as given by their positions, orientations, dipole moments, and free enthalpies of partition (ΔG part ) and penetration (ΔG pen ) are reported.

	Compound	Position MD (nm)	Position COSMOmic (nm)	Orientation of longest axis (degree)	Dipole moment (Debye)	ΔG part (kcal mol -1 )	ΔG pen (kcal mol -1 )
	catechin	1.50 ± 0.25	1.85	89 ± 24	4.38	-5.6	13.2
	cyanidin	1.30 ± 0.20	1.76	95 ± 25	4.60	-11.4	21.6
	cyanidin-3-O-glucoside	1.68 ± 0.34	1.76	73 ± 27	5.46	-10.8	27.7
	galangin	1.37 ± 0.25	1.44	140 ± 20	2.53	-7.9	4.3
	luteolin	1.53 ± 0.32	2.26	89 ± 24	9.26	-7.1	10.7
	quercetin	1.50 ± 0.30	1.60	95 ± 31	3.58	-7.5	8.8
	resveratrol	1.46 ± 0.21	1.93	85 ± 22	3.94	-5.7	10.1
	rosmarinate	2.20 ± 0.25	2.34	74 ± 33	15.65	-9.8	36.3
	rosmarinic acid	-	2.26	-	-	-10.2	16.8
	rutin	1.40 ± 0.30	2.50	76 ± 19	7.93	-6.8	26.3
	secoisolariciresinol	1.28 ± 0.21	1.44	92 ± 20	5.28	-6.8	5.5
	silybin	1.43 ± 0.22	1.60	50 ± 24	4.49	-8.0	7.4
	taxifolin	1.32 ± 0.15	2.26	98 ± 29	4.62	-7.2	17.0

Table 9 :

 9 Association energies and enthalpies (kcal.mol

-1

3 Liposome formation and fluorescence DOPC

  , vitE and DMSO were purchased from Sigma. Quercetin dihydrate was from Merck Millipore. HBS solution (HEPES 20 mM pH 7.4, NaCl 150 mM) was prepared in Milli-Q water (resistivity higher than 18.2 M.cm) and filtered (0.22 µM). All chemicals were analytical grade reagents.

	780, has been
	validated by high-level SCS-MP2 calculations and experimental values on polyphenol non-
	covalent complexes [295,334]. The COSMO implicit solvent was used to model benzene
	(non-polar) and water (polar) solvation. Bond dissociation enthalpy calculations were
	performed with B3P86/6-311+G(d,p) as previously reported [159].
	IV.4.3.

  1,2 , Anna Itkin 3 , Evgeniy Salnikov 3 , Michal Otyepka 2 , Norbert A. Dencher4 , The currently available medications are symptomatic and their efficacy is questioned. One of the characteristic histopathological markers of AD is the presence of Aβ peptides that consist of 40 or 42 amino acids (Aβ 40 and Aβ 42 ). Recent research suggests that soluble oligomers of Aβ are responsible for AD symptoms[337]. Aβ peptides originate from the cleavage of the amyloid precursor protein (APP) by γ-secretase, a large transmembrane enzyme complex[338]. They are continuously produced but harmless in the healthy brain. In case of concentration increase, they acquire the capacity to aggregate and form plaques. An increase in either total levels of Aβ or relative concentrations of Aβ 40 and Aβ 42 have been implicated in AD pathogenesis. One of the promising strategies for AD therapy is modulation of the γ-secretase activity to control the ratio of Aβ fragments without affecting other activities of this enzyme complex (e.g., cleavage of Notch). As for most of drugs, γ-secretase inhibitors must be sufficiently selective to avoid interactions with other enzymes.

	Thomas Hauß 5 , Patrick Trouillas* ,2,6 , Burkhard Bechinger* ,3
	1 LCSN EA1069, Faculté de Pharmacie, Université de Limoges, 2 rue de Dr. Marcland, 87025
	Limoges, France
	2 Regional Centre of Advanced Technologies and Materials, Department of Physical Chemistry,
	Faculty of Science, Palacký University, tř. 17 listopadu 12, 771 46 Olomouc, Czech Republic
	3 Université de Strasbourg / CNRS UMR 7177, Institut de chimie, 1, rue Blaise Pascal, 67070
	Strasbourg, France
	IV.5.1. Introduction
	Alzheimer's disease (AD) affects ca. 20 million people worldwide [336] and it is expected to
	double or triple by 2030 or 2050, respectively, if no significant progresses are made in the
	research for new efficient treatments.

Table 10 :

 10 Positions of COM and deuterated moieities of benzyl-and sulfonylcarprofen in POPC:SM:Chol and DOPC at 288 K and 310 K, from neutron diffraction and MD simulations. Free energy differences for membrane affinity (ΔG wat ) and crossing (ΔG pen ) under these various conditions are also reported.

			288 K			310 K
		POPC:SM:Chol	DOPC	POPC:SM:Chol	DOPC
		benzylca	sulfonylc	benzylca	sulfonylc	benzylca	sulfonylc	benzylca	Sulfonylc
		rprofen	arprofen	rprofen	arprofen	rprofen	arprofen	rprofen	arprofen
	Deuterated moiety
	position from neutron	2.64	1.99	-	-	-	1.58	-	-
	diffraction (nm)
	Deuterated moiety position from free MD simulations (nm)	2.5 ± 0.4 2.7 ± 0.5 2.2 ± 0.3	3.0 ± 0.4 1.2 ± 0.2	1.6 ± 0.3 2.7 ± 0.4	2.9 ± 0.4 1.6 ± 0.4	1.4 ± 0.3 1.1 ± 0.2
	COM position from free MD simulations	2.6 ± 0.4 2.6 ± 0.2 2.3 ± 0.5	2.7 ± 0.2 1.7 ± 0.2	1.7 ± 0.3 2.5 ± 0.2	2.7 ± 0.3 1.6 ± 0.3 1.8 ± 0.2
	COM position from ΔG profiles (nm)	2.6	2.4	-	-	1.8	3.0	2.4	1.6
	Membrane affinity ΔGwat (kcal mol -1 )	-7.5 ± 1.4	-8.6 ± 2.8	-	-	-9.0 ± 0.7	-4.7 ± 0.3	-6.7 ± 0.1	-10.4 ± 0.1
	Barrier for membrane crossing ΔGpen (kcal mol -1 )	46.7 ± 9.0	53.8 ± 6.4	-	-	18.7 ± 2.3	39.2 ± 3.3	17.8 ± 0.1	20.1 ± 0.1

Table 11 :

 11 Experimental and theoretical order parameters for the C-D bonds of both carprofen derivatives.

	C-D bond	Experimental Δν Q (kHz)	Experimental |S CD |	Theoretical S CD * (COM < 2 nm)	Theoretical S CD (COM > 2 nm)
	benzylcarprofen C β -D 1	20 ± 10	0.09	-0.39	-0.06
	benzylcarprofen C β -D 2	20 ± 10	0.09	-0.26	0.09
	benzylcarprofen C δ -D	46 ± 5	0.19	-0.18	0.01
	benzylcarprofen C ε -D	46 ± 5	0.19	-0.18	0.01
	benzylcarprofen C ξ -D	66 ± 3	0.28	-0.27	0.12
	sulfonylcarprofen C-D 3	56	0.23	-0.25	-0.04

*

These data are taken only from the last 150 ns of the MD simulation, for one molecule which orientation matches with experiment.

Table S4 :

 S4 Penetration barriers ΔG pen calculated by considered FF and COSMOmic. The mean differences and mean absolute differences are calculated with respect to values from Slipids FF.

	Molecule				ΔG pen (kcal/mol)
		Berger	Slipids	CHARMM36	GAFFlipids	GROMOS 43A1-S3	COSMOmic
	glycerol	6.58	5.81	6.92	6.04	9.14	5.03
	methanol	5.05	3.09	3.53	3.46	2.11	3.57
	acetone	1.80	1.22	1.31	9.38	0.91	0.90
	1-butanol	4.73	1.66	2.46	2.77	3.12	2.16
	benzylalcohol	1.70	3.60	3.20	2.56	2.60	2.12
	aniline	2.30	3.03	2.84	2.47	2.43	1.93
	2-nitrotoluene	2.23	2.03	1.12	2.26	1.08	1.99
	xylene	0.00	0.00	0.94	0.29	0.34	0.03
	4-chloro-3-methylphenol	2.31	2.90	1.00	4.66	0.58	0.72
	2,4,5-trichloroaniline	1.51	3.18	0.90	3.59	0.58	0.24
	hexachlorobenzene	0.00	0.17	0.85	0.68	0.00	0.00
	Statistics
	Mean difference	0.14		-0.15	1.04 (0.33)	-0.35 (-0.29)	-0.73
	Mean absolute difference	1.06		0.89	1.33 (0.65)	1.28 (1.31)	0.91

Table S5 :

 S5 Water/lipids barriers ΔG wat calculated by all FF and COSMOmic. The mean differences and mean absolute differences are calculated with respect to values from Slipids FF.

	Molecule				ΔG wat (kcal/mol)
		Berger	Slipids	CHARMM36	GAFFlipids	GROMOS 43A1-S3	COSMOmic
	glycerol	0.42	0.77	0.04	1.10	0.67	0.93
	methanol	0.25	0.09	0.10	0.40	0.34	0.73
	acetone	1.58	1.09	0.00	0.00	1.08	0.92
	1-butanol	2.54	3.02	2.71	3.57	0.91	2.62
	benzylalcohol	6.65	2.24	2.47	2.23	1.76	3.11
	aniline	4.32	2.19	1.91	3.03	1.02	3.24
	2-nitrotoluene	9.54	2.94	4.43	4.16	9.89	4.89
	xylene	5.16	4.95	5.70	4.29	3.02	5.92
	4-chloro-3-methylphenol	5.79	5.54	4.33	6.38	4.02	4.50
	2,4,5-trichloroaniline	8.72	6.96	6.65	5.90	4.78	4.89
	hexachlorobenzene	14.95	8.73	7.20	7.72	7.26	8.05
	Statistics
	Mean difference	1.94		-0.27	0.02 (0.14)	-0.34 (-1.07)	0.12
	Mean absolute difference	2.09		0.72	0.72 (0.68)	1.65 (1.12)	0.91

Table S6 :

 S6 Positions of free energy minima calculated by all FF and COSMOmic. The mean differences and mean absolute differences are calculated with respect to values from Slipids.

	Molecule			Position of free energy minimum (nm)
		Berger Slipids CHARMM36 GAFFlipids GROMOS 43A1-S3 COSMOmic
	glycerol	1.80	1.30	2.89	1.40	2.40	2.11
	methanol	1.21	2.51	2.71	1.30	2.11	1.90
	acetone	1.19	0.99	3.10	3.19	0.89	1.90
	1-butanol	0.89	0.99	1.01	1.10	1.98	1.11
	benzylalcohol	0.91	1.00	1.01	1.11	0.89	1.11
	aniline	1.19	1.20	1.01	1.29	1.80	1.11
	2-nitrotoluene	0.91	1.00	0.81	1.30	0.90	0.89
	xylene	0.00	0.01	1.00	0.71	0.61	0.75
	4-chloro-3-methylphenol	1.20	1.00	1.00	1.11	0.59	1.04
	2,4,5-trichloroaniline	0.90	1.00	0.90	1.10	0.58	1.04
	hexachlorobenzene	0.00	0.49	0.71	0.80	0.00	0.00
	Statistics
	Mean difference	-0.12		0.42	0.26 (0.07)	0.11 (0.14)	0.13
	Mean absolute difference	0.28		0.51	0.49 (0.31)	0.48 (0.52)	0.37

Table S7 :

 S7 Diffusion coefficients of antioxidants in x,y plane and on z-axis.

		Lateral diffusion coefficient (x,y-	Diffusion coefficient on z-axis (10 -8
		plane) in 10 -8 cm 2 .s -1	cm 2 .s -1 )
	vitC	17 ± 2	2 ± 2
	quercetin	17 ± 2	1 ± 1
	vitE	22 ± 5	3 ± 2

Table S8 :

 S8 Characterization of antioxidant pairs: intermolecular distance, position in the membrane, and ratio of the time spent in contact.

	Pairs	Average distance between the centers of mass of each antioxidant (nm)	Average distance between the centers of mass of the pair (nm) and bilayer center	Ratio for the pair of time spent closer than the 10 Å cutoff
	vitC:vitC	0.66 ± 0.12	1.94 ± 0.26	96 ± 5 %
	vitC vitE	0.74 ± 0.16	1.77 ± 0.13 1.28 ± 0.12	90 ± 12 %
	vitC quercetin	0.73 ± 0.18	1.88 ± 0.12 1.77 ± 0.36	95 ± 9 %
	vitE:vitE	0.77 ± 0.28	1.21 ± 0.06	87 ± 14%
	vitE quercetin	0.72 ± 0.18	1.15 ± 0.07 1.49 ± 0.27	93 ± 13 %
	quercetin:quercetin	0.75	2.09 ± 0.14	99%

Table S9 :

 S9 Occurrence of antioxidant pairs involving quercetin and vitE.

	Pairs	Number observed (ratio) of	complexes	Number expected for a random of complexes distribution
	vitE:vitE	45 (30%)		40
	quercetin:vitE	65 (43%)		75
	quercetin:quercetin	40 (27%)

Department of Chemistry, Physical Biochemistry, Technical University Darmstadt, Darmstadt,

Acknowledgements Acknowledgement

This work was supported by the Grant Agency of the Czech Republic (P208/12/G016), the Operational Program Research and Development for Innovations-European Regional Development Fund (CZ.1.05/2.1.00/03.0058), the Operational Program Education for Competitiveness-European Social Fund (CZ.1.07/2.3.00/20.0017), and a student project of Palacký University (IGA_PrF_2014023). PT and GF thank the "Conseil Régional du Limousin" for financial support and CALI (CAlcul en LImousin).

Acknowledgements

The authors thank Cali (Calcul en Limousin) for computer facilities. A.H. gratefully acknowledges the Fonds der Chemischen Industrie for a doctoral fellowship. P.T. and R.S. thank the COST CM1004 support. The project was supported the Cluster of Excellence "Unifying Concepts in Catalysis" coordinated by the Technische Universität Berlin. Financial support from the Czech Science Foundation (P208/12/G016), the Ministry of Education, Youth and Sports of the Czech Republic (project LO1305), the Operational Program Education for Competitiveness-European Social Fund (project CZ.1.07/2.3.00/20.0058 of the Ministry of Education, Youth and Sports of the Czech Republic), and IGA_PrF_2015_027 of Palacky University is also gratefully acknowledged.

Acknowledgements

The authors thank CALI (CAlcul en LImousin) for computational resources. Financial support from the Czech Science Foundation (P208/12/G016), the Ministry of Education, Youth and Sports of the Czech Republic (project LO1305) and the Operational Program Education for Competitiveness-European Social Fund (project CZ.1.07/2.3.00/20.0058 of the Ministry of Education, Youth and Sports of the Czech Republic) is also gratefully acknowledged. Acknowledgements CALI, Czech Science Foundation (P208/12/G016), Ministry of Education, Youth and Sports of the Czech Republic (project LO1305), Operational Program Education for Competitiveness-European Social Fund (project CZ.1.07/2.3.00/20.0058 of the Ministry of Education, Youth and Sports of the Czech Republic), Marie Curie Research Training Network CHEBANA, FP7 ITN 2010-264772, ADF.

Acknowledgements

The authors thank Région Limousin, COST Chemistry project and RCPTM for funding and Cali (Calcul en Limousin) for computational facilities.

Abbreviations