remercier Eulalie, notre dernière venue, pour son aide, et son incroyable bonne humeur. Grâce à vous, nous avons lancé une dynamique dans le labo et avons maintenant un groupe incroyable d'amis.

Une thèse est aussi un travail d'équipe. En ce qui me concerne, je remercie mon équipe Mcube (M 3 ) pour toute leur aide durant ma thèse mais aussi pour cet incroyable ambiance de travail qu'on retrouve rarement de nos jours.

Ainsi je remercie Carlos, notre directeur d'équipe. Je tiens aussi à remercier Pierre pour sa disponibilité et sa bonne humeur et surtout à sa capacité à trouver une réponse à chacune de mes questions même les plus imprévisibles. Je remercie également Tatiana pour son aide et ses conseils scientifiques et sportifs notamment par rapport aux étirements. Comment ne pas remercier Thierry, qui comme Fabrice a été mon enseignant et que je considère comme un modèle. Je remercie également André pour ses remarques constructives et les blagues qu'il nous faisait durant les « pauses-café mathématiques ». Je remercie aussi Jean Wolff pour les quelques discussions que nous avons pu avoir. Je remercie Marc pour son aide durant ma thèse que ce soit pour les problèmes matériels, pour nous avoir préparé régulièrement des petits gâteaux ou encore pour avoir géré les jeunes stagiaires avec moi. Enfin, je remercie Antonio pour ses conseils, son aide et le dynamisme qu'il apporte avec lui. Je remercie également Jean Farago et Hendrik Meyer de l'équipe TSP pour leur aide et conseils notamment concernants les mathématiques pour la physique. J'ai aussi eu l'occasion de participer à plusieurs évènements avec l'institut Charles Sadron tel que la fête de la science ou encore « Mature your PhD », j'aimerais donc remercier toutes les personnes avec qui j'aurais pu travailler sur ces projets, je pense notamment à François Schosseler, Olivier Felix, Wiebke Drenckhan, Marc Schmutz, Laure Biniek, Antonio Stocco et Amparo Ruiz Carretero avec qui nous avons travaillé sans relâche pour le projet du château d'eau mais aussi aux designers Raphaël et Juliette. Enfin je remercie notre directeur Christian Gauthier pour toutes ses recommandations et son aide pour ces évènements.

Je remercie Marc Schmutz de m'avoir fourni des images de vésicules de DPPC mais aussi pour sa curiosité qui nous a permis de visualiser une plume de volant de badminton au microscope électronique.

Je remercie notre personnel administratif qui est toujours disponible pour résoudre nos problèmes. Je pense notamment à Odile, Nadia et Katia qui m'ont souvent permis de résoudre mes problèmes mais aussi à Magali avec qui ma vie de doctorant a été rendue très agréable.

Je souhaiterais remercier notre « Lunch group » (Jean-Mulhouse, Vaibhav-G

J'aimerais aussi remercier Olivier Benzerara pour son aide précieuse concernant les problèmes de simulation et surtout pour son aide concernant le monde mystérieux de l'informatique. Parfois, on ne sait pas pourquoi nos scripts refusent de faire ce qu'on veut et c'est là qu'Olivier intervient en ajoutant la virgule ou les deux points magiques qui font que tout fonctionne.

Durant ma thèse j'ai eu la chance d'avoir différents compagnons de bureau avec qui j'ai partagé cette aventure. Tout d'abord, je remercie Tiago Espinosa de Oliveira pour tout le temps qu'il a pris pour m'aider sur GROMACS mais aussi pour nos différents voyages à travers « Google map » et évidemment nos petites pauses vidéo. « Muchas gracìas, Pepito. Vamos comer ».

Comment ne pas remercier Adrien pour toute son aide lors de ma fin de thèse. Il possède un sens scientifique très aiguisé avec lequel il m'a permis d'améliorer mon travail, chaque détail compte. Je le remercie également de m'avoir appis à utiliser certains logiciels ainsi pour le temps qu'il a pris pour m'aider pour améliorer mes graphiques.

Enfin, je souhaiterais remercier Vaibhav et Elise pour la bonne humeur qu'ils ont apportée dans notre bureau et pour nos nombreuses rigolades. Et enfin, J'aimerais [START_REF] Periole | Chapter 20 the martini coarse-grained force field in[END_REF]. The types Q a , Q 0 , N a and C 1 are identical as for the mapping of DSPC. . . . . . . . . . [START_REF] Falk | Interleaflet sliding in lipidic bilayers under shear flow: comparison of the gel and fluid phases using reversed non-equilibrium molecular dynamics simulations[END_REF] When opposing forces are exerted on each leaflet (case a) the resulting stationary state sees the two leaflets sliding at constant relative velocity, surrounded by a uniform solvent velocity gradient profile, as emphasized in the sub-picture (b) where the simulation box boundary has been purposely shifted to sit exactly at the mid-plane of the bilayer. Case (a-b) are subsequently referred as a Couette situation. When a uniform force is exerted on both leaflets and an opposing force on the solvent beads (case c), a symmetric parabolic velocity profile builds up in the solvent, assuming sticking boundary conditions at the interface with the bilayer (d). Case (c-d) will be subsequently referred as Poiseuille situation. In all cases, the total momentum of the system is constant and vanishes. . . . . . . . . . . . . . . . . . . . . . . 3.5 Constant force pulling experiments in the fluid state: (A) single leaflet COM displacement X(t) (a) , starting from the simulation box center (ca 6.6 nm) and (B) averaged displacement hXi(t) ' 1/50 P 50 a=1 X (a) (t) vs time. A bootstrap procedure (C) estimates the dispersion s b (X u (t)) caused by the finiteness of the sample {a}. Vertical bars represent the confidence interval of 10 selected points from the second half of the trajectory (5000 < t < 10000 ps), taken as twice the estimated bootstrap standard deviation. The vertical bars are used to provide a confidence interval for the drift velocity (slope of the averaged displacement curve). /TrajectoryAndBootstrapFluidPull . . . . . . . . . . 3.6 Normalized averaged displacements (upper leaflet) X u (t)/F for a set of increasing pulling forces 10, . . . , 2000 kJ.mol 1 .nm 1 (equivalently stresses t = 0.1 . . . 190 bars). The displacement for F = 2000 lies clearly beyond the linear regime and the force F = 10 is competing with thermal agitation. Inset: close-up focus on the normalized displacement for stresses larger than 0.955 bars. . . . . . . . . . . . . . . . 3.7 Averaged drift velocities vs applied forces (lower horizontal scale bar)

List of Figures

Examples of lipid mapping to CG Martini beads

or stresses (higher horizontal scale bar). A shear thinning deviation is seen at t > 100 bars. Inset: focus on the linear regime region. . . . . . . xi 3.8 Ten different bootstrap realizations of the displacement hX u i(t) for an initial kick impulsion V 0 = 0.09 nm.ps 1 (upper leaflet). Inset: closeup look of the first peak. . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.9 Normalized averaged displacements hX u i(t)/V 0 for a set of increasing velocities 0.01 and 0.03, . . . , 0.1 nm.ps 1 . . . . . . . . . . . . . . . . . 3.10 Normalized averaged velocities hV u i(t) = hdX u /d ti(t)/V 0 . The velocity starts at an initial value of 1, decreases fast to 0 (coinciding with the sharp peak in the displacement curve) reaches a negative minimum and finally slowly regresses to 0 from below, coinciding with the slowly decreasing approach of the displacement plateau value. Inset: close-up look at the minimum. . . . . . . . . . . . . . . . . . . . . . . . 3.11 Ten different bootstrap realizations of the displacement hX u i(t) for an initial kick impulsion V 0 = 0.09 nm.ps 1 in the gel phase (upper leaflet). Inset: close-up look at the first peak. . . . . . . . . . . . . . . . 3.12 Normalized averaged displacements hX u i(t)/V 0 for a set of increasing impulsions 0.01 and 0.03, . . . , 0.1 nm.ps 1 . The plateau value is clearly increasing with the initial applied velocity, and the normalized displacements do not appear to collapse onto a master curve, pointing to an absence of linear response. . . . . . . . . . . . . . . . . . . . . . . 3.13 Normalized averaged displacements hX u i(t)/V 0 for a set of increasing impulsions 0.01 and 0.03, . . . , 0.1 and sample size 1  a < 200. . . . 3.14 Normalized averaged velocity hdX h /d ti(t)/V 0 . Starting from 1, the normalized velocity crosses 0 fast to reach a minimum, and then relaxes slowly to 0. Curves do not collapse to a master curve. The shape of the normalized relaxation curve is strinkingly different from the fluid case state. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.15 Velocity V u -applied shear stress t characteristics in the gel and fluid state. If the fluid state display a linear characteristics up to t=100 bars (linear regime) the characteristics seems nowhere linear in the gel state. 3.16 Velocity V u -applied shear stress t characteristics in the gel and fluid states, using a logarithmic representation. The cross-over from linear to shear-thinning regime is visible in the fluid state, while the gel regime seems consistent with a power-low relation with exponent V u ⇠ t 1.44 as a function of the applied stress t. . . . . . . . . . . . . . . 3.17 Apparent friction coefficient b + h/L w from impulsion simulations as a function of the initial induced velocity in the fluid phase. A vertical error bar is inferred from a bootstrap procedure with 10 synthetic averaged displacement curves. The error bar is very large for the initial velocities smaller than 0.7 nm.ps 1 , and becomes small for initial velocities larger than 0.2 nm.ps 1 , for which a departure from linear response can bee seen in Fig. ferent parts of the system (Top/Bot leaflets). Also shown are the top Water, bottom Water and the motionless substrate coating as reference. On this graph, the system center of mass is free to move (no center of mass velocity removal). Each leaflet was pulled in opposite directions along the x direction with a force equal to 250 kJ.mol 1 .nm 1 . 4.8 Normalized averaged displacement of the center of mass of the different parts of the system (top and bottom leaflets, free and confined water above and frozen substrate). The removal of center of mass is activated. Each leaflet was pulled in the opposite direction along the x axis with a force equal to 250 kJ.mol 1 .nm 1 . . . . . . . . . . . . . . 4.9 Normalized averaged displacement of the center of mass of the different system subsets (top and bottom leaflets, top and bottom water above and motionless substrate. The global center of mass is free.

The entire bilayer is pulled in one direction and the substrate in the opposite direction, parallel to the x axis with a force of 250 kJ.mol 1 .nm 1 . Again the substrate is steady. . . . . . . . . . . . . . . . . . . . . . . . . 4.10 Normalized averaged displacement of the center of mass of the different subsystems (top and bottom leaflets, top and bottom water, and motionless substrate. The center of mass velocity removal is activated. The entire bilayer is pulled in one direction and the substrate in the opposite direction parallel to the x axis with a force equal to 250 kJ.mol 1 .nm 1 . The substrate stays at rest. . R HEOLOGY is the branch of physics that deals with the flow and deformations of matter under the effect of forces. Its etymology comes from Greek rhéo ( < r 'ew) which means "flow" and logia (log' ia) "study of" . Let's imagine now two plates separated by a material; what is happening if shear is applied, in other words, what is the behavior of the material if the plates are pulled in the opposite direction? The answer will depend not only on the physical state of the material but also on the applied force.

All these universal phenomena occur all around us, indeed, it is also a wellknown fact that nature produces water-based lubricant systems which are outclassing most of the man-made devices [START_REF] Urbakh | The nonlinear nature of friction[END_REF]. Biological contacts (knees, hips...) undergoes severe conditions (i.e low speed, high load) and could be subject to injury or breakage, that is why the most important thing results in reducing frictions and wears. This is true not only for biological systems but also for the majority of technological devices which often works in extreme conditions (MPa-GPa range) [START_REF] Yamada | Friction and Adhesion Hysteresis of Fluorocarbon Surfactant Monolayer-Coated Surfaces Measured with the Surface Forces Apparatus[END_REF]. Animals joints (Fig. 1) are formed by cartilage surfaces separated by a pocket of synovial fluid composed by phospholipids (multi-layers, vesicles), hyaluronic acid and lubricin. This complex structure plays a central role in bio-lubrication. It presents a very low macroscopic coefficient of friction (µ = F T /F N = 0.005 0.02), F T and F N being the transverse and normal components of the force applied to the surface [START_REF] Schroeder | Surface Active Phospholipids as Cartilage Lubricants en[END_REF]. Chapter 0. General introduction

The understanding of such biophysical systems implicates both fundamental and applied challenges, for instance in regards to osteoarthritis and the manufacture of joint prostheses [START_REF] Klein | Repair or Replacement-A Joint Perspective[END_REF]. The important point is undoubtedly the comprehension from a microscopic point of view the role of phospholipids in the bio-lubrication. It is assuredly of great biochemical, biophysical, and technical interest. Moreover, common tribological and lubrication techniques used for macroscopic objects can be ineffective at nanometric length scale [START_REF] Urbakh | The nonlinear nature of friction[END_REF].

This work aims towards a better understanding of the role of phospholipids bilayers in bio-lubrication by developing numerical models at molecular length scale with a coarse-grained model [START_REF] Marrink | The MARTINI Force Field: Coarse Grained Model for Biomolecular Simulations[END_REF] to deeply understand the mechanical properties.

Foremost, chapter 1 State-of-the-art introduces the reader to the world of lipids by giving an overview of what lipids are and why they are so important in living system. An introduction to molecular dynamics is given. Some important results about the domain are presented.

Then, chapter 2 concerns the methodology and describes how the systems are prepared and analyzed. Chapter 3 describes the different results obtained for an ordinary free lipid bilayer subject to different types of solicitations. Then, chapter 4 addresses the case of a supported membrane. In this situation, the simulation approaches of chapter 3 are adapted to account for the presence of a substrate. Chapter 5 is a brief account of a method that we developed in order to equilibrate very small vesicles. This method is certainly promising as far as equilibrating and relaxing such systems is concerned. Finally chapter 6 concludes our work and unveils general perspectives.

Chapter 1

State-of-the-art

"Do not grow old, no matter how long you live. Never cease to stand like curious children before the Great Mystery into which we were born."

Albert Einstein (1879-1955)

Lipids

Introduction

L IPIDS refer to a large group of molecules that are soluble in nonpolar solvents.

Biochemists define lipids as hydrophobic or amphiphilic molecules. Oils and fats belong to this large family, differing in their physical states. Oils are usually found in a liquid state, as for instance olive oil, whereas fats like butter and waxes are usually solid. However, it is possible to melt butter or freeze olive oil until it becomes first highly viscous and finally solid. The state that fats adopt depends mainly on the nature of the long aliphatic chains which form an essential part of many of these compounds [START_REF] Mouritsen | Lipids in a Membrane Biophysics Perspective 2[END_REF]. Amphiphilic lipid molecules comprise one hydrophilic group associated to one or more hydrophobic group. In many cases the hydrophobic group consists in fatty acid chains, that is to say a skeleton of carbon atoms linked to each other by single bonds (saturated) or by one or several cis-double bonds (unsaturated), with a carboxylic acid group that can participate in esterification reactions with a number of hydrophilic groups. Such fatty acids are depicted in fig. 1.1.

Lipids are divided into eight categories depending on their chemical structure (fatty acids, glycerolipids, glycerophospholipids, sphingolipids, saccharolipids, polyketides, sterol lipids and prenol lipids (fig. 1.2) [START_REF] Fahy | Lipid classification, structures and tools[END_REF]. This work focuses mainly only on phospholipids to understand their role in bio-lubrication.

A brief history of fat

Since the beginning, the nature of life has fascinated humanity. However, it is not until recently that the importance of lipids for life has been fully recognized.

The first isolated lipid (cholesterol) was discovered in bile and gallstones around 1789 by François Poulletier de La Salle [START_REF] Olson | Discovery of the Lipoproteins, Their Role in Fat Transport and Their Significance as Risk Factors[END_REF]. Later on, Michel-Eugène Chevreul exposed the concept of fatty acids and identified the saponification process [START_REF] Chevreul | Recherches chimiques sur les corps gras d'origine animale fr[END_REF]. It is only after the contribution of William Prout in 1827 that fat started to be considered as important to the same extent as sugars and proteins [START_REF] Prout | On the Ultimate Composition of Simple Alimentary Substances; with Some Preliminary Remarks on the Analysis of Organized Bodies in General eng[END_REF]. Finally came the discovery of lecithin in 1847 by Theodore Nicolas Gobley. He was the first to identify the phospholipids class, which is central to this study.

More recently (1917), Irving Langmuir developed an apparatus where lipids can be spread onto an air-water interface, to form a monolayer. He suggested that in such monolayers, hydrophobic tails were directed towards the air region and the hydrophilic head groups were embedded into the water region [START_REF] Heimburg | Digital encyclopedia of applied physics[END_REF]. Langmuir troughs are nowadays routinely used to study air-water lipid monolayers, and to transfer monolayer films onto solid surfaces. Evert Gorter and François Grendel demonstrated in 1925 that the total area of red blood cells was half of the area of the monolayer created with the lipids extracted from the same cells. This pioneering observation was key to solving the puzzle of phospholipid assembly, and they suggested immediately that lipids were naturally forming bilayer like structures [START_REF] Gorter | On Bimolecular Layers of Lipoids on the Chromocytes of the Blood[END_REF]. At this time, scientist started to unravel the structure and the complexity of the cell biomembranes. In 1972, Singer and Nicolson proposed a new picture of cell membranes which became known as the Fluid Mosaic model [START_REF] Singer | The Fluid Mosaic Model of the Structure of Cell Membranes[END_REF]. These authors considers that a membrane is an organized fluid in which some proteins are inserted, either into one single monolayer, as transmembrane components or as peripheral bounded units (fig. 1.3). FIGURE 1.3 -A popular cell membrane representation, according to the fluid mosaic picture [START_REF] Singer | The Fluid Mosaic Model of the Structure of Cell Membranes[END_REF].

Yet in the fluid mosaic approach, lipids were only viewed as a sophisticated passive solvent, playing essentially a structural and compartimentation role. It was suggested later than lipid membranes were capable of acting on their proteic components, by means of hydrophilic-hydrophobic interactions [START_REF] Mouritsen | Mattress model of lipid-protein interactions in membranes[END_REF]. The membrane was found to be asymmetric in composition and laterally heterogeneous [START_REF] Devaux | Transmembrane Asymmetry and Lateral Domains in Biological Membranes[END_REF][START_REF] Simons | Functional rafts in cell membranes[END_REF]. Some lipids are today considered as being actively involved in some signaling chain processes (functional lipids) while others were found to be necessary to facilitate lipid membrane remodeling (vesicle scission/fusion). Finally, lipids are in charge of controlling the biological membrane fluidity, as it was said above [START_REF] Mouritsen | Lipids in a Membrane Biophysics Perspective 2[END_REF].

Lipids are everywhere and essential to your body and constitute a major part of the brain. Our lipids are metabolized in the liver, which distribute fatty acids and address them to the other organs. Fatty acids are used as energy storage (insoluble triglyceride fat) or in a more noble way, incorporated in various proportion into phospholipid compounds. All cells are capable of lipid metabolism (nuclear membrane and endoplasmic reticulum) and involved in active lipid traffic processes between organelles [START_REF] Van Meer | Membrane lipids: where they are and how they behave[END_REF]. Lipids are the building scaffold of living matter and allow the cells to create barriers between interior and exterior, to compartment the activity of within cells and organelles (see fig. 1.4). Many lipids are hormones precursors (steroids, vitamins) and act as bio-receptors (gangliosides) [START_REF] Mouritsen | Lipids in a Membrane Biophysics Perspective 2[END_REF]. One can perfectly claim that "Fat is life" . In relation, lipids have found numerous applications in the fields of cosmetics, food industry and medicine [START_REF] Mashaghi | Lipid Nanotechnology[END_REF]. More specifically related to our research domain are the animal joints (fig. 1). These fascinating mechanisms are composed by two cartilage surfaces separated by a pocket of synovial fluid which together realize a really efficient bio-lubricant device with very low coefficient of friction. This viscous fluid contains the following components: hyaluronic acid, lubricin, phospholipids (phosphatidylcholines (41%), phosphatidylethanolamine (27%) and sphingomyelin (32%)) [START_REF] Sarma | Phospholipid composition of articular cartilage boundary lubricant[END_REF]. In particular, phospholipid compounds are arranged as three to seven stacked lipids bilayers, trapping physiological serum in between [START_REF] Schwarz | Surface-active phospholipid as the lubricating component of lubricin[END_REF] and thus yielding a remarkable bio-lubricant. Bio-lubrication is not restricted to joints, but everywhere tissues are sliding relative to each other: cardiac muscle beats, eyes motion. . .

Phospholipids

Phospholipids can be divided into two main classes (fig. 1.5) depending on their backbones [START_REF] Bailey | Food lipids: chemistry, nutrition and biotechnology (food science and technology series/88) edited by casimir c. akoh (university of georgia) and david b. min (university of ohio). marcel dekker: new york and basel[END_REF].

• Glycerophospholipids: based on a glycerol group • Sphingolipids: built on a sphingosine backbone Glycerophospholipids are amphiphilic molecules composed of one bulky zwitterionic hydrophilic head group (e.g. phosphocholine, or phosphatydil-choline aka PC) and two hydrophobic aliphatic hydrocarbon chains with a chain length comprising between 12 and 20 carbons, esterified around a glycerol group (fig. 1.6). The phosphocholine group associates a negative phosphate charge (PO 4 ) and a positive ternary amine N(CH 3 ) + , liked by two carbons. When the both chains consists of stearic acids, with 18 carbons long linear unsaturated carbon chains, the phospholipid is called di-stearoyl-phosphatydilcholine, abbreviated as DSPC. Another common lipid is the di-palmitoyl-phosphatydilcholine, abridged as DPPC, consisting of two palmitic fatty acids with 16 carbons along their unsaturated chains. DPPC bears its name from palm oil, an important agricultural resource whose environmental cost is a matter of strong controversy.

Glycerol Sphingosine

Phosphate FIGURE 1.5 -Molecular structure of two important phospholipid classes: glycerophospholipids, built around a glycerol with two variable fatty acids, and sphingolipids (built around a sphingosine with a variable fatty acid chain). The phosphate group is common to both.

Note that this lipid nomenclature is common but not unique. IUPAC lipid names and other denominations are detailed in biochemistry handbooks, such as the Handbook of Lipid Bilayers [START_REF] Marsh | Handbook of Lipid Bilayers en[END_REF]. Unlike fats, phospholipid are not completely insoluble in water, but do not disperse well either. These molecules self-assemble due to weak interactions as flat bilayers, made of two leaflets with tail to tail opposing lipid molecules. The driving force of lipid self-assembly is, as in the case of protein folding, the hydrophobic effect [START_REF] Israelachvili | Intermolecular and surface forces[END_REF]. As hydrocarbon chain disrupt the preexisting network of hydrogen bonds characteristic of liquid water, an unfavorable free-energy is associated to the hydrocarbon-water interface. The self-assembly mechanism consists in bunching together the aliphatic part of the molecules in order to to minimize the free-energy of water solvation. An intuitive and correct picture amounts to saying that fatty chains tries to shield themselves from water.

The phenomenon is common to all surfactant molecules (soaps, detergents). Other possible configurations are then observed such as micelles or worm-like micelles (fig. 1.6) depending on geometrical packing parameter considerations [START_REF] Israelachvili | Intermolecular and surface forces[END_REF]. PC lipids almost always arrange themselves as bilayers or membranes. The bilayer thickness is of the order of 4-5 nm, while its lateral dimensions reaches easily tens of micrometers. This bilayer thickness was first determined experimentally by X-rays crystallography by V. Luzzati and coworkers at the Institut Charles Sadron [START_REF] Luzzati | The Structure of the Liquid-Crystalline Phases of Lipid-Water Systems[END_REF].

There is in principle no thermodynamic limit to the bilayer lateral dimensions. In most cases, both leaflets have a symmetric composition and the bilayer do not show any spontaneous curvature. Nevertheless, the presence of a free bilayer boundary is unfavorable beyond a given critical size, and the bilayer closes itself spontaneously to form a vesicle. The size of vesicles is not prescribed by the thermodynamics or the lipid composition, but depends on the detail of the lipid preparation protocol and the age of the solution. Small vesicles, with radius of the order of 50 to 100 nm are called liposomes. They are easy to prepare in large quantities and have found many application in pharmaceuticals and cosmetics. Liposomes are also an essential model systems for lipid studies. When comprising a unique bilayer, their are called Small Unilamellar Vesicles (SUV).

Vesicles of larger sizes are also very interesting. Giant Unilamellar Vesicles (GUV) are large enough to be observed by optical microscopy and manipulated using tweezers or micropipettes [START_REF] Dimova | The giant vesicle book 1st[END_REF]. An alternative to lipid vesicles are the Supported Lipid Bilayers (SLB). They consists in a single bilayer in close contact with a flat hydrophilic surface. The cost of the free boundary in there compensated by the favorable interaction with the surface. Langmuir film transfer or liposomes deposition can be used to form SLB. One can compare these membranes to an ocean of lipid with molecules free to move relative to each other in the plane. Biological membrane fluidity depends on three main factors:

• Temperature: it is a measure of the average kinetic energy of the systems.

As most condensed bodies, membrane dilate under the effect of temperature: lipids are tightly packed at low temperature whereas they expand at high temperature. As a matter of fact, membrane fluidity is believed to increase with temperature according to the Arrhenius law, i.e. thermal activation of elementary lipid moves is easier and faster at high temperatures.

• Presence of saturation in the fatty acid chain: cis-unsaturations are associated to kicks in the chain skeletons and brings disorder into the lipid packing, thus decreasing the mutual cohesion of the assembly. The insaturated to saturated ratio is one key parameter that membranes controls in order to maintain the membrane fluidity to the desired level. It is also distinctive of the various vegetal oils properties.

• Cholesterol: This molecule interact very specifically with saturated PC molecules, to create a very ordered with still fluid phase called Liquid ordered state. It also change the mechanical and permeation properties of the membrane.

Pure lipid bilayers can be found in different physical states, and are subject to a sequence of phase transitions as temperature increases. The most important transition is the main, or melting transition, associated to a sharp calorimetric signature. Below the melting transition T m , the bilayer is in a highly viscous state called "gel" phase. Lipid chains are tightly packed, with consecutive CC bonds in all trans configurations. Strictly speaking PC molecules are observed not in one but two different gel phases. At a temperature lower than the pretransition temperature, the bilayer is laterally homogeneous with tilted packed chains, the L 0 b state. Above the pretransition, an instability takes place which modulates laterally the membrane thickness, in the so-called ripple phase P b . Above the melting transition temperature T m the bilayer is in a fluid phase. Lipid chains are melted, with a finite fraction of gauche isomerization kinks along the chains. The bilayer is thinner and more extended in the lateral dimension. The melting temperature depends on the lipid headgroup nature, increases with the number of carbons in the chains and decreases when cis insaturations are present [START_REF] Heimburg | Thermal biophysics of membranes[END_REF].

The lipid phase influences strongly the formation of membranar objects. Only a fluid phase can grow as giant vesicles while a gel state facilitates the transfer of Langmuir films onto the glass, silicon or mica surfaces.

The particular case of DSPC

This work is essentially based on the 1,2-distearoyl-sn-glycero-3-phosphocholine (fig. 1.7) which is a saturated phospholipid. It is composed of two (D) stearoyl (S) chains groups that is to say, two apolar chains of 18 carbons. The polar head group is composed by a phosphatidylcholine (PC) which is an association of a N,N,N-trimethylethanolammonium cation and a phosphate.

Common phospholipids (e.g. dipalmitoylphosphatidylcholine DPPC, distearoylphosphatidylcholine DSPC) do not interdigitate under standard conditions, and the leaflets are relatively weakly bound together [START_REF] Marsh | Structural and thermodynamic determinants of chain-melting transition temperatures for phospholipid and glycolipids membranes[END_REF]. Bilayer fluidity depends significantly on temperature [START_REF] Cevc | Phospholipid bilayers. physical principles and models[END_REF]. Pure DSPC bilayers melt at T m = 55 C while DPPC bilayers melt at T m = 41 C [START_REF] Marsh | Handbook of Lipid Bilayers en[END_REF]. Above melting, lipid tails are isomerically disordered, weakly cohesive, conferring fluidity to the bilayer with Arrhenius dependence in temperature. Below melting, lipid tails adopt all-trans conformations, are subject to stronger cohesion, displaying solid type dynamics at short time scales, while remaining a viscous fluid on longer scales. Cryo transmission electron microscopy (Cryo-TEM) pictures of DPPC liposomes quenched from the gel and fluid phases are shown in fig 1.9. Vesicles in the gel state display characteristic polyhedral shapes revealing a rigid structure. Vesicles in the fluid state are smooth and spherical and the double leaflet structure is clearly visible. FIGURE 1.9 -Cryo-TEM pictures of DPPC vesicles in fluid phase (at 46 C) and gel phase (at room temperature RT). The picture was provided through the courtesy of Dr. Marc Schmutz. The preparation protocol is exposed in ref [START_REF] Krafft | Fluorinated Vesicles Allow Intrabilayer Polymerization of a Hydrophobic Monomer, Yielding Polymerized Microcapsules[END_REF] .

An important geometric characteristics of membranes is the area per lipid, APL (eq.1.1) determined as the total area divided by the number of lipids present in a leaflet. This parameter is not easy to determine experimentally but is well known, around 63.8 Ų for DSPC at 60°C [START_REF] Kučerka | Fluid phase lipid areas and bilayer thicknesses of commonly used phosphatidylcholines as a function of temperature[END_REF].

APL =

Total area Number of lipids (1.1)

Many reasons led us us to choosing DSPC for studying the role of phospholipids in the bio-lubrication:

• Phosphatidylcholine headgroups are known to form a major component in lubricant systems (41% for joints).

• The presence of experimental set-up in our group dedicated to the measure of the friction and mobility of supported lipid bilayers (TribologyFluorescence recovery after patterned photobleaching FRAPP ) [START_REF] Fu | A new tribological experimental setup to study confined and sheared monolayers[END_REF].

• The preparation of mono or trilayers of DSPC at the solid-air interface by Langmuir trough deposition is facilitated by the gel state at room temperature.

Molecular Dynamics 1.2.1 Introduction

Molecular dynamics (MD) is a computational method that focusses on the time dependent properties of molecular systems, such as proteins, lipids, nucleic acids. . . It allows to determine important molecular properties based on statistical physics principles, in realistic situations where analytical approaches are completely impossible. It informs us on the structure, stability and dynamics of a given system. It can also be used to determine certain thermodynamical properties. In the present approach, MD will be our main tool to investigate the behavior of phospholipid bilayers subject to shear stresses. Equilibrium statistical physics was formulated in terms of ensemble averages by Gibbs. Only few cases are analytically tractable or can be treated by means of an accurate perturbative approach. In the majority of the cases of interest, the desired ensemble average can only be achieved by sampling the phase space (configuration and velocity space) with a numerical strategy.

The principle of molecular dynamics is to calculate realistic molecular trajectories in order to sample as accurately as possible the desired statistical ensemble. The main idea is to solve numerically the Newton's equations of motion (eq:1.2) of a mechanical system represented by N beads (or atoms) subject to interaction forces, on a given time interval t. In other words, it amounts to solving numerically a large number N of coupled second order differential equations in three dimensions.

X i ! F i = m i ⇥ ! a i , i = 1...N (1.2)
• ! F i : external force acting on the point-like center of force i, usually due to neighboring bodies.

• m i : mass associated to the center of force.

• ! a : acceleration of the center of force, or center of mass of the bead/atom considered

The integration of the equations of motion results in a time discrete trajectory which is recorded into a file. The trajectory contains enough information to compute time averages of physical observables, i.e. functions of the phase space variables that have a special and physically relevant meaning. Meanwhile, parameters of thermodynamic relevance are recorded in the course of the simulation, such as pressure, temperature, internal energy, etc. . . It is expected on general grounds that such observable time averages converge to the equivalent ensemble average, in a process known as ergodicity. Equilibrium is said to be reached when time averages are stabilized and do not longer evolve when increasing the simulation time. A compromise must be found in terms of trajectory size, system size and simulation time to obtain the most accurate result with reasonable computational effort.

Molecular dynamics simulations share some similarity with real experiments [START_REF] Frenkel | Understanding Molecular Simulation: From Algorithms to Applications en[END_REF]. It deals with "samples" (the prepared simulation box) and "instruments" (the observables) to determine say, the temperature. Data acquisition takes place during finite time (the simulated time), and the quality of the measure depends on the sampling time. The quality of the recorded data increases with the acquisition time and the number of repetitions of the experiments. Putting a thermometer inside a water bath for 1 second is inefficient because the recorded signal will be noisy whereas if we let it for 5 minutes a more realistic and representative average value of the temperature will be obtained. Obviously, statistics and quality of sampling plays an essential role in MD. On the other hand, there are also significant differences. The time and length scales are typically so small that direct comparison with experiments is not always possible. Systems with slow dynamics or large-scale organisation are in practice either difficult or impossible to directly simulate. Finally, the advantage of MD is that all the parameters can be controlled and tuned, and almost all the parameters of physical interest can be computed. By contrast, real experiments are always subject to uncertainties and unknown factors, for instance in terms of composition or external parameter control. MD simulation cannot reproduce exactly experiments, but may give insight onto the underlying physical mechanisms [START_REF] Hamming | Numerical methods for scientists and engineers 2nd[END_REF].

A brief history of computer simulation

World War II was undoubtedly a powerful motor for the development of computing machines for military purposes (atomic bomb, code breaking etc...). From the 1950s, the computational methods started to be applied to civil academic research [START_REF] Israelachvili | Intermolecular and surface forces[END_REF]. In 1953, Metropolis et al. invented the Monte Carlo method to simulate the behavior of 2D rigid spheres, using the famous machine "MANIAC" at Los Alamos [START_REF] Metropolis | Equation of State Calculations by Fast Computing Machines[END_REF]. This method is based on repeated random moves for sampling the configuration space (coordinates) but gives no information about the dynamics. A few years after, Alder and Wainwright returned to a system of 3 dimensional hard spheres but applying molecular dynamics [START_REF] Alder | Phase Transition for a Hard Sphere System[END_REF]. In 1964, Rahman performed the first simulation using a realistic potential for liquid argon [START_REF] Rahman | Correlations in the Motion of Atoms in Liquid Argon[END_REF]. While in hard sphere systems the dynamics is dominated by collisions, with particles moving along straight lines in between, the motion of argon atoms require the evaluation of the forces at every step of the computation, as in the modern MD studies of realistic molecular systems. The complexity and realism of the MD studies have only increased since this time.

In 1977, McCammon carried out the first protein simulation with the simulation of the bovine pancreatic trypsin inhibitor (BPTI) [START_REF] Mccammon | Dynamics of folded proteins[END_REF]. MD started to become more and more popular in the studies of solvated proteins, protein-DNA, or lipid membranes. Indeed, in 1982, P.van der Ploeg and H. J. C. Berendsen carried out the first MD simulation of an amphiphilic bilayer membrane (decanoate molecules) [START_REF] Ploeg | Molecular dynamics simulation of a bilayer membrane[END_REF]. In 1993, Heller et al. observed the phospholipid melting transition [START_REF] Heller | Molecular dynamics simulation of a bilayer of 200 lipids in the gel and in the liquid-crystal phases[END_REF].

MD simulation is nowadays an indispensable tool to investigate the thermodynamical, dynamical and structural properties of complex molecular systems. Its success is due to the availability of powerful simulation packages that can exploit the full computational power of parallel CPU (computation processor unit) and GPU (graphic processor unit) clusters, and also to the constantly improving quality of the available force fields (i.e. the mathematical parameterization of the intermolecular interactions, cf below). Thanks to these tools, it is now easier to sample the equilibrium properties of a system with a MD code than to write a dedicated Monte-Carlo software, which should be in principle more efficient. An example of use of MD is the screening of drug candidates in pharmaceutical science.

Principle of Molecular dynamics 1.2.3.1 Statistical mechanics

Molecular dynamics simulations produce a large amount of data at the molecular scale, such as atoms instantaneous positions and velocities. This raw data are unintelligible, and a conversion of the microscopic level information to the macroscopic level observable (e.g. pressure, temperature, energy . . . ) is necessary [START_REF] Ngo | Physique statistique -Introduction DUNOD[END_REF]. Statistical Physics is the branch of physics that offers these interpretations tools. It aims at explaining the system macroscopic properties starting from a set of molecular interaction rules, and provides a link between the microscopic and macroscopic worlds. Statistical physics applies as soon as the number N of interacting centers is large enough to provide averages representative of the whole physical system.

Statistical physics relates the thermodynamic state of system with N components to an ensemble a microscopic configurations. There are different possible ensembles of microscopic states, each with its own characteristics:

• Microcanonical ensemble (NVE): The microcanonical ensemble is characterized by a fixed number of particles N contained in a fixed volume V and where the total energy E is conserved. It is often described as an isolated system. This ensemble is the most simple to conceive but is not a very realistic model, and is not easy to implement, due to the difficulty of discrete algorithms to conserve the energy.

• Canonical ensemble (NVT): The canonical ensemble is a closed system, like the microcanonical ensemble, which can exchange heat with an external reservoir at constant temperature called thermostat. The amount of particles N, volume V and temperature T are therefore fixed. As a consequence the total energy of the system fluctuates. This implementation of the ensemble with molecular dynamics is relatively convenient.

• Isobaric-Isothermal ensemble (NPT): The Isobaric-Isothermal ensemble is characterized by a constant number of N particles, a fixed temperature T and pressure P. The pressure is maintained by a barostat, the temperature by a thermostat. In this situation, the system is closed, with varying volume V and total energy. Simulation of lipid membranes require special versions of the NPT ensemble.

• Grand Canonical ensemble (µVT): The grand Canonical ensemble is characterized by a constant Volume V, temperature T but a variable number of particles N. It corresponds to an open system. Though theoretically very interesting, its implementation in the case of complex molecular systems is difficult.

For each ensemble, a long MD trajectory samples a representative set of microscopic configurations. Time averages of observable built from this long trajectory should converge to the corresponding ensemble average, according to the ergodic principle. This is only true if the system dynamics is not too slow and if the sampling is sufficiently large.

Classical mechanics in the NVE ensemble

Molecular dynamics simulation integrates the Newton's second law (eq. 1.2) for each bead, or center of force. Classical molecular dynamics do not consider the state nor the motion of the electrons around the atoms. It focuses on the motion of the atom nuclei, taken as representative of the whole atoms, and treated classically according to the Born-Oppenheimer approximation. Atoms are considered as point-like and act on each other by means of classical potentials. In this approximation, forces derive from the gradient of a potential energy function V (eq: 1.3) which supposedly accounts for all the electronic contributions.

! F i = ! r i V (1.3)
By combining and rewriting eq. 1.2 and eq. 1.3, we can express the relation between the derivative of the potential energy V and the positionsr i (eq: 1.4).

rV = m i d 2 ri dt 2 (1.4)
Forces are evaluated once per computation cycle and are therefore constant during every MD single step. Let us therefore suppose that the accelerationã is constant for some time, we have:

ã = dṽ dt (1.5)
and integrating 1.5:

ṽ = ãt + ṽ0 (1.6)
repeating the same steps, position is obtained

ṽ = dr dt (1.7) r = ṽt +r 0 (1.8)
Finally, the position x at time t as function of the force rV(r 0 ), initial position v 0 and initial velocity v 0 reads: r = rV

2m i t 2 + ṽ0 t +r 0 (1.9)
The above equation (eq: 1.9) expresses the position and velocity at time t in terms of the initial positions at time 0. The process can be iterated step by step, keeping the integration time t short enough for the forces to remain approximately constant. Integration in the presence of a stepwise constant force is a discrete deterministic dynamical system. However, each step introduces a small position error of the order of t 2 ||ã(t) ã(0)|| ⇠ t 3 . As a result, the total energy is not conserved as it should if Newton equations were exactly solved.

Prior to integrating the equation of motion, an initial configuration must be determined. Sometimes, experimental data suggests the nature of the starting point (X-ray crystallographic structure for instance). In general, one starts from either an ordered or random conformation, which is built with a sequence of gentle energy minimization and thermalization schemes, describe below. The initial distribution velocities (eq: 1.12) is exactly known, due to the peculiar form of the classical energy function describing the system

H = X i ✓ m i v 2 i 2 ◆ + V (1.10)
The velocity distribution follows a Maxwell distribution law (gaussian). In addition, it is often required to set the total momentum P of the system to zero. Equations (eq: 1.11) and (1.13) enable us to sample velocities from the exact velocity distribution at temperature T.

P = N X i=1 m i ṽi = 0 (1.11) r(ṽ i )dṽ i = ✓ m i 2pk b T ◆ 3/2 ⇥ exp  1 2
m i v i 2 k b T dṽ i (1.12)
The total kinetic energy defines the instantaneous kinetic temperature.

T = 1 3k b N N X i=1 m i v 2 i 2 (1.13)

Integration algorithms

The system composed of N particles generates 6N first order differential equations that must be numerically integrated. The discrete dynamics consists in finding the new coordinates at a time t + dt knowing the past positions and velocities. The knowledge of the position enables to compute the forces, according to a model detailed below (section 1.2.3.2).

In mathematics, many algorithms exist that can solve this system of differential equations. Approximation for integrating of a first order equation dy/dt = f (y, t) usually starts by writing the equivalent integral equation y(t) = y(0) + R t 0 dt 0 f (y(t 0 ), t 0 ) and repeatedly substituting the left hand side into the right hand side, reducing the error at each step. A verification of the resulting numerical scheme can be done by comparing the resulting solution to a Taylor expansion of the solution around the initial time 0 [START_REF] Frenkel | Understanding Molecular Simulation: From Algorithms to Applications en[END_REF].

In this work the Leap-Frog algorithm was used. This integration method is based in the Verlet positions algorithm developed in 1967 by Loup Verlet [START_REF] Verlet | Computer "Experiments" on Classical Fluids. I. Thermodynamical Properties of Lennard-Jones Molecules[END_REF]. This algorithm is simple but efficient, and presents the best performance for a minimum of integration errors. It preserves well the phase space volume and has good reversibility properties, as do the exact Newton's equations. The discrete trajectories are not exact, but they are qualitatively comparable to real trajectories. The leap-frog algorithm executes the following steps:

r(t + dt) =r(t) + ṽ ✓ t + 1 2 dt ◆ dt (1.14) ṽ(t + 1 2 dt) = ṽ ✓ t 1 2 dt ◆ +ã(t)dt (1.15)
In this picture, the velocities at time t + 1 2 dt (1.15) are calculated first, then the positions at the time t + dt are obtained (1.14). Therefore, the velocities and positions are not calculated for the same time, each one leaps over the other (fig. 1.11 hence giving the name of the method. The advantage is that velocities are explicitly known on a half integer space grid. When required, the velocities at t can be approximated by the following half sum relation [START_REF] Frenkel | Understanding Molecular Simulation: From Algorithms to Applications en[END_REF]:

ṽ(t) = 1 2  ṽ ✓ t 1 2 dt ◆ + ṽ ✓ t + 1 2 dt ◆ (1.16)

Controlling temperature

In order to work in NVT or NPT ensemble, the temperature needs to be constant over the simulation. Working in these ensembles cannot achieved by simply integrating Newton's equations. It is necessary to devise a system of differential equation whose solutions are known to sample the canonical or NPT ensemble instead of the constant energy microcanonical ensemble. In the mean time, the equations of motion must only differ slightly from the original Newton's equation, otherwise the dynamical properties of the system would not be realistic. The idea is to modify the original system of equations with a term mimicking a slow of exchange of energy with the external thermostat. The assumption is that details regarding the energy exchange matter little and are in fine all equivalent.

In practice, one needs to control the kinetic temperature of the system of interest. In this work, we mostly use the "V-rescale" thermostat, which derives from on the Nose-Hoover scheme [START_REF] Bussi | Canonical sampling through velocity rescaling[END_REF][START_REF] Van Der Spoel | GROMACS: fast, flexible, and free[END_REF]. This is achived by introducing one (or more) extra degree of freedom Q coupled to the kinetic energy term. This coupling pumps kinetic energy in and out of the system, and the whole scheme was demonstrated to sample ergodically the canonical ensemble.

A simpler and more robust scheme to enforce a prescribed value for the kinetic energy of the system is the Berendsen scheme. This consists in setting a target for the kinetic energy, and rescale the velocities according to the difference between the target and instantaneous kinetic energy term. The Berendsen scheme is very efficient to thermalize a system, i.e. to obtain structural and kinetic configurations typical from the desired temperature, which can serve later as a starting point for MD NVT simulations. The Berendsen scheme does not sample the canonical ensemble, as the fluctuations of the total energy are significantly smaller that thoses expected in the NVT ensemble, as predicted from statistical physics arguments. The Berendsen thermostat is therefore a favorite choice for preparing in a quick and robust way a system at a given prescribed temperature, but not for analysing it.

Controlling pressure

In the same manner as the thermostat, a barostat is an alteration of the dynamical equations of motion that reproduces the effect of external "pistons" acting on the system size, compressing and expanding it according to well defined rules. System size rescaling is done globally, by rescaling the particle positions with the simulation box sizes (L x , L y , L z ) and finding a dynamical law for the vector of box sizes. Changing the box dimensions amounts to perform uniform axis dilations, acting simultaneously on all particles. Changing the box size during the course of a MD simulation is a crucial feature in some cases, as for instance when a phase transition is accompanied by a change in density of the physical system. The isothermal-isobaric NPT ensemble adds an extra bias to the canonical probability of a system configuration proportional to exp( P 0 V/k B T) with P 0 the fixed pressure, V the volume of the system, k B the Boltzmann constant and T the temperature. Barostats work by confronting the instantaneous system pressure to a target value. In practice the instantaneous pressure of a system can be obtained from the forces, velocities and positions of the particles.

P = 2 3V (E k VIR) VIR = 1 2 X i<j rij • Fij (1.17)
with i, j running over the possible particle indices,r ij = ri rj , E k the kinetic energy and VIR is the internal virial for pairwise additive potentials. We used two barostats in this study:

• Berendsen barostat: the Berendsen barostat is a robust barostat usually used during the initial preparation stage, to reach quickly the correct system size [START_REF] Berendsen | Molecular dynamics with coupling to an external bath[END_REF]. Let us assume now that the instantaneous pressure P defined in eq: 1.17 is known. As it differs from the set pressure P 0 , a first order objective is assigned to the pressure variation

✓ dP dt ◆ obj = P 0 P t p (1.18)
It is not possible to act directly on P, but one can force a variation of volume that would match the corresponding pressure variation if the system was at equilibrium, with a compressibility coefficient c. Therefore, the relative system volume variation reads

dV V = c P 0 P t p (1.19)
where t p and c are a set of user-defined parameters. The Berendsen barostat does not sample the NPT ensemble but reaches the correct system size average value.

• Parrinello-Rhaman barostat: the Parrinello-Rahman is another type of barostat, which allows more fluctuations of the size of the simulation box, and is expected to sample the NPT isothermal-isobaric ensemble [START_REF] Parrinello | Crystal Structure and Pair Potentials: A Molecular-Dynamics Study[END_REF]. It also relies on rescaling the system volume, and introducing a set of supplementary dynamical variables conjugated to the system box sizes.

In lipid membrane simulations, it is necessary to allow the system box size to fluctuate independently in the xy and z directions. This is particularly clear in the case of the bilayer lipid melting, where the system simultaneously expands in the lateral xy and contracts in the z direction. These are called semi-isotropic barostats.

Periodicity

In order to avoid problems due to the edges of the simulation box, periodic boundary conditions (PBC) are applied. This method consist into surrounding the simulation box by copies of itself, called images (fig: 1.12).

The introduction of images of the simulated box introduces the risk that a particle interact with its own image, i.e. with itself. The risk is real in the case of electrostatic interactions which decreases slowly in space. Special summation techniques are required to handle this situation. Short range interactions are usually cut-off beyond a critical radius. The use of periodic boundary conditions requires that the interaction cut-off to be sufficiently short (see Fig. 1.12).

In the presence of PBC, the distance involved in the computation of the mutual interactions between particles is the given by the minimal image convention.

r ij ⌘ min k[j] ||r i rk || (1.20)
with the index k running over all the images of particle j.

Concept of Force Field

Molecular dynamics is based on the integration of the mutual atomic interactions.

These forces derive from a potential energy function V which one must now specify: this is the notion of Force Field (FF) or potential of interaction. A force field is an empirical description of the forces originating from the quantum dynamics of the electrons around their nuclei, and from electrostatic interactions between static charges. The force field shapes the molecules. The accuracy of the descriptions relies largely on it. This potential energy is split into two parts: • Intramolecular interactions: these interactions describes the chemical bonds that keep together the atoms in a molecule. They are referred as bonded interactions. The list of connections between atoms within a molecule is referred as the molecular topology.

• Intermolecular interactions: it concerns the interactions between atoms in different molecules, or atoms far apart in a molecule. They include Van der Waals and electrostatic contributions, and are referred as non-bonded interactions.

The non bonded interactions actually shape the atoms by assigning a Van der Waals radius to them, representing approximately the size of the cloud of electrons orbiting around the atom centers.

The total potential energy V(r) can be described by equation (1.21). One part correspond to the total energy from the bonded interactions V bonded and the other one to the non_bonded ones V non_bonded .

V tot = V bonded + V non_bonded (1.21)
The bonded interactions interactions are based on a fixed list of atoms (topology file). It describes the 2-bodies chemical bonding between pairs of adjacent atoms: stretching stiffness and equilibrium distance. It also includes 3-bodies bending potentials defining the rest angle between three consecutive atoms along a molecular skeleton, and the bending rigidity. Such bending potentials can describe for instance for the kink angle linked to a cis unsaturated double bond. In addition, the bonded potential includes 4-bodies torsion potentials that enforce the possible values of the dihedral angles defined between 4 consecutive atoms along a chain skeleton. Dihedral potentials describe for instance the izomerisation potential within aliphatic chains, which are involved in the lipid main transition.

V bonded = V stretching + V bending + V torsion (1.22)
• Bond stretching: it is related to the elongation or compression of a bond (fig:1.13).

Usually compared to a Hookean spring, it refers to a shifted harmonic potential and is expressed mathematically by the following relation:

V stretching (i, j) = 1 2 X ij K bond (r ij r 0 ) 2 (1.23)
where r 0 is the equilibrium distance between the atoms i and j, r ij is the actual distance or bond length and K bond the stiffness constant. Such bonds are unbreakable. Other potentials could be used if one considered the possibility of breaking bonds, for instance Morse potentials.

V bending (i, j, k) = 1 2 X ijk K angle (q ijk q 0 ) 2 (1.24)
where K angle is the spring constant of the angular deformation, q i,j,k the actual angle between atoms i, j, k and q 0 the equilibrium angle.

V dihedral (i, j, k, l) = M X n=0 K n 2 ⇥ 1 + cos(nf ijkl d) ⇤ (1.25) FIGURE 1.

-Bond rotation or torsion

where K n is the constant force of torsion, n an integer linked to the periodicity, f ijkl the angle of torsion (dihedral angle) and d the phase of the angle.

The non bonded interactions are classified into two types (eq: 1.26). The first concerns the short-ranged, or Van der Waals (VdW) interactions. The second concerns the so called long-ranged interactions (of electrostatic nature). These interactions act on all pairs of atoms, excluding the pairs that are too close neighbors along the molecular skeleton, in order not to interfere with the bonded 2, 3 and 4-bodies potentials. Non bonded interaction confers a "volume" to the atoms and are therefore crucial for predicting the actual density of a close-packed system.

In the practical implementation of the force fields, it is crucial to cut-off potentials at a finite distance, in order to keep the computational time manageable. Indeed, the computational time increase very fast with the cut-off distance, and becomes dominated by the evaluation of the forces which increase quadratically with the number of centers of forces, and in particular grows like r 3 cut , with r cut the interaction range.

V Non bonded = V VdW + V Electrostatics (1.26)
• Van der Waals interactions: they are considered as fast decreasing ( A/r 6 ) and weak interactions between atoms. The Van der Waals interactions are of three types [START_REF] Parsegian | Van der Waals Forces by V. Adrian Parsegian en[END_REF] with different coefficients A. It is a combination of dipoledipole interactions.

-Keesom interactions (permanent dipole / permanent dipole) represent the attractive interactions of two molecules presenting a permanent dipole moment such as the water molecule H 2 O. The interaction arise from the thermal fluctuations of orientations among the dipoles [START_REF] Parsegian | Van der Waals Forces by V. Adrian Parsegian en[END_REF]. As the dipoles are associated with the atom partial charges, the Keesom interaction is as a matter of fact accounted for by the long range electrostatic contribution.

C Keesom r 6 (1.27)

-Debye interactions (permanent dipole / induced dipole) represent the attractive interaction of a permanent dipole (A) on an apolar atom (B) which is polarized the permanent dipole (A). For instance, if an H 2 O molecule get closer to an H 2 , it will polarize the H 2 and create an induced dipole [START_REF] Parsegian | Van der Waals Forces by V. Adrian Parsegian en[END_REF]. If some force fields are able to account for the electronic polarization (polarizable force fields) this is not the case of the force fields used in the present approach.

C Debye r 6 (1.28)

-London interactions (induced dipole / induced dipole) take their origins from the energy reduction induced by the fluctuating dipoles associated with the electrons. This leads to a correlation of the directions of their dipoles [START_REF] Parsegian | Van der Waals Forces by V. Adrian Parsegian en[END_REF]. Such interactions are accounted for on a phenomenological basis by the following mathematical expression:

C London r 6 (1.29)
-Short range repulsive interaction. As the distance between two nuclei atoms decreases, the electronic cloud around them starts to overlap. The combination of electron confinement, electrostatic repulsion and Pauli exclusion principle results in a strongly short range repulsive interaction between nuclei.

A common mathematical expression that combines short range repulsion and Van der Waals attraction is the Lennard-Jones potential (fig. 1.16). This intermolecular potential is the sum of two contributions. The first part corresponds to a really short-ranged repulsive interaction (A/r 12 ) (the Pauli repulsion) and accounts for the energetically costly interpenetration of electronic clouds. The second part is London-Van der Waals interactions acting on longer distances.

The Lennard-Jones (1.30) is parameterized by two numbers, a depth e that measure how much two particles attract each other, and a distance s which represents the Van der Waals radius, e.g. the distance where the electronic clouds starts to overlap significantly. The Lennard-Jones potential has a single minimum equal to 1.12s, which represents the preferred distance between pairs of atoms. This distance is commensurate to the typical atomic packing distances in solid or fluid condensed phases. As a matter of fact, the Van der Waals attraction is responsible for liquid-vapor coexistence, and contributes significantly to the cohesion of condensed phases.

V(r) = 4e  ⇣ s r ⌘ 12 ⇣ s r ⌘ 6 (1.30) FIGURE 1.

-Representation of a Lennard-Jones potential

In practice, Lennard-Jones interactions are cut-off on a finite distance R cut , owing to the fast decrease of the r 6 Van der Waals contributions. This approximation is considered as minor in isotropic condensed systems.

• Electrostatic interactions: also known as Coulombic potential, they account for the long-ranged interactions between charged particles. Its mathematical expression is given by:

V elec = X i X j6 =i q i q j 4pe r e 0 r ij (1.31)
where e 0 is the permittivity of vacuum, e r the relative permittivity of the medium, related to the electronic polarisability (optical index), q the charge of the particle and r ij the distance between the two charged particles. Electrostatic interaction cannot be cut-off without altering deeply the properties of the system. To get around this obstacle, efficient computation methods have been introduced. The interactions are divided into a short-range and long-range part. The long range part is efficiently computed using Fourier transform techniques (Particle-Mesh-Ewald method). The short range part can be cut-off and treated in the same way as the Lennard-Jones interactions [49]

Different force field models

There are many different force fields, depending on the desired accuracy in the description of the molecular system. One can distinguish three different types of force fields.

• All atom force field: the atomistic simulations associates a center of force to all the atoms present in the system, including hydrogens. Sometimes, for the purpose of introducing polarization effects, more than one single center of force can be assigned to each given atom. Atomistic simulations are very realistic, but requires imprtant computational ressources. This is due to the large number of degrees of freedom, and to the importance of the electrostatic interactions associated with atomistic partial charges (i.e. the unequal electronic distribution across each chemical bond). The computation time is especially important when hydrogen atoms are present and their dynamics taken into account.

• United atoms description: the united atom simulations are quite similar to the all atoms ones excepted that the non-polar hydrogens are not explicitly taken into account. This implies a gain in the simulation time, but a loss of accuracy in the description and the kinetic properties.

• Coarse grained description: the coarse-grained simulations are another type of modelling where more than one atom are grouped together and assigned to a single center of force. In this case the simplification comes from the drastic reduction in the number of degrees of freedom.

In the Martini lipid force field (cf chapter 2.2), one single water bead represents four real water molecules. This brings about a huge gain in simulation time and allows to consider large systems simulated during long times. This gains come at the expense of a loss in accuracy and non quantitative dynamical properties. The Martini force field is used extensively in our approach (see section 2). A force field can be built from first principles (quantum chemistry calculations) or empirical considerations (structural and thermodynamic data). None of them is fully predictive and much work is currently devoted to improving them.

Existing methods of determination of friction

Coating a solid surface with a dense phospholipid monolayer modifies the sliding friction properties significantly. Experiments report a significant decrease in the sliding friction coefficient when both surfaces are covered with lipids in a dense, or gel conformation [START_REF] Briscoe | Boundary lubrication under water[END_REF]. This issue is relevant in the field of biolubrication, such as for instance the mechanism of articular joints. As a matter of fact, synovial fluid combines lipid and biopolymer molecules for optimal lubrication, the role of each component being still a topic of investigation.

It is difficult to relate the macroscopic friction between a pair of surfaces with microscopic mechanisms involved at the molecular scale [START_REF] Persson | Sliding friction, physical principles and applications[END_REF]. In the case of hydrated lipid bilayers, a lateral shear displacement involves the solvent viscosity, the sliding leaflet friction and possibly some sliding of the solvent on top of the hydrophilic bilayer surface. In the framework of linear response, sheared lipid bilayers display a viscous response, characterized by an interleaflet friction coefficient b, a Newtonian transverse viscosity h for the solvent, and a solvent-bilayer friction b 0 which quantifies the importance of the sliding of the fluid at the bilayer interface.

Interleaflet friction

When a couple of transverse shear stresses t is applied to the bilayer leaflets, the latter acquire after a transient response a finite relative velocity Dv (fig. 1.18). The consequence is an interleaflet friction stress defined by the following linear relation [START_REF] Zgorski | Surface Shear Viscosity and Interleaflet Friction from Nonequilibrium Simulations of Lipid Bilayers[END_REF]: The friction coefficient may not exists if the stress-velocity characteristics does not show any linear regime. As a rule, the linear regime extends only over a finite range of applied stress of velocity. When the solicitation is large, a velocity dependent coefficient b(v) is expected. A shear-thinning behavior corresponds to a decreasing t/Dv ratio.

t = b ⇥ Dv (1.32)
The relation between the linear response friction coefficient b and the solid friction coefficient µ is by no means obvious. The latter is a macroscopic concept while the former is a local microscopic one. The transition from lubricated (or viscous) friction to dry friction is yet not completely understood.

Finally, one must keep in mind that the dynamical predictions from coarse-grained simulations cannot be quantitatively compared to real physical quantities. In coarsegrained simulations, the dynamics of the system is comparatively faster compared to all-atom simulations. Various factors contribute to make this happen. CG beads may not scatter momentum as the original molecular group. They are smooth and lack the steric hindrance caused by small atomic groups (methyles, hydrogen). Finally, in the case of Martini, the mass repartition is not identical to the one of the original molecules. Coarse-graining dynamical results must therefore be interpreted qualitatively, or comparatively.

Existing method

The experimental determination of b (and b 0 ) is difficult. Evans and Yeung suggested that b dominates the resistance of a bilayer when pulling a lipid nanotube from a giant vesicle, with a micropipette or an optical tweezer device [START_REF] Evans | Hidden dynamics in rapid changes of bilayer shape[END_REF]. Tube pulling experiments have since become a standard protocol for probing membrane physical properties, including the case of living cells [START_REF] Gauthier | Mechanical feedback between membrane tension and dynamics[END_REF]. Leroy et al. were able to estimate the dissipation induced by the friction of the interfacial water beneath a supported lipid bilayer deposited onto a mica surface using a surface force apparatus (SFA) [START_REF] Leroy | Probing biolubrication with a nanoscale flow[END_REF]. More recently, simulations by Schlaich et al. [START_REF] Schlaich | Hydration friction in nanoconfinement: from bulk via interfacial to dry friction[END_REF] investigated in details the nature of the friction between amphiphilic surfaces separated by a variable amount of interfacial water using atomistic molecular dynamics simulations.

Seifer and Langer [START_REF] Seifert | Viscous Modes of Fluid Bilayer Membranes[END_REF] showed how the relaxation dynamics of the transverse membrane undulation modes depend on h and b, and interpreted in this way experimental data from inelastic neutron scattering [START_REF] Pfeiffer | Neutron spin echo study of membrane undulations in lipid multibilayers[END_REF]. This formalism was successfully used by den Otter and Shkulipa for estimating b for various numerical model of lipids, using equilibrium molecular dynamics (MD) [START_REF] Otter | Intermonolayer Friction and Surface Shear Viscosity of Lipid Bilayer Membranes[END_REF]. Müller and Müller-Plathe showed how the bilayer friction and viscosity parameters could be obtained from reverse non-equilibrium molecular dynamics (RNEMD) simulations [START_REF] Müller | Determining the Local Shear Viscosity of a Lipid Bilayer System by Reverse Non-Equilibrium Molecular Dynamics Simulations[END_REF]. Falk et al. managed to determine b for a coarse-grained bilayer in both the fluid and gel states using RNEMD [START_REF] Falk | Interleaflet sliding in lipidic bilayers under shear flow: comparison of the gel and fluid phases using reversed non-equilibrium molecular dynamics simulations[END_REF] for shearing the solvent on both sides across the bilayer. In particular, the authors reached the conclusion that there was only minor sliding velocity effects at the solvent-lipid interface. Using a similar method, Zgorski et al. determined b and the membrane transverse viscosity for DPPC Martini models [START_REF] Zgorski | Surface Shear Viscosity and Interleaflet Friction from Nonequilibrium Simulations of Lipid Bilayers[END_REF]. The approaches of den Otter and Shkulipa (fig. 1.19), or Falk et al. cannot easily be generalized to supported bilayers in close interaction with a flat solid surface. It is known for instance that a proximal solid surface influences the lipid diffusion dynamics, as shown in Scomparin et al. [START_REF] Scomparin | Diffusion in supported lipid bilayers: Influence of substrate and preparation technique on the internal dynamics[END_REF]. There is therefore a need for simple approaches for determining the friction properties of lipid bilayers interacting with solid surfaces.

A natural idea consist in pulling directly on various system components (lipid or water layers) and measuring the resulting velocity profile. Alternatively, information can be obtained by observing how a system initially prepared with mutual nonvanishing relative sliding velocities relaxes to its equilibrium state. When linear response from the system holds, it is expected on general grounds that both approaches give consistent results. In the present work we show how a constant pull force and momentum relaxation methods can be used to determine the interleaflet friction coefficient in the simple case of a coarse-grained lipid bilayer in water. T HE present chapter focuses on the methodology employed during this thesis, that is to say how the systems are prepared and analysed. In particular the molecular dynamics engine GROMACS [START_REF] Van Der Spoel | GROMACS: fast, flexible, and free[END_REF] will be presented as well as the coarsegrained lipid force field Martini [START_REF] Marrink | The MARTINI Force Field: Coarse Grained Model for Biomolecular Simulations[END_REF].

Molecular Dynamics software

GROningen MAchine for Chemical Simulations (GROMACS) is the molecular dynamics software used during this thesis (fig. 2.1). It is primarily designed for studying biochemical molecules (proteins, lipids and nucleic acids) but can also be used for polymers. This is a fast and popular and free software (open-source software: http://www.gromacs.org/) [START_REF] Van Der Spoel | GROMACS: fast, flexible, and free[END_REF]. Due to popularity, there are many available tutorials dedicated to GROMACS simulations. A nice one can be found on the following web site: (http://www.mdtutorials. com/gmx/) [START_REF] Lemkul | From Proteins to Perturbed Hamiltonians: A Suite of Tutorials for the GROMACS-2018 Molecular Simulation Package [Article v1.0[END_REF]. The list of strong points in favor of GROMACS are the following ones:

1. It is a well distributed, well tested software, with an active community of users.

2. It is the software for which the Martini force field was originally designed and implemented.

3. It is user friendly and allows for a rational and rigorous organization of the user data files.

4. It comes with a complete integrated suite of dedicated analysis tools, e.g. for the measuring thermodynamic parameters, diffusion coefficients, density distributions. It also provides system configurations building and editing tools.

5. It works on parallel clusters, and is optimized for simulations of molecules in water solutions.

However, there are also shortcomings, or weaker points such as:

1. A lack of flexibility and fine control on the molecular dynamics conditions. It is not possible to exert anything but constant or harmonic forces on the system. Continuous temperature annealing or quenching of the system is not possible.

2. Optimization was sometimes achieved at the expense of performance. In particular the conservation of energy is not very good.

3. It offers only a limited choice of thermostat options.

Martini force field 2.2.1 Overview of Martini

We use in this study the Martini force field (fig. 2.2) [START_REF] Marrink | The MARTINI Force Field: Coarse Grained Model for Biomolecular Simulations[END_REF][START_REF] Marrink | Coarse grained model for semiquantitative lipid simulations[END_REF]. As in the case of GRO-MACS, the software was initially developed in the city of Groningen. Its name is in reality the city's nickname where a Martini tower can be found. The main web site is http://cgmartini.nl/, where the force field parameters can be downloaded and the latest developments are regularly presented.

Martini is a coarse-grained force field originally designed for phospholipids, and since extended to describe the other components of a lipid membranes: cholesterol, proteins, glycolipids, synthetic polymers and amphiphiles. . . The Martini model is based on a four-to-one mapping (four-to-two in the case of ring-like molecules, and other special cases). In phospholipid systems, the water bead represents four water molecules, 1 bead represents the charged choline group, 1 bead the charged phosphate group, 2 beads for the glycerol, 1 bead for about 4 carbons in the fatty acid tail (e.g. butane) (fig. 2.3). Each type is also distinguished by a number from 1 to 5 which indicates the degree of polarity (ex: C1, C2, . . . , C5) or by a letter in connection with the hydrogenbonding capabilities (d = donor, a = acceptor, da = both and 0 = none). These are the subtypes. [START_REF] Periole | Chapter 20 the martini coarse-grained force field in[END_REF].

The mass of the coarse-grained beads is identical and fixed to 72 amu (atomic mass units or g.mol 1 ). This corresponds to the weight of 4 water molecules. (see appendix A for more information about the DSPC parameters).

Interactions potentials

The Non-bonded interactions are based on a shifted Lennard-Jones (LJ) 12-6 potential function (eq. 2.1) where 1.12s ij is the equilibrium mutual distance between two particles (ij), and e ij the strength of the interaction between these particles [START_REF] Marrink | The MARTINI Force Field: Coarse Grained Model for Biomolecular Simulations[END_REF]. It is necessary to define a matrix of Lennard-Jones parameter for every possible pair in the Martini model. As the number of bead types and subtypes is already quite large, it is necessary to restrict the number of entries in this matrix. The rationalization of the list of possible Lennard-Jones parameters consists in defining an attractivityrepulsivity unidimensional scale. In the case of standard Martini bead types, the value of s is always set equal to 0.47 nm except but in the super-repulsive situation, involving the Q-type and the more hydrophobic beads (C1, C2), where s is set equal to 0.62 nm.

V LJ = 4e ij " ✓ s ij r ◆ 12 ✓ s ij r ◆ 6 # (2.1)
The attractivity-repulsivity scale is divided into 10 different levels, each of which characterized by a strength e (kJ/mol) and radius s (nm).

In the case of CG Martini, the choice of a Lennard-Jones potential is a matter of simplicity and convenience. If this functional form is well accepted when describing atomistic models, its relevance for a coarse-grained models is not obvious. design choice for Martini was to stick to simple standard Lennard-Jones functions instead of parameterizing from scratch an effective inter-bead potential. As a result, the r 12 short-range repulsion might appear slightly too sharp and unphysical. If it does not affect the thermodynamic properties, which Martini reproduces well, it certainly has consequences as far as the dynamical properties are concerned: collisions between two Lennard-Jones beads are different from collisions between soft potential beads, and the transfer of momentum and the related transport coefficients are poised to be different. Nevertheless, Martini has become a standard and a study of its dynamical predictions is certainly worth considering.

The second part of the non-bonded interactions is the shifted Coulombic potential energy (eq. 2.2). As Martini is meant to simulate physiological situations in the presence of ionic buffers, long range electrostatic interactions are expected to be screened on rather short distances. As a result, and for computational efficiency, it was decided that the Martini coulombic interactions would be truncated also at 1.2 nm. Moreover, the Martini water beads do not have any electric dipole, and the solvent cannot reproduce the large value of the static water dielectric constant e r = 70 80. As the dielectric constant of a lipid bilayer interior is much smaller (e r ' 2) an effective intermediate value of the dielectric constant e r = 15 is used in the Martini force-field [START_REF] Marrink | The MARTINI Force Field: Coarse Grained Model for Biomolecular Simulations[END_REF].

U elec (r) = q i q j 4pe 0 e r r (2.2)
Bonded interactions in the Martini description of phospholipids are quite simple. The bond spring constant is set to K bond = 1250 kJ.mol 1 .nm 2 , and the angle bending constant to K angle = 45 kJ mol 1 . Equilibrium bending angles are usually 180 except for the phosphate-glycerol which is 120 and the glycerol-carbon link which is not subject to a bending potential. The equilibrium distances between consecutive beads (bonding distance) are 0.47 nm except the glycerol beads separated by 0.37 nm. Dihedral potentials are not necessary and do not appear in the DSPC lipid parameterization [START_REF] Marrink | The MARTINI Force Field: Coarse Grained Model for Biomolecular Simulations[END_REF] .

In the Martini version used in this work, DPPC chains contains 4 apolar beads while DSPC chains have 5 apolar identical beads. The two models are therefore differents. There is a degree of appreciation as far as associating a number of carbons to a number of apolar beads. Standard lipids have an even number of beads in their chain (L (lauric acid) 12, M (myristic acid) 14, P (palmitic acid) 16, S (stearic acid) 18). It is impossible with a 4:1 mapping to enter into this degree of detail. As a result, the Martini model cannot distinguish for instance a DMPC from a DPPC molecule, while experimentally there are clear structural and thermodynamic differences.

Cut-off distance

Lennard-Jones interactions are shifted and truncated in order to vanish exactly at some cut-off distance value r cut = 1.2 nm. When using the force-shift directive, the interaction potential starts to be modified at r shift = 0.9 nm and vanishes completely for r 1.2 nm [START_REF] Marrink | The MARTINI Force Field: Coarse Grained Model for Biomolecular Simulations[END_REF][START_REF] Van Der Spoel | GROMACS: fast, flexible, and free[END_REF][START_REF] Lindahl | GROMACS 3.0: a package for molecular simulation and trajectory analysis[END_REF]. Both the force and the potential are continuous at the cut-off radius in order to maintain an acceptable degree of conservation of energy from the Verlet leap-frog integrator. Using a short cut-off makes it possible for the algorithm to ignore all the particle pairs separated by a distance larger than some augmented cut-off value (say r 0 cut = 1.4 nm) which are truly non interacting, resulting in shorter computations. Latest versions of Martini tend to use a smaller cut-off to reduce the computational effort. However, changing the cut-off is not harmless, as it is fully part of the parameterization of the model. Changing the cut-off radius without changing e and s leads to a different force-field.

Organization of a simulation

A simulation is planned as follows. First the initial configuration is prepared by putting together the various components and submitted the system to a sequence of short energy minimization and thermal equilibration runs. This is a crucial step because poorly prepared configurations may result in a crash of the software, due to excessive built-in stresses. Quite often, the effort required for designing a starting configuration is the most demanding of all.

Once the system is prepared, production runs are generated, on the scale that depends on the question one wishes to answer. This step is most efficiently carried out on High Performance Clusters, shared by many users and ruled by a system of queues.

The production runs create some amount of data, stored for analysis purposes. The analysis can be on the fly (simple thermodynamic quantities, or if the software allows it) or a posteriori. The most emblematic kind of data produced by MD softwares are the trajectory files. They corresponds to the periodic recording of all beads positions and velocities (sometimes only positions) into a usually large file. The size of the produced data increases with the simulation time, the dumping frequency and the system size. In some cases the trajectory size is the limiting factor of a MD study.

Initial system set-up

The very first step is to define the "topology" of the system, namely to list the number of molecules (or residues) and to enumerate the nature of the beads and the bonding geometry of these molecules. The topology file refers to the force field parameters ff.itp that are common to all the possible simulations using a given version. For instance DSPC molecules were already defined and available. The topology file must also states explicitly how many molecules of each type are to be included in the simulation.

Our system is typically made of a bilayer containing 512 DSPC molecules for 10 times more water beads, i.e. 5120 CG water beads. Once the topology is defined, one must create a numerical table of all the beads positions (3 coordinates per bead), called configuration file (configuration.gro in the GROMACS world). The box size must also be specified, which in the orthogonal case requires three numbers {L x , L y , L z }, (L x , 0, 0), (0, L y , 0) and (0, 0, L z ) being the three primitive vectors of the periodic images lattice. There are infinitely many different configurations sharing the same topology.

Building the initial configuration require the assistance of home-made scripts or configuration builders. One of them is the software named Packmol. This tool builds an initial configuration based on the concept of packing optimization. The user must provide a single molecule configuration template for each species present in the system and instruction file detailing targets and restrictions regarding the position of the molecules. Packmol then inserts molecules according to the template model by inserting, translating and rotating while accounting for the geometrical constraints (e.g. size of the box, distance between the two monolayers etc. . . ) [START_REF] Martínez | PACKMOL: a package for building initial configurations for molecular dynamics simulations[END_REF]. Another way consists in using the Martini Maker section of the on-line website CHARMM-GUI (http://www.charmm-gui.org/?doc=input/martini.bilayer). The server offers the possibility to choose the lipids of interest, assemble the bilayer (or other structures) and solvate it with the desired water amount. In addition it minimizes the energy of the resulting structure [START_REF] Jo | CHARMM-GUI: A web-based graphical user interface for CHARMM[END_REF][START_REF] Qi | CHARMM-GUI Martini Maker for Coarse-Grained Simulations with the Martini Force Field[END_REF].

Energy minimization

After the initial system preparation, the simulation box contains molecules packed by geometrical constraints (especially by using Packmol). Such configuration state is in general unsuitable for starting a MD simulation. In many cases, the random or periodic location of the beads causes an important local stress, when two beads interacting via a Lennard-Jones hard core repulsion are situated too close from each other. At the opposite the system may open an empty cavity. Starting a simulation from a internally stressed configuration releases a large amount of energy and generate forces that thermostats find it difficult to control, eventually leading to a crash. In order to overcome the problem, it is required to relax the inner stresses by minimizing the energy. This is equivalent to a submitting the system to a viscous dynamics at vanishing temperature.

With GROMACS, the input file file.mdp (molecular dynamics parameter) summarizes all the parameters of a simulation. It specifies the integration time step (dt), the type of integrator (md for molecular dynamics, em for energy minimization,. . . ), the thermostat, the barostat and so on.

In this work, the Steepest Descent energy minimisation is used. It attempts to minimize the multidimensional energy function by displacing the particles along the direction of the gradient rV. However, due to the large dimensionality and the complexity of the energy landscape, the length of the discrete gradient descent steps must be tuned in permanence. The algorithm starts with a displacement of all atoms equal to h 0 . Then the potential energy and the forces F are computed, and the new atom virtual positions are calculated following (eq.2.

3).

r(n+1

) i = r(n) i + F(n) i max i (|| F(n) i ||) h (n) (2.
3)

It is checked whether the potential energy had actually decreased following this move, and the move is accepted according to the rule n) : the new positions are accepted and h (n+1) = 1.2h (n) • if V (n+1) > V (n) : the new positions are rejected and h n+1 = 0.2h n and the step length modified in accordance. This algorithm stops after a certain number of steps has been reached or when the potential energy is no longer evolving.

• if V (n+1) < V (

Equilibration

Following the minimization step, the system can be considered as relaxed at zero temperature. However it may not have the required density nor thermal energy. Even if the velocity distribution is known, it is not enough to draw at random the velocities values. It is also necessary to raise the potential energy to its canonical average value, including the thermalization of all the harmonic degrees of freedom. It is therefore necessary to simulate the system with a thermostat. If in addition a fixed pressure must be set, a second thermalization run in the presence of a barostat is required. The system is therefore subject to a sequence of NVT and NPT simulations. It may happen that Nose-Hoover (v-rescale) and Parinello-Rahman thermostats and barostats fail to bring the system to the desired state. It is then necessary to change the strength of the coupling with the barostat/thermostat (a less physical but more efficient situation) or to use the Berendsen barostat/thermostat which is very stable.

Both of them have a fixed integration time equal to 0.020 ps and are inside a semi-isotropic box ((x, y) and (z)).

Determination of the optimal system size in the NPT ensemble

While the above considerations are generic and apply to all the molecular simulations, we now detail a more specific procedure that only concerns the approach used in this manuscript.

We first enforce the temperature and pressure with a V-rescale thermostat with time coupling constant t t = 1.0 ps and a Berendsen barostat with time constant t p = 12.0 ps as suggested by the GROMACS developers. As soon as the pressure stabilizes, the barostat is changed to a semi-isotropic Parrinello-Rahman to simulate a true NP x P z T ensemble. The system is simulated for 40 ns with P x = P z = 1 bar and the average box size is determined.

At the end of the NPT step, the simulation box size (L x , L y , L z ) is edited by using the built-in command (gmx editconf ) of GROMACS and is set to the desired averaged values (hL x i, hL y i, hL z i).

Verification of the bilayer state

It is possible to check the validity of the structure obtained at the end of the previous steps. Such structural characterization include for instance parameters such as the area per lipid (1.1) and the order parameter.

The area per lipid is given by the ratio of the total area to the number of lipids of one leaflet, provided undulations of the membranes are neglected. For system comprising 512 lipids, the bilayer patch is almost flat. It is also possible to use a software named FATSLIM [START_REF] Buchoux | FATSLiM: a fast and robust software to analyze MD simulations of membranes[END_REF] which determines the parameters more accurately, by accounting for the height variations of the membrane. The APL depends on temperature, it changes discontinuously at the melting transition, and increases steadily otherwise.

The chain order parameter measures the alignment of the lipid tails with the bilayer normal direction. The concept originates from deuterium quadrupolar splitting nuclear magnetic resonance (D-NMR) analysis, a spectroscopic technique sensitive to the average angular orientation of the C-H (C-D) bonds with respect to the direction of the magnetic field. It was shown that this order parameter was equivalent to the average angular distribution of the C-C bonds in the aliphatic chains. In the Martini model, single carbons atoms are not resolved, but it is assumed that the vector linking two consecutive beads along the chain is representative of this C-C orientation. The order parameter is therefore defined as

P = 1 2 (3hcos 2 (q)i 1) (2.4)
with q the angular tilt of the carbon bead-carbon bead bond with respect to the bilayer normal direction z. This definition is similar to the order parameter introduced in the field of nematic liquid crystals where one has:

• P = 1: total alignment (q = 0) between the bonds and the normal bilayer

• P = 0.5: total anti-alignment (q = p/2).

• P = 0: random isotropic orientation of q Tails are expected to be more elongated, stretched and aligned in the gel phase than in the fluid phase (fig. 2.6). At the melting transition, the order parameter decreases discontinuously.

The Martini website provides a script do_order to carry out this order parameter analysis. The GROMACS tool gmx order is an alternative. According to the figure (2.6), we can admit that our lipid are more aligned in the gel phase than the fluid phase due to a higher order parameter. Snapshots in gel and fluid phase are provide in the figure (3.1) and (3.2).

Production runs

Production runs are the one where the material used in the analysis of the simulations is produced. It can consist in a single very long run, or a sequence or equivalent short runs.

In our case, we used many (50-150) repeated moderately long simulation to acquire the necessary data.

During production, the GROMACS engine gmx mdrun records periodically the thermodynamic data associated to the simulation (bow size, temperature, energy, pressure,. . . ) as time series (edr files. It also dumps periodically the bead trajectories: positions and/or velocities in a compressed binary format (trr files).

Trajectories are often required during the simulation post-treatment, or analysis. They represent a massive amount of data, of the order of many gigabytes, and may limit the size of the investigated systems due to storage limits.

In our case, despite the moderate simulation time we had to dump the trajectories at a high rate. We therefore decided to erase systematically the trajectories following analysis. In this precise situation, it is cheaper and more convenient to regenerate the trajectories rather than to archive them.

Prior to performing this sequences of runs, we first created a reference NVT simulation (50 ps) from which where extracted up to 500 initial starting configuration files, using the conversion command gmx trjconv.

Analysis tools

Visualization programs

The spatial representation of a molecular system and its organization is a very important operation. More insight is obtain from a visual representation of complex objects than from graphs of numbers. Quite often, abnormal MD simulation results find an immediate explanation when a picture of the simulated system becomes available. Vizualisation tools enable to convert rows and lines of number into human friendly representations. These pictures convey at a glance an large number of essential information. In our case, it allows to look at the structure of the system and the global organisation of the lipids as far as (fig. 2.7 and 2.8).

Among the several available tools, we used VMD (Visual Molecular Dynamics) to monitor the evolution of the system during the simulation [START_REF] Humphrey | VMD -Visual Molecular Dynamics[END_REF]. We also generated appealing snapshots of bilayer membranes with Ovito (Open Visualization Tool) [START_REF] Stukowski | Visualization and analysis of atomistic simulation data with OVITO the Open Visualization Tool[END_REF].

Home made tools

During this work, many script were written in order to automatize the simulation cycles, to analyse the trajectory files and treat the data of interest. They were preferentially written in Python 3 or in C. Python 3 was mainly used to control the simulation cycles, to copy and create folders where temporary and auxiliary files were stored for the simulations and analysis. These scripts were also used to call C executables for some analysis parts requiring intensive fast computations.

FATSLiM

FATSLiM (Fast Analysis Toolbox for Simulations of Lipid Membranes) is convenient tool to analyse MD configations of membranes. For instance, when providing index files for the lipid headgroup beads, it may calculate the thickness of the membrane, identify the leaflets or estimate the area per lipid ( flat or curved membranes either) [START_REF] Buchoux | FATSLiM: a fast and robust software to analyze MD simulations of membranes[END_REF].

Bootstrap analysis

The bootstrap is an empirical statistical method that provides a quantitative estimate for the confidence interval of an average sampled quantity [START_REF] Press | Numerical Recipes in C (2Nd Ed.)[END_REF]. It is often qualified as a re-sampling method [START_REF] Efron | Bootstrap Methods: Another Look at the Jackknife[END_REF]. In the absence of extra information regarding the nature of the statistical process under investigation, the bootstrap approach uses only available sample values to build this estimate.

Considering a set of N sampled values S (0) = {x i }, i = 1 . . . N as main input information, one can generate an number M of synthetic samples S (a) = {x

h f i (a) = P N i=1 f (x (a) i ) N (2.5)
as a function on the synthetic sample S (a) , provides us with a confidence interval 2s for the sampled average, using the following estimator

s 2 ⌘ 1 M 1 M X a=1 0 @ h f i (a) 1 M ( X b h f i (b) ) 1 A 2 (2.6)
with M large enough, M = 10 in our case. As discussed in [START_REF] Press | Numerical Recipes in C (2Nd Ed.)[END_REF], the bootstrap approach makes optimal use of the sole available information contained in S (0) . The bootstrap analyses proceeded as follows. In each case, a number N s of realisations X

hX µ i = 1/N s P a X (a)
µ (t). The relative variation of the quantities of interest, deduced from hX µ i(t) such as a plateau value, or an average velocity, provides a confidence interval for the quantity of interest. Bootstrap amounts to randomly selecting subsets of the working sample in order to infer its intrinsic variability. The whole procedure is a kind of Monte-Carlo estimate of an average value, using the working sample as integration space. For large and independent enough samples, the bootstrap approach should converge to the true variability of the desired average value.

The figure (2.9) presents an example of a set of ten re-sampled displacements from a pool of 50 simulations for a pulling force equal to 50 kJ.mol 1 .nm 1. This statistical method provides an interesting point of view of the dispersion of the variable of interest around the average.

Rheological properties of free lipid bilayers

"It is by intuition that we discover and by logic we prove."

Henri Pointcaré (1854-1912)

Introduction

In the present work we show how a constant pull force and momentum relaxation methods in the linear regime can be used to determine the interleaflet friction coefficient in the simple case of a coarse-grained lipid bilayer in water.

We chose to study DSPC molecules, parameterized using the version v2.0 of the Martini model [START_REF] Marrink | The MARTINI Force Field: Coarse Grained Model for Biomolecular Simulations[END_REF]. DSPC lipids possess two saturated 18 carbons chains. This choice was driven by experimental considerations, as DSPC supported lipid bilayers obtained by Langmuir deposition constitute a robust and well studied model systems [START_REF] Scomparin | Diffusion in supported lipid bilayers: Influence of substrate and preparation technique on the internal dynamics[END_REF][START_REF] Hemmerle | Controlling interactions in supported bilayers from weak electrostatic repulsion to high osmotic pressure[END_REF]

Methodology

Relaxation and forced sheared experiments.

Our purpose is to characterize the response of a supported bilayer sheared parallel to its longitudinal xy directions, as it may provide clue on the experimentally observed drag friction reduction upon coating surfaces with deposited lipid mono or multilayers. The determination of the interleaflet friction bilayer in water solution is therefore a first step towards the desired answer, which will be later extended to lipid layers deposited onto solid surfaces.

Two strategies were used in the present work. Both were implemented using the Gromacs molecular dynamics simulation tool [START_REF] Hess | GROMACS 4:Algorithms for highly efficient, Load-Balanced, and scalable molecular simulation[END_REF]. In the first approach, here referred as constant pull force (CPF), a non-equilibrium stationary pull of each membrane leaflets was set up, resulting in a constant drift velocity of the bilayer. The pulling force-velocity ratio gives access into the value of the interleaflet friction coefficient b. In the second approach, referred as force kick relaxation (FKR), the relaxation stage of a leaflet consecutive to an initial step increase in its center of mass (COM) velocity was measured. The displacement response curve of the leaflet gives another estimate of the interleaflet coefficient b. It provides in addition a direct picture of the transient bilayer response following a sudden shear force kick.

A linear response regime is expected provided the pulling forces (CPF) and initial velocities (FKR) remain below their respective threshold values. In the linear regime, both drift velocities and displacements compete with random equilibrium fluctuations, a situation corresponding to a small Peclet number. The extraction of the signal (drift displacement and velocity) out of the noise (equilibrium fluctuations) requires averaging over many independent simulation runs. The statistical significance of the bilayer response curves was estimated by means of a bootstrap statistical procedure. In our case, for every simulation condition (external constant force, or initial force kick), a sample of ca N s ⇠ 50 150 independent runs was subject to random reweighting, in order to infer a reliable value of the statistical uncertainty associated with sample averaging. Details on our numerical simulation procedure and the associated statistical analysis are deferred to the appendix.

Standard hydrodynamic description

A natural interpretation frame for our numerical simulations is the classical hydrodynamics model. In this framework, both lipid leaflets are described as rigid solid slabs (thickness L b , area A), surrounded by a water layer considered as a Newtonian fluid (thickness L w , viscosity h). Inertia of lipids (leaflet mass M) and fluid (volumetric mass density r) components are accounted for. The upper and lower leaflets move with respective velocities V u , V d along the horizontal x direction. Water is described by a Eulerian velocity field v(z)ẽ x , where the vertical coordinate z, normal to the bilayer, varies in the interval z > L b /2; z < L b /2 with periodic boundary conditions v(z + L) = v(z) (PBC), and x is one of the horizontal direction, without loss of generality (Fig. 3.3). The fluid is subject to a Newtonian shear stress t zx (z), abbreviated as t(z). Sticking boundary conditions at the lipid water interface z = ±L b /2 are assumed (or equivalently an infinite lipid-fluid friction b 0 = •). We assume that leaflets experience a friction proportional to their mutual relative sliding velocity V u V d , leading to an interleaflet shear stress t(0) obeying

t(0) = b(V u V d ), (3.1)
with b the inter-layer friction coefficient. An average fluid velocity can be defined as:

V w = 1 L w Z L b /2+L⌘ L b /2 L b /2 dz v(z). (3.2)
In addition, we consider the possibility to act upon each leaflet, and the water layer, by means of a uniform force acting on the center of mass of the corresponding subsystem. Such forces are respectively denoted F u , F d , F w , and directed along x. For convenience, one introduces the corresponding stresses f µ (F µ = Af µ ) with µ = u (upper leaflet), µ = d (lower leaflet) and µ = w (water region). One restricts ourselves to the physical case of a vanishing total force F u + F d + F w = 0, henceforth preserving the total momentum rL w V w + MV u + MV d of the hydrodynamic system.

The stationary solution of the hydrodynamic problem corresponds to a parabolic flow. Two stationary velocity profiles are of particular interest. The linear Couette profile corresponds to

f u = f d , f w = 0, V u = V d , V w = 0 and 2 ✓ b + h L w ◆ V u = f u . (3.3)
The Poiseuille flow profile corresponds to

f u = f d = f w /2, V u = V d and 6 h L w (V u V w ) = f u = f w 2 . ( 3.4)
The relation above can be further simplified as the total momentum is assumed to vanish 2Mf u + rL w f w = 0. Both flows are represented in Fig. 3.4.

Viscoelastic relaxation model.

As the Results section demonstrates, the hydrodynamic model is useful but does not accurately represent the observed numerical behavior. We therefore introduce here a more general viscoelastic model. We assume that a transient linear response of a bilayer subject to a suddenly applied external force exists, that can be expressed by means of a retarded memory function. Using the same notations as above, but now with time dependent velocity fields V u (t), V d (t), V w (t) one has:

M A Vu (t) = Z t • ds [g bu (t s)(V d (s) V u (s)) + g wu (t s)(V w (s) V u (s))] + f u (t); (3.5) M A Vd (t) = Z t • ds [g bu (t s)(V u (s) V d (s)) + g wu (t s)(V w (s) V d (s))] + f d (t); (3.6) rL w Vw (t) = Z t • ds [g wu (t s)(2V w (s) V u (s) V d (s))] + f w (t). (3.7)
The retarded response involves two memory functions. A first kernel g bu (t) accounts for the interleaflet interaction, including interleaflet dynamic friction b, lipid inertia as well as viscoelastic lipid elastic tilt and stretch modes. A second kernel g wu (t) accounts for all the water leaflet interactions, which possibly includes solvent sliding friction, retardation of the fluid motion due to inertia, and again viscoelasticity arising from lipid tilt and stretch. The same kernel is used for both leaflets, as a consequence of the up-down z symmetry of the flow. External stresses f u (t), f d (t), f w (t) are arbitrary functions of time.

We now restrict ourselves to two main situation of interests, namely Couette

V w = 0, V u (t) = V d (t), f u (t) = f d (t), f w = 0 and Poiseuille V u (t) = V d (t) = V w (t)rL w /2M, f u (t) = f d (t) = f w (t)/2 (see Fig. 3.4
). The retarded motion equations are in the Couette case:

M A Vu = Z t • ds (2g ud + g wu )(t s)V u (s) + f u (t); V w = 0, (3.8)
and in the Poiseuille case:

M A Vu = Z t • ds g wu (t s) ✓ 1 + 2M Ar w L w ◆ V u (s) + f u (t); rL w Vw = Z t • ds g wu (t s) ✓ 2 + Ar w L w M ◆ V w (s) + f w (t). (3.9)
Of particular importance in the present study is the response to a couple of force kicks (Couette case)

f u = f d = M A V 0 d(t), (3.10)
that confers instantly a momentum MV 0 ẽx to the upper leaflet, and MV 0 ẽx to the lower leaflet. Velocity profiles can be inversed by Laplace transforms of the velocity, stress and memory functions, e.g.

Vu (p) = Z • 0 dt e pt V u (t), (3.11)
leading to

✓ M A p + 2 ĝud + ĝwu ◆ Vu (p) = M A V u (t = 0). (3.12)
In particular, the impulsional displacement

DX u = R • 0 dt V u (t) = V(p = 0) obeys the relation DX u = M A V u (0) 2 ĝud (0) + ĝwu (0) . (3.13)
By contrast, a stationary stress f u should result in an asymptotically constant velocity V u and Vu = 0:

Z • 0 ds (2g ud (t) + g wu (t))V u = (2 ĝud (0) + ĝwu (0))V u = f u (3.14)
One recovers the hydrodynamic limit f u /V u = 2(b + h/L w ) and therefore

(2 ĝud + ĝwu )(p = 0) = 2 ✓ b + h L w ◆ (3.15)
In conclusion, one obtains a useful relation between the impulsional displacement DX u and the interleaflet friction coefficient.

b + h L w = M A V u (0) 2DX u (3.16)
This viscoelastic model assumes a linear relation between forces (the cause) and displacement or velocity (the effect). A master curve X (f) (t) can be introduced to represent the normalized drift displacement (X u (t) X u (0))/f u associated to a step stress H(t) = 1 for t 0 and H(t) = 0 for t < 0 (Heaviside function). This master curve obeys to

X (f) (t) = 0 for t < 0;
(3.17)

M A Ẍ(f) (t) = Z t • ds (2g ud + g wu )(t s) Ẋ(f) (s) + 1 for t 0. (3.18)
In the mean time, a master curve for the normalized displacement X (V) (t) = (X u (t) X u (0))/V 0 can be introduced for the impulsion case, which obeys:

X (V) (t) = 0 for t < 0; (3.19) M A Ẍ(V) (t) = Z t • ds (2g ud + g wu )(t s) Ẋ(V) (s) + M A d(t) for t 0. (3.20)
Both master curves can be related to the memory function 2g ud + g wu in Laplace space. FIGURE 3.4 -When opposing forces are exerted on each leaflet (case a) the resulting stationary state sees the two leaflets sliding at constant relative velocity, surrounded by a uniform solvent velocity gradient profile, as emphasized in the sub-picture (b) where the simulation box boundary has been purposely shifted to sit exactly at the midplane of the bilayer. Case (a-b) are subsequently referred as a Couette situation. When a uniform force is exerted on both leaflets and an opposing force on the solvent beads (case c), a symmetric parabolic velocity profile builds up in the solvent, assuming sticking boundary conditions at the interface with the bilayer (d). Case (c-d) will be subsequently referred as Poiseuille situation. In all cases, the total momentum of the system is constant and vanishes.

Diffusion of the lipids and water centers of mass

Simulations deal with finite size systems, and thermal fluctuations are always present. In our case, the center of mass of each of the three main components of the simulated system (upper (u) and lower (d) leaflets, water (w)) is subject to brownian motion, while the global center of mass is fixed, as required by weak coupling to Nose-Hoover or v-rescale thermostats. It results that the instantaneous kinetic energy of the upper, lower leaflets and water does is not given by the usual equipartition of energy theorem. However, the order of magnitude of the instantaneous kinetic energies M µ V 2 µ /2, µ = {d, u, w} remains of the order of k B T/2. We therefore distinguish the average, non fluctuating hydrodynamic displacements X µ (t), µ = {d, u, w} from the sampled, brownian trajectories X µ (t) (a) , with a an index relative to a given center of mass trajectory realization, or simply X µ (t), µ = {d, u, w} when referring to a generic trajectory. Similarly, one introduces the brownian instance of the velocity response V (a) µ (t) or generically V µ (t). In order to quantify the magnitude of the brownian fluctuations acting on the positions X µ (t), one naturally defines the diffusion coefficient D COM,µ of the center of mass of the subcomponent µ (not to be confused with the molecular diffusion coefficient), based on the mean quadratic displacements h(X µ (t) X µ (0)) 2 i. Hydrodynamic and brownian displacements are related by canonical ensemble averages X µ (t) = hX µ (t)i. So are the velocities V µ (t) = hV µ (t)i.

Expression (3.16) relates the dissipation b + h L w to the normalized displacement

DX u /V u (0). Noting that DX u /V u (0) = R • 0 dt V u (t)/V u (0), one can write ADX u MV u (0) = A MV u (0) 2 Z • 0 dt V u (t)V u (0) = 2 b + h/L w . ( 3.21)
By analogy with brownian motion, where the diffusion coefficient is linked to the velocity autocorrelation function, one has 2D COM,u t ' 2t R • 0 dt hV u (t)V u (0)i, and obtain from (3.21) an heuristic "Stokes-Einstein" relation:

D COM,u ⇠ 2hMV(0) 2 i A(b + h/L w ) ⇠ k B T A(b + h/L w ) . ( 3

.22)

The precise relation between the relative quadratic displacements matrix of the various system subcomponents (leaflets, water. . . ) and the hydrodynamic friction coefficients (b, h. . . ) when the global center of mass is fixed is non trivial. Eq. 3.22 provides however an order of magnitude for D COM,u .

Constant pulling force simulations.

A direct estimate of the asymptotic stationary drift velocity hV µ (t)i obtained as a result of a piecewise constant step increase of the external applied stresses f u = f d , f w = 0 can be obtained by pulling directly on the leaflets. Even though the out of equilibrium features of the molecular dynamics software that we use are somewhat limited, it is possible to exert a constant force to the upper leaflet while exerting the opposite force on the lower leaflet (section B.1). This features comes as part of tools available to perform biased, constrained umbrella sampling simulation schemes. The displacement hX u (t)i can then be read directly from the trajectory and fitted to an affine time function x 0 + tV u .

The possibility of imposing a pulling force for long times enables a quite precise determination of the relative stationary drift velocity of the leaflets.

Force kick relaxation simulations.

Starting from an equilibrium trajectory configuration (reference NVT run), an initial condition C a is prepared by adding an identical V 0 ẽx constant velocity to all the beads pertaining to the upper leaflet, and the opposite velocity to all the beads in the lower leaflet. In the Martini model, all beads possesses the same mass (72 a.m.u., 1008 Da for a DSPC molecule), and the upper leaflet center of mass acquires a finite momentum MV 0 ẽx as a result, with M the mass of all beads in a leaflet. The veloc- ity of the water beads is unaltered. Physically, this corresponds to an instantaneous force torque (MV u ẽx d(t), MV u ẽx d(t)) applied to the bilayer, and the total momentum of the system is preserved. In particular, the system center of mass remains fixed, as required when using a Nose-Hoover or velocity-rescale thermostat. Following the force kick, the kinetic energy of the bilayer is increased by an amount

N b X i=1 m 2 (ṽ i ± V 0 ẽx ) 2 = N b X i=1 m 2 ṽ2 i + N b m 2 V 2 0 +V 0 • ( N b X i=1 ±ṽ i .ẽ x ), (3.23)
where N b stands for the number of beads (center of forces) present in the simulated system, and m the associated (here identical) bead masses. The third term is a statistical O( p N) fluctuation. The kinetic energy term is therefore increased by a relative amount

MV 2 0 3N b k b T . (3.24)
This sets an upper bound V max for the velocity shift V 0 that can be applied without requiring the thermostat to pump energy out of the system, of the order of V max = (N b k b T/M) 1/2 ' 0.2 nm.ps 1 , using N b = 256 ⇥ 14 = 3584, M = N b ⇥ 72 amu and T = 340 K.

Assigning to each leaflet a too small initial velocity value results in lowering the signal to noise ratio, the signal being the forward displacement and the noise the brownian displacement of the leaflet center of mass. Assuming it takes a characteristic time t relax for the leaflets to return to equilibrium, and that a given initial drift velocity V 0 drives the leaflet over a distance DX u , the ratio between ballistic and random displacement reads DX u / p D COM,u t relax at the end of the relaxation stage. If in addition, the simple and naive scaling DX u = V 0 t relax holds, the ballistic to random displacement ratio assumes a familiar Peclet number expression Pe 1/2 with Pe = V 0 DX u /D COM,u .

The displacement DX u (t) is monitored as a function of time t. As each run provides a noisy brownian response DX (a) u (t), the procedure must be repeated many times, until a significant displacement hDX u (t = •)i emerges from the thermal noise. Meaningful information can only be obtained in the linear response regime, i.e. when the ratio DX u /V u (0) is constant up to some uncertainty. Too large velocity kicks V u (0)

V lr deviate from the linear regime and cannot be described within the framework of retarded linear response functions. The velocity scale V lr until which the linear regime is expected to hold must be empirically determined and is expected to be smaller than V max determined above. In the opposite limit, a too low kick V 0 does not give any useful result as the signal to noise ratio becomes too large. Again, to estimate a confidence interval for hDX u i, one resorts to a statistical bootstrap procedure.

Preparation of the initial configurations.

The system was equilibrated first at 340 K (fluid phase) and 280 K (gel phase) using a thermostat and a semi-isotropic barostat (see section B.1). This thermalization stage makes it possible to determine the average system size in the absence of external stress, or equivalently vanishing surface tension, respectively in the fluid and the gel phases. Out of equilibrium simulations were then run for some times, using a thermostat and constant box size conditions (L x , L z ), where L x , L z were the result of the previous step. Coupling to a thermostat was however still required to preserve the mechanical energy of the system. For each phase, configurations from a reference canonical, constant volume (NVT) runs were then periodically recorded and stored, providing a set of up to 150 initial conditions, in relation with the bootstrap and ensemble averaging procedures. The resulting equilibrium lipid bilayer geometrical characteristics are summarized in Table 3

Results

Fluid phase constant pull force (CPF) simulations

The bilayer was submitted to a sequence of increasing pulling stresses f u , resulting in an average displacement curve hX u i. Each external pulling force condition was repeated about 50 times (Table B.1), resulting in a sample set of raw displacement curves Fig. 3.5(A). As seen in this figure, a typical pulling experiment generates a brownian displacement of the leaflet center of mass superimposed with a constant velocity horizontal translation. An example of bootstrap averaging of the trajectories is shown in Fig 3 .5(C). Displacements curves start with a short transient regime, dominated by inertial and viscoelastic contributions. It is followed by a linear regime associated with stationary hydrodynamic dissipation and constant velocity translation V u . The bootstrap analysis shows a dispersion among synthetic displacement curves, only slowly decreasing with the size of the set of trajectories, and inversely proportional to the applied stress f u . Flat averages of the normalized displacement curves hX u (t)i/f u are shown in Fig 3 .6. In the framework of linear response, the averaged normalized displacements should converge to a master curve X (f) (t). This is indeed the case for a set of applied stresses within an interval 4.8 ⇥ 10 5  f u  20 ⇥ 10 5 Pa (applied forces in the range 50  F  200 kJ.mol 1 .nm 1 ). A too small applied stress f u = 10 5 Pa (force F = 10 kJ.mol 1 .nm 1 ) departs from the master curve due to strong brownian fluctuations. Large applied stresses clearly brings about strong deviations from linear response, associated with shear-thinning behavior. In the present case, a value A = 174 nm 2 was used for the area in the force-stress conversion, with f u = 9600F, f u in Pa and F in kJ.mol 1 .nm 1 . Taking the bilayer thickness L b = 4.8 nm as a characteristic length, the upper limit of validity of the linear response regime (50 bars) can be turned into a surface tension f u L b of magnitude 25 mN.m 1 , typical of the oil-water surface tension (35 mN.m 1 ).

The determination of DX u using CPF leads to a value for b + h/L w , following eq. (3.16) equal to 2.75 ± 0.08 ⇥ 10 6 Pa.s.m 1 . This value was further confirmed by using a larger sample of 1024 lipids with the same hydration of 10 water beads per lipid. (A) single leaflet COM displacement X(t) (a) , starting from the simulation box center (ca 6.6 nm) and (B) averaged displacement hXi(t) ' 1/50 P 50 a=1 X (a) (t) vs time. A bootstrap procedure (C) estimates the dispersion s b (X u (t)) caused by the finiteness of the sample {a}. Vertical bars represent the confidence interval of 10 selected points from the second half of the trajectory (5000 < t < 10000 ps), taken as twice the estimated bootstrap standard deviation. The vertical bars are used to provide a confidence interval for the drift velocity (slope of the averaged displacement curve

Fluid phase force kick relaxation (FKR) Couette simulations

Repeated kicks were there applied, starting from 150 different configurations. A bootstrap sample of the upper leaflet displacements is shown in Fig. 3.8. The typical averaged displacement curve hX u (t)i increases first linearly, as a natural consequence of the initial force kick that confers a uniform translation velocity to the leaflet (Fig. 3.9). The initial impulsion dissipates fast and vanishes within 5 ps. The displacement curve starts to decrease, or equivalently the leaflet velocity becomes negative. This peak is followed by a much slower relaxation to an apparent plateau value, also associated with a negative velocity, which extends on a few hundred

i(t) = hdX u /d ti(t)/V 0 .
The velocity starts at an initial value of 1, decreases fast to 0 (coinciding with the sharp peak in the displacement curve) reaches a negative minimum and finally slowly regresses to 0 from below, coinciding with the slowly decreasing approach of the displacement plateau value. Inset: close-up look at the minimum.

ps. The apparent plateau value is associated to a relaxation time t relax such that DX u ' hX u (t relax ) X u (0)i, with t relax of the order of 500 ps. The striking main feature of the impulsion relaxation curve is the non monotonic behavior of the displacement X u (t) (Fig. 3.9) and the velocity V u (t) (Fig. 3.10). It is not possible to account for such a behavior without a conservative, elastic contribution to the membrane relaxation. Fig. 3.10 therefore suggests that the mechanical response of a sheared bilayer is viscoelastic in the time scale of t vel ⇠ t relax , with t vel a bilayer internal viscoelastic relaxation time.

As in the constant pulling force experiments, it is possible to define a linear response regime, by plotting the displacement normalized with the initial velocity hX u (t)i/V 0 as a function of time. A master curve X (V) (t) is expected to describe this averaged, normalized displacements in the short and intermediate time regime t  t relax . The normalized displacement velocity Ẋu /V 0 = V u /V 0 is dimensionless, and can be interpreted a velocity autocorrelation linked to the momentum scattering efficiency of the mutual interleaflet molecular interactions.

Figure 3.10 describes the normalized velocity relaxations V u /V 0 for a set of increasing V 0 , and show a deviation of the relaxation from the master curve X (V) (t) at V 0 larger than 0.1-0.2 nm.s 1 . Correspondingly, the effective normalized translation shift (plateau) DX u (t relax )/V 0 starts to increase, pointing again to a shear-thinning behavior. The empirical upper bound V l.r. of the linear response regime is therefore found to be of the same magnitude as the maximal velocity V max deduced from equation 3.24.

While the convergence to a finite plateau value is a reasonable expectation for the averaged displacement curve, simulated trajectories are subject to the thermal motion of the leaflet center of mass, which is expected to be asymptotically dominant at large times. Given a sample size N s , the thermal motion of the sample averaged displacement curve is set to scale as (D COM,u /N s ) 1/2 t 1/2 . The determination of DX u from MD sampling is therefore empirical to a certain extent, as any finite sample average eventually departs from the plateau value. The sample size must be large enough to keep the combination (D COM,u /N s ) 1/2 t 1/2 relax smaller than DX u . Equation (3.22) provides a theoretical estimate of the accuracy of hDX u i. The bootstrap estimate of the variance of DX u (eq. 2.6) is another independent path to estimating the sample dependence in DX u .

It turns out that the condition DX u ⇠ V 0 t relax is not met. Instead DX u /V 0 is of the order of 0.5 ps (see Discussion) and the long relaxation time t relax DX u /V 0 enhances the effect of brownian fluctuations.

Fluid phase Poiseuille flow geometry

Constant pulling rate experiments can be performed in the Poiseuille geometry, when both leaflets are pulled in one direction and the solvent homogeneously pulled in the reverse direction. Assuming that the solvent does not slip at the lipid-solvent interface, the average relative drift velocities obey relation (3.4). Using L w = 3.5 nm in the fluid state (T = 340 K), one finds a value of the coarse grained Martini water viscosity h = 8 ⇥ 10 4 Pa.s (eq. 3.4). Repeating the simulation with a larger number of water beads (10240 solvent beads for 512 lipids, L w = 7.2 nm), the resulting water viscosity changes to h = 7 ⇥ 10 4 Pa.s. Independent simulations using reverse non-equilibrium molecular dynamics [START_REF] Müller-Plathe | Reversing the perturbation in nonequilibrium molecular dynamics: An easy way to calculate the shear viscosity of fluids[END_REF] (with LAMMPS, using an equivalent fluid of truncated Lennard-Jones particles at the same temperature) confirms that the solvent viscosity lies close to h = 7 ⇥ 10 4 Pa.s. The slightly larger value obtained in the presence of a thin water layer is likely to be due to lipid water interfacial effects. The Poiseuille flow simulation design described above can therefore be considered as a viable route to estimate the viscosity of a solvent, provided interfacial effects are small. It is worth noting that the Martini water viscosity value lies quite close to the true water viscosity value, a feature hardly expected from a coarse grained unrealistic water model.

Gel phase CPF and FKR simulations

An ordered phase of the lipid bilayer is obtained at low temperature T = 280 K. A number of bootstrap realizations of the displacement X u (t) corresponding to an initial velocity step of V 0 = 0.9 nm.ps 1 is shown in Fig. 3.11. Normalized averaged displacements curves X u (t)/V 0 are represented in Fig. 3.12, for increasing initial velocities ranging from 0.01 to 0.5 nm.ps 1 . The normalized displacements do not superimpose well, even in the low velocity regime, and a master curve X (V) (t) may not exist at low temperatures. This is especially visible in Fig. 3.13, where the displacements DX u are plotted as a function of the initial velocity V 0 . Unlike the fluid phase, the gel phase curve does not display any established linear regime.

Normalized velocities in the low temperature phase are shown in Fig. 3.14. It is distinctly different from the equivalent fluid counterpart Fig. 3.10. Correspondingly, the initial displacement peak X u (t) X u (0) (inset of Fig. 3.11) is smoother than in the fluid situation. However, the leaflet velocity change of sign during the relaxation stage is seen in both at high and low temperatures.

The CPF result in the gel phase is summarized in Fig. 3.15. It represents the average drift velocity as a function of the pull force (velocity-force characteristics). Unlike the fluid phase, the gel phase do not show any linear regime. A log-scale representation of the velocity-force characteristics seems to indicate a power-low behavior over almost two decades, with apparent exponent hVi ⇠ t 1.44 FIGURE 3.12 -Normalized averaged displacements hX u i(t)/V 0 for a set of increasing impulsions 0.01 and 0.03, . . . , 0.1 nm.ps 1 . The plateau value is clearly increasing with the initial applied velocity, and the normalized displacements do not appear to collapse onto a master curve, pointing to an absence of linear response.

Lipid tilt modes

The non monotonic velocity relaxation curve consecutive to an external force kick at t = 0 cannot be accounted for by a simple hydrodynamic model. Instead, it suggests that some elastic response is involved in the leaflet translational relaxation. All the numerical evidence suggest that the bilayer remain flat, with negligible out-of-plane bending strain. On the other hand, the simulations are held at constant volume, ruling out standard membrane stretching (or compressibility) contribution. We therefore checked whether lipid tilt modes were activated as a result of the interleaflet friction. We estimated the average lipid tilt angle, defined as a polarization vector linking the first and last carbon beads in each chain (Fig. B.1). Figure 3.19 shows, on an enlarged scale, that the average bilayer tilt angle is less than 0.5 at equilibrium. When . τ

0 0.
1.44 FIGURE 3.16 -Velocity V u -applied shear stress t characteristics in the gel and fluid states, using a logarithmic representation. The crossover from linear to shear-thinning regime is visible in the fluid state, while the gel regime seems consistent with a power-low relation with exponent V u ⇠ t 1.44 as a function of the applied stress t.

the bilayer is submitted to a CPF, the angle deviates from its vanishing average, proportionally to the applied force (in the limit of linear response and small angles) as shown in Fig. 3.20. The tilt angle in the fluid phase reaches a well defined asymptotic stationary value, while in the gel phase, the angle seems to be still evolving on the figure time scale (5 ns). In addition, the tilt angle in the gel phase has a larger magnitude than in the fluid phase. The ratio between the average tilt angle and the applied stress is of the order of q/t ' 2.5/24 ' 0.1 .bar 1 or 1.8 ⇥ 10 3 bar 1 with q expressed in radians. Impulsional FKR tilt angle results are shown in Fig. 3 velocity V 0 = 0.09 nm.ps 1 , for which the agreement between the impulsional and constant force approaches holds. For comparison, we also represent equilibrium curves, in the absence of bilayer sollicitation. If the fluid phase tilt angle relaxation tends to relax to its vanishing equilibrium value within the time interval considered (600 ps), the gel phase tilt angle does not recede to zero on the same scale. A residual tilt, reminiscent of a slowly relaxing strain, seems to hold in the low temperature phase, accompanying the observed non-linear rheological behavior.

Discussion

The interleaflet sliding kinetics in the fluid phase display an extended linear regime, both in the CPF and the FKR regimes. Deviations from the linear behavior emerge as the pulling stress exceeds 100 bars (Fig. 3.7) or the initial velocity exceeds 0.2 nm.ps 1 (Fig. 3.10). Inversely, small Peclet number considerations make it impractical to use both approaches for too small initial velocities or pull stresses. This lower limit is not intrinsically related to the physical system considered, but to simulation box size limitations: increasing the sample size amounts to decreasing the collective center of mass diffusion coefficient and the sensitivity of the method. Unfortunately, unconfined large bilayers systems are subject to strong unfavorable undulation fluctuations, and do not constitute a viable option.

The CPF linear regime yields a consistent estimate for b + h/L w , provided one neglects the sliding velocity of the solvent. Substracting off the viscous contribution h/L w , our estimate for b is 2.54 ± 0.10 ⇥ 10 6 Pa.s.m 1 , with h/L w = 2.0 ± 0.12 ⇥ 10 5 Pa.s.m 1 and its confidence interval obtained from our stationary Poiseuille flow pulling simulations. Different other approaches of the Martini water viscosity (7 ⇥ 10 4 Pa.s.m 1 in [START_REF] Otter | Intermonolayer Friction and Surface Shear Viscosity of Lipid Bilayer Membranes[END_REF], or using reverse non-equilibrium molecular dynamics with LAMMPS [START_REF] Plimpton | Fast parallel algorithms for short-range molecular dynamics[END_REF]) provides fully consistent estimates. This value compares well with the A 55 model of den Otter and Shkulipa obtained using a completely different method (RNEMD shear of the surrounding solvent), for which the quoted value for b is between 2.7 ⇥ 10 6 and 2.8 ⇥ 10 6 Pa.s.m 1 . The A 55 is a similar lipid with 5 beads in each chain (as our DSPC) parameterized using the values of the Martini model, and which was simulated at 323 K. The agreement between both models is very good, given the difference between the approaches and also the 13 K temperature difference.

Falk et al. [START_REF] Falk | Interleaflet sliding in lipidic bilayers under shear flow: comparison of the gel and fluid phases using reversed non-equilibrium molecular dynamics simulations[END_REF] simulated a different coarse-grained model (SDK [START_REF] Shinoda | Zwitterionic lipid assemblies: molecular dynamics studies of monolayers, bilayers, and vesicles using a new coarse grain force field[END_REF], see also [START_REF] Seo | SPICA force field for lipid membranes: domain formation induced by cholesterol[END_REF]) obtained a b value of 1.4 ⇥ 10 6 Pa.s.m 1 . The difference may be attributed to a difference of parameterization between the SDK and Martini model. This difference is significant enough to change qualitatively the nature of the gel phase. In the SDK model, the low temperature state is a L 0 b tilted chain phase. It results that the SDK phase displays anisotropic friction properties, with the direction parallel to the tilt direction displaying a b coefficient close to the fluid case (1.3 ⇥ 10 6 Pa.s.m 1 ) and a yield force in the direction perpendicular to the tilt. In our case also, the apparent b value is similar in the gel and fluid case (Fig. 3.18). However, due to the absence of linear regime, we cannot provide anything but a qualitative behavior of the coefficient b.

Zgorski et al. [START_REF] Zgorski | Surface Shear Viscosity and Interleaflet Friction from Nonequilibrium Simulations of Lipid Bilayers[END_REF] performed RNEMD simulations to shear the solvent and the bilayer and obtain b, a similar approach as Falk et. al. They compare the old and new version of DPPC Martini lipids (4 beads chains) and obtain a value in the range of 4.3 ⇥ 10 6 to 5.5 ⇥ 10 6 Pa.s.m 1 . These value are significantly larger than ours (even though not strictly comparable) and also than den Otter and Shkulipa (2.4 ⇥ 10 6 Pa.s.m 1 for the 4 beads chain model A 44 ). Interestingly, Zgorski et al. have determined b for the atomistic CHARMM 36 model, reaching values of the order of 1.1 ⇥ 10 7 Pa.s.m 1 , still an order of magnitude smaller than the experimental estimates of Evans and Yeung or Pfeiffer et al. [START_REF] Evans | Hidden dynamics in rapid changes of bilayer shape[END_REF][START_REF] Pfeiffer | Neutron spin echo study of membrane undulations in lipid multibilayers[END_REF]. More work is therefore needed, both on the experimental and simulation sides, to determine how accurately current atomistic simulations reproduce the local interlayer friction phenomenon.

The FKR predictions for b + h/L w are summarized in Fig. 3.17 and the only numerical values in the linear regime with reasonable error bars are those with V 0 = 0.08, 0.09 and 0.1 nm.ps 1 . They are respectively 2.22 ± 0.57 ⇥ 10 6 , 2.22 ± 0.42 ⇥ 10 6 , 2.67 ± 0.51 ⇥ 10 6 Pa.s.m 1 with confidence intervals inferred from the bootstrap analysis. All three values are slightly smaller than the CPF values, although almost consistent. The estimation of DX µ used in eq. (3.16) was obtained by computing the average stationary value of the relaxation curves featured in Fig. 3.9. The position of the plateau may have been underestimated as the displacements X µ (t) relaxes slowly to its asymptotic limit. Extending the analysis to longer time scales does not improve much the determination of the displacement because the brownian random diffusion increases, and the signal to noise decreases with the elapsed time. It is therefore necessary to both simulate for longer times and to increase in parallel the number of independent trajectories. We therefore conclude that there is a rough agreement between the CPF and FKR methods. Such an agreement is expected based on linear response considerations, which is only seen in the fluid phase. The numerically observed upper limit of validity of the linear response regime V lr = 0.1 nm.ps 1 is remarkably similar to the velocity V max deduced from the system kinetic energy argument. This does not directly prove that the excess of kinetic energy is responsible for the breakdown of linear response, but it indicates that not other process limiting linear response occurs until the V max limit is reached.

The FKR approach gives insight on the transient mechanical response of the bilayer, and predicts a sign inversion of the leaflet COM velocity following the positive impulsional initial velocity. We interpret this phenomenon as the consequence of a slowly relaxing lipid chain tilt angle, causing a reactive (non dissipative) stress contribution. Following the initial velocity kick, an elastic stress builds up, and is further dissipated.

We note that a different transient regime would occur if the initial force kick was applied non uniformly to the bilayer leaflets, for instance on the lipid headgroups only. Linear response arguments suggests that the macroscopic hydrodynamic coefficient b + h/L w must not depend on the location of the applied pulling force or force kick. However, the transient response is expected to depend on the way forces are exerted. Further work is needed to compare the current procedure to other possibilities, that would more closely mimic a real shear force pulling experiment. The uniform pulling force used in the current approach corresponds to a uniform body force applied on each leaflet, due to the fact that all Martini beads have an identical mass.

A transient shear stress response can be inferred from the retarded memory function formalism exposed in the methodology section. This response can be probed by any spectroscopic shear force experiment, using electromagnetic [START_REF] Canale | Nanorheology of interfacial water during ice gliding[END_REF] or piezoelectric vibrations (dissipative quartz-crystal microbalance QCM-D [START_REF] Rodahl | A simple setup to simultaneously measure the resonant frequency and the absolute dissipation factor of a quartz crystal microbalance[END_REF][START_REF] Johannsmann | Springer ser chem sens biosens[END_REF]). So far, none of these techniques reaches the frequency domain of the observed viscoelastic regime. The characteristic ("Maxwell") relaxation time predicted by the Martini model is about 100-1000 ps (Figs. 3.9,3.12)). The connection between Martini coarsegrained and atomistic kinetic properties is quite loose. At room temperature, the Martini lipid diffusion coefficients (ca 70 µm 2 .s 1 in the DSPC fluid phase at 340 K) are predicted to exceed by a factor 10 the actual values (ca 15 µm 2 .s 1 for DSPC at 60 C [START_REF] Marsh | Handbook of Lipid Bilayers en[END_REF]). On the other hand, the predicted Martini water viscosity (0.7 mPa.s) is in reasonable agreement with the real value (1 mPa.s). These examples show that the difference between the coarse-grained and atomistic time scales may stretch from 1 to 10, depending on the phenomenon considered. Assuming that the actual relaxation dynamics associated with the leaflet viscoelastic response falls between 1 and 10 times the corresponding numerical prediction, one may estimate the real Maxwell relaxation time to be of order of 1 ns, and a frequency response possibly in the GHz range.

In addition to the intrinsic membrane elastic response, the water gap probed by the sliding leaflets (Fig. 3.4) is also expected to respond according to a viscoelastic memory pattern. Stokes hydrodynamics predicts that rigid slabs cannot drag the interstitial fluid instantaneously. The stress-velocity response function can be computed analytically for sticking boundary conditions, using for instance Duhamel's principle [START_REF] Carlslaw | Conduction of heat in solids[END_REF]. However, if there were no elastic contribution, the viscous memory function alone would not lead to a reversal of the COM velocity.

The transient response is characterized by a sharp initial increase. We attribute it to the fast loading of the bond springs connecting the beads in the interleaflet area. A characteristic time scale can be obtained as the period t fast ⇠ 5 ⇥ 2p/ p m/K bond of a chain of 5 harmonic springs of stiffness K bond = 1250. kJ.mol 1 .nm 2 and bead mass m = 72 a.m.u (g.mol 1 or atomic mass unit), typical from the Martini force field used in this approach. One finds t fast ⇠ 7 ps, in reasonable agreement with the observed initial peak dynamics in Fig. 3.9.

Let us now estimate the hydrodynamic damping time t hyd resulting from balancing inertia with interleaflet friction. One has t hyd = MN l /(2Ab). With M ' 1000 g.mol 1 , N l = 256, b = 2.5 ⇥ 10 6 Pa.s.m 1 , A = (13.2 nm) 2 , we obtain t relax = 0.5 ps. This time scale est extremely short. We note that this time scale is (it is not a coincidence) of the order of magnitude of the normalized plateau value X (V) (•) = DX u /V 0 ' 0.5 ps in Fig. 3.9. If the displacement curve following the initial force kick was a single exponential dominated by a balance between friction and inertia, one would see a very fast asymptotic approach to the plateau value, on the same time scale as the first peak. Such a fast relaxation would indeed describe the hydrodynamic response of an incompressible rigid slab subject to solvent and interlayer friction. On the other hand, taking c t ' 1000 km.s 1 as the celerity of the transverse sound waves in the bilayer (a typical magnitude for a fluid sound wave celerity) it would take at least 5 ps for the sudden shear stress wave following the force kick to establish itself across a 5 nm thick membrane. This proves that the ideal incompressible solid relaxation result cannot describe the current situation. It also provides an alternative estimate of the characteristic time scale of the initial displacement peak position.

A Poiseuille characteristic time scale t Poiseuille can be defined as the slowest relaxation time of the Stokes hydrodynamic flow in a flat slab rL 2 w /(hp 2 ) involving the channel gap L w = 3.5 nm and the water kinematic viscosity h/r ' 7 ⇥ 10 7 m 2 .s 1 .

Its value is t Poiseuille = 1.6 ps, and also much shorter than the observed relaxation time. We therefore conclude that the sliding leaflets relaxation time t relax has a viscoelastic origin, and we denote it t vel . We therefore conclude that the relaxation seen on Fig. 3.9 results from slow membrane internal relaxation dynamics and is not limited to the interfacial sliding region. Slow lipid tilt modes relaxation, such as depicted in Fig. 3.20, certainly contribute to the observed slow viscoelastic response of the bilayer FKR.

It turns out that the condition DX u ⇠ V 0 t relax is not met. Instead DX u /V 0 is of the order of t hyd ⇠ 0.5 ps and the long relaxation time t relax ⇠ t vel DX u /V 0 enhances the effect of brownian fluctuations. Following eq. (3.22) one expects a COM diffusion coefficient of the order of 10 µm 2 .s 1 . A numerical estimate based on the COM mean squared displacement yields a value D COM,u ' 3.4 µm 2 .s 1 (Figures 3.22 and 3.23).

6.2 6.3 6.4 6.5 6.6 6.7 6.8 6.9 6.2 6.3 6.4 6.5 6.6 6.7 6.8 6.9 As a consequence, the ballistic to brownian displacement ratio r bal/br equals

r bal/br = DX u p D COM,u t vel = V 2 0 t hyd D COM,u ! 1/2 ✓ t hyd t vel ◆ 1/2 (3.25)
For an initial velocity jump V 0 = 0.09 nm.ps 1 , it gives a ratio r bal/br ' 1. According to the above expression, brownian displacement and ballistic drift are for each single run X (a) (t) of the same order of magnitude. If a number N s of independent runs is used to reduce the brownian uncertainty, the effective ballistic to brownian ratio is increased by a factor p 150 ' 12. Such a large (at least 10%) relative uncertainty is consistent with the empirical bootstrap estimates as represented in Fig. 3.17.

We finally note that that the friction coefficient b = 2.54 ⇥ 10 6 Pa.s.m 1 can be interpreted as a Newtonian fluid sheared between two infinitely thin parallel planes separated by a 5 nm thick gap, with an equivalent dynamic viscosity h equiv = 13 mPa.s.

Let us now consider the friction properties of the bilayer in the gel state. The most prominent characteristics is the absence of visible linear response regime. This is particularly clear from Figures 3.13, 3.15 and 3.18. The effective b + h/L w coefficients decreases with the external pulling stress (CPF) and the initial force kick (FKR), atypical shear-thinning behavior. As the solvent viscosity does not change at the transition, the interleaflet friction coefficient b is responsible for the observed behavior. If it is not possible to affirm for sure that no linear regime exists at lower pulling stresses, such a linear regime clearly lies beyond our current simulations capacities.

Shear thinning behavior is the hallmark of complex fluids dynamics. In the CPF regime, the effective friction b appears to follow an approximate power-law regime V u ⇠ t 1.44 , or equivalently b = t/V u ⇠ V 0.44 u . Outside linear response, one does not expect equivalence between CPF and FKR measurements in the gel phase.

The tilt relaxation dynamics (Figs 3.20 and 3.21) suggest that the lipid tilt relaxation occurs slowly in the gel phase. A possible explanation would be that irreversible or slowly reversible plastic deformations are involved in the gel sheared bilayer. Fig. 3.24 illustrates a possible mechanism. Assuming that straight lipid chains are sufficiently ordered in the gel phase, a stationary tilt angle arises as a consequence of the lateral stress. In the absence of lattice defect, Hookian (linear) elasticity is expected from the strained structure. By introducing lattice defects and rotating the structure, it is possible to tilt the lipid chains without, or with smaller elastic stress. Such defects should not be be present in bulk crystalline structures, but become much more likely in a slab of finite thickness. One expects in this case an overall softer elastic response t(q) of the bilayer which could account for the nonlinear friction response. It would also account for the remnant tilt angle following the suppression of the applied pulling stress. However, we have not yet a quantitative explanation for the apparent power-law exponent of the velocity-force characteristics. The box has a periodicity in the x and y directions. In the z direction, there is no periodicity but semi-infinite walls.

The bottom wall is polar (P4) and in the top wall is apolar (C4). These walls are represented by a 9-3 Lennard-Jones potential which has a longer range and is smoother than the 12-6 interparticle Lennard-Jones potential (shorter range and sharper).

V wall (z) = 2p#rs 3 3

⇢ 2 15 ⇣ s z ⌘ 9 ⇣ s z ⌘ 3 (4.1)
The potential V wall (z) acts only on z and preserves the xy translation invariance.

The exponent -3 corresponds to the Van der Waals interaction of a polarizable particle with a three dimensional homogeneous semi-infinite half-space [START_REF] Parsegian | Van der Waals Forces by V. Adrian Parsegian en[END_REF]. The parameters e and s are related to 12-6 Lennard-Jones potential involving the "wall bead type" P4 and a free bead in the bulk. Eq: 4.1 is nothing but the integration of the 12-6 potential with a volume density r of P4 (or C4) beads other a semi-infinite continuous space. At the bottom, there is also an explicit 0.5 nm thick coating layer composed of 717 beads of polar type (P4) providing some roughness to the bottom wall. These beads are "frozen" in the sense that they do not move irrespective of the force exerted on them. As a result, the frozen layer breaks slightly the xy translation invariance.

Above the substrate (wall and coating), one finds a layer of water molecules (P4 type) whose thickness is L w = 1.8 nm estimated from the histogram fig. 4.2, then a bilayer comprising 512 DSPC distributed into two leaflets of 256 lipids each. The DSPC molecules are identical to the ones used in the previous chapter.

At last, there is water layer showing liquid-vapor coexistence in the NVT ensemble. In experimental samples, a strong attraction maintains the bilayer in close contact with the surface. The difficulty here is to reproduce this physical interaction while the Van der Waals attraction between beads (but not walls) is truncated at 1.2 nm. When more than two species interact with a surface, the effective attraction by the wall results from differences in the mutual interaction body 1-body 2, body 1 wall, body 2 wall. . . In order to have a preferential interaction of the lipid heads with the substrate compared with water and substrate, the force field file was modified by setting the wall/water (WUP-P4) and the coating/water interaction (SSUP/P4) to "attractive", such that the water-surface (wall and coating) interaction became globally "attractive". Between the lipid heads (Q 0 , Q a ) and the substrate (SSUP/WUP), the interaction was set to "super-attractive" in order to increase the affinity of the lipid heads with the bottom surface and favor the binding of the membrane (table 4.1). If the water-substrate interaction was set to super-attractive, they would be no reason for the lipids to remain close to the membrane, and the confined water layer could grow without limit, provided that water is able to cross the bilayer. The thickness of the confined water was adjusted to match available neutron scattering data (around 2-3 nm).

Interaction of interest

Methods

The goal of this analysis is to investigate the influence of a fixed substrate on the friction dynamics of a bilayer (interleaflet friction in the presence of a substrate and confined water). Two methods relying on a pulling force strategy inspired from chapter 3 were used.

• The first method consists in pulling the center of mass of each leaflet in opposite directions in the presence of a surface. It is analogue to the Couette shear geometry in chapter 3.

• The second experiment is a bit more subtle. The center of mass of the bilayer was pulled in one direction and the substrate (though fixed) was pulled opposite to it. To be more precise, two groups of atoms were pulled relative to each other, one being the lipids, the other the coating layer of the substrate. As the coating layer is frozen (GROMACS directives freezegrps and freezedim (Y Y Y)), it did not move. The second method is therefore equivalent to applying a uniform non vanishing external force.

The most significant difference compared with the free bilayer case is that the frozen surface coating breaks the translation invariance of the system and can exert forces onto the rest of the system. While a constant external force would not be admissible in the free bilayer case, it is admitted in the present situation.

The first practical step was to create an index file. This latter contains the list of the bead indices (numbers) of the upper (DSPC TOP), lower leaflet (DSPC BOT), entire bilayer (DSPC), water (W), confined water (W BOT), free water on top (W TOP) and finally substrate (su = coating). The index file was used to select which group one wishes to pull or analyse.

Then, the directive "Pull" from GROMACS was used in order to shear with a constant pull force equal to 250 kJ.mol 1 .nm 1 . The choice is based on the linear regime found previously in the free bilayer case. The temperature was set to 340 K by using a V-rescale thermostat. The DSPC bilayer is in the fluid phase at this temperature.

The most important practical detail was the removal of the usual directive concerning the fixation of the center of mass of the system. While in the free-bilayer case it was important to keep the system COM at rest in conjunction with the use of a Nose-Hoover or V-rescale thermostat, it is important here to let the system COM free to move. One must therefore deactivate the center of mass removal (directive nstcomm =0 and comm-mode=none). To stress upon this point, a comparison between activated and deactivated center of mass removal is presented. It is safe to keep the center of mass free if the system is not translation invariant thanks to the presence of the frozen group of "su" beads. 50 NVT pull simulations of 10 ns were performed in order to accumulate statistics. To analyse the results, home made programs were written in order to read the binary trajectory file (GROMACS format trr). For instance, beads selection can be done using the index file or using the spatial localization ( z slices) and the displacement of the center of mass of each subgroup of interest computed. Moreover, one can use GROMACS tools to determine the mean square lateral displacements (MSD), based on index files.

Results

The results are based on 50 independent runs, e.g. 50 initial configurations. The first part presents the MSD results in the case of supported bilayer. Then, results of the dynamics of interleaflet pull experiments are exposed, followed by the substratebilayer pull experiments. Finally, the results are finally quantitatively interpreted.

Mean square displacements

We start by investigating the mean squared displacement of the water beads and lipids on a single NVT simulation (10 ns). The interaction with the substrate, mediated by the confined water layer is susceptible of changing the numerical value of the 2D-diffusion coefficient. Unlike the free bilayer case where complex hydrodynamical corrections were shown to be important [START_REF] Camley | Strong influence of periodic boundary conditions on lateral diffusion in lipid bilayer membranes[END_REF], the presence of the coating layer should screen all the possible long range hydrodynamic corrections. The 2D lateral diffusion should be well defined in this case. The MSD was calculated by the relation (4.2) where r(t) is the position at time t, r(0) the initial position, n the dimension, D the diffusion coefficient and Dt the time.

h(r(t) r(0))

2 i = 2nDt (4.2)
2D mean squared displacement What can be clearly seen in this figure is the difference between the upper and lower leaflet MSDs. The lower leaflet diffuses slowly, as expected due to its interaction with the confined water layer. Moreover, the confined water diffuses also much more slowly than the free water.

3D mean squared displacement

In the same manner, we investigated the 3D mean squared displacements. Results are presented in figure 4.4. Again, there is a clear difference between the upper and lower leaflets of the DSPC bilayer, and between the confined and the free water. Table (4.3) presents the diffusion coefficients (in nm 2 .ps 1 ) which is the slope of the graph representing the MSD with respect to time divided by 4 (or 6) respectively for 2D (and 3D) (eq. 4.2). The slope was estimated from the linear part of the curves situated between t = 2 and 8 ns through linear regression. As the z direction is confined for all the components, the 2D and 3D MSDs naturally follow the same linear slope. 0.0 2.0×10 The graphs 4.5 and 4.6 show respectively the MSD in 2D and 3D of the DSPC bilayer (upper and lower leaflets) and the water (below and above the membrane). The slope in 2D and 3D of the DSPC between 2 and 8 ns is clearly the same. The contribution of the vertical coordinate is simply to shift the MSD by a constant (Dz) 2 of the order of magnitude of the vertical confinement. The lipids are self-assembled into a bilayer. The lateral diffusion is simple for the lipids, they can diffuse freely between each other. However, the diffusion along the z axis (which is the normal of the membrane) is more difficult due to the organization of the membrane. The water do not have this difference of diffusion in 2D or 3D.

Inter-leaflet pull

This simulation consists in pulling two leaflets of a bilayer in opposite directions with a constant force F in the presence of a substrate. Figure 4.7 represents the displacement of the center of mass of different parts of the system. In contrast to earlier findings concerning the free bilayer and where each leaflet were symmetrically moving in opposite directions with relative velocity DV, the lower leaflet here seems to be stuck. One sees in figure 4.7 that DSPC Bot does not move significantly despite the external pulling force. In figure 4.7 there is no removal of the center of mass velocity, and the center of mass can freely move. A second experiment was performed under the same conditions but keeping the periodic (every 100 steps) center of mass velocity removal of the whole system. The results are indeed completely different in this case compared to the previous one. There, one observes the lower leaflet (DSPC Bot) moving in the direction of decreasing x, as well as the dragged confined water (Water Bot). There is no stationary leaflet in this case contrary to the free COM case. shown are the top Water, bottom Water and the motionless substrate coating as reference. On this graph, the system center of mass is free to move (no center of mass velocity removal). Each leaflet was pulled in opposite directions along the x direction with a force equal to 250 kJ.mol 1 .nm 1 . The normalized velocities of the centers of mass of each system components are presented in the two cases in table 4.4. One notices that when the COM removal is activated (fig. 4.8), the center of mass of the upper and lower leaflets moves at similar speed. Hence, as soon as one deactivates the COM removal, the lower leaflet stays at rest compared to the upper leaflet which moves with double velocity.

COM removal DSPC

Substrate-bilayer pull

The second experiment consists in pulling the whole membrane in one direction parallel to the x axis and the substrate in the opposite direction. However, the substrate has a position restrained and cannot move. The force on the substrate has no effect on it. Again, experiments with both activated and deactivated removal COM option are compared. Figure (4.9) present the results of the simulation with deactivated COM removal. It depicts the displacement of the COM of the whole bilayer as well as of each leaflet. It is also presents the displacement of the confined water (Water Bot) and free water above the membrane (Water Top). The substrate is constrained and so remains motionless. The top leaflet moves slightly faster than the bottom leaflet. The free water (Water Top) has the same velocity as the upper leaflet and the confined water (Water Bot) has roughly half the velocity of the lower leaflet (DSPC Bot).

Figure 4.10 presents the case where the COM velocity removal is active. Again, one sees the influence of the COM velocity removal in the results. The full bilayer (lower and upper leaflets) moves with the same velocity (see table 4.5). The confined water have a larger velocity than the free water in this situation.

Discussion

We describe now the hydrodynamic treatment of a stack of four layers on top of a steady, frozen surface. The hydrodynamic treatment is expected to describe accurately the stationary regime where each layer moves at constant velocity under the effect of the external applied stresses.

The system comprises a stack of five layers. The bottom surface S contains frozen beads. These beads can stand any arbitrary external force without moving. In reaction the surface layer can exert external forces to the system depending on the motion of the beads pertaining to the second layer. The second layer is the bottom water layer (BW) and its center of mass is animated with a velocity V BW . Similarly the bottom leaflet BL, top layer TL and upper water TW center of masses move with velocities V BL , V TL , V TW . The relative motion of the various layers creates a dissipation per unit of area and time, characterized by 3 coefficients g 1 , b and g 2 . The total power dissipated is proportional to

P dis = Ag 1 2 (V TW V TL ) 2 + Ab 2 (V TL V BL ) 2 + Ag 2 2 (V BL V BW ) 2 + Ag 2 2 V 2 BW . (4.3)
By symmetry, the same coefficient g 2 is used for the two layers in contact with the bottom water layer. Let us consider first the case of a couple of opposing forces (F, F) exerted on the leaflets. Dissipation balances the power injected by the externally applied stresses. Minimizing the sum of the dissipation function and the injected power FV TL + FV BL with respect to the velocities gives the following set thickness. Therefore AhV BL = FL w and V TL = (1 + h/(2bL w ))V BL .

Thus, by knowing the relative velocities of the leaflets and the thickness of the water we can calculate the viscosity of the water h and then determine the friction coefficient b from the following relations:

b = V BL ⇥ h 2L w (V TL V BL ) , (4.8)
with

h = F ⇥ L w A ⇥ V BL . (4.9)
The viscosity of the Martini water determined in the literature at 323 K is equal to 7 ⇥ 10 4 Pa.s. Our calculation gives us a viscosity for the confined water equal to h 340K = 11 ⇥ 10 4 Pa.s (eq. 4.9) with L w = 1.8 nm and an interleaflet friction coefficient b = 2.63 ⇥ 10 6 Pa.s.m 1 (eq. 4.8) with data from table 4.5. When using the data from table 4.4 and the relation V UL = F/bA from eq. (4.5), one finds b = 2.68 ⇥ 10 6 Pa.s.m 1 .

The viscosity obtained from the thin confined bilayer motion (11 ⇥ 10 4 Pa.s) is higher than the one obtained in chapter 3 (7 8 ⇥ 10 4 Pa.s). We observe that the viscosity seems to increase when the water thickness decreases. A possibility is that the water beads that interacts directly (within the cut-off distance) with external beads (lipid headgroups or substrate) contribute to make the value of the viscosity value larger, and their contribution is all the more important than the water layer becomes thinner.

Looking carefully at the velocities in table 4.5, one observes that the ratio V BW /V BL = 0.46 is slightly inferior to 1/2. This could indicate that the friction properties of the water with the lipid heads and with the coating layer are not strictly equivalent, and that the up-down symmetry of the confined water does not hold exactly. The strong interaction between water and the coating could induce a stronger friction coefficient at the lower interface of the confined water. The 9-3 Lennard-Jones bottom potential may also play a significant role in confining locally the confined water near the coating, increasing the effective viscosity and the associated friction.

The interlayer friction coefficients obtained from pulling the leaflets in the same and opposite directions, are in good agreement (2.63 ⇥ 10 6 and 2.68 ⇥ 10 6 Pa.s.m 1 ). We conclude that the hydrodynamic interpretation of the leaflet and water layers drift velocities is convincing.

Chapter 5

Artificial Flip-Flop for Small Unilamellar Vesicle equilibration "There exists everywhere a medium in things, determined by equilibrium."

Dmitri Mendeleev (1834-1907)

Introduction

T HE previous works dealt only with planar membrane samples, for which were studied the rheological properties. There are however circumstances where the lipid membrane is highly curved. It was shown experimentally in our group that silica nanoparticles were in some circumstances strongly interacting with fluid dioleoylphosphatidylcholine (DOPC) up to the point of wrapping themselves as a single bilayer around the nanoparticles.

Given the size of the nanoparticles (radius smaller than 10 nm), one conceives that the degree of curvature of the bilayer is very large, and much larger than the smallest SUV obtained by lipid tip sonication or extrusion (diameter of 50 nm). A natural question arises to determine whether the bilayer in such extreme conditions differs both thermodynamically and dynamically from flat bilayers. To some extent, the wrapped silica nanoparticles are similar to supported lipid bilayers in the sense that they are stabilized by a solid inorganic surface.

The previous section showed how a supported lipid bilayer on a silica surface can be modelled in the Martini world. It would therefore be possible to reconstitute entirely the silica NPs covered with bilayers. However, there are first important questions to address before trying to simulate this complex situation:

1. Are highly curved vesicles stable in Martini simulations ? 2. How can we equilibrate the lipid populations of the two leaflets ? 3. Does it exist a minimal number of lipid molecules required to form a stable (or metastable) closed spherical vesicle, and what is its value ?

4. Is curvature changing significantly the dynamical properties of the highly curved vesicles ?

In real vesicles, there is always some imbalance between the outer and inner leaflet populations. Since the seminal work by McConnell and Kornberg [START_REF] Mcconnell | Inside-outside transitions of phospholipids in vesicle membranes[END_REF] one knows that lipid exchange between leaflets is very slow, of the order of an hour. In simulations also, one does not see spontaneous leaflet exchange of lipid molecules despite the time acceleration of the coarse-grained molecular dynamics. On the other hand geometrical packing considerations suggest that very small vesicles count many more lipids in their outer leaflet than in their inner leaflet. As dynamical properties, such as the lipid diffusion constant, are strongly dependent on the local compression, or area per lipid of the leaflet, it is paramount to be able to balance properly the populations of the two leaflets before studying their dynamics.

This chapter deals with an attempt to perform such leaflet population balance for very small vesicles (r < 12 nm).

Composition of the system

The systems (vesicles) are composed of 650 DSPC lipids surrounded by around 27000 beads of water (depending on solvatation) in a cubic box of 16 nm lateral size. The systems differ in the initial number of lipids in each leaflet, or ratio R between the outer and inner leaflet.

Method

This method is a little bit tricky. A new type of lipid was created and added to the force field. This new type of lipid was named HLPC (Hydrophobic Lipid Phosphatidylcholine). The parameters of HLPC are the same as DSPC except that the hydrophilic headgroups in the force field were muted into hydrophobic beads of "C1" type (see fig. 2.3). Therefore the new HLPC is totally hydrophobic, only composed of beads of C1 type.

The idea is to force the flip-flop rate by alternatively changing a number of DSPC into HLPC, and then back into DSPC, and repeat the cycle many times. We call this "Activation" and "Deactivation" of the head-groups. Figure 5.1 presents each situation during one cycle.

• (A) The first step consist into selecting randomly the lipids one wants to mutate as HLPC

• (B) After a short minimisation stage, one simulates for 100 ns. As the HLPC are hydrophobic, they dive into the membrane core (in the tail hydrophobic region).

• (C) The nature of the HLPC is changed back to the original lipid (DSPC). Again, after a short minimisation, the system is simulated for 100 ns.

However, the number of mutated lipids matters significantly. If one mutates a large number of lipids in a row, the membrane can be destroyed, or during the reactivation step (C) long lived inverse micelles may form. Thus, for the equilibration part, only 1% of the total number of lipids were mutated during a cycle. In our case, this corresponds to 8 mutated lipids each time (4 lipids by leaflet). Moreover, a randomization step was important in order to avoid to always mutate the same lipids, and also to give an equal chance for each lipid to "choose" which side they want to go.

This method was first tried on a flat bilayer in order to set-up the procedure, then it was applied to vesicles (fig. 5.2). The vesicles were prepared using Packmol for the spherical vesicle template and the "solvate" tool from GROMACS in order to fill the box with water. They are localized in between the tails due to the fact that they became strongly hydrophobic.

Results

The first part of the results presents the equilibration of a vesicle made of 650 DSPC molecules. The initial number of lipids in each leaflet was fixed to 400 DSPC outside and 250 DSPC inside. Moreover, there were around 800 water beads inside. Thus the initial ratio (eq: 5.1) was 1.6 (400/250).

Ratio =

Lipid number from the outer leaflet Lipid number from the inner leaflet = OUT IN (5.1) Figure 5.3 presents the evolution of the ratio with respect to the cycles of activation and deactivation. One notice that the ratio increases. This corresponds to an increase in the number of lipids of the outer leaflet.

Figure 5.4 shows the evolution of the area per lipid during the equilibration. Again, one can see that the APL in the outer leaflet decreases while the APL of the inner leaflet increases. During the process, some lipids prefer to move outside rather than to come back inside. The area per lipid is calculated by considering the ratio of the surface composed by the inner/outer leaflet (surface of a sphere) and the total number of lipids of the leaflet of interest. The radii of the two spheres of interest are respectively around 6.0 nm and 2.7 nm. The outer APL is larger due to the high curvature of the membrane. Moreover, we investigated the number of water beads inside the vesicle during the process of equilibration. The results are shown in figure 5.5. First, the water quantity decreases, then it reaches a plateau around 380 beads.

The next step of the experiment consisted in trying to equilibrate different initial ratios, using the same method from 1:1 to 6:1 (where 1:1 is the same amount of lipid outside and inside and 6:1 is 6 times more lipids outside than inside). The hope is that those curves would converge from above and below to the optimal outer/inner lipid ratio.

Figure 5.6 presents the result of the evolution for different initial ratios. The ratio is defined a bit differently than previously by the following relation:

Ratio = OUT IN Total number of lipids (5.2)
The amount of water is approximately equal to 400 beads for all the ratios (except the initial ratio 1.6 ( around 800 water beads), 5:1 (around 1075 water beads) and 6:1 (around 1100 water beads)).

The systems 1:1 to 4:1 seem to converge to a plateau between 0.5 and 0.6 which is the probable range of stability of a bilayer of this size. The vesicles are spherical and not deformed until the ratio 6:1 (fig. 5.7). However, we notice that the equilibration occurs much faster when the ratio is small. For the systems 5:1 and 6:1, due to the large differences between the number of lipids in each leaflet, the vesicles became strongly deformed (fig. 5.7). In the case 6:1, the vesicle is deformed and has a fold, or bulge, at its side.

Discussion

The principle of the method is to deactivate and reactivate cyclically the hydrophilic character of the lipid head-groups of a few lipids. These lipids migrate to the membrane interior. It is hoped that when they are reactivated, the difference in chemical potential within the two leaflets makes one of the two choices more attractive than the other, and therefore accelerate the migration from the less favorable to the most favorable place. The mutated lipids are given the choice of which side to go. This choice is influenced by many parameters such as the number of water beads inside and the initial number of lipids. Through many cycles of activation and deactivation, the randomly selected lipids can decide somehow where they want to go.

It turns out that the ratio between the outer and inner number of lipids is influenced by the number of water beads inside the vesicle. This number co-evolve with the number of lipids. This is likely to be due to an increase in the porosity of the membrane in the course of the activation-deactivation cycle, which is associated to a flux of water.

As water is necessary to stabilize the amphiphilic assembly, a decrease in the amount of water inside the vesicle is unfavorable to the inner leaflet. A runaway cycle may start, with lipid leaving the inner leaflet, and water leaving the interior as a result, reinforcing the unbalance. Controlling the amount of water as well as the lipid leaflet imbalance seems crucial. However, after a certain number of simulation cycles, the amount of water seems to reach a stable value.

In the most extreme imbalance case, a deformation and the creation of a fold (a kind of of membrane bilayer patch grafted to the vesicle side) appears (fig. 5.7). It looks like the strong excess of lipids in the outer leaflet needs to be expelled as a small flat bilayer patch.

We therefore think that there is a range of stability for a membrane of 650 lipids. It should contain about 380-400 water beads (1500 -1600 water molecules according to the 4:1 mapping in the Martini model). The outer/inner ratio should range between 0.55 and 0.65. Beyond this ratio, an instability of the outer leaflet takes place.

We are convinced that the effective thermodynamic potential controlling the energy cost of creating a leaflet population imbalance is quite flat, and that the restoring thermodynamic forces, which would drive the excess lipids from one leaflet to the other are small. As a result, the convergence of a vesicle towards equilibrium is slow, and the equilibrium population range rather undetermined. In any case, the exchange of lipids cannot be decorrelated from the exchanged of water.

The activation-deactivation method may prove useful as an alternative to pores for fast lipid equilibration between leaflets. Even though there are still open questions regarding the stability of nanoscopic vesicles, we learnt a number of relevant facts on very small unilamellar vesicles. The scripts written for vesicle equilibration by activation-deactivation formed a very helpful starting point for the Force Kick Relaxation studies of the previous chapter Chapter 6

Conclusion and perspective

"In science one tries to tell people, in such a way as to be understood by everyone, something that no one ever knew before. But in poetry, it's the exact opposite."

Paul Dirac W E developed an original approach for studying the lipid bilayer friction through molecular dynamics simulations. In the free bilayer case, the constant pull force (CPF) method gave a friction coefficient in agreement with literature. A linear regime was found to hold, and the force kick relaxation (FKR) and constant pull force (CPF) methods gave consistent results. We clearly put into evidence a complex viscoelastic behavior, that we ascribe to lipid tilt and stretch elasticity, lasting for a time of the order of a ns. In the gel state, we discovered a non-linear power law regime, corresponding to a shear-thinning behavior, but no linear response. The apparent friction coefficient b was found to decrease as a function of the stress with an exponent -0.44 (relatively close to 1/2) when stress increases. The bilayer in the gel state is therefore subject to a slower and more complex dynamics than in the fluid state. For each regime, the validity range and the magnitude of the brownian error are reasonably well understood Our approach has the advantage that it can be generalized to supported membranes. We showed that an interacting bilayer with a substrate shows expected interesting features e.g. a difference in lipid diffusion coefficient between the upper and lower leaflets. The CPF methods provides us a reasonable friction coefficient by two independent means and the confined Martini water was found to display a slightly larger viscosity.

We also developed an original method based on an artificial lipid Flip-Flop in a vesicle of very small radius. This method is based on cycles of activation and deactivation of the polar nature of the lipid heads and accelerates the lipid exchange between leaflets, by letting them choose freely in which leaflet they prefer to go. However, we found that the exchange of lipids cannot be decorrelated from the exchange of water, and we found it difficult to determine the optimal water and lipid population content. We suggest however our method as an alternative to pores for leaflet equilibration.

Exploration of other types of lipid bilayers such as unsaturated lipids like POPC (1 -Palmitoyl-2-oleoylphosphatidylcholine) or charged lipids like DPPS (1,2-dipalmitoyl -sn-glycero-3-phospho-L-serine) with our methods could bring new element to the understanding of lipid coating rheology.

The analysis of a tri-layer DSPC membrane will be one of our near future objective. Indeed, it is inspired by local experimental work in which a tri-layer membrane of DSPC is prepared with the help of a Langmuir trough, and then studied Chapter 6. Conclusion and perspective with a FRAPP (Fluorescence recovery after patterned photo-bleaching) technique. Therefore, our simulations will provide the possibility of investigating this specific system in which one can pull a monolayer and observe the behavior of the others in the presence of a substrate.

Regarding the Flip-Flop method, there are also perspectives. First, one could improve the time efficiency of the simulation. Indeed, 200 ns for a cycle is quite long, maybe 10 ns per step could be enough. The second perspective concerns the exploration of increasing size of vesicles in order to see how the equilibrium ratio between the number of lipids in each leaflet changes as the total number of lipids increases. Conversely, one might attempt to obtain the smallest possible vesicle in this way.

We would like to determine the evolution of the lipids angular mean-squared displacement and the angular diffusion coefficient of an equilibrated lipid vesicle. Finally it would be interesting to estimate the friction in the case of curved lipid membranes. This could be achieved by initiating a rotational motion of the vesicle and monitoring its relaxation to equilibrium, by analogy with force kick relaxation experiments. i j funct length force.c. Center of mass (COM) fixation (nstcomm) deserves a special attention. It is required to fix the system COM to a constant position as soon as the system in invariant by translation conjugated with the use of a Nose-Hoover or V-rescale thermostat. In the Couette flow situation, the bilayer and water groups have a separately vanishing linear momentum. In the Poiseuille flow, only the system COM is stationary. One must therefore apply the constraint on the system center of mass (which would otherwise not be perfectly steady due to the approximate treatment of intermolecular forces), and not separately to the subsystems. Constant force pulling was implemented using umbrella sampling control parameters, such as in the following example for pulling in the Couette geometry with a constant force of 250 kJ.mol 1 .nm 1 :

Pull =yes pull_ngroups = 2 pull_ncoords = 1 pull_group1_name = up pull_group2_name = down pull_coord1_type = constant-force pull-coord1-vec = 1 0 0 pull_coord1_geometry = direction-periodic pull_coord1_groups = 1 2 pull_coord1_dim = Y N N pull_coord1_k = 250 pull_coord1_start = yes

In the CPF analysis, 50 trajectories of 10 ns were used and combined for each pulling stress condition.

Force-kick relaxation simulations were realised by changing with a python script the x components of the velocities in the initial configuration file (gro when using Gromacs 5.1) as suggested in eq. 3.23 and using the new velocities as a starting configuration. The FKR relaxation dynamics is unusual in terms of short characteristic relaxation times, of the order of 1 ps. To perform our statistical analysis, trajectory frames were dumped every 10 time-steps (0.2 ps) and 150 trajectories of 25000 steps (500 ps) were used and combined for each initial velocity condition.

Tables B ; MARTINI and TEMPERATURE/PRESSURE ; normal temperature and pressure coupling schemes can be used.

; It is recommended to couple individual groups in your system separately.

; Good temperature control can be achieved with the velocity rescale (V-rescale) ; thermostat using a coupling constant of the order of 1 ps. Even better ; temperature control can be achieved by reducing the temperature coupling ; constant to 0.1 ps, although with such tight coupling (approaching ; the time step) one can no longer speak of a weak-coupling scheme.

; We therefore recommend a coupling time constant of at least 0.5 ps.

; The Berendsen thermostat is less suited since it does not give ; a well described thermodynamic ensemble. ; ; Pressure can be controlled with the Parrinello-Rahman barostat, ; with a coupling constant in the range 4-8 ps and typical compressibility ; in the order of 10e-4 -10e-5 bar-1. Note that, for equilibration purposes, ; the Berendsen barostat probably gives better results, as the Parrinello-; Rahman is prone to oscillating behaviour. For bilayer systems the pressure ; coupling should be done semiisotropic. Notre contribution a consisté à développer une approche numérique appliquée à des systèmes modèles (bicouches lipidiques avec et sans substrat), visant à déterminer la nature de la friction inter-feuillet dans le cas d'une phase fluide et d'une phase gel (qui dans le contexte de la biophysique des membranes désigne un état de chaînes ordonnées de nature très visqueuse). Notre approche théorique découle directement de l'expérience « NanoTribo-FRAPP » développée au sein de notre laboratoire (ICS), qui permet d'appliquer un cisaillement à des multicouches lipides déposés sur lame de verre, tout en mesurant la vitesse d'écoulement locale de chaque couche moléculaire [START_REF] Falk | Interleaflet sliding in lipidic bilayers under shear flow: comparison of the gel and fluid phases using reversed non-equilibrium molecular dynamics simulations[END_REF].

Avant tout, le chapitre 1 intitulé état de l'art présente au lecteur le monde des lipides en donnant un aperçu de ce que sont les lipides et pourquoi ils sont si importants dans les systèmes du vivant. Par ailleurs, une introduction à la dynamique moléculaire est donnée. Enfin, quelques résultats importants sur le domaine sont présentés dans ce chapitre.

Le chapitre 2 concerne la méthodologie et décrit comment les systèmes sont préparés et analysés en donnant les détails spécifiques et originaux de notre approche, outils d'analyse et étude statistique de bootstrap [START_REF] Press | Numerical Recipes in C (2Nd Ed.)[END_REF].

Le chapitre 3 décrit les différents résultats obtenus pour une bicouche lipidique libre ordinaire soumise à différents types de sollicitations. Ensuite, le chapitre 4 traite du cas d'une membrane supportée. Dans cette situation, les approches de simulation du chapitre 3 sont adaptées pour tenir compte de la présence d'un substrat.

Le chapitre 5 est un bref compte rendu d'une méthode que nous avons développée afin d'équilibrer de très petites vésicules. Cette méthode est certainement prometteuse en ce qui concerne l'équilibrage et la relaxation de tels systèmes.

Enfin le chapitre 6 conclut notre travail et dévoile des perspectives générales. Notre étude porte sur les propriétés rhéologiques d'une bicouche de DSPC (1,2distearoyl-sn-glycerol-3-phosphocholine) formée de 512 molécules de DSPC entourée de 5120 billes d'eau. Le champ de force est une description à gros grains « Martini » du lipide [START_REF] Marrink | The MARTINI Force Field: Coarse Grained Model for Biomolecular Simulations[END_REF]. Les systèmes sont simulés à des températures de 340K et 280K afin d'obtenir respectivement un état fluide et un état gel, sachant que la température de coexistence Tm du DSPC est de 328K.

Nous simulons à l'aide du moteur de dynamique moléculaire Gromacs [START_REF] Bekker | Gromacs: A parallel computer for molecular dynamics simulations in[END_REF] une bicouche lipide subissant une contrainte appliquée de différentes manières, et observons la réponse de la membrane en fonction de l'état (fluide ou gel) dans laquelle elle se trouve. Nous définissons un protocole de simulation compatible avec les possibilités de l'environnement de simulation Gromacs. La première méthode (« Constant Pull Force CPF ») consiste à appliquer un couple de forces constantes à deux sous-systèmes, et induire par conséquent une déplacement relatif stationnaire de ces deux sous-groupes. La seconde méthode (« Force Kick Relaxation FKR ») consiste à préparer les deux sous-systèmes avec une vitesse relative finie opposée, de façon à garantir l'immobilité du centre de masse du système complet, et à suivre au cours du temps la relaxation mécanique de ces systèmes vers l'équilibre.

Les deux méthodes sont potentiellement équivalentes lorsqu'un régime de réponse linéaire existe. La première approche, typique des méthodes de simulations hors équilibre, place le système dans un régime stationnaire légèrement dissipatif, et donne les coefficients frictions recherchés directement. La seconde méthode permet de découvrir la réponse transitoire, potentiellement viscoélastique du système. Cette approche originale est transposable aux systèmes de bicouches supportées en contact avec des substrats solides lisses (fig. Les résultats moyennés sur 50 simulations sont en accord avec la littérature et donne un coefficient de friction proche de 2,7⇥10 6 Pa.s.m 1 [START_REF] Otter | Intermonolayer Friction and Surface Shear Viscosity of Lipid Bilayer Membranes[END_REF] pour la phase fluide. Dans la situation de Poiseuille où la bicouche est tirée d'un côté, alors que l'eau est envoyée de l'autre côté (Fig. A.3), il est possible de déterminer la viscosité de l'eau « Martini » ; Celle-ci vaut environ 7.01⇥10 4 Pa.s pour des forces comprises entre 150 et 250 kJ.mol 1 nm 1 , en accord avec les résultats publiés [START_REF] Falk | Interleaflet sliding in lipidic bilayers under shear flow: comparison of the gel and fluid phases using reversed non-equilibrium molecular dynamics simulations[END_REF].

A.3 Conclusion

Nous avons développé une méthodologie originale pour l'étude de la friction des bicouches lipides. Notre approche a le mérite de pouvoir se généraliser à des membranes supportées. Nous avons clairement mis en évidence un régime viscoélastique que nous attribuons à l'élasticité d'inclinaison des lipides. Dans l'état gel, nous avons découvert un régime non linéaire, correspondant à une rhéo-fluidification. Le coefficient de friction apparent b semble diminuer avec un exposant -0,48 (proche de ½) quand la force augmente. La bicouche dans l'état gel est donc sujette à une dynamique lente et plus complexe que dans l'état fluide.

Nous avons également développé une méthode originale basée sur des "Flip-Flop" artificiels dans des vésicules de très petit rayon. Cette méthode est basée sur des cycles d'activation et de désactivation de la nature polaire des têtes lipidiques et accélère ainsi l'échange lipidique entre les feuillets, en leur laissant choisir librement dans quelle feuillet ils préfèrent aller. Cependant, nous avons constaté que l'échange de lipides ne peut pas être entiérement décorrélé de l'échange d'eau. En effet, le ratio optimal entre lipide et molecule d'eau reste difficile à determiner avec précision.

Rhéologie des couches de lipides à l'échelle moléculaire Résumé

Les phospholipides ont sans aucun doute un rôle capital dans les propriétés de lubrification de biosystèmes. L'exemple le plus probant est le cas des contacts biologiques tels que les articulations ; celles-ci présentent un coefficient de friction très faible. La compréhension de tels systèmes biophysiques présente des enjeux aussi bien fondamentaux qu'appliqués, en lien par exemple avec l'arthrose et la confection de prothèses articulaires. Nous simulons à l'aide du moteur de dynamique moléculaire GROMACS une bicouche de lipide (DSPC) subissant une contrainte appliquée de différentes manières, et observons la réponse de la membrane en fonction de l'état (fluide ou gel) dans laquelle elle se trouve. La première méthode (« CPF : constant pull force ») consiste à appliquer un couple de forces constantes à deux sous-systèmes, et induire par conséquent une déplacement relatif stationnaire de ces deux sous-groupes. La seconde méthode («FKR : force kick relaxation ») consiste à préparer les deux sous-systèmes avec une vitesse relative finie opposée, de façon à garantir l'immobilité du centre de masse du système complet, et à suivre au cours du temps la relaxation mécanique de ces systèmes vers l'équilibre.

Nous avons clairement mis en évidence un régime viscoélastique que nous attribuons à l'élasticité d'inclinaison des lipides. Dans l'état gel, nous avons découvert un régime non linéaire, correspondant à une rhéo-fluidification. Le coefficient de friction apparent b semble diminuer quand la force augmente. La bicouche dans l'état gel est donc sujette à une dynamique lente et plus complexe que dans l'état fluide. Par ailleurs notre approche a le mérite de se généraliser aux membranes supportées, pour lesquelles nous obtenons des résultats pour la diffusion et la friction des différentes couches.

Mots clefs : lipides, bicouche, friction, simulation numérique, GROMACS, MARTINI, membranes supportées

Résumé en anglais

Phospholipids undoubtedly have a crucial role in the lubrication properties of biosystems. The most convincing example is the case of biological contacts such as human joints, which have a very low coefficient of friction. The comprehension of such biophysical systems implicates both fundamental as well as applied challenges, for instance in regard to osteoarthritis and the manufacture of joint prostheses. Using the molecular dynamics engine GROMACS, we simulate a lipid bilayer (DSPC) undergoing stress applied in different ways and observe the response of the membrane depending on the physical state (fluid or gel). The first method (CPF: constant pull force) consists in applying a pair of constant forces to two subsystems, and therefore induce a stationary relative displacement of these two subgroups. The second method (FKR: force kick relaxation) consists in preparing the two subsystems with an opposite finite relative velocity, so as to guarantee the immobility of the center of mass of the complete system, and to follow over time the mechanical relaxation of these systems towards equilibrium.

The results clearly show the presence of a viscoelastic regime that we attribute to the elasticity of lipids tilt. In the gel state, we observe a nonlinear regime, corresponding to a shear thinning. The apparent friction coefficient b tends to decrease when the force increases. Therefore, the bilayer in the gel state is subject to a slow and more complex dynamics than in the fluid state. Moreover, our approach can be generalized to supported bilayers for which we obtain results on the diffusion and friction of the different layers.

Keyword: lipids, bilayer, friction, numerical simulation, GROMACS, MARTINI, supported membranes

  Scheme of a synovial joint. (B) Scheme of a contact point with confined lipid bilayers, stacked in layers and vesicles. [4] . . . . . . . . 1.1 Representation of fatty acids. They are composed by a -COOH carboxyl function attached to a hydrocarbon chain. The cis unsaturations induce a 120 kink of the axis of the oligomeric chain. . . . . . . . . . . 1.2 Lipids classification and examples of molecules [8]. . . . . . . . . . . . 1.3 A popular cell membrane representation, according to the fluid mosaic picture [14]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4 Schematic representation of a cell [20] with organelles delimited or made of lipid membranes. . . . . . . . . . . . . . . . . . . . . . . . . . . 1.5 Molecular structure of two important phospholipid classes: glycerophospholipids, built around a glycerol with two variable fatty acids, and sphingolipids (built around a sphingosine with a variable fatty acid chain). The phosphate group is common to both. . . . . . . . . . . . . . 1.6 Representation of a DSPC lipid with its molecular structure (left) and a cartoon depicting the hydrophilic head and the two hydrophobic tail groups. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.7 Example of phospholipid organisations in aqueous media . . . . . . . 1.8 DSPC molecule formed by Stearoyl (C18) aliphatic chains and phosphocholine head. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.9 Cryo-TEM pictures of DPPC vesicles in fluid phase (at 46 C) and gel phase (at room temperature RT). The picture was provided through the courtesy of Dr. Marc Schmutz. The preparation protocol is exposed in ref [31] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.10 Representation of two configurations of a bilayer (only upper leaflet) from the work of P. van der Ploeg and H. J. C. Berendsen [40]. . . . . . 1.11 Diagram depicting the leap-frog scheme. Time difference between position and velocity update is dt/2. . . . . . . . . . . . . . . . . . . . . . 1.12 Periodic boundary conditions (PBC) in 2D[45] . . . . . . . . . . . . . . 1.13 Bond stretching . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.14 Angle bending . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.15 Bond rotation or torsion . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.16 Representation of a Lennard-Jones potential . . . . . . . . . . . . . . . . 1.17 Different representation of lipids from all-atom (left) to coarse grained (right). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.18 Draw of a membrane bilayer under normal shear stress . . . . . . . . . 1.19 Shear solvent method by out of equilibrium dynamics [59] . . . . . . . 2.1 Logo of GROMACS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Logo of Martini force field . . . . . . . . . . . . . . . . . . . . . . . . . .

3 . 10 .

 310 The horizontal dashed line stands for the CPF value b + h/L w = 2.75 ⇥ 10 6 . . . . . . . . . . . . . . . . . . 3.18 Apparent friction coefficient b + h/L w from impulsion simulations as a function of the initial induced velocity in the fluid and gel phases. The nonlinear shear thinning is much stronger in the gel phase. . . . . 3.19 Reference tilt angle hqi at equilibrium. . . . . . . . . . . . . . . . . . . . 3.20 Evolution of the average tilt angle hqi during a constant pull force experiment with a force F = 250 kJ.mol 1 .nm 1 (stress t = 24 bars). . 3.21 Average tilt angle curves hqi during an impulsion experiment (V 0 = 0.09 nm.ps 1 in the gel and fluid states respectively). . . . . . . . . . . xii 3.22 Brownian xy trajectory of the upper leaflet center of mass, observed during 50 ns in the fluid phase. . . . . . . . . . . . . . . . . . . . . . . . 3.23 Mean square displacement curve of the upper leaflet center of mass at 340 K and a linear adjustment with 2D = 6.8 ⇥ 10 5 nm 2 .ps 1 . . . . . . 3.24 (a) Schematic representation of unstrained lipid chains ordered phase. (b) As a consequence of sliding friction, strain builds up in the gel phase. (c) strain can be decreased by shifting the crystalline network (pairs of edge dislocations) and rotating the resulting "staircase" shape. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1 Snapshot of a simulation of a supported membrane (in collaboration with Claire Loison ILM, Lyon). From bottom to top: substrate, thin layer of water, bilayer, water, water vapor. . . . . . . . . . . . . . . . . . 4.2 Density profiles of the substrate, lipid and water in the simulation box. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Lateral mean squared displacement of a supported DSPC bilayer. . . . 4.4 Volume diffusion of a supported DSPC bilayer in water (3D) . . . . . . 4.5 Comparison of the 2D and 3D lipid MSD in a supported DSPC bilayer. 4.6 Comparison of the 2D and 3D water MSD in a supported DSPC bilayer. 4.7 Normalized averaged displacements of the center of mass of the dif-

  . . . . . . . . . . . . . . . 4.11 Schematic of the different layers of the supported bilayer. From the bottom to the top : S:Substrate ; BW:Bottom Water ; BL:Bottom Layer ; TP:Top Layer and TW:Top Water. . . . . . . . . . . . . . . . . . . . . . 5.1 Procedure of activation/deactivation of the hydrophilic head-groups. (A) one randomly picks up a certain number of lipids, (B) one changes their nature into HLPC, (C) One proceeds to the reactivation of the true nature of the lipids. . . . . . . . . . . . . . . . . . . . . . . . . . . . xiii 5.2 Half-vesicle with 8 HLPC lipids with deactivated hydrophilic groups. The beads in red correspond to the mutated lipids. They are localized in between the tails due to the fact that they became strongly hydrophobic. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3 Evolution of the ratio between the number of lipids from the outer and inner leaflet with respect to the number of simulations cycles of activation and deactivation. The initial ratio corresponds to 400 DSPC molecules outside for 250 DSPC molecules inside. . . . . . . . . . . . . 5.4 Evolution of the area per lipid (of the outer and inner leaflet) with respect to the cycles of simulations. . . . . . . . . . . . . . . . . . . . . . 5.5 Evolution of the number of water beads inside the vesicle during the process of equilibration. . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.6 Evolution of the ratio between the inner and outer lipids of a vesicle during the process of equilibration for different initial ratios. . . . . . 5.7 Snapshot of a slice of a deformed vesicle from an initial ratio 6:1 after more than 100 equilibration cycles. . . . . . . . . . . . . . . . . . . . . . A.1 Examples of lipid mapping to CG Martini beads [64]. The types Q a , Q 0 , N a and C 1 are identical as for the mapping of DSPC. The DSPC has one more bead of type C1 in the tails. . . . . . . . . . . . . . . . . . B.1 Martini CG representation of a DSPC molecule. Bead 1: choline, bead 2: phosphate, beads 3-4: glycerol, beads: 5-9 and 10-14 hydrophobic chains. A vector linking the first and last carbons of each chain is used for defining the average lipid tilt angle q. . . . . . . . . . . . . . . B.2 Evolution of the horizontal and vertical box sizes during the NPT simulation, used for determining the average box size, in the fluid phase. /boxFluid . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B.3 Evolution of the horizontal and vertical box sizes during the NPT simulation, used for determining the average box size, in the gel phase. /boxGel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . D.1 Displacement of COM of one leaflet vs time in fluid phase for respectively simulation number 1, 2, 3, 42, 43 and 44 . . . . . . . . . . . . . . . D.2 Displacement of COM of one leaflet vs time in gel phase for respectively simulation number 1, 2, 3, 42, 43 and 44 . . . . . . . . . . . . . . . D.3 Displacement of COM of one leaflet vs time in fluid phase for respectively simulation number 1, 2, 3, 42, 143, 144 and 150 . . . . . . . . . . . D.4 Displacement of COM of one leaflet vs time in gel phase for respectively simulation number 1, 2, 3, 42, 143, 144 and 150 . . . . . . . . . . . A.1 (A) Schéma d'une articulation synoviale. (B) Schéma d'un point de contact avec des bicouches lipidique confinées, assemblée en piles et vésicule. [4] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.2 Modélisation d'une bicouche de DSPC dans de l'eau . . . . . . . . . . . A.3 Différents types de cisaillement . . . . . . . . . . . . . . . . . . . . . . . xiv A.4 Déplacements moyens normalisés (feuillet supérieur) X u (t)/F pour un ensemble de forces de traction croissantes 10, . . . , 2000 kJ.mol 1 .nm (de manière équivalente en contrainte t = 0.1 . . . 190 bars). Le déplacement pour F = 2000 se situe clairement au-delà du régime linéaire et la force F = 10 est en concurrence avec l'agitation thermique. Encart: Gros plan sur le déplacement normalisé pour les contraintes supérieures à 0,955 bar. . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.5 Dérive des vitesses moyennes en fonction des forces appliquées (barre d'échelle horizontale inférieure) ou contraintes (barre d'échelle horizontale supérieure

FIGURE 1 -

 1 FIGURE 1 -(A) Scheme of a synovial joint. (B) Scheme of a contact point with confined lipid bilayers, stacked in layers and vesicles.[START_REF] Falk | Interleaflet sliding in lipidic bilayers under shear flow: comparison of the gel and fluid phases using reversed non-equilibrium molecular dynamics simulations[END_REF]

FIGURE 1 . 1 -

 11 FIGURE 1.1 -Representation of fatty acids. They are composed by a -COOH carboxyl function attached to a hydrocarbon chain. The cis unsaturations induce a 120 kink of the axis of the oligomeric chain.

FIGURE 1 . 2 -

 12 FIGURE 1.2 -Lipids classification and examples of molecules [8].

FIGURE 1 . 4 -

 14 FIGURE 1.4 -Schematic representation of a cell [20] with organelles delimited or made of lipid membranes.

FIGURE 1 . 6 -

 16 FIGURE 1.6 -Representation of a DSPC lipid with its molecular structure (left) and a cartoon depicting the hydrophilic head and the two hydrophobic tail groups.

FIGURE 1 . 7 -

 17 FIGURE 1.7 -Example of phospholipid organisations in aqueous media

FIGURE 1 . 8 -

 18 FIGURE 1.8 -DSPC molecule formed by Stearoyl (C18) aliphatic chains and phosphocholine head.

FIGURE 1 . 10 -

 110 FIGURE 1.10 -Representation of two configurations of a bilayer (only upper leaflet) from the work of P. van der Ploeg and H. J. C. Berendsen [40].

FIGURE 1 . 11 -

 111 FIGURE 1.11 -Diagram depicting the leap-frog scheme. Time difference between position and velocity update is dt/2.

FIGURE 1 . 12 -

 112 FIGURE 1.12 -Periodic boundary conditions (PBC) in 2D[START_REF] Van Der Spoel | GROMACS: fast, flexible, and free[END_REF]

FIGURE 1 .

 1 FIGURE 1.13 -Bond stretching

FIGURE 1 .

 1 FIGURE 1.14 -Angle bending

FIGURE 1 . 17 -

 117 FIGURE 1.17 -Different representation of lipids from all-atom (left) to coarse grained (right).

FIGURE 1 . 18 -

 118 FIGURE 1.18 -Draw of a membrane bilayer under normal shear stress

FIGURE 1 . 19 -

 119 FIGURE 1.19 -Shear solvent method by out of equilibrium dynamics[START_REF] Otter | Intermonolayer Friction and Surface Shear Viscosity of Lipid Bilayer Membranes[END_REF]

Chapter 2 Methodology"

 2 Science is what scientists do, and there are as many scientific methods as there are individual scientists."Percy Williams Bridgman

FIGURE 2 . 1 -

 21 FIGURE 2.1 -Logo of GROMACS The version 5.1.4 of GROMACS was used for this work. Many other and more recent versions are available, which differ little from the one we used. We decided for simplicity and stability to keep the same version of the software during the thesis.Due to popularity, there are many available tutorials dedicated to GROMACS simulations. A nice one can be found on the following web site: (http://www.mdtutorials. com/gmx/)[START_REF] Lemkul | From Proteins to Perturbed Hamiltonians: A Suite of Tutorials for the GROMACS-2018 Molecular Simulation Package [Article v1.0[END_REF]. The list of strong points in favor of GROMACS are the following ones:

FIGURE 2 . 2 -

 22 FIGURE 2.2 -Logo of Martini force field

FIGURE 2 . 3 -

 23 FIGURE 2.3 -Examples of lipid mapping to CG Martini beads [64]. The types Q a , Q 0 , N a and C 1 are identical as for the mapping of DSPC.

FIGURE 2 . 4 -

 24 FIGURE 2.4 -Topology file: it includes the force field informations, the structures and total number of molecules

FIGURE 2 . 5 -

 25 FIGURE 2.5 -Example of a plot representing the evolution of the potential energy during a minimization

FIGURE 2 . 6 -

 26 FIGURE 2.6 -Order parameter of our membrane bilayer in the gel and fluid phase.

FIGURE 2 . 7 -

 27 FIGURE 2.7 -Snapshot of a coarse-grained lipid DSPC bilayer in water with VMD (VDW representation style).

FIGURE 2 . 8 -

 28 FIGURE 2.8 -Snapshot of a bonded lipid DSPC bilayer in water using the cg_bonds.tcl script.

  i = 1 . . . N, a = 1 . . . M by drawing with repetition, at random, N elements of S (0) . The variability of the average

µb = 1 )µ

 1 (t) of given procedure (CPF, FKR . . . with different input parameters) is taken as working sample. Prior to analysing, a collection of weight vectors w (a) b ; b = 1 . . . M; a = 1 . . . N s was drawn at random, where for each givenb, N s independent draws of integers I 2 [1, N s ] were performed and w (I) (b) was set equal to the number of times I was drawn (with repetition) during the process, and divided by N s . In this way w . The flat sample average corresponds to the special vector w 0 (a) = 1/N s . Each bootstrap realisation corresponds e.g. to a contraction X b,µ (t) = (t) of the working sample. Each function X b,µ (t) represents a randomly resampled average of the original working sample, close to the flat average

FIGURE 2 . 9 -

 29 FIGURE 2.9 -Example of a set of 10 synthetic graphs from a bootstrap procedure for the upper and lower leaflet.

  which we intend to simulate in a near future. Our simulated systems comprise a single bilayer alongside a single water slab, with periodic boundary conditions in the three dimensions. Two representative snapshots are shown in Figures. 3.1 and 3.2.

FIGURE 3 . 1 -

 31 FIGURE 3.1 -Snapshot of a configuration of a coarse-grained bilayer containing 256 DSPC lipids per leaflet, with 2560 water beads molecules on both sides, in the high temperature fluid state at 340 K.

FIGURE 3 . 2 -

 32 FIGURE 3.2 -Snapshot of a configuration in the low temperature gel state at 280 K. Compared to the fluid case, the bilayer is less extended in the xy direction and thicker. No appreciable lipid chain tilt angle is visible.

FIGURE 3 . 3 -

 33 FIGURE 3.3 -Geometric parameterization of the system used in the present study, with L = L b + L w .

FIGURE 3 . 5 -

 35 FIGURE 3.5 -Constant force pulling experiments in the fluid state:(A) single leaflet COM displacement X(t)(a) , starting from the simulation box center (ca 6.6 nm) and (B) averaged displacement hXi(t) ' 1/50 P 50 a=1 X (a) (t) vs time. A bootstrap procedure (C) estimates the dispersion s b (X u (t)) caused by the finiteness of the sample {a}. Vertical bars represent the confidence interval of 10 selected points from the second half of the trajectory (5000 < t < 10000 ps), taken as twice the estimated bootstrap standard deviation. The vertical bars are used to provide a confidence interval for the drift velocity (slope of the averaged displacement curve). /TrajectoryAndBootstrapFluid-Pull

FIGURE 3 . 15 -

 315 FIGURE 3.15 -Velocity V u -applied shear stress t characteristics in the gel and fluid state. If the fluid state display a linear characteristics up to t=100 bars (linear regime) the characteristics seems nowhere linear in the gel state.

FIGURE 3 . 18 -

 318 FIGURE 3.18 -Apparent friction coefficient b + h/L w from impulsion simulations as a function of the initial induced velocity in the fluid and gel phases. The nonlinear shear thinning is much stronger in the gel phase.

FIGURE 3 . 19 -FIGURE 3 . 20 -

 319320 FIGURE 3.19 -Reference tilt angle hqi at equilibrium.

FIGURE 3 . 21 -

 321 FIGURE 3.21 -Average tilt angle curves hqi during an impulsion experiment (V 0 = 0.09 nm.ps 1 in the gel and fluid states respectively).

FIGURE 3 . 22 -

 322 FIGURE 3.22 -Brownian xy trajectory of the upper leaflet center of mass, observed during 50 ns in the fluid phase.

2 )FIGURE 3 . 23 -

 2323 FIGURE 3.23 -Mean square displacement curve of the upper leaflet center of mass at 340 K and a linear adjustment with 2D = 6.8 ⇥ 10 5 nm 2 .ps 1 .

FIGURE 3 .

 3 FIGURE 3.24 -(a) Schematic representation of unstrained lipid chains ordered phase. (b) As a consequence of sliding friction, strain builds up in the gel phase. (c) strain can be decreased by shifting the crystalline network (pairs of edge dislocations) and rotating the resulting "staircase" shape.

Figure 4 .

 4 2 represents the density profile of the various components in the simulation box.

FIGURE 4 . 2 -

 42 FIGURE 4.2 -Density profiles of the substrate, lipid and water in the simulation box.

Figure 4 .

 4 Figure 4.3 shows the MSD of four different subsets of the system (confined water, free water, upper and lower leaflet of the bilayer).

FIGURE 4 . 3 -

 43 FIGURE 4.3 -Lateral mean squared displacement of a supported DSPC bilayer.

FIGURE 4 . 5 -

 45 FIGURE 4.5 -Comparison of the 2D and 3D lipid MSD in a supported DSPC bilayer.

FIGURE 4 . 6 -

 46 FIGURE 4.6 -Comparison of the 2D and 3D water MSD in a supported DSPC bilayer.

FIGURE 4 . 7 -

 47 FIGURE 4.7 -Normalized averaged displacements of the center of mass of the different parts of the system (Top/Bot leaflets). Also shown are the top Water, bottom Water and the motionless substrate coating as reference. On this graph, the system center of mass is free to move (no center of mass velocity removal). Each leaflet was pulled in opposite directions along the x direction with a force equal to 250 kJ.mol 1 .nm 1 .

FIGURE 4 . 8 -

 48 FIGURE 4.8 -Normalized averaged displacement of the center of mass of the different parts of the system (top and bottom leaflets, free and confined water above and frozen substrate). The removal of center of mass is activated. Each leaflet was pulled in the opposite direction along the x axis with a force equal to 250 kJ.mol1 .nm1

FIGURE 4 . 11 -

 411 FIGURE 4.11 -Schematic of the different layers of the supported bilayer. From the bottom to the top : S:Substrate ; BW:Bottom Water ; BL:Bottom Layer ; TP:Top Layer and TW:Top Water.

FIGURE 5 . 1 -

 51 FIGURE 5.1 -Procedure of activation/deactivation of the hydrophilic head-groups. (A) one randomly picks up a certain number of lipids, (B) one changes their nature into HLPC, (C) One proceeds to the reactivation of the true nature of the lipids.

FIGURE 5 . 2 -

 52 FIGURE 5.2 -Half-vesicle with 8 HLPC lipids with deactivated hydrophilic groups. The beads in red correspond to the mutated lipids.They are localized in between the tails due to the fact that they became strongly hydrophobic.

4FIGURE 5 . 3 -

 53 FIGURE 5.3 -Evolution of the ratio between the number of lipids from the outer and inner leaflet with respect to the number of simulations cycles of activation and deactivation. The initial ratio corresponds to 400 DSPC molecules outside for 250 DSPC molecules inside.

FIGURE 5 . 4 -

 54 FIGURE 5.4 -Evolution of the area per lipid (of the outer and inner leaflet) with respect to the cycles of simulations.

FIGURE 5 . 5 -

 55 FIGURE 5.5 -Evolution of the number of water beads inside the vesicle during the process of equilibration.

FIGURE 5 . 6 -FIGURE 5 . 7 -

 5657 FIGURE 5.6 -Evolution of the ratio between the inner and outer lipids of a vesicle during the process of equilibration for different initial ratios.

FIGURE B. 2 -FIGURE B. 3 -

 23 FIGURE B.2 -Evolution of the horizontal and vertical box sizes during the NPT simulation, used for determining the average box size, in the fluid phase. /boxFluid

FIGURE D. 2 - 44 D. 2 FKRFIGURE D. 3 -FIGURE D. 4 -

 244234 FIGURE D.2 -Displacement of COM of one leaflet vs time in gel phase for respectively simulation number 1, 2, 3, 42, 43 and 44

  A.11 ci-dessous). Lorsque la paire de force est appliquée sur les deux feuillets opposés (Fig.A.3 a) l'écoulement résultant dans le solvant est un écoulement de Couette plan (Fig.A.3 b). Les simulations nous donnent alors accès au coefficient de friction inter-feuillet. Lorsque la paire de forces est appliquée respectivement à l'eau et aux lipides (Fig.A.3 c) l'écoulement résultant est un écoulement de Poiseuille (Fig.A.3 d), et l'information porte sur la viscosité du solvant.L'analyse des résultats se fait à l'aide de nombreux scripts et programmes (python ou C) développés spécifiquement pour nos simulations, et permettant une analyse de type « bootstrap » (ré-échantillonnage) pour estimer l'intervalle de confiance des quantités mesurées.

FIGURE A. 2 -FIGURE A. 5 -

 25 FIGURE A.2 -Modélisation d'une bicouche de DSPC dans de l'eau

FIGURE A. 11 -

 11 FIGURE A.11 -Capture d'une simulation d'une membrane supportée (collaboration Claire Loison, ILM Lyon). De bas en haut : le substrat, une couche d'eau, une bicouche, une couche d'eau, de la vapeur d'eau

  2.4 Topology file: it includes the force field informations, the structures and total number of molecules . . . . . . . . . . . . . . . . . . . . . . . 2.5 Example of a plot representing the evolution of the potential energy during a minimization . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.6 Order parameter of our membrane bilayer in the gel and fluid phase. . 2.7 Snapshot of a coarse-grained lipid DSPC bilayer in water with VMD (VDW representation style). . . . . . . . . . . . . . . . . . . . . . . . . . 2.8 Snapshot of a bonded lipid DSPC bilayer in water using the cg_bonds.tcl script. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.9 Example of a set of 10 synthetic graphs from a bootstrap procedure for the upper and lower leaflet. . . . . . . . . . . . . . . . . . . . . . . . 3.1 Snapshot of a configuration of a coarse-grained bilayer containing 256 DSPC lipids per leaflet, with 2560 water beads molecules on both sides, in the high temperature fluid state at 340 K. . . . . . . . . . . . . 3.2 Snapshot of a configuration in the low temperature gel state at 280 K. Compared to the fluid case, the bilayer is less extended in the xy direction and thicker. No appreciable lipid chain tilt angle is visible. . . . 3.3 Geometric parameterization of the system used in the present study, with L = L b + L w . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.

Membranes FF Force Field FKR Force Kick Relaxation FRAPP Fluorescence Recovery After Patterned Photo-bleaching GPU Graphics Processing Unit GROMACS GROningen MAchine Chemical Simulations GUV Giant Unilamellar Vesicles HLPC Hydrophobic Lipid IUPAC International Union of Pure and Applied Chemistry LAMMPS Largescale Atomic Molecular Massively Parrallel Simulator LJ

  ). Un comportement de "shear thinning" est observé à partir de t > 100 bars. Encart: zoom sur la région du régime linéaire. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.6 Déplacements moyens normalisés hX u i(t)/V 0 pour un ensemble de vitesses croissantes allant de 0, 01 à 0, 09 nm.ps 1 et de 0, 10.5 nm.ps 1 .117 A.7 Trajectoire individuelle (A), échantillon de 50 trajectoires moyenné Table of the different LJ interactions [6] . . . . . . . . . . . . . . . . . . . 3.1 Geometric characteristics of the simulated systems in the fluid and gel regimes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1 Interaction types for substrate-lipid and substrate-water interaction (see also figure 2.3 and table 2.1). . . . . . . . . . . . . . . . . . . . . . . 4.2 Summary of the parameters of the simulations . . . . . . . . . . . . . . 4.3 Lipid diffusion coefficient of a supported DSPC bilayer. DSPC Top and Bot represent the upper and lower leaflets with respect to the sub-List of simulations used in the constant pull force statistics. . . . . . . . A.1 Valeurs des coefficients de friction en fonction de la force appliquée, en régime de réponse linéaire . . . . . . . . . . . . . . . . . . . . . . . . A.2 Quelques coefficients de friction en fonction de l'impulsion initiale . . = 1.380 648 52 ⇥ 10 23 J K 1 Permittivity of vacuum e 0 = 8.854 187 812 8 ⇥ 10 12 F m 1

		xviii		xvii xix
		SFA	Surface Force Apparatus
		SLB	Supported Lipid Bilayers
	List of Abbreviations 2D Two dimensions 3D Three dimensions APL SUV Small Unilamellar Vesicles UP Upper Layer UW Upper Water VdW Van derWaals VMD Visual Molecular Dynamics V-rescale Velocity rescale Physical Constants Boltzmann constant A Area Chapter 0 General introduction K b List of Symbols a Acceleration Area PerLipid BL Bottom Layer BPTI b Coefficient of friction D Diffusion coefficient "Ever-newer waters for those who step m 2 m/s 2 Pa s m 1 cm 2 /ps into the same rivers: Panta Rei" Bovine Pancreatic Trypsin Inhibitor BW Bottom Water CG CHARMM COM CPF Constant Pull Force CPU DOPC k Box height (z) nm 1,2-dioleoyl-sn-glycero-3-phosphocholine In Number of lipids in the inner leaflet dimensionless Computation Processor Unit F Pulling force in GROMACS kJ mol 1 nm 1 F w Pulling force of the water N Center Of Mass F d Pulling force lower leaflet N Chemistry at HARvard Macromolecular Mechanics F u Pulling force upper leaflet N Coarse Grained F T Tangential force N F N Normal force N Heraclitus of Ephesus (c.535-c.475 BC)
	k		Box height (z)	nm
	L		Box height (z)	nm
	L b		Bilayer thickness	nm
	Water thickness m in gromacs Mass u (unified atomic mass unit) L w M Leaflet mass	nm u=1.660 538 921 ⇥ 10 27 kg u
	OUT	Number of lipids in the outer leaflet	dimensionless
	P		Order parameter	dimensionless
	Pe		Peclet number	dimensionless
	r		vector position (x,y,z)	nm
	t		Time	s or ps in gromacs
	T		Temperature	K
	V x		Velocity of the corresponding part of the system nm ps 1
	h		Lennard-Jones viscosity	Pa s
	MD t(z)	Molecular Dynamics Shear stress	Pa
	f x	MSD	Mean-Square Displacement Shear stress of the corresponding force F x	Pa
		NMR	Nuclear Magnetic Resonance
		NP	Nanoparticle
		NPT	Isobaric-Isothermal ensemble
		NVE	Microcanonical ensemble
		NVT	Canonical ensemble
		Ovito	Open Visualization Tool
		PBC	Periodic Boundary Conditions
		PC	PhosphoCholine
		PDB	Protein Data Bank
		POPC	1-palmitoyl-2-oleoyl-glycero-3-phosphocholine
		QCM-D	Dissipative Quartz-crystal microbalance
		RNEMD	Reverse Non-Equilibrium Molecular Dynamics
		S	S ubstrate
		SDK	Shinoda-DeVane-Klein model

(B), dispersion du déplacement moyen après ré-échantillonnage (bootstrap)(C) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.8 Snapshot d'une configuration à basse temperature (280 K). . . . . . . . A.9 Courbe caractéristique du déplacement du centre de masse des feuillets en fonction de la contrainte appliquée. Une loi de puissance (cercles, état gel) est clairement visible. . . . . . . . . . . . . . . . . . . . . . A.10 Schéma des differentes couches d'une membrane supportée. De bas en haut : S:Substrate ; BW:Bottom Water ; BL:Bottom Layer ; TP:Top Layer and TW:Top Water. . . . . . . . . . . . . . . . . . . . . . . . . . . A.11 Capture d'une simulation d'une membrane supportée (collaboration Claire Loison, ILM Lyon). De bas en haut : le substrat, une couche d'eau, une bicouche, une couche d'eau, de la vapeur d'eau . . . . . . . xv List of Tables 2.1 strate or z axis. Water Bot stands for the confined water and Water Top for the free water above the membrane. . . . . . . . . . . . . . . . . . . 4.4 Normalized velocities hVi/F of the center of mass of the different parts in a supported DSPC bilayer system, with both activated and deactivated COM removal. The normalized velocities are given in nm.ps 1 /(kJ.mol 1 .nm 1 ) . . . . . . . . . . . . . . . . . . . . . . . . . . 4.5 Normalized relative velocities when both leaflets are pulled in the same direction and the substrate in the opposite. The normalized velocities are expressed in nm.ps 1 /(kJ.mol 1 .nm 1 ) . . . . . . . . . . . . A.1 Atoms of DSPC molecule . . . . . . . . . . . . . . . . . . . . . . . . . . A.2 Bonds of DSPC molecule . . . . . . . . . . . . . . . . . . . . . . . . . . . A.3 Angles of DSPC molecule . . . . . . . . . . . . . . . . . . . . . . . . . . B.1 List of simulations used for the force kick relaxation statistics. The notation 512-10W stands for 512 lipids and 10 water beads per lipid. . B.2 DPPC 1,2-dipalmitoyl-sn-glycero-3-phosphocholine DPPS 1,2-dipalmitoyl -sn-glycero-3-phospho-L-serine DSPC 1,2-distearoyl-sn-glycero-3-phosphocholine EM Energy Minimization FATSLiM Fast Analysis Toolbox for Simulations of Lipid xxiii Dedicated to my parents and my little brother. . .

  The

	Name	Description	e (kJ/mol) s (nm)
	O	Super attractive	5.6	0.47
	I	Attractive	5.0	0.47
	II	Almost attractive	4.5	0.47
	III	Semi attractive	4.0	0.47
	IV	Intermediate	3.5	0.47
	V	Almost intermediate	3.1	0.47
	VI	Semi repulsive	2.7	0.47
	VII	Almost repulsive	2.3	0.47
	VIII	Repulsive	2.0	0.47
	IX	Super repulsive	2.0	0.62

TABLE 2 .

 2 1 -Table of the different LJ interactions[START_REF] Marrink | The MARTINI Force Field: Coarse Grained Model for Biomolecular Simulations[END_REF]

  .1.

	State L x (nm) L z (nm) A (nm 2 ) L b (nm) L w (nm)
	Fluid	13.2	8.2	174.	4.6	3.6
	Gel	11.1	10.6	124.	5.6	5.1

TABLE 3 .

 3

1 -Geometric characteristics of the simulated systems in the fluid and gel regimes.

  Normalized averaged displacements (upper leaflet) X u (t)/F for a set of increasing pulling forces 10, . . . , 2000 kJ.mol 1 .nm 1 (equivalently stresses t = 0.1 . . . 190 bars). The displacement for F = 2000 lies clearly beyond the linear regime and the force F = 10 is competing with thermal agitation. Inset: close-up focus on the normalized displacement for stresses larger than 0.955 bars. Ten different bootstrap realizations of the displacement hX u i(t) for an initial kick impulsion V 0 = 0.09 nm.ps 1 (upper leaflet). Inset: close-up look of the first peak.

	Displacement Bootstrap (V0=0.09 nm/ps)	0 -0.02 0 0.02 0.04 0.015 7000 7500 8000 8500 9000 9500 10000 2000 4000 6000 Time (ps) Time (ps) 0.01 0.02 0.025 0 Normalized displacement ({x(t) in nm} /F)) 100 200 300 6.63 6.66 6.69 6.72 6.75 0 10 6.7 6.72 6.74 6.76	400 20	8000	10000 0.0955 0.955 4.78 9.55 14.3 19.1 23.9 47.8 95.5 191 287 500 30
	FIGURE 3.6 -0 0 0.001 0 0.002 0.004 0.006 0.008 COM velocity (nm/ps) FIGURE 3.8 -0 100 0 0.5 1 1.5 Normalized displacement (x(t)/impulsion)	50	200 100 Applied pressure ( bar) 400 600 150 200 200 300 400	250 500 imp 0.01 imp 0.02 imp 0.03 imp 0.04 imp 0.05 imp 0.06 imp 0.07 imp 0.08 imp 0.09 imp 0.1 imp 0.2 imp 0.3 imp 0.4 imp 0.5	300
					Time (ps)
	0 0 FIGURE 3.9 -Normalized averaged displacements hX u i(t)/V 0 for a 500 1000 1500 2000 2500 3000 set of increasing velocities 0.01 and 0.03, . . . , 0.1 nm.ps 1 .
					Applied Force (kJ/mol/nm)
			FIGURE 3.7 -Averaged drift velocities vs applied forces (lower hor-
	Pull	izontal scale bar) or stresses (higher horizontal scale bar). A shear ). /TrajectoryAndBootstrapFluid-thinning deviation is seen at t > 100 bars. Inset: focus on the linear regime region.

  .FIGURE 3.11 -Ten different bootstrap realizations of the displacement hX u i(t) for an initial kick impulsion V 0 = 0.09 nm.ps 1 in the gel phase (upper leaflet). Inset: close-up look at the first peak.

	Displacement Bootstrap (V0=0.09 nm/ps)	0 5.6 5.65 5.7 5.75	100	200 0 5.55 5.6 5.65 5.7 5.75	300 10	400 20		500 30
					Time (ps)
	Normalized displacement (x(t)/impulsion)	0 0 0.5 1 1.5 2 2.5	100	200	300	400	500	600 imp 0.01 imp 0.02 imp 0.03 imp 0.04 imp 0.05 imp 0.06 imp 0.07 imp 0.08 imp 0.09 imp 0.1 imp 0.2 imp 0.3 imp 0.4 imp 0.5
				Time (ps)

  FIGURE 3.17 -Apparent friction coefficient b + h/L w from impulsion simulations as a function of the initial induced velocity in the fluid phase. A vertical error bar is inferred from a bootstrap procedure with 10 synthetic averaged displacement curves. The error bar is very large for the initial velocities smaller than 0.7 nm.ps 1 , and becomes small for initial velocities larger than 0.2 nm.ps 1 , for which a departure from linear response can bee seen in Fig.3.10. The horizontal dashed line stands for the CPF value b + h/L w = 2.75 ⇥ 10 6 .

		4.0×10	6
	Friction Coefficient (Pa.s.m ) -1	1.5×10 6 2.0×10 6 2.5×10 6 3.0×10 6 3.5×10 6			Over all simulations (fluid) From Boostrap method (fluid)
		1.0×10 6	0.01		0.1	1
					Initial velocity V 0 (nm.ps)	-1
	) -1 Friction Coefficient (Pa.s.m	5.0×10 1.0×10 6 7			Over all simulations (fluid) From Boostrap method (fluid) Over all simulations (gel) From Bootstrap method (gel)
		0 0.0	0.1	0.2	0.3	0.4	0.5
					Initial velocity V 0 (nm.ps)	-1
							.21, associated to an initial

TABLE 4 .

 4

1 -Interaction types for substrate-lipid and substrate-water interaction (see also figure 2.3 and table 2.1).

TABLE 4 .

 4 2 -Summary of the parameters of the simulations

	Sim type State (K)	Box size (nm)	Lipids Other beads Stress(bar) Nbr Sim
	Pull	340	x, y, z = 13.2 nm	512	5120 Water 770 Substrate	F = 250 t = 23.9	50

Comparison between 2D and 3D: Diffusion coefficients DSPC Top DSPC Bot WATER Top WATER Bot 2D (cm 2 .ps 1 ) 1

  .12 ⇥ 10 6 8.54 ⇥ 10 7 8.19 ⇥ 10 4 5.08 ⇥ 10 5

	2
	Water Bot
	Water Top
	FIGURE 4.4 -Volume diffusion of a supported DSPC bilayer in water
	(3D)

TABLE 4 .

 4 3 -Lipid diffusion coefficient of a supported DSPC bilayer. DSPC Top and Bot represent the upper and lower leaflets with respect to the substrate or z axis. Water Bot stands for the confined water and Water Top for the free water above the membrane.

TABLE 4 .

 4 [START_REF] Falk | Interleaflet sliding in lipidic bilayers under shear flow: comparison of the gel and fluid phases using reversed non-equilibrium molecular dynamics simulations[END_REF] -Normalized velocities hVi/F of the center of mass of the different parts in a supported DSPC bilayer system, with both activated and deactivated COM removal. The normalized velocities are given in nm.ps 1 /(kJ.mol 1 .nm1 )

TABLE 4 .

 4 Normalized averaged displacement of the center of mass of the different system subsets (top and bottom leaflets, top and bottom water above and motionless substrate. The global center of mass is free. The entire bilayer is pulled in one direction and the substrate in the opposite direction, parallel to the x axis with a force of 250 kJ.mol 1 .nm1 . Again the substrate is steady. The entire bilayer is pulled in one direction and the substrate in the opposite direction parallel to the x axis with a force equal to 250 kJ.mol 1 .nm1 . The substrate stays at rest. ⇥ 10 6 2.00 ⇥ 10 6 1.93 ⇥ 10 6 1.67 ⇥ 10 7 6.57 ⇥ 10 7 Off 1.64 ⇥ 10 5 1.73 ⇥ 10 5 1.55 ⇥ 10 5 1.72 ⇥ 10 5 7.15 ⇥ 10 6 5 -Normalized relative velocities when both leaflets are pulled in the same direction and the substrate in the opposite. The normalized velocities are expressed in nm.ps 1 /(kJ.mol 1 .nm1 )

		0.2
				DSPC bilayer
	Normalized displacement (x(t)/F)	0.05 0.1 0.15	W_Bot W_Top DSPC_bot Substrat DSPC_top
		0	0	2000	4000	6000	8000	10000
					Time (ps)
		FIGURE 4.9 -

TABLE A .

 A

3 -Angles of DSPC molecule

  .1 and B.2 summarizes the characteristics of the trajectories used in the present study.

	Simulation type State (K)	Box size (nm)	nbr Lipids Velocity (nm.ps 1 ) nbr Runs
	FKR	fluid (340)	x-y: 13.18; z: 8.17	512-10W	0.01:0.1 step=0.01 0.1-0.5 step=0.1	150
			x-y: 13.15 ; z: 11.88	512-20W	0.05:0.1 step 0.01 0.2
			x-y: 18.55 ; z: 8.24	1024-10W	0.04:0.1 step=0.01 0.2
	FKR-Tilt		x-y: 13.18 ; z: 8.17	512-10W	0.09 0.00	500
	FKR	gel (280)	x-y: 11.14 ; z: 10.64	512-10W	0.01:0.1 step=0.01 0.1-0.5 step=0.1	150
	FKR-Tilt		x-y: 11.14 ; z: 10.64	512-10W	0.09 0.00	500

TABLE B .

 B 1 -List of simulations used for the force kick relaxation statistics. The notation 512-10W stands for 512 lipids and 10 water beads per lipid.

	Simulation type State (K)	Box size (nm)	Lipids	Stress (bar)	Number of runs
					F=10 ; t = 0.955
					F=50 ; t = 4.78
					F=100 ; t = 9.55
					F=150 ; t = 14.3
	CPF	fluid (340) xy:13.18 ; z:8.17 512-10W	F=200 ; t = 19.1 F=250 ; t = 23.9	50
					F=500 ; t = 47.8
					F=1000 ; t = 95.5
					F=2000 ; t = 191
					F=3000 ; t = 2.87
					F=200; t = 9.66
	CPF	fluid (340) xy:18.55; z:8.24 1024-10W	F=250; t = 12.1	50
					F=400; t = 19.3
					F=10 ; t = 1.34
					F=50 ; t = 6.69
					F=100 ; t = 13.4
					F=150 ; t = 20.1
	CPF	gel (280)	xy:11.14; z:10.6 512-10W	F=200 ; t = 26.8 F=250 ; t = 33.4	50
					F=500 ; t = 66.9
					F=1000 ; t = 134.
					F=2000 ; t = 268
					F=3000 ; t = 401.

TABLE B .

 B 2 -List of simulations used in the constant pull force statistics.; non-bonded cutoff (1.2 nm). However, to improve energy conservation ; or to avoid local heating/cooling, you may increase the update frequency ; and/or enlarge the neighbourlist cut-off (to 1.4 nm). The latter option ; is computationally less expensive and leads to improved energy conservation Standard cut-off schemes are used for the non-bonded interactions ; in the Martini model: LJ interactions are shifted to zero in the ; range 0.9-1.2 nm, and electrostatic interactions in the range 0.0-1.2 nm. ; The treatment of the non-bonded cut-offs is considered to be part of ; the force field parameterization, so we recommend not to touch these ; values as they will alter the overall balance of the force field. ; In principle you can include long range electrostatics through the use ; of PME, which could be more realistic in certain applications ; Please realize that electrostatic interactions in the Martini model are ; not considered to be very accurate to begin with, especially as the ; screening in the system is set to be uniform across the system with ; a screening constant of 15. When using PME, please make sure your For use in combination with the Verlet-pairlist algorithm implemented ; in Gromacs 4.6 a straight cutoff in combination with the potential ; modifiers can be used. Although this will change the potential shape, ; preliminary results indicate that forcefield properties do not change a lot ; when the LJ cutoff is reduced to 1.1 nm. Be sure to test the effects for ; your particular system. The advantage is a gain of speed of 50-100%.

	;verlet-buffer-drift	= 0.005
	; MARTINI and CONSTRAINTS
	; for ring systems and stiff bonds constraints are defined
	nstlist ; which are best handled using Lincs. = 10
	ns_type	= grid
	pbc	= xyz
	rlist	= 1.4
	; MARTINI and NONBONDED
	; ; system properties are still reasonable.
	;
	; With the polarizable water model, the relative electrostatic screening
	; (epsilon_r) should have a value of 2.5, representative of a low-dielectric
	; apolar solvent. The polarizable water itself will perform the explicit screening
	; in aqueous environment. In this case, the use of PME is more realistic.
	;
	; coulombtype	= cut-off ;
	rcoulomb_switch	= 0.0
	rcoulomb	= 1.2
	epsilon_r	= 15 ; 2.5 (with polarizable water)
	vdw_type	= Cut-off ;
	rvdw_switch	= 0.9
	rvdw	= 1.2 ;1.1 (for use with Verlet-pairlist)
	;cutoff-scheme	= verlet
	;coulomb-modifier	= Potential-shift
	;vdw-modifier	= Potential-shift
	;epsilon_rf	= 0	; epsilon_rf = 0 really means epsilon_rf = infinity

2 Détermination du paramètre de friction inter-feuillet.

  FIGURE A.6 -Déplacements moyens normalisés hX u i(t)/V 0 pour un ensemble de vitesses croissantes allant de 0, 01 à 0, 09 nm.ps 1 et de 0, 10.5 nm.ps 1 .

	Normalized displacement (x(t)/impulsion)	imp 0.01 imp 0.02 imp 0.03 imp 0.04 imp 0.05 imp 0.06 imp 0.07 imp 0.08 imp 0.09 imp 0.1 imp 0.2 imp 0.3 imp 0.4 imp 0.5 FIGURE A.7 -Trajectoire individuelle (A), échantillon de 50 tra-jectoires moyenné (B), dispersion du déplacement moyen après ré-échantillonnage (bootstrap)(C) Force (kJ.mol-1.nm-1) b (Pa.s.m-1) Db 100 2,67⇥10 6 9,79⇥10 4 150 2,78⇥10 6 1,01⇥10 5 200 2,75⇥10 6 7,69⇥10 4 250 2,76⇥10 6 6,59⇥10 4
		500	2,75⇥10 6	2,78⇥10 4

A.2.

A.2.2.1 Contrainte appliquée stationnaire « CPF »

L'hydrodynamique de Stokes permet de relier le coefficient de friction b à la vitesse de déplacement V h du centre de masse du feuillet, et aux différents paramètres tels que la viscosité de l'eau h, l'épaisseur d'eau L w , et la contrainte F h (eq.A.1).

TABLE A .

 A 1 -Valeurs des coefficients de friction en fonction de la force appliquée, en régime de réponse linéaire

Acknowledgements

La réalisation de cette thèse a été possible grâce au concours de plusieurs personnes à qui je souhaiterais témoigner toute ma gratitude.

Je tiens à remercier en tout premier lieu mon directeur de thèse, Fabrice Thalmann. Je le connais maintenant depuis mes premières années à l'université. Il a d'abord été pour un moi un professeur en licence avant de devenir mon directeur lors d'un formidable stage en licence 2 qui par ailleurs m'a introduit au monde des phospholipides et de la fluorescence. Par la suite, il m'a enseigné en master pour enfin devenir mon directeur de thèse à l'ICS. Il a énormément contribué à ce que ma thèse se déroule dans d'excellentes conditions professionnelles et humaines notamment en étant toujours à l'écoute et disponible, même quand il s'agit de discuter de notre passion commune qui est le badminton. J'adresse tous mes remerciements au Pr. Joly ainsi qu'au Pr. Sommer, de l'honneur qu'ils m'ont fait en acceptant d'être rapporteur de cette thèse et d'avoir pris le temps de lire et commenter mon travail.

Je remercie également Pr. Le Houérou de l'honneur qu'il m'a fait d'avoir accepté de présider ma défense ainsi que d'avoir examiné mon travail.

Je tiens particulièrement à remercier Dr. Loison d'avoir accepté d'être examinatrice ainsi que d'avoir partagé son savoir avec moi concernant la simulation des lipides durant ma thèse.

Je pense que je ne serais probablement pas là où je suis sans l'aide de Vivien Walter. Le « stagiaire » comme tu disais, à maintenant grandi et pense à verrouiller son ordinateur pour que personne ne modifie son bash_profile quand il est absent et ceci grâce à toi. Je ne pourrais jamais te remercier assez pour toute ton aide durant mon stage de master 2 et ma thèse notamment pour ta patience incroyable, ta pédagogie extraordinaire et surtout ton amitié.

Chapter 4

Supported Bilayers

"Quand l'hypothèse est soumise à la méthode expérimentale, elle devient une théorie; tandis que si elle est soumise à la logique seule, elle devient un système."

Claude Bernard (1813-1878)

Introduction

T HE previous chapter was about a free lipid bilayer in water. In this one we intro- duce a solid substrate close and parallel to the bilayer, interacting with it. This system is a more realistic approximation of the lubricated surfaces, it introduces a confined water layer.

Composition of the system

The system is built inside a cubic box of side 13.2 nm (fig. 4.1). The z direction is perpendicular to the lipid membrane plane. This system was initially set-up by Li Fu, Florian Benedetti and Claire Loison (Institut Lumière Matière, Lyon). of equations:

which are easily solved as:

and gives in the end V BL = V BW = 0 and V TL = F/(Ab). The bottom water and leaflet layers stay at rest with vanishing velocities. This may seem counterintuitive at first, as one pulls on the bottom leaflet also. In fact, the union of the top four layers (TW, TL, BL, BW) taken as a single mechanical system is subject to no resulting force (the pair of forces cancels). If there were no substrate, the four top layers could move together at any arbitrary average velocity. However, the presence of a fixed surface layer breaks the translation invariance of the whole system. Out of all possible solutions, V BL = V BW = 0 is the one which minimizes the dissipation in the bottom layer.

If one now considers the case of a total force F exerted on both leaflets, one needs to minimize P dis F(V BL + V TL )/2 with respect to the velocities. This gives

leading to Ag 2 V BL = 2F, and V TL = (1 + g 2 /(4b))V BL . The coefficient g 2 /2 can easily be identified to h/L w with h the solvent viscosity and L w the water layer We used the Martini lipid version v2.0 and Gromacs 5.1. The representation of a DSPC lipid is described in Fig. B.1. It consists in 14 beads located at various levels on an hydrophilicity scale, interacting with Lennard-Jones interactions with radius r 0 = 0.47 nm, connected with harmonic springs of stiffness k 0 = 1250 kJ.mol 1 [START_REF] Marrink | The MARTINI Force Field: Coarse Grained Model for Biomolecular Simulations[END_REF][START_REF] Marrink | Coarse grained model for semiquantitative lipid simulations[END_REF].

In all the simulations, the standard Gromacs md leap-frog molecular dynamics integrator was used, with a time step of 20 fs. The velocity rescale [START_REF] Bussi | Canonical sampling through velocity rescaling[END_REF] was used to keep the internal energy constant in the simulation. This barostat is an alternative to the Nose-Hoover thermostat based on the same principle of extended phase space, with a single supplementary stochastic coordinate Q, ensuring canonical ensemble ergodicity for the simulated system. Lipid and solvent groups of molecules were separately coupled to two V-rescale thermostats, with a coupling time constant of 1 ps. For constant pressure simulations, we used a semi-isotropic Parinello-Rahman barostat with a time coupling constant of 12 ps and a compressibility 3 ⇥ 10 4 bar 1 in the xy and z directions.

Appendix C

Molecular dynamics parameters

This section presents some examples of molecular dynamics parameter used during this thesis. These latter are inspired from the one found in the Martini website.

C.1 Energy minimization: EM

Resumé en Français

"A different language is a different vision of life"

Federico Fellini (1920-1993)

A.1 Introduction

L A nature produit des systèmes lubrifiants à base d'eau qui surclassent la majorité des dispositifs créés par l'homme [START_REF] Urbakh | The nonlinear nature of friction[END_REF]. L'exemple le plus probant est le cas des contacts biologiques tels que les articulations (Fig. A.1). Ces systèmes présentent deux surfaces cartilagineuses séparées par une poche de synovie. On retrouve dans la composition du fluide synovial des multicouches de phospholipides qui jouent un rôle majeur dans l'obtention d'un coefficient de friction très faible (m=Force transverse/Force normale = 0,005-0,02) [START_REF] Schroeder | Surface Active Phospholipids as Cartilage Lubricants en[END_REF]. La compréhension de tels systèmes biophysiques présente des enjeux aussi bien fondamentaux qu'appliqués, en lien par exemple avec l'arthrose et la confection de prothèses articulaires [START_REF] Klein | Repair or Replacement-A Joint Perspective[END_REF].Le but de notre étude est de développer des systèmes numériques à l'échelle moléculaire (approche à gros grains [START_REF] Marrink | The MARTINI Force Field: Coarse Grained Model for Biomolecular Simulations[END_REF]), au contact d'une surface solide, rigide et d'en étudier les propriétés mécaniques. Nous envisageons aussi le cas de membranes lipides placées dans des situations non conventionnelles, telles que des vésicules de très forte courbure, entourant des nanoparticules de silice par exemple. « b », associé au mouvement de glissement relatif des feuillets. Aux temps intermédiaires, la relation stress t(t)-vitesse de déplacement V(t) requiert des fonctions de mémoire généralisant la notion de module viscoélastique.

A.2.1 Obtention du régime linéaire dans l'état fluide

Un régime linéaire se manifeste par la validité du principe de superposition entre causes et effets. L'existence d'un régime linéaire suppose un temps de relaxation du système suffisamment court, et l'absence de changement structural induit par la sollicitation extérieure. Il se traduit en pratique par l'existence de courbes maitresses décrivant le rapport entre l'effet et la cause. Le régime linéaire a été déterminé de façon empirique pour chacune des deux méthodes précédentes. Dans le premier cas, le système a été soumis à une gamme de contrainte appliquée comprise entre 0,1 et 300 bars, au moyen d'un couple de forces appliqué en chaque centre de masse des feuillets. La représentation du déplacement rapporté à la force en fonction du temps permet bien de déterminer un régime décrit par une courbe maitresse (fig.A.4, fig.A.5). De la même manière, une courbe maitresse peut être obtenue dans le cas impulsionnel, en imposant un incrément de vitesse allant de 0,01 à 0,5 nm.ps-1 (fig.

A.6). Un régime linéaire est visible pour des contraintes allant de 5 à 50 bars pour la méthode du « pull » et des vitesses initiales allant de 0,05 à 0,1 nm.ps-1 pour les impulsions. Pour des forces ou des vitesses trop faibles, le déplacement est dominé par l'agitation thermique. Au contraire, si les sollicitations extérieures sont trop élevées, le système sort du régime linéaire. Le systeme cette fois se trouve dans un état qualifié de gel dans lequel les phospholipides possèdent une plus grande organisation structurelle (fig.A.8). A basse température la caractéristique force-appliquée-vitesse de déplacement n'est nulle part linéaire. Au lieu de cela, un régime de rhéo-amincissement (décroissance du coefficient de friction avec la force exercée est clairement visible, donnant lieu à une loi de puissance apparente V u ⇠ t 1.44

A.2.3 Absence de régime linéaire dans l'état gel

A.2.4 Généralisation à la géométrie d'une membrane supportée

Nous avons montré que toutes les méthodes développées ci-dessus pouvaient se généraliser à des membranes au contact de substrat solides tels que celle représentée sur la figure A.11. Le système comprend un empilement de cinq couches. La couche inférieur S contient des billes gelées. Celle-ci peuvent supporter n'importe quelle force externe arbitraire sans se déplacer.

Nous suggérons cependant que notre méthode soit une alternative aux pores pour l'équilibration des feuillets de petites vésicules.