Daniel Abergel

Mark Johnson

Keywords: subdi↵usion, simulations, molecular dynamics, lipids, membranes

A PhD thesis might appear as a solitary undertaking, but in reality there are many people who contributed in various ways to the research in this work and created opportunities for me to gain invaluable experience and knowledge.

First and foremost, I would like to thank Gerald Kneller for being my supervisor and mentor. I appreciate all his help, ideas and constructive critiques, which taught me professionalism and encouraged me to immerse myself in scientific work. During our conversations we often talked beyond academic matters and also these discussions were important life lessons, which I enjoyed greatly. Furthermore, I am grateful to Konrad Hinsen for his valuable tips

Résumé

De nombreuses études récentes, expérimentales et in-silico, montrent que la di↵usion latérale des molécules dans des membranes biologiques présentent des anomalies, dans le sens que les déplacements carré moyen moléculaires évoluent de façon sublinéaire au lieu de linéaire avec le temps. Mathématiquement, ce type de di↵usion peut être modélisé par des équations de di↵usion généralisées, dans lequelles une dérivée fractionnaire du temps s'ajoute à l'équation de di↵usion normale correspondante.

Le but de cette thèse est d'obtenir un aperçu plus physique des processus de di↵usion dans des membranes biologiques. A cet e↵et, des simulations de dynamique moléculaire d'une bicouche lipidique POPC sont analysées en utilisant les concepts de la physique statistique des liquides. La queue aux temps longs de la fonction d'autocorrélation de vitesses des centres de masse, qui réflète le régime de di↵usion des molécules en question, est mise en relation avec leur dynamique aux basses fréquences et avec la structure dynamique autour de chaque molécule. La dernière est caractérisée par la fonction de corrélation de paires dynamique de van Hove pour leurs centres de masse. Il est en particulier montré que la première couche de voisins d'une molécule lipidique qui di↵use ne se désintègre que très lentement avec t , où 0 < < 1. Ce résultat est en accord avec l'observation faite par d'autres auteurs, que les molécules dans une bicouche lipidique ont la tendance de se déplacer d'une manière concertée.

Afin d'évaluer l'impact du champ de force sur la nature des processus de di↵usion observés, la bicouche lipidique POPC a été simulée avec le champ de force tout atome OPLS et avec le champ de force à gros grains MARTINI. Dans le second, quatre atomes lourds forment un seul "grain". Dans les deux cas, on observe une sous-di↵usion latérale, avec des exposants similaires, mais la di↵usion des lipides obtenue avec le champ de force à gros grains est d'environ trois fois plus rapide et elle apparaît aussi plus rapide que dans des expériences correspondantes. Ce résultat est confirmé par la dynamique des molécules POPC à basse fréquence et par la structure dynamique de leur environnement local.

Mots clés : sous-di↵usion, simulations, dynamique moleculaire, lipides, membrane

Résumé Substantiel I. INTRODUCTION

Une membrane biologique est une enveloppe continue qui joue le rôle de barrière physique délimitant la cellule: elle sépare le milieu intracellulaire (tout ce qui compose la cellule) du milieu extracellulaire (l'environnement de la cellule). Il existe également des membranes biologiques intracellulaires, qui délimitent les organites (par exemple les membranes mitochondriale, nucléaire, lysosomiale, etc.) [Berg et al. [2005]]. Toutes les membranes sont composées principalement de lipides, dont les plus abondants sont les phospholipides (Fig. 1). Ces lipides forment des mono-ou bicouches fermées. Une membrane biologique contient également de nombreuses autres molécules, comme des protéines membranaires ou des glucides [START_REF] Alberts | Essential cell biology[END_REF]]. Leur présence permet le transport bi-directionnel de petites molécules à travers la membrane, ce qui intervient dans l'adhésion des cellules et dans des nombreux processus enzymatiques. Tous ces processus sont fortement influencés par la di↵usion latérale (dans le plan de la bicouche) des lipides, qui peut être caractérisée par l'évolution dans le temps du déplacement carré moyen (DCM) d'une molécule, qui, dans le régime asymptotique t ! 1, prend la forme

W (t) = h(x(t) x(0)) 2 i t!1 ⇠ 2D ↵ t ↵ , 0 < ↵ < 2. ( 1
)
Ici x signifie la position de la molécule , D ↵ est la constante fractionnaire de di↵usion et le paramètre ↵ décrit le type de di↵usion. Pour ↵ = 1, nous avons une di↵usion v dite normale, ce qui caractérise les liquides simples (par exemple l'eau). Le régime 0 < ↵ < 1 est typique pour les systèmes "crowded". On l'appelle la sous-di↵usion. Pour 1 < ↵ < 2, nous avons la super-di↵usion, ce qui signifie que les mouvements des molécules sont accélérés. Les expériences et les simulations par dynamique moléculaire montrent des résultats contradictoires pour la di↵usion latérale des lipides. Les expériences de di↵usion ]. La sous-di↵usion peut être comprise par un "e↵et de cage" qui s'exprime dans la fonction d'autocorrélation des vitesses C vv = hv(0) • v(t)i, où v signifie la vitesse de la molécule, par une valeur négative à des temps longs [START_REF] Kneller | Generalized Kubo relations and conditions for anomalous di↵usion: physical insights from a mathematical theorem[END_REF]]. Dans cette thèse, cette "cage", formée par les voisins les plus proches autour de chaque lipide, est étudiée par la fonction de corrélation de paire vue comme fonction du temps [START_REF] Rahman | Intermediate Scattering Function in Slow Neutron Scattering[END_REF]; Van Hove [1954]]:

G D (r, t) = V N n(r, t) V . (2)
Elle montre les fluctuations de la structure formée par les voisins lipidiques à la distance r de la molécule entourée. Ici V et N sont le volume total du système et le nombre de molécules, respectivement, et n(r, t) est le nombre de molécules présentes dans l'élément de volume/surface V à la distance r.

L'objectif de cette thèse est d'obtenir une description plus physique des processus de di↵usion des lipides, en s'appuyant sur la simulation par dynamique moléculaire (MD) et la théorie des liquides. Les trajectoires MD ont été analysées pour en extraire le comportement à des temps longs des fonctions de corrélation.

II. SIMULATIONS DE DYNAMIQUE MOLECULAIRE

Une série de simulations par dynamique moléculaire d'une bicouche lipidique de 1-palmitoyl-2-oleoylphosphatidylcholine (POPC) a été calculée pour deux modèles moléculaires : toutatome OPLS [START_REF] Jorgensen | The OPLS [optimized potentials for liquid simulations] potential functions for proteins, energy minimizations for crystals of cyclic peptides and crambin[END_REF]; Jorgensen et al. [1996]] et gros-grain MARTINI [START_REF] Marrink | Coarse Grained Model for Semiquantitative Lipid Simulations[END_REF][START_REF] Marrink | The MARTINI Force Field: Coarse Grained Model for Biomolecular Simulations[END_REF]] (Fig. 2).

Une membrane tout-atome consistant de 274 molécules de POPCa été simulée penvii AA CG dant 15 ns et 150 ns, dans l'ensemble isobare-isotherme [START_REF] Stachura | Anomalous lateral di↵usion in lipid bilayers observed by molecular dynamics simulations with atomistic and coarse-grained force fields[END_REF]]. La température de 310 K a été contrôlée par un thermostat de Nosé-Hoover [START_REF] Hoover | Canonical dynamics: Equilibrium phase-space distributions[END_REF];

Nosé [1984]] et la pression de 1 atm par un barostat de Parrinello-Rahman [START_REF] Parrinello | Polymorphic transitions in single crystals: A new molecular dynamics method[END_REF]].

Afin d'avoir un accès à des temps de simulation plus longs pour une plus grande bicouche de POPC, des simulations MARTINI ont été e↵ectuées pour un système de 2033 lipides pendant 600 ns, dans deux conditions thermodynamiques: NV T et N Ap z T (pression ajustée uniquement dans la direction perpendiculaire à la bicouche) [Stachura & Kneller viii [2013]]. La température de 320 K a été fixée par le thermostat de Berendsen et la pression de 1 atm pour l'ensemble N Ap z T par la barostat de Berendsen [START_REF] Berendsen | Molecular dynamics with coupling to an external bath[END_REF]. Ces calculs ont été aussi servi pour tester la capacité du modèle gros-grains de reproduire les propriétés dynamiques des systèmes tout-atome.

III. RESULTATS

Les analyses de l'évolution temporelle du déplacement carré moyen (DCM) du centre de masse des lipides ont montré des devations de la forme linéaire caractéristique pour la di↵usion normale. Les fits de l' éq. 1 (Fig. 3) ont donné des paramètres ↵ ⇡ 0.68 et

Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡
Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê ÊÊ ÊÊÊÊ Ê ÊÊÊ Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê 0.1 0.2 0.
= 0 Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê ÊÊ Ê ÊÊÊ ÊÊÊ ÊÊ Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê 0.

LogHt @nsDL

LogHGDHr

= 0 Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê ÊÊ ÊÊÊ ÊÊÊ ÊÊ Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê 0.

LogHt @nsDL

LogHGDHr,tLL

NApZT CG POPC

Fit of t -b (d)

t @psD

G D Hr,tL

Data points

Fit of exp@têtD

Figure 6: La fonction de corrélation de paire distincte dépendantes du temps, G D (r, t), pour la simulation de l'eau SPC/E. L'e↵et de cage est préservé pendant des temps allant jusqu'à 10% de la longueur totale des trajectoires MD. En revanche, le G D (r, t) pour la décroissance de l'eau est beaucoup plus rapide (3 ordres de grandeur) et aussi exponentielle (Fig. 6). Ce résultat confirme xii l'existence d'une forte corrélation entre les lipides et leurs voisins, l'e↵et de cage, qui est aussi une explication plausible pour la sous-di↵usion des lipides POPC dans la bicouche.

Dans le cas des simulations gros-grains, la décroissance de la première couche de voisins est aussi algébrique, mais la décorrélation complète des fonctions G D (r, t) est visible pour des temps inférieurs à 1% de la longueur totale des trajectoires MD. Ceci pourrait indiquer les limites du modèle gros-grains MARTINI pour reproduire complètement les propriétés dynamiques des systèmes lipidiques.

IV. CONCLUSION

Dans cette thèse, on a été montré que les mouvements latéraux des lipides POPC sont sous-di↵usifs pour deux modèles: le modèle tout-atoms OPLS et le modèle gros-grains MARTINI. L'e↵et de cage associé à la sous-di↵usion a été mis en évidence par la fonction de corrélation de paire dépendante du temps, G d (r, t). Pour les simulations des lipides POPC, G d (r, t) a montré une couche de voisins les plus proches très stable dans le temps, qui se désintègre avec le temps de façon algébrique.

Results of one of the first molecular dynamics simulation published by

Alder [START_REF] Alder | Phase Transition for a Hard Sphere System[END_REF][START_REF] Alder | Studies in Molecular Dynamics. I. General Method[END_REF] Lipids are a very diverse group of chemical compounds that are either hydrophobic or amphipathic (soluble in both -water and fats) [Berg et al. [2005]]. The lipids can refer to amassed fatty acids in a volume, more complex steroids, or phospholipids. Table 1.1 shows that cholesterol, cardiolipin and phospoholipids are the most abundant lipids in animal cell membranes. The percentage of the total masses in di↵erent tissues indicates that the most abundant are phosphatidylcholine and phosphatidylethanolamine (phospholipids)

and cholesterol for erythrocytes.

.4 ± 0.2 - 1.4 ± 0.3 2.0 ± 1.0 1.8 Cardiolipin 2.2 ± 0.3 - 14.8 ± 1.2 1.6 ± 1.3 - 3.
Tail(s) consisting of hydrocarbon chains that are derivatives of fatty acids. They can be either saturated (all the carbon atoms in the chains are connected by single bonds) or mono-and poli-unsaturated (at least one double bond).

The lipid composition of membranes presented in Table 1.2 for various organelles shows again that the most common ones are phosphatidylcholines and phosphatidylethanolamines.

The abundance of the first does not fall below 40% of the total weight of the membrane for all the presented organelles and is the lowest for mitochondria, which exhibits on the other hand the highest concentration of cardiolipin and phosphatidylethanolamine. The latter is the second most abundant type of lipid, with the lowest abundance of 20% in a plasma membrane. According to these findings, phosphatidylethanolamine and in particular phosphatidylcholine are the basic lipids constituting biological membranes in cells.

They stand for 51% of the total weight of all membranes in a tissue.

- - - - 5 8
The comparison above referred only to the properties of the headgroup moieties. The other di↵erences concern the lengths of the hydrocarbon chains (tails) and their level of saturation (number of double bonds). In Table 1.3 the content of diverse fatty acid tails is presented for six animal species: rat, pig, duck, horse, herring and seal. The nomenclature in the first column should be read as -the total number of carbon atoms in a chain : number of double bonds ("20 : 1" means for example a tail consisting of twenty carbon atoms with one double bond). As one can easily see, the most abundant tails are the saturated palmitoleic fatty acid rests (sixteen carbon atoms) and oleic acid rests with one double bond (eighteen carbon atoms). This result suggests that the most common lipid would be a phosphocholine with palmitoleic and oleic tails, which would be a 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC). An interesting observation is that higher order lipids with the unsaturated hydrocarbon chains of twenty and more carbon atoms are found only in aquatic animals leaving in cold environment, here herring and seal. The surface area of a membrane is huge in comparison to macromolecules, so the probability of finding a specific membrane protein is low. Many processes in cells involve several proteins, and for that reason macromolecules taking part in the same reaction must stay localized and close to each other in biological membranes [Berg et al. [2005]].

To raise the probability of finding the necessary proteins in the reaction volume, nature developed several strategies. One of them is a confinement of proteins within physical barriers, which are formed by the cytoskeleton attached to the membrane. Another strat-egy, more interesting from the perspective of this work, is to slow down the di↵usion of the lipids and proteins, keeping them more localized in the bilayer. In this case, there is no need for physical barriers, since the very nature of the di↵usion of molecules would already raise the probability of keeping them in the same place. One of the examples is the self-organization of lipids into domains (so-called rafts), which seem to attract proteins [START_REF] Risselada | The molecular face of lipid rafts in model membranes[END_REF]; [START_REF] Silva | Ceramide-Domain Formation and Collapse in Lipid Rafts: Membrane Reorganization by an Apoptotic Lipid[END_REF]]. These proteins di↵use with the lipids in the raft.

Transport of molecules in cells and membranes

The di↵usion of molecules observed in living cells deviates in many cases from normal di↵usion (the definition and its generalization is described in more details in the next section), where normality is defined by Einstein's di↵usion law [START_REF] Einstein | Uber die von der molekularkinetischen Theorie der Warme gefordete Bewegung von in ruhenden Flussigkeit suspendierten Teilchen[END_REF]], which predicts that the mean square displacement of the di↵using particle grows linearly with time. In biological membranes one observes instead an anomalous sublinear growth. As biological bilayers are vivid mosaics containing not only lipids but also proteins, the scope of experimentalists also focused on the di↵usion of the latter. FCS studies of the Golgi resident membrane proteins [START_REF] Weiss | Anomalous Protein Di↵usion in Living Cells as Seen by Fluorescence Correlation Spectroscopy[END_REF]] tracked the subdi↵usion of these molecules in the endoplasmic reticulum. Other findings using SPT on potassium membrane channels, imply anomalous di↵usion of these proteins [START_REF] Weigel | Ergodic and nonergodic processes coexist in the plasma membrane as observed by single-molecule tracking[END_REF]].

Anomalous di↵usion has been observed in cells by

Movements which are supposedly a↵ected by the actin cytoskeleton network.

Modeling anomalous di↵usion

In 1855 Adolf Fick published one of the first theoretical papers on di↵usion, which was inspired by Fourier's work on heat transfer. The basic assumption in his work is that the direction of the di↵usive flux, J, is from higher to lower concentrations of a substance [START_REF] Fick | Uber Di↵usion[END_REF]],

J(x, t) = DrC(x, t), (1.1)
where D is the di↵usion coe cient. Expressing the conservation of the total number of particles in form of en equation of continuity, @ t C(x, t) + r • J(x, t) = 0, it follows then that C(x, t) fulfills the well-known di↵usion equation,

@ @t C(x, t) = D C(x, t).
(1.2)

General Introduction

The solution of Eq. (1.2) is a Gaussian function. Later Einstein has given a probabilistic interpretation for Fick's model, in which he has shown that this Gaussian function is the mean squared displacement of the di↵using particles [START_REF] Einstein | Uber die von der molekularkinetischen Theorie der Warme gefordete Bewegung von in ruhenden Flussigkeit suspendierten Teilchen[END_REF]]:

W (t) = h(x(t) x(0)) 2 i = 2D t. ( 1
W (t) = 2D ↵ t ↵ , 0 < ↵ < 2, (1.4)
where the di↵usion coe cient D is replaced by its fractional counterpart. The ↵ parameter is defined in the range ↵ 2 (0, 2) and characterizes the type of di↵usion. The regime 0 < ↵ < 1 is typical for crowded systems and is called subdi↵usion. For 1 < ↵ < 2 we have superdi↵usion, suggesting the presence of forces accelerating the motions of the particle under consideration. Anomalous di↵usion can be modelled by considering probability density functions of waiting times, t , and lengths of jumps, j , which is the base of the continuous time random walk model [START_REF] Montroll | Random walks on lattices. IV. Continuoustime walks and influence of absorbing boundaries[END_REF]; [START_REF] Montroll | Random Walks on Lattices[END_REF]; [START_REF] Shlesinger | Asymptotic solutions of continuous-time random walks[END_REF]]. In case of subdi↵usive processes, waiting times can be very long and the characteristic waiting time, R 1 0 dt t (t), diverges. The corresponding di↵usion equation has the form of a 1. General Introduction fractional di↵erential equation [START_REF] Metzler | Subdi↵usive transport close to thermal equilibrium: From the Langevin equation to fractional di↵usion[END_REF]; [START_REF] Wyss | The fractional di↵usion equation[END_REF]]

@P (x, t) @t = 0 @ 1 ↵ t D ↵ P (x, t), (1.5)
where 0 @ 1 ↵ t denotes the factional Rieman-Liouville derivative of order (↵ 1) [START_REF] Oldham | The Fractional Calculus[END_REF]]. In general

0 @ m t f (t) = d m /dt m R t 0 ( ) 1 (t t 0 ) 1 f (t 0 ) is the fractional derivative of order m
, where m = 0, 1, 2 and 0. The convolution expresses non-Markovian memory e↵ects and for ↵ = 1 the normal di↵usion equation is retrieved.

In contrast to Eq. (1.2), P (x, t) ⌘ P (x, t|x 0 , t 0 ) is, however, the conditional transition probability of a stochastic process and not a concentration.

In the framework of the classical theory of liquids, the MSD can be expressed as [Boon & Yip [1991]]

W (t) = h(x(t) x(0)) 2 i = 2 Z t 0 dt 0 (t t 0 )C vv (t 0 ), (1.6)
where C vv = hv(t) • v(0)i is the velocity autocorrelation function of the di↵using particle.

As shown in Fig. 1.8, the MSD exhibits for short times a so-called ballistic behavior, W (t) = h|v| 2 it 2 , which corresponds to movements without any collisions. For long times, the asymptotic forms of the MSDs show either normal di↵usion, subdi↵usion or superdiffusion [START_REF] Kneller | Generalized Kubo relations and conditions for anomalous di↵usion: physical insights from a mathematical theorem[END_REF]]. Equation (1.4) is true only in asymptotic regime for t ! 1,

W (t) t!1 ⇠ 2D ↵ t ↵ , 0 < ↵ < 2.
(1.7)

Recent simulation work

Molecular dynamics simulation is a numerical method allowing to simulate the dynamics of complex liquids and solids at atomic resolution [START_REF] Allen | Computer Simulation of Liquids[END_REF]]. This computer One can in particular address biological relevant problems as protein folding and ligand binding.

Classical Molecular Dynamics Simulations

Concept and potential function

Molecular dynamics simulations are conducted by solving Newton's equations of motion for a system of N atoms on a discrete time scale,

m i d 2 r i dt 2 = F i , i = 1, 2, ..., N. (2.1)
Here F i is the acting force on the i th atom, where m i is its mass. The force is computed as the negative gradient of the potential energy V (r 1 , r 2 , ...r N ),

F i = @V i @r i (2.2)
The potential energy describes all interactions in the simulated system, which are treated in terms of classical mechanics. The equations of motions are solved for the atoms in the system and the electrons are not considered explicitly. The dynamics of electrons being much faster than one of the atomic nuclei, so one can consider that the movements of electrons instantly adjust to the movements of atom. This is the Born-Oppenheimer approximation. Moreover electrons are implicitly always in the ground states, so MD can not model processes of electron transfers or systems with electrons in excited states. It is one of the reasons why it is not possible to treat properly chemical reactions in classical MD. This is also true for all the quantum e↵ects that can appear in the system, e.g.

tunneling of protons or existence of noncovalent intermediates. Such phenomena cannot be simulated by classical molecular dynamics simulations. Despite these limitations, classical MD simulations have proven to be very successful in reproducing experimental results and in predicting new properties of the systems under consideration.

The potential energy describes so-called bonded and non-bonded interactions,

V (r 1 , r 2 , ..., r N ) = X bonds K b (b b 0 ) 2 (2.3) + X angles K ✓ (✓ ✓ 0 ) 2 (2.4) + X dihedrals V n 2 (1 + cos(n 0 )) (2.5) + X i<j A r 12 ij B r 6 ij (2.6) + X i<j q i q j ✏r ij . (2.7)
Bonded interactions result from covalent bonds between atoms and include bond-stretchings with a harmonic force with force constant K b around the equilibrium bond length, b 0 (2.3), harmonic bending of valence angles (2.4) with harmonic constants K ✓ , and displacements of dihedral angles (2.5). Here V n is the height of the energy barrier and n-th number of minima in the range (0, 2⇡). Non-bonded interactions include van der Waals (shortrange) and Coulomb interactions (long-range), which depend on the relative distances r ij between interacting atoms. Van der Waals interactions describe repulsion between atoms due to volume exclusion, as well as dipole-induced attractions between atoms. These interactions are described by the Lennard-Jones function (2.6). The Coulomb term in the potential function 2.7, represents the electrostatic interactions between atom charges q i and q j , and ✏ is the permittivity constant for the medium in which the charges are embedded.

The values for the aforementioned physical quantities (e.g K b , b 0 , ✓ 0 , q i ...) define the force-field (↵). They are fitted to experimental data and/or quantum mechanical calculations. The calculation of the potential energy with all-atom force-fields is very timeconsuming, since the computational e↵ort for the two-body interactions between free N-particles (not bonded covalently) grows as N 2 . One of the approaches to reduce the computational e↵ort is therefore to reduce the number of interaction centres. Three categories of force-fields are used today: all-atom (AA), united atom (UA) and coarse-grained (CG) ↵s. The UA and CG ↵s group several atoms into one particle. In the united atom model, one particle describes one heavy carbon atom with all aliphatic hydrogen atoms bound to it. In coarse-grained force-fields, several heavy atoms (with all hydrogen atoms bound to them) are represented as one interaction center, a so-called "bead", allowing for simulations on much larger time and length scales.

As already mentioned, Newton's equations of motion are numerically solved by a discretization scheme. One of these schemes is the leap-frog algorithm. This integrator uses the positions r i at time t and velocities v i at t 1 2 t to calculate the positions at time t + t [START_REF] Frenkel | Understanding Molecular Simulations[END_REF]] :

v i (t + 1 2 t) = v i (t 1 2 t) + t m i F i (t), (2.8)
r i (t + t) = r i (t) + t v i (t + 1 2 t).
(2.9)

The new coordinates at t + t are calculated after a leap over velocities at t + 1 2 t (therefore the name leap-frog). The algorithm is robust and time-reversible, and the integration error is of order 3 in t. The total energy can, however, not be calculated directly at a given time step, as the potential and kinetic energies are not computed at the same time. For this reason the the so-called Velocity-Verlet algorithm has been proposed.

As the Leap-Frog integrator, it integrates explicitly positions and velocities, but both are available at the same time [START_REF] Frenkel | Understanding Molecular Simulations[END_REF]] :

r i (t + t) = r i (t) + tv i (t) + t 2 2m i F i (t),
(2.10)

v i (t + t) = v i (t) + t 2m i (F i (t + t) + F i (t)).
(2.11)

In both algorithms the explicit integration of velocities can be used to simulate systems in contact with an external heat bath, which requires the scaling of the velocities.

Periodic boundary conditions and treatment of electrostatic interactions

Molecular systems in MD simulations are many orders of magnitude smaller than macroscopic systems studied by experimental methods. A direct simulation would produce unwanted artefacts, especially at the edges of the simulated systems. In order to avoid such artefacts, one can apply periodic boundary conditions (PBC).

This means an infinite replication of the simulation box in all the directions. Assuming a cubic simulation box, the images of the i-th atom have the positions

r 0 i = r i + nL = (x i + n x L, y i + n y L, z i + n z L) ( 2 . 1 2 )
where L is the dimension of the box and n = (n x , n y , n z ), with n x,y,z 2 N. An illustration id given in Fig. 2.3. During MD simulations the use of PBC just means that molecules can leave the box and enter it on the opposite side. To avoid jumps in the forces the minimum image convention is used; interactions are computed for the closest site (either a "real" atom or its image). This convention implies also use of a cut-o↵ radius for the short-range interactions, which is generally defined as half the length of the shortest linear box dimension (r cutof f  L). This approach is, however, not correct for the long-ranged Coulomb interactions, which exceed the box dimensions. Here the simulated system is The Coulomb energy of a system with periodic boundary conditions reads [START_REF] Frenkel | Understanding Molecular Simulations[END_REF]]

V electrostatic (r i ) = 1 X |n|=0 0 N X i=1 N X j=1 q i q j |r ij + nL| . (2.13)
The series above does not include the interactions between identical atoms in the central box (|n| = 0) which is indicated by the prime on the first sum. Series (2.13) does not converge and depends on the order of summation. The most common choice is spherical shell summation, which means summation over the lattice of simulation cells within a finite radius R >> L, taking the limit R ! 1. Since the system under consideration is periodic, the electrostatic potential and the charge distribution are written as Fourier series. The Fourier series for the charge distribution is split into a component which is made convergent by a spherical Gaussian damping factor and the complement which is evaluated in real space [START_REF] Ballenegger | Communication: On the origin of the surface term in the Ewald formula[END_REF]; [START_REF] Ewald | Die Berechnung optischer und elektrostatischer Gitterpotentiale[END_REF]]. An e ciency gain for the computation of the first component can be obtained by evaluating charge distribution on a mesh (grid) and computing its Fourier transform by the Fast Fourier Transform algorithm (FFT). The is the so-called particle-mesh Ewald (PME) summation [START_REF] Darden | Particle mesh Ewald: An NlogN method for Ewald sums in large systems[END_REF]; [START_REF] Deserno | How to mesh up Ewald sums. I. A theoretical and numerical comparison of various particle mesh routines[END_REF]].

Treatment of thermodynamic conditions

With the conventional molecular dynamics method systems are simulated in microcanonical conditions (NV E ensemble), i.e. with a constant number of atoms, N , a constant volume, V , and a constant total energy E. Most processes in nature happen, however, in conditions where the temperature, T , and/or pressure, p, are constant. More realistic for simulating biological systems is the canonical ensemble, where the total number of particles, the volume and the average temperature are kept constant (N, V, T = const), or the isothermal-isobaric ensemble (N, p, T = const), where the total number of particles, the average pressure and the average temperature are fixed (N, p, T = const).

In order to perform molecular dynamics simulations of a system in the NpT ensemble, Andersen and Nosé introduced the extended system method [START_REF] Andersen | Molecular dynamics simulations at constant pressure and/or temperature[END_REF]; Hoover [1985]; [START_REF] Nosé | A unified formulation of the constant temperature molecular dynamics methods[END_REF]]. In this approach scaling variables are introduced, which are additional degrees of freedom of the dynamical system and which represent, respectively, the thermostat and the barostat. By scaling the positions and the velocities of the system with these these additional variables, one can obtain the desired ensemble. The equations of motions for the isothermal-isobaric ensemble read

ṙi = V 3V r i + p i m i , (2.14) ṗi = F i V 3V p i ⇠p i , (2.15) V = p V W V , (2.16) ṗV = 1 3V n N X i=1 p 2 i m i N X i=1 r i • F i o p ⇠p V , (2.17) ⇠ = 1 W S n N X i=1 p 2 i m i + p 2 V 2W V (3N + 1)k B T o .
(2.18)

Here r i are the atomic positions, p i the associated momenta, m i is the mass of atom i, and V is the (fluctuating) volume of the simulation box. The variable ⇠ plays the role of a "friction constant", which can, however, take positive and negative values. The constants W S and W V are pseudo masses describing the inertia of the thermostat and barostat, respectively, which determine the reaction time of the system to adapt to the imposed values of pressure and temperature. As usual, k B denotes the Boltzmann constant.

Another method to regulate the temperature and pressure in the simulated system is the Berendsen coupling scheme [START_REF] Berendsen | Molecular dynamics with coupling to an external bath[END_REF]]. The algorithm imposes corrections to deviations from the reference temperature or/and pressure at each step of the simulation. The Berendsen method is very robust, but does not correspond to the usual thermodynamic NV T or NpT ensembles.

An generalization of the extended system method, which allows for pressure control in anisotropic systems, was proposed by Parrinello and Rahman [START_REF] Parrinello | Polymorphic transitions in single crystals: A new molecular dynamics method[END_REF]]. Here the components of the stress tensor instead of the (scalar) pressure are adjusted to desired values.

Chapter 3

Simulated systems and simulations

System setup

In this thesis two types of force-fields were used for molecular dynamics simulations of POPC bilayers -the all-atom Optimized Potentials for Liquid Simulations (OPLS) [Jor-

gensen & Tirado-Rives [1988]; Jorgensen et al. [1996]] and the coarse-grained MARTINI force field [START_REF] Marrink | Coarse Grained Model for Semiquantitative Lipid Simulations[END_REF][START_REF] Marrink | The MARTINI Force Field: Coarse Grained Model for Biomolecular Simulations[END_REF]]. The longest production run using the OPLS-AA ↵ was of 150 ns, given the limitations of computer power. In order to attain longer time scales and larger systems, additional simulations were performed with a coarse-grained model.

All-atom OPLS force-field

The OPLS-AA force-field was originally developed for the simulations of proteins in liquid

Coarse-grained MARTINI force-field

The MARTINI force-field [START_REF] Marrink | Coarse Grained Model for Semiquantitative Lipid Simulations[END_REF][START_REF] Marrink | The MARTINI Force Field: Coarse Grained Model for Biomolecular Simulations[END_REF]] was chosen for the extended simulations with a coarse-grained force-field. In this ↵ each molecule is built from beads, which represent each four heavy atoms with all hydrogen atoms bonded to them. Therefore one water bead equals four water molecules. (amu means atomic mass units). To compensate for the higher masses the simulations were conducted at a temperature higher by 10 K compared to the temperature for the

System equilibration

The equilibration of the system was performed according to the method proposed by [START_REF] Berger | Molecular dynamics simulations of a fluid bilayer of dipalmitoylphosphatidylcholine at full hydration, constant pressure, and constant temperature[END_REF]. The process has several stages:

1. A short 100 ps simulation in the canonical ensemble. The starting structure being far from equilibrium, a constant temperature was adjusted with the Berendsen thermostat, using a relatively short relaxation constant of 90 fs. The initial velocities of the atoms were generated according to the Maxwell-Boltzmann distribution at a temperature of 310 K.

2. Afterwards the system was coupled for 10 ps to the Berendsen barostat at 1 atm applying the pressure semi-isotropically (adjusting the pressure independently in perpendicular and vertical directions to the membrane). The relaxation time constants were set to 50 fs and 20 fs for pressure and thermal baths, respectively.

3. In the last stage, much longer isothermal-isobaric (NpT) simulations were conducted, where the system was coupled to the Nosé-Hoover thermostat and the Parrinello-Rahman barostat with the time constants of 200 fs and 500 fs, respectively. Again a pressure of 1 atm was adjusted semi-isotropically.

The last part of the equilibration process was performed until the total energy and surface area per lipid attained stable fluctuations around the desired mean values. Three consecutive simulations were conducted -50 ns, then 30 ns and lastly 150 ns (230 ns in total), before the system was considered to be in equilibrium.

Short production run -15 ns

After the equilibration phase a 15 ns production simulation was performed for further analyses. The total energy of the POPC membrane does not show any signs of decaythe energy fluctuations have the same amplitude around the mean value (Fig. 3.4). The same can be said about the box fluctuations in the X, Y and Z directions (see Fig. 3.5a).

The pressure applied during the simulation was treated isotropically in the XY-plane.

The fluctuations in the X and Y directions are thus identical. The plots show fluctuations of : a) the simulation box in all three dimensions; b) the area surface per lipid. Both are stable with the mean value of the latter of 0.638 ± 0.004 nm 2 , which indicates that the system is in equilibrium.

Extended simulation -150 ns

For comparison, a longer 150 ns simulation was performed. Here the starting structure for the calculations was the last frame of the 15 ns simulation. Fig. 3.6 shows that the total energy of the membrane is stable. This is also true for the fluctuations of the simulation box (Fig. 3.7a). The calculated mean surface area per POPC lipid, 0.637 ± 0.006 nm 2 , is very close to the value obtained for the 15 ns simulation. With the di↵erence that the mean fluctuations are nearly twice larger. Nevertheless, the system can still be considered stable, as this value is much smaller than the di↵erence of 0.03 nm 2 between the experimental surface areas obtained for 20 mN/m and 30 nM/m (Fig. 3.7b).

MARTINI force-field

The MARTINI force-field parameters chosen in this thesis are very close to the ones published by the developers of the force-field [START_REF] Marrink | Coarse Grained Model for Semiquantitative Lipid Simulations[END_REF][START_REF] Marrink | The MARTINI Force Field: Coarse Grained Model for Biomolecular Simulations[END_REF]]. The reduction of interaction centres in this ↵ allows for longer integration time-steps. In this study the integration time-step was was chosen to be 30 fs for all simulations (the proposed values are between 20 fs and 40 fs). Coulomb and van der Waals interactions were smoothly switched to zero in a distance range, r min  r  r max [START_REF] Marrink | The MARTINI Force Field: Coarse Grained Model for Biomolecular Simulations[END_REF]]. For short-range interactions r min = 0.9 nm and r max = 1.2 nm were used, and for long-range interactions r min = 0 nm and rmax = 1.2 nm. The simulations were conducted in two thermodynamic conditions -canonical (NVT) and pseudo-isothermal-isobaric (N Ap z T ) [START_REF] Marrink | The MARTINI Force Field: Coarse Grained Model for Biomolecular Simulations[END_REF]]. A constant temperature of 320 K was kept using the Berendsen thermal bath with a relaxation time of ⌧ T = 1.5 ps. In case of N Ap z T simulations the surface area of the bilayer was constant and the Berendsen barostat adjusted only the pressure of 1 atm in direction perpendicular to the surface of the membrane. The coupling time constant was ⌧ p = 3 ps and the compressibility = 4.5 ⇥ 10 5 bar 1 .

System equilibration

In comparison to the all-atom calculations the procedure for the MARTINI system equilibration was much simpler, as the starting structure used for the simulations was already equilibrated. The system equilibration process was conducted for both systems for 1.2µs in the respective ensembles, until drifts were su ciently small. The mean surface area per POPC lipid is fixed at 0.666 nm 2 , which is very close to the experimentally estimated value of 0.673 ± 0.013 nm 2 at a temperature of 323 K [START_REF] Kučerka | Fluid phase lipid areas and bilayer thicknesses of commonly used phosphatidylcholines as a function of temperature[END_REF]].

600 ns production run in the canonical ensemble

After equilibration the system was simulated for additional 600 ns. The total energy of the bilayer, which is shown in Fig. 3

600 ns production run in the N Ap z T ensemble

The simulations in the N Ap z T ensemble were performed for 600 ns, as well. The total energy of the bilayer is also reasonably stable, which is shown in Fig. 3.9.

Bulk water reference simulation

For comparison, a 100 ps simulation of SPC/E bulk water was performed in the canonical ensemble. The cubic simulation box of water of 5 nm edge length was created using GROMACS package [START_REF] Hess | GROMACS 4: Algorithms for Highly E cient, Load-Balanced, and Scalable Molecular Simulation[END_REF]] and energetically equilibrated for 1 ns prior to the production run. The starting velocities of atoms were generated according to the Maxwell-Boltzmann distribution for the temperature of 310 K. The system was coupled to the Berendsen thermostat to fix the temperature at 310 K, using a relaxation time of 90 fs. In the Fig.

C(t) = hA(0) A(t)i, (4.1)
where h• • • i denotes a time average and A is the physical quantity under consideration.

As an example, A(t) can be the velocity of the molecule,

C vv (t) = hv(0) • v(t)i. (4.2)
Here v is the velocity vector, which implies that C vv (t) is averaged over the three directions, X, Y, Z. In this thesis, all autocorrelation functions are calculated in the plane of the model membranes (lateral autocorrelation functions). It was mentioned earlier, the dynamics of the POPC lipids within the bilayers, is of main interest for this work.

The presented analyses of the OPLS all-atom and MARTINI coarse-grained simulations

Mean Square Displacements

The lateral mean square displacements, W (t) = h(x(t) x(0)) 2 i, were calculated for the centers of mass (CM) of the all-atom and coarse-grained POPC lipids, using an e cient FFT implementation [START_REF] Kneller | nMOLDYN -a program package for a neutron scattering oriented analysis of molecular dynamics simulations[END_REF]] of

W i (n) ⇡ 1 N f n N f n 1 X k=0 (x i (k + n) x i (n)) 2 . (4.3)
Here N f is the total number of frames (time-steps) of the analysed trajectories, x(n) ⌘

x(n t), and t is the integration time-step. In order to raise statistical accuracy the center-of-mass MSDs were averaged over the total number of all lipids,

W m = 1 N N X j=1 W i (n) ( 4 . 4 )
The maximum accessible lag-time of the analyses was chosen to be shorter than 10%

of the total length of the MD trajectories. This choice was made in order to guarantee the statistical relevance of the results. For longer lag times the time averaging of the calculated functions might prove insu cient to obtain reliable results. Therefore the maximal lag-time for the 15 ns OPLS-AA simulation was 1.5 ns, 15 ns for the 150 ns OPLS-AA simulation and 50 ns for the MARTINI calculations. [ 2011]]. This comparison shows that the characteristics of lateral di↵usion in these biological membranes is the same. The fitted fractional di↵usion constants D ↵ for POPC have values of 0.016 nm 2 ns ↵ (150 ns) and 0.018 nm 2 ns ↵ (15 ns). In order to compare obtained fractional di↵usion constants with the value for DOPC (D ↵ = 0.101 nm 2 ns ↵ ), one can calculate the typical time scale ⌧ vacf , which quantifies the asymptotic regime for VACF and depends on D ↵ and ↵ parameter (more details are given in next section). The typical time scale ⌧ vacf is related to friction in a system as ⌧ vacf = 1 . Therefore the larger ⌧ vacf the faster observed di↵usion. For the POPC simulations calculated ⌧ vacf = 0.065 ps and ⌧ vacf = 0.046 ps for 15 ns and 150 ns simulations, respectively, are smaller than for the DOPC ⌧ vacf = 0.345 ps, which shows the time scale of subdi↵usive motions for the both systems is clearly di↵erent. The lateral di↵usion of DOPC lipids in the bilayer is thus several times faster than for POPC membrane. This can be explained by the fact that the latter is less fluid at the same temperature -its main phase transition from gel to liquid, occurs at 270 K, while the transition for DOPC happens already

Coarse-grained force-field results

I recall that the mass of one lipid is M AA = 760.08 amu, while the total mass of coarsegrained POPC lipid, M CG , is 936 amu ("amu" : means atomic mass units). Obviously motions of heavier lipid molecules would be slower in comparison to lighter lipids and hence have a reduced fractional di↵usion constant, D ↵ , which is proportional to the amplitude of the velocity autocorrelation function, C vv [START_REF] Kneller | Generalized Kubo relations and conditions for anomalous di↵usion: physical insights from a mathematical theorem[END_REF]],

D ↵ = 1 (↵ + 1) Z 1 0 dt 0 @ ↵ 1 t C vv (t), (4.5)
where

0 @ ↵ 1 t C vv (t) = d/dt R t 0 dt 0 (2 ↵) 1 (t t 0 ) 1 ↵ C vv (t 0
) is the factional Rieman-Liouville derivative of order (↵ 1). As C vv (0) is the mean squared velocity, h|v| 2 i = 2k B T /m, the reduction of the D ↵ is of around 20% ( M AA M CG ⇡ 0.8). In order to compensate for that, the simulations with the MARTINI force-field were conducted at a higher Table 4.2: Experimental values of the di↵usion constant of POPC lipids for di↵erent temperatures [START_REF] Köchy | Lateral di↵usion coe cients of phospholipids in spherical bilayers on a solid support measured by 2H-nuclear-magnetic-resonance relaxation[END_REF]]. 4.2 [START_REF] Köchy | Lateral di↵usion coe cients of phospholipids in spherical bilayers on a solid support measured by 2H-nuclear-magnetic-resonance relaxation[END_REF]], the rise of 10 K (10 C), from 310 K to 320 K, results in estimated increase of the di↵usion constant by ⇡ 23%. Therefore MARTINI simulations were performed at a temperature of 320 K. The lateral mean square displacements for both MARTINI simulations (with canonical, NVT, and pseudo-isothermal-isobaric, N Ap z T , thermodynamic conditions) were calculated for maximal time-lags of 50 ns, which is less than 10 % of the total trajectory length In comparison to the results from the fits for 50 ns, the fractional di↵usion constants obtained from the fits for a maximum time lag of 1.5 ps are slightly higher, but the values of ↵ are smaller, meaning that for the longer time-lags the motions of molecules are "less normal". The calculated mean ⌧ vacf = 0.44 ps and ⌧ vacf = 0.30 ps for the shorter 1.5 ns) and longer (50 ns) maximal time-lag, respectively, are very similar, which shows that analyses for shorter time-lags as short as 1.5 ns are already representative. seems to be slightly more linear than for the coarse-grained case. This observation is reflected in a higher value of ↵ parameter (see Table 4.4), implying more "normal" diffusional motions of all-atom POPC lipids. In all cases the fractional di↵usion constants (D ↵ ⇡ 0.05 nm 2 ns ↵ ) and ⌧ vacf = 0.44 ps and ⌧ vacf = 0.30 ps for the shorter (1.5 ns) and longer (50 ns) maximal time-lags are several times higher than the values for the all-atom simulations (D ↵ = 0.018 nm 2 ns ↵ ) with ⌧ vacf = 0.065 ps. This result is consistent with other findings showing that calculated di↵usion constants from the MARTINI simulations are in most cases several times higher than experiments predict [START_REF] Marrink | Coarse Grained Model for Semiquantitative Lipid Simulations[END_REF][START_REF] Marrink | The MARTINI Force Field: Coarse Grained Model for Biomolecular Simulations[END_REF]].

Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡
Moreover the scaling value of 4 has been mentioned in these works in case of motions of MARTINI water beads, where it was shown that di↵usion constant is 4 times higher than the experimental value. This di↵erence can be attributed to less detailed structure and therefore smaller number of individual interactions between them. This could lead to less hindered movements, so faster lateral movements of POPC lipids.

The above results show, however, also that the dynamics obtained by the coarse-grained simulations with the MARTINI force-field, preserves the subdi↵usive character of the lipid motions.

It is worth noting that a transition from anomalous to normal di↵usion was reported by Flenner et al. for the all-atom and the 600 ns coarse-grained POPC lipids simulations, respectively. A possible explanation for these discrepancies is the statistical error, which is induced by using appropriate maximal lag-times.

Bulk water comparison with membrane hydration water

In order to verify if the simulated system of bulk water reproduces the behavior of a real water system, in particular normal di↵usion, it is compulsory to compare the di↵usion constant obtained from the simulation with an experimental value. Although this subject was not in the focus of this thesis, some short analysis can be given on the water motions in studied POPC systems in this thesis. The hydration level of atom detailed OPLS POPC bilayer and coarse-grained MARTNI membrane is ⇡ 38 water molecules and ⇡ 28 water beads per POPC lipid, respectively. The reported number of water molecules forming hydrogen bridges vary with the type of lipids, but in case of fully hydrated POPC it may include up to average ⇡ 9.4 molecules per one POPC lipid [START_REF] Gawrisch | Hydration of POPC bilayers studied by 1H-PFG-MAS-NOESY and neutron di↵raction[END_REF]]. This means that ⇡ 25% of all water molecules surrounding simulated bilayers are potentially involved in forming water bridges. This should have an influence on the average di↵usion constant of the simulated water molecules (alternatively water beads). Figure 4.5 shows plots of the mean square displacements for the centre of mass of SPC/E water from 15 ns OPLS-AA simulations (Fig. 4.5a) and for water beads of both, NVT and N Ap z T 600 ns MARTINI simulations (Fig. 4.5b). The calculated di↵usion constants of 1.721•10 9 m 2 s (1.721 nm 2 ns ) (OPLS-AA) and 1.399•10 9 m 2 s (1.399 nm 2 ns ) (MARTINI) are indeed smaller than for the bulk SPC/E water -3.528•10 9 m 2 s (3.528 nm 2 ns ), i.e. by a factor of ⇡ 2 (Table 4.5). Table 4.5: Comparison of the di↵usion constants of a bulk and surrounding lipid bilayer SPC/E water, water in the MARTINI POPC simulations and the experimental value [START_REF] Mills | Self-di↵usion in normal and heavy water in the range 1-45 C[END_REF]]. WATER bulk SPC/E SPC/E (mem.) CG (mem.) exper. D (10 9 m 2 /s) 3.528 1.721 1.399 2.919

Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 10
This indicates that in the simulated systems presented in this thesis, similarly to earlier findings, exist strong water interactions involving water molecules/beads and POPC lipids in model membranes. Very interesting is the fact that di↵usion constants for water beads in MARTINI systems (which are exactly the same for both thermodynamic conditions) is smaller than for all-atom case. This would suggest a much stronger e↵ect of waterbridging (especially when the reported di↵usion constant of bulk MARTINI water is ⇡ 4 higher than experimental result [START_REF] Marrink | The MARTINI Force Field: Coarse Grained Model for Biomolecular Simulations[END_REF]]). This matter would need more further investigation to give a more coherent view of the problem.

Velocity autocorrelation function analysis

For the following considerations we need that the MSD and VACF are related via [Boon & Yip [1991]]

W (t) = 2 Z t 0 (t t 0 )C vv (t 0 )dt 0 , (4.6)
It follows now from Eq. (1.7) that its Laplace transform behaves as [START_REF] Kneller | Generalized Kubo relations and conditions for anomalous di↵usion: physical insights from a mathematical theorem[END_REF]]

Ŵ (s) s!0 ⇠ 2D ↵ (↵ + 1) s ↵+1 , (4.7)
and on the other hand it follows from Eq. (4.6) :

Ŵ (s) = 2 Ĉvv (s) s 2 , (4.8)
For an arbitrary function f (t) the Laplace transform is defined as f (s) = R 1 0 dt exp( st)f (t) (<{s} > 0). The relation between asymptotic form of MSD and VACF can be obtained by comparing Eqs. (4.7) and (4.8), which leads to [START_REF] Kneller | Generalized Kubo relations and conditions for anomalous di↵usion: physical insights from a mathematical theorem[END_REF]]:

Ĉvv (s) s s!0 ⇠ D ↵ (↵ + 1)s ↵ , (4.9)
where Ĉvv (s)/s is the Laplace transform of the integral

f c (t) = Z t 0 dt 0 C vv (t 0 ). (4.10)
Since Eq. (4.9) is singular for any ↵ > 0, as s approaches zero, one obtains the equivalence fc (s)

s!0 ⇠ D ↵ (↵ + 1)s ↵ , f c (t) t!1 ⇠ ↵D ↵ 1 t ↵ 1 (4.11)
and finally [START_REF] Kneller | Generalized Kubo relations and conditions for anomalous di↵usion: physical insights from a mathematical theorem[END_REF]]

W (t) t!1 ⇠ 2D ↵ t ↵ , f c (t) t!1 ⇠ ↵D ↵ 1 t ↵ 1 . (4.12)
Since C vv (t) = df c (t)/dt, one obtains from Expression (4.12) the necessary condition

W (t) t!1 ⇠ 2D ↵ t ↵ ) C vv (t) t!1 ⇠ D ↵ ↵(↵ 1)t ↵ 2 (4.13)
for the long-time tail of the VACF. I mention here that Expression (4.5) for the fractional di↵usion constant follows from Eq. (4.9).

The asymptotic regime can be quantified as t ⌧ vacf , where [START_REF] Kneller | Communication: Consistent picture of lateral subdi↵usion in lipid bilayers: Molecular dynamics simulation and exact results[END_REF]]:

⌧ vacf = ⇣ D ↵ hv 2 i ⌘ 1 2 ↵ . (4.14)
To facilitate comparisons in the time domain, the VACFs in this thesis are presented normalized,

(t) = C vv (t) C vv (0) (4.15)
and the long-time tail of the integral (4.10) behaves then:

f (t) t!1 ⇠ D ↵ ↵ hv 2 i t ↵ 1 . (4.16)
The velocity autocorrelation functions were computed for the center of masses (CM) of the POPC lipids according to

C i (n) ⇡ 1 N f n N f n 1 X k=0 (v i (k) • v i (k + n)), (4.17)
where N f is the total number of frames (time-steps) of the analysed trajectories. Here

x(n) ⌘ x(n t) and t is the integration time-step. For reasons of e ciency, Expression (4.17) is computed by a FFT-based Fast Correlation Algorithm [START_REF] Kneller | nMOLDYN -a program package for a neutron scattering oriented analysis of molecular dynamics simulations[END_REF]].

Similarly to the computation of the MSDs, Eq. (4.4), all the VACFs were averaged over the total number of all the lipids, N :

C m = 1 N N X j=1 C i (n) ( 4 . 1 8 )
The integration of the numerical normalized VACFs of the timestep, t, was performed by using the central di↵erence scheme [START_REF] Stachura | Anomalous lateral di↵usion in lipid bilayers observed by molecular dynamics simulations with atomistic and coarse-grained force fields[END_REF]],

f (n) ⇡ 0<n<Nt 1 X k=0 t 12 (8 (k) + 5 (k 1) (k + 1)). (4.19)
The velocity autocorrelation functions for the 15 ns all-atom OPLS and 600 ns coarsegrained MARTINI (NV T and N Ap z T ) simulations are presented in Fig. 4.6. The VACF curves for the MARTINI case are identical for both ensembles. This is consistent with the observation of the MSDs for short time-lags of maximum 1.5 ns (section 4.1.1.2), which were also indistinguishable for the NV T and N Ap z T conditions. One distinguishes two regimes in the VACF: a short-time regime, containing the minimum of the VACF, which extends up to ⇡ 3ps, and a long-time regime. In case of the OPLS-AA ↵ the minimum is deeper and starts earlier than for the MARTINI ↵. This suggests a weaker "cage e↵ect"

for the latter and faster lateral motions of the coarse-grained lipids. This is true in light of the MSD calculations, for which the fractional di↵usion constants are higher than for the atom detailed simulations. In addition the VACFs of MARTINI simulations have a "bump" around 0.7 ps. The origins of it have not been discerned and do not have any theoretical explanation. More important is the long-time tail of the VACFs, which for both, the all-atom and the coarse-grained calculations, get closer to each other for longer time-lags and both approach asymptotically zero from negative values. Negative velocity autoccorelations of the molecules can be expressed in physical terms, namely as the tendency of reversing the direction of motion and staying localized for longer times ("cage e↵ect"). The calculated typical time scales ⌧ vacf , Eq. (4.14), of 0.065 ps for OPLS-AA and 0.293 ps for MARTINI (the mean squared velocities hv 2 i = 0.0057 nm 2 ps 2 and hv 2 i = 0.0068 nm 2 ps 2 for the coarse-grained and all-atom lipids, respectively) define the asymptotic regime, which stands for t

Ê Ê Ê Ê Ê Ê Ê ÊÊ Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê ÊÊÊÊ ÊÊ Ê Ê Ê Ê ÊÊ ÊÊ Ê Ê Ê Ê Ê ÊÊ Ê Ê
C vv (t) = lim ⌧ !1 1 ⌧ Z ⌧ 2 ⌧ 2 dt v(t) v(t + t 0 ) ( 4 . 2 2 )
According to the Wiener-Khintchine theorem [Boon & Yip [1991]], the Fourier transformed VACF, Cvv (!), in limit of ⌧ ! 1, equals the mean squared velocity in !-domain, which is proportional to the frequency spectrum of the kinetic energy,

Cvv (!) = lim ⌧ !1 |ṽ(!)| 2 ⌧ / E k (!). (4.23)
In order to obtain an analytical expression for the DOS in the low frequency region, Eq. (4.9), and Cvv (!) = <{ Ĉvv (i!)} are used, where Cvv ⌘ g(!) is the cosine transform of C vv (t) and Ĉvv (s) its Laplace transform,

g(!) = 1 2 Cvv (!) !!0 ⇠ D ↵ 2 (↵ + 1) ! 1 ↵ sin ⇣ ⇡ ↵ 2 ⌘ . (4.24)
In case of normal di↵usion, ↵ = 1, the model derived above, tends for ! ! 0 to the di↵usion constant. In general:

D ↵ = lim !!0 g(!) ! (1+↵) (↵ + 1) sin ⇣ ⇡ ↵ 2 ⌘ . (4.25)
This alternative expression for D ↵ is certainly easier to use in practice than Expression (4.5). The Fourier spectra of ensemble averaged MARTINI and OPLS-AA VACFs, Eq. (4.18) were computed only for the 15 ns all-atom and 600 ns N Ap z T coarse-grained simulations (corresponding VACF are nearly identical for both NV T and N Ap z T ensembles).

In order to retrieve the fractional di↵usion constant from g(!), one has to work with the unnormalized VACFs, where C vv (0) = hv 2 i. Here for the MARTINI POPC lipids hv 2 i = 0.00568 nm 2 ps 2 and hv 2 i = 0.00678 nm 2 ps 2 for the all-atom lipids.

For numerical calculations, the continuous Fourier transform

Cvv (!) = Z +1 1 dt exp( i!t)C vv (t), (4.26)
may be computed via the discrete Fourier transform (DFT) [START_REF] Brigham | Fast Fourier Transform[END_REF]]: where Ĉvv (k) ⌘ C vv (k !) and ! = 2⇡/( t(2N 1)). With increasing number of points, for N ! 1, the Fourier transform of VACF becomes more and more noisy. For this reason the Cvv (!) is smoothed by a window function w(t) and then the Fourier spectrum is calculated according to the following integral:

Ĉvv (k) = N 1 X n= (N 1) exp ✓ 2⇡i kn 2N 1 ◆ C vv (n). ( 4
C(w) vv (!) = Z +1 1 dt exp( i!t)w(t)C vv (t). (4.29)
The smoothing function w(t) is positive and fulfills the conditions,

w(0) = 1, (4.30)
w(t) = w( t), (4.31)
where the first condition conserves the integral of the spectrum

Z +1 1 d! Cw vv (!) = Z +1 1 d! Cvv (!). (4.32)
According to the convolution theorem of the Fourier transform [START_REF] Brigham | Fast Fourier Transform[END_REF]] the smoothed Fourier transform is then

Cw vv (!) = 1 2⇡ Z +1 1 d! 0 w(! ! 0 ) Cvv (! 0 ). (4.33)
In this thesis the Fourier transform is smoothed with a Gaussian function [START_REF] Harris | On the use of windows for harmonic analysis with the discrete Fourier transform[END_REF]]:

w(t) = exp 1 2  t 2 ! ! w(!) = p 2⇡ exp ✓ 1 2 [ !] 2 ◆ , ( 4
Ĉw vv (k) = N 1 X n= (N 1) exp ✓ 2⇡i kn 2N 1 ◆ w(n)C vv (n). (4.36)
In this thesis the densities of states were calculated with = 10% of the total length of the MD trajectories.

Figure 4.8 shows the comparison between the computed g(!) for the coarse-grained (main plot) and the all-atom simulations (inset). In case of the latter the peak is less pronounced and occurs at ! ⇡ 3.5 T Hz, whereas for the MARTINI simulation the peak is at ! ⇡ 1.5 T Hz. Moreover, g(!) is broader for the all-atom simulation.

The fits of Relation (4.24) to the calculated g(!) for small angular frequencies are presented in the Fig. 4.9. The choice of the region that would be considered as for small ! is not clearly defined. The final choice of the maximal value of ! defining this region was taken to be smaller than 10% of the calculated maximal frequency (! = 20 T Hz), with a condition that the small frequencies part refers only to ! smaller than the frequency for which the peak in g(!) occurs (! 2 (0, 1.5 T Hz) for the coarse-grained calculations and ! 2 (0, 3.5 T Hz) for all-atom case). Following this reasoning the fitted parts of the g(!) were up to ! = 0.6 T Hz and ! = 2.0 T Hz for the MARTINI and OPLS-AA ↵s, respectively. The fractional di↵usion constants obtained from the fits have values of 0.059 nm 2 ns ↵ for the MARTINI and 0.020 nm 2 ns ↵ for the OPLS-AA simulations, which are slightly larger than the results from the MSDs fits (Tab. 4.6). The obtained ↵ parameters of 0.437 and 0.467 for the all-atom and the coarse-grained DOS, respectively, are considerably smaller than their MSD counterparts. Especially in case of the former, the ↵ parameter is ⇡ 30% lower, which suggests a stronger subdi↵usional behavior of the lipids than the MSD results indicate. This observation is also true for the MARTINI calculations, for which the ↵ parameter has a very similar value to the ↵ for the OPLS-AA. The calculated ⌧ vacf = 0.29 ps for the OPLS-AA simulation is much higher than the ones computed for the MSD D ↵ and ↵ parameters (⌧ vacf = 0.065 ps), which implies observation of faster di↵usion in the DOS. But, as expected, is smaller than for the ⌧ vacf = 0.50 ps calculated for the MARTINI simulation, which is also higher than the mean ⌧ vacf = 0.30 ps obtained from the MSD. The higher ⌧ vacf for the coarse-grained model than for the all-atom case means smaller friction in the system, according to = (⌧ vacf ) 1 . The values of = 3.46 ps 1 for the OPLS-AA simulation is higher than for the MARTINI force-field ( = 2.01 ps 1 ), which corresponds to broader peak for the DOS in case of the former.

Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê 0.0 0.
Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê

Density of the States from the Autoregressive Model

A discrete approximation of the continuous Fourier transform (4.26), which does not assume periodicity in time, is the infinite sum [START_REF] Papoulis | Signal Analysis[END_REF][START_REF] Papoulis | Probability, Random Variables, and Stochastic Processes[END_REF]]

Cvv (!) = t +1 X n= 1 C vv (n) exp( i!n t), (4.37)
where C vv (n) ⌘ C vv (n t) and t again the sampling interval of the VACF. Then the z-transform of C vv (n), defined as [START_REF] Papoulis | Signal Analysis[END_REF][START_REF] Papoulis | Probability, Random Variables, and Stochastic Processes[END_REF]]

C z vv (z) = +1 X n= 1 C vv (n)z n , (4.38)
is related as follows to Expression (4.37)

Cvv (!) = t C z vv (exp(i! t)). (4.39)
Concurrently, the velocity autocorrelation function can be considered in terms of a sample

v N (n) of v(t), C vv (n) ⌘ lim N !1 1 2N 1 N 1 X k= (N 1) v N (n + k)v N (k), (4.40)
where the length of the sample is 2N 1 and it follows from the convolution theorem of the z-transform, (the asterisk denotes a complex conjugate)

(f g)(n) ! F (z)G ⇤ (1/z ⇤ ) ( 4 . 4 1 ) that C z vv (z) = lim M !1 1 2N 1 V N (z)V ⇤ N (1/z ⇤ ). (4.42)
The z-transform of the an autocorrelation function -here the VACF -can be easily calculated if one assumes an autoregressive model for the time series of the underlying time series [START_REF] Papoulis | Signal Analysis[END_REF][START_REF] Papoulis | Probability, Random Variables, and Stochastic Processes[END_REF]]

v(n) = M X k=1 a k v(n k) + R(n). (4.43)
Here M is the order of the AR process, a k are constant coe cients and R(n) is white noise with zero mean and autocorrelation hR(n)R(k)i = 2 M nk . The amplitude of the noise R(n) can be determined by using the fact that R(n) is not correlated with past values of

v(n), hR(n)v(n k)i = 0, 2 M = C vv (0) M X k=0 a k C vv (k), (4.44)
Since R(n) is white noise, it is not correlated with the velocity signal for di↵erent times.

This property can be used to derive the following set of linear equations for the coe cients Performing now the z-transform of Eq. (4.43) and using the parameters for the autoregressive model, one obtains

V N (z) = R(z) 1 P M k=1 a k z k , (4.46)
where R(z) is the z-transform of R(n). In the next step the expression above is inserted to Equation (4.42). Using that lim

N !1 1 2N 1 R(z)R ⇤ (1/z ⇤ ) = C rr (z) = 2 . (4.47)
the z-transform of the VACF is found to be [START_REF] Kneller | Computing memory functions from molecular dynamics simulations[END_REF]]:

C z vv (z) = 2 M (1 P M k=0 a k z k )( P M l=0 a k z l ) (4.48)
Setting z = exp(i! t), the AR density of states becomes

g AR (!) = t 2 C z vv exp(i! t), (4.49)
The coe cients of the autoregressive process of order M for the estimation of the densities of states, g AR (!), were obtained by solving the linear Yule-Walker equations (4.45) [START_REF] Walker | On Periodicity in Series of Related Terms[END_REF]; [START_REF] Yule | On a Method of Investigating Periodicities in Disturbed Series, with Special Reference to Wolfer's Sunspot Numbers[END_REF]] with order M = 666 for the analysis OPLS-AA simulations and with M = 3333 for MARTINI simulations.

The estimated densities of states obtained from the AR model are very similar to those calculated by discrete Fourier transform (see Fig. 4.10). They show in particular the same di↵erences between the DOS of the coarse-grained and the all-atom simulations (Fig. 4.11).

In case of the latter the peak is also shifted to ⇡ 3.5 T Hz (for CG ⇡ 1.5 T Hz). In transform method. Nevertheless, ⌧ vacf = 0.33 ps and = 3.04 ps 1 for the OPLS-AA force-field and ⌧ vacf = 0.50 ps and = 1.99 ps 1 for the MARTINI force-field obtained from the AR model are very close to the FT calculations, which shows that both methods give nearly identical results.

Ê Ê Ê Ê Ê Ê Ê Ê ÊÊ Ê ÊÊ Ê ÊÊ Ê Ê Ê Ê Ê Ê Ê Ê ÊÊ Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê ÊÊ Ê Ê ÊÊ ÊÊ ÊÊ Ê Ê Ê 0.0 0.
Ê ÊÊÊ Ê Ê Ê Ê Ê Ê Ê Ê Ê ÊÊ Ê Ê Ê Ê Ê Ê Ê ÊÊ Ê Ê Ê Ê Ê Ê ÊÊ Ê Ê Ê ÊÊ Ê ÊÊ Ê Ê Ê Ê ÊÊÊ ÊÊÊ ÊÊ ÊÊ Ê Ê Ê D a =0

Conclusions

In this chapter a single molecular dynamics of the OPLS all-atom and MARTINI coarsegrained POPC lipids was presented. The MD trajectories were analysed in terms of the mean square displacements, velocity autocorrelation functions and the densities of states.

All the proposed models clearly show that the movements of POPC lipids within model membrane exhibit a subdi↵usive characteristics, within the accuracy of the analysis. The analyses of the MSDs predict an ↵ parameter of ⇡ 0.7 and ⇡ 0.6 for the OPLS-AA and MARTINI simulation, respectively.The analyses of the DOS indicate a stronger subdi↵usive characteristics, with ↵ parameter of ⇡ 0.4. The independent estimation of DOS using autoregressive process confirms the results for g(!), obtained by the Fourier transform of the velocity autocorrelation function.

Collective motions of lipids

In contrast to the earlier chapter, which covered the single molecular dynamics, this chapter concentrates solely on the collective motions of POPC lipids in the plane of model membrane. The main goal is to visualize the cage of nearest neighbors of a tagged di↵using particle. This subject will be presented in terms of spatio-temporal density fluctuations within the simulated POPC membranes.

Pair Correlation Function

Structural and dynamical correlations in condensed matter systems are described by the van Hove correlation function [Van Hove [1954]]

g (V H) (r, t) = 1 N N X ↵, =1 h (r [R (t) R ↵ (0)])i (4.50)
where R ↵ (t) is the position of particle ↵, which is the center of mass of a lipid molecule in the context of this thesis. The van Hove correlation function plays a fundamental role in the interpretation of neutron and light scattering experiments on liquids and colloids [START_REF] Bée | Quasielastic Neutron Scattering: Principles and Applications in Solid State Chemistry[END_REF]; [START_REF] Berne | Dynamic Light Scattering[END_REF]; [START_REF] Lovesey | Theory of Neutron Scattering from Condensed Matter[END_REF]]. It splits into a self-and a distinct part,

g (V H) (r, t) = g (V H) S (r, t) + g (V H) D (r, t), (4.51)
where

g (V H) S (r, t) = 1 N N X ↵=1 h (r [R ↵ (t) R ↵ (0)])i, (4.52) g (V H) D (r, t) = 1 N N X ↵6 = =1 h (r [R (t) R ↵ (0)])i. (4.53)
The distinct part describes is of interest in the following since it quantifies the structural correlations in the system under consideration. One verifies easily that its spatial average over the volume V of the system yields the mean particle density, ⇢ = N/V , for large N , This motivates the definition of the time-dependent pair correlation function through the first peak vanishes already after 1.5 ps (dotted line) -but also the character of this decay is di↵erent than for the POPC lipids, which was already visible in Fig. 4.17. It is worth noticing that the computed pair distribution function for t = 0 (solid line), matches very well known molecular dynamics results [START_REF] Rahman | Molecular dynamics study of liquid water[END_REF]; Sciortino et al. correlation function for a bulk water [START_REF] Narten | Liquid Water: Molecular Correlation Functions from X-Ray Di↵raction[END_REF]].

1 V Z V dV g (V H) D (r, t) = N 1 V ⇡ ⇢. ( 4
G D (r, t) = 1 ⇢ 1 V r Z Vr dV g (V H) D (r, t), ( 4
Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê ÊÊ ÊÊÊ ÊÊ Ê ÊÊ Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê 0.1 0.2 0. 5
In order to further investigate the nature of motions in the all-atom POPC bilayer, the time decay of the first shell of neighbors was analysed, by taking into account only their long-time behavior. The insets in Fig. respectively. The parameters are, however, nearly twice smaller than for the all-atom simulations and similar for both thermodynamic ensembles. The choice of the time range for the fits of Expression (4.58) is, of course, debatable. Figure 4.21 shows the dependence of the -parameter on the initial time t start for the fit interval. The fitted values are relatively stable for t start up to 70 % and to 90 % of the maximal analysed time-lags for the OPLS-AA and MARTINI simulations, respectively.

Conclusions

In this chapter the analysis of the local structure dynamical structure around each lipid molecule in a POPC membrane was analyzed by analyzing the behavior of the lateral time-dependent pair correlation functions. They reveal a persisting cage of the first shell of neighboring lipids (for over 15 ns in case of the 150 ns OPLS-AA simulation), which are consistent with the observed negative long-time tails of the calculated VACFs. The analyses of the long-time tails of the first peaks of the G D (r, t), exhibit a power-law decay, with ⇡ 0.55 and ⇡ 0.25 for the all-atom OPLS and the coarse-grained MARTINI simulations, respectively. These results contrast strongly with findings for the bulk SPC/E water, for which the first shell of neighbors vanishes rapidly (already after 1.5 ps) and, moreover, exponentially.

The existence of the cage e↵ect in POPC lipid membrane for long correlation times is in good agreement with molecular dynamics simulations findings of concerted movements in lipid bilayers, forming so-called "collective flow patterns" [START_REF] Falck | Lateral Di↵usion in Lipid Membranes through Collective Flows[END_REF]].

The connection between the ↵-parameter describing the long-time behavior of the MSD, Eq. (1.4) and the -parameter describing the slow decay in the first peak of G D (r, t) is not yet established.

Chapter 5 Conclusions

The main goal of this thesis was to characterize the nature of the lateral movements of lipids within a model membrane. For this purpose a series of molecular dynamics simulations of a 1-palmitoyl-2-oleoylphosphatidylcholine (POPC) lipid bilayer was conducted with the all-atom (OPLS) and the coarse-grained (MARTINI) force fields. An all-atom energetically equilibrated POPC membrane of 274 lipids in total was simulated for 15 ns and 150 ns in full hydration. The equilibration process had several intermediate steps, which account for the total simulation time of 230 ns, and the 150 ns simulation was an extension to the 15 ns one. The systems were considered to be in equilibrium when the fluctuations of the total energy and the average area per lipid were stable, with the value of the latter being close to the experimental ones. In order to have access to longer time scales and larger system sizes, simulations of 600 ns for a system of 2033 lipid molecules were conducted with the MARTINI force field, using the NV T and the N Ap z T ensemble.

In the latter ensemble the pressure was adjusted only in the direction perpendicular to the bilayer. The production runs for both ensembles were performed after energy equilibration simulations of 1.2 µs. The MARTINI force-field calculations were also used to test the ability of coarse-grained models to reproduce dynamical properties of all-atom simulations.

The lateral di↵usion of the lipid molecules was analyzed using di↵erent physical quantities and numerical methods. For all simulations, the time evolution of the computed molecular mean square displacements (MSDs) showed clear deviations from the asymptotic linear form which characterizes normal di↵usion. The fits of the relation for the MSD gave average ↵ parameters of ⇡ 0.69 and ⇡ 0.55 for the OPLS-AA and MARTINI simulations, respectively. The observed di↵usion of POPC lipids was 6-8 times faster for the MARTINI simulations than for the OPLS-AA, which is also coherent with publications showing that in general the di↵usion constants obtained from the MARTINI simulations are several times higher than the experimental values. The MSD results are consistent with the analyses of the VACF power spectra, which were calculated with two di↵erent methods: (a) by a windowed discrete Fourier transform of the VACF, and (b) by fitting an autoregressive model to the underlying velocity time series. For both OPLS-AA and MARTINI, the corresponding ↵-parameters were, however, slightly smaller than those obtained from the analyses of the MSDs. Similarly to the MSD analysis, the di↵usion of the coarse-grained lipids was faster than the di↵usion of the all-atom ones. The so-called cage e↵ect, which manifests itself by negative algebraic long-time tails in the VACFs, was further analyzed by examining the time-dependent local structure around each lipid molecule by means of the time-dependent pair correlation function. The first peak of the latter exhibits a slow power-law decay, indicating that the local cage formed by the shell of nearest neighbors decays very slowly. This is a direct physical manifestation of the cage e↵ect, which is suggested by the negative long-time tails of the VACF. The persistence of the local structure is in agreement with flow-like dynamical patterns in membranes, which have been observed recently by other authors [START_REF] Falck | Lateral Di↵usion in Lipid Membranes through Collective Flows[END_REF]]. The comparison of the simulations with the all-atom OPLS and the coarse-grained MARTINI force field showed clearly that both exhibit subdi↵usion, with similar ↵-parameters, but the latter is several times faster. In agreement with these findings, the molecular VACFs exhibit a more shallow minimum and the cage of nearest neighbors of the lipid molecules decays faster. These observations can be attributed to the fact that the coarse-grained representation of the lipid molecules in the MARTINI force field leads to less interactions with the neighbors and thus to less entanglement.

Figure 1 :

 1 Figure 1: La complexité d'une membrane biologique 1 .

Figure 2 :

 2 Figure 2: Comparaison entre les modèles tout-atome OPLS (AA) et gros-grain MARTINI (CG) de membrane de POPC.

Figure 3 :Figure 4 :

 34 Figure 3: Les fits de l' éq. 1 de déplacement carré moyen pour les simulations: a) 15 ns OPLS-AA, b) 150 ns OPLS-AA, c) NV T et N Ap z T MARTINI pour des temps maximales de 50 ns et d) NV T et N Ap z T MARTINI pour des temps maximales de 1.5 ns.

Figure 5 :

 5 Figure 5: Les fonctions de corrélation de paire distincte dépendantes du temps, G D (r, t) pour les simulations: a) 15 ns OPLS-AA, b) 150 ns OPLS-AA, c) NV T MARTINI et d) et N Ap z T MARTINI. Dans les encarts sont montrés les fits de des temps longs.

4 1. 3

 43 ns OPLS-AA, b) 150 ns OPLS-AA, c) NV T et N Ap z T MARTINI pour des temps maximales de 50 ns et d) NV T et N Ap z T MARTINI pour des temps maximales de 1.5 ns. . . . . . . . . . . . . . . . . . . . . . . . . . . ix 4 Comparaison entre les fonctions d'autocorrélation des vitesses pour les simulations de OPLS-AA et MARTINI. . . . . . . . . . . . . . . . . . . . . . . x 5 Les fonctions de corrélation de paire distincte dépendantes du temps, G D (r, t) pour les simulations: a) 15 ns OPLS-AA, b) 150 ns OPLS-AA, c) NV T MARTINI et d) et N Ap z T MARTINI. Dans les encarts sont montrés les fits de des temps longs. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xi 6 La fonction de corrélation de paire distincte dépendantes du temps, G D (r, t), pour la simulation de l'eau SPC/E. . . . . . . . . . . . . . . . . . . . . . . xii 1.1 Illustration 1 showing the complexity of a membrane in cells and the structure of a lipid molecule. [1.3cm] . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2 Schematic 2 of the phosphatidylcholine chemical structure. . . . . . . . . . Sketch of the POPC chemical structure 3 . . . . . . . . . . . . . . . . . . . 6 xvii 1.4 Crowded cell compartments visualized by cryoelectron microscopy 4 Actin filaments are shown in orange, ribosomes and other macromolecular assemblies in grey, and membrane structures in blue. . . . . . . . . . . . . . . . . 1.5 Illustration of a chromosome 5 , depicting its parts; the endings of each of the four arms(3,4), called telomeres (1) and the centromere (2) connecting the two chromatides. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.6 Visualization of ion tra cking in cells involving Golgi apparatus and rectifier potassium channels (inward and delayed) 6 . The ion channels are synthesized in endoplasmic reticulum and through Golgi are transported to the plasma membrane. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.7 Recorded di↵usion of putty grains published by Perrin [1909]. . . . . . . . 1.8 MSD of a model system for various values of ↵ [Kneller [2011]]. . . . . . . 1.9 Plot representing the displacements of the DPPC lipids in x and y planes for the upper leaflet (a,b,c) and lower leaflet (d); (a) shows a 50 ps interval from 5.00 to 5.05 ns of the simulation, (b) a 500 ps interval from 5.0 to 5.5 ns, (c) a 5 ns interval from 5 to 10 ns, and (d) a 30 ns interval from 10 to 40 ns [Falck et al. [2008]]. . . . . . . . . . . . . . . . . . . . . . . . . . . .

4. 18

 18 Plots show the calculated G D (r, t) of the OPLS-AA 15 ns (a) and 150 ns (b) simulations for three time-lags. Visible slow time decay of the first shell of neighbours in long time-tail analysis, the log-log plots in the insets, indicate power-law behavior. . . . . . . . . . . . . . . . . . . . . . . . . . . 4.19 Pair correlation function for the centre of mass of SPC/E water molecules for three time slices : t = 0 (solid line), t = 300 fs (dashed line) and t = 1.5 ps (dotted line). The first shell of neighbours, in log plot shown on the inset, decays rapidly and exponentially. . . . . . . . . . . . . . . . . . . . . 4.20 The G D (r, t) of the MARTINI NV T (a) and N Ap z T (b) simulations for three time lags : t = 0 (solid line), t = 500 ps (dashed line) and t = 1.5 ns (dotted line). Similarly to the all-atom case, the long-time tails (the insets) decay as a power-law. . . . . . . . . . . . . . . . . . . . . . . . . . . xxi 4.21 Figures show the parameter dependence on the time-lag, defining the beginning of the G D (r, t) long-time tail for the simulations : (a) 15ns OPLS-AA, (b) 150 ns OPLS-AA, (c) 600 ns NV T MARTINI and (d) 600 ns N Ap z T MARTINI. The maximal time-lag for which the G D (r, t) functions were computed is denoted by t MAX . . . . . . . . . . . . . . . . . . . . . . membranesBiological membranes are one of the most important constituents of prokaryotic and eukaryotic cells. The very beginning of life can be related to the formation of a biological compartments protecting the very first organisms from their environment. This was possible only with the existence of biological membranes -physical barriers between cells and their exterior. A following step in the evolution of cells was a further compartmentalization of cells, leading to the formation of organelles. In this way, biological membranes allowed various microenviornments to coexist within cells, which increased their ability to adapt to di↵erent environments. A biological membrane consists of lipid molecules forming mono-and bi-layers. These layers, in cellular conditions, always form closed structures because of hydrophobic e↵ects resulting from interactions between lipids and water molecules[Berg et al. [2005]]. The reason is that a lipid molecule has a hydrophobic part -a water-repulsive tail, and a hydrophilic one -a water-attractive headgroup. In order to shield hydrophobic regions from water, lipids organize themselves in an aggregated phase. A biological membrane contains also many other diverse types of molecules and actually is a very complex and heterogeneous system[START_REF] Alberts | Essential cell biology[END_REF]]. This is 1 1. General Introduction why biological membranes play also an important role in the bidirectional transport of molecules across the bilayer (or monolayer), in relaying signals, in the adhesion of cells, and in many enzymatic processes.

Figure 1 . 1 :

 11 Figure 1.1: Illustration 1 showing the complexity of a membrane in cells and the structure of a lipid molecule.

Figure 1 . 2 :

 12 Figure 1.2: Schematic 2 of the phosphatidylcholine chemical structure.

Figure 1 . 3 :

 13 Figure 1.3: Sketch of the POPC chemical structure 3 .

Figure 1 . 4 :Figure 1 . 5 :

 1415 Figure 1.4: Crowded cell compartments visualized by cryoelectron microscopy 4 Actin filaments are shown in orange, ribosomes and other macromolecular assemblies in grey, and membrane structures in blue.

Figure 1 . 6 :

 16 Figure 1.6: Visualization of ion tra cking in cells involving Golgi apparatus and rectifier potassium channels (inward and delayed) 6 . The ion channels are synthesized in endoplasmic reticulum and through Golgi are transported to the plasma membrane.

. 3 )

 3 Here x is the position of the particle, D is the di↵usion constant and h• • • i is an average over all realizations of a random walk. Experimental illustrations of Brownian movements were published by Jean Baptiste Perrin [1909], where he recorded the random walk trajectories of putty grains and reconstructed the Gaussian solution of the di↵usion equation by histograms (Fig. 1.7).

Figure 1. 7 :

 7 Figure 1.7: Recorded di↵usion of putty grains published by Perrin [1909].

Figure 1 . 8 :

 18 Figure 1.8: MSD of a model system for various values of ↵ [Kneller [2011]].

  Figure 1.9: Plot representing the displacements of the DPPC lipids in x and y planes for the upper leaflet (a,b,c) and lower leaflet (d); (a) shows a 50 ps interval from 5.00 to 5.05 ns of the simulation, (b) a 500 ps interval from 5.0 to 5.5 ns, (c) a 5 ns interval from 5 to 10 ns, and (d) a 30 ns interval from 10 to 40 ns [Falck et al. [2008]].

Figure 2 . 1 :

 21 Figure 2.1: Results of one of the first molecular dynamics simulation published by Alder [Alder & Wainwright [1957, 1959]] -the figure shows traces of the simulated hard sphere particles for about 3000 collisions.

Figure 2 . 2 :

 22 Figure 2.2: Structure of the first simulated protein -bovine pancreatic trypsin inhibitor (BPTI) 10 . On the bottom of the figure is written full amino acid sequence of this protein.

Figure 2 . 3 :

 23 Figure 2.3: Ilustration of the periodic boundary conditions (PBC) method with the minimum image convention. The box in the middle is replicated in all directions and the non-bonded interactions are calculated only for the closest atoms or their images.

  solutions[START_REF] Jorgensen | The OPLS [optimized potentials for liquid simulations] potential functions for proteins, energy minimizations for crystals of cyclic peptides and crambin[END_REF]; Jorgensen et al. [1996]] and has been extended later for the simulation of macromolecules [Kony et al. [2002]; McDonald & Jorgensen [1998]; Takeuchi [2012]; Watkins & Jorgensen [2001]], including lipids [Pasenkiewicz-Gierula et al. [1999]; Róg et al. [2002]; Shinoda & Okazaki [2001]; Takaoka et al. [2000]].

Figure 3 . 1 :

 31 Figure 3.1: Simulated box of POPC all-atom model membrane. The system consists of total 274 lipids embedded in 10 471 water molecules.

Figure 3 . 2 :Figure 3 . 3 :

 3233 Figure 3.2: Simulated box of POPC coarse-grained model membrane. The system consists of 2033 lipids hydrated by 231 808 water molecules.

Figure 3 . 4 :

 34 Figure 3.4: Total energy of the POPC bilayer during 15 ns simulation with OPLS-AA force-field. The system is in equilibrium.

Figure 3 . 5 :

 35 Figure 3.5: Plots representing the behavior of the system during the 15 ns production run.The plots show fluctuations of : a) the simulation box in all three dimensions; b) the area surface per lipid. Both are stable with the mean value of the latter of 0.638 ± 0.004 nm 2 , which indicates that the system is in equilibrium.

Figure 3 . 6 :

 36 Figure 3.6: Total energy of the POPC bilayer of 150 ns all-atom simulation -the system is in energetic equilibrium.

Figure 3 . 7 :

 37 Figure 3.7: Plots show the fluctuations of the box (a) and the surface area per POPC (b) for 150 ns simulation. The mean area per lipid of 0.637 ± 0.006 nm 2 is close to the value obtained for the 15 ns simulation.

. 8 ,Figure 3 . 8 :

 838 Figure 3.8: Total energy of 600 ns MARTINI simulation in NVT conditions.

Figure 3 . 9 :

 39 Figure 3.9: Total energy fluctuations of 600 ns MARTINI simulation in the N Ap z T ensemble.

Figure 3 . 10 :

 310 Figure 3.10: Fluctuations of the total energy of the reference box of bulk SPC/E water simulated for 100 ps.

  molecule dynamics in a lipid bilayer This chapter focuses on movements of a single lipid molecule. The single particle dynamics is described by autocorrelation functions [McQuarrie [2000]],

  were performed with the nMoldyn program[START_REF] Hinsen | nMoldyn 3: Using task farming for a parallel spectroscopy-oriented analysis of molecular dynamics simulations[END_REF];[START_REF] Kneller | nMOLDYN -a program package for a neutron scattering oriented analysis of molecular dynamics simulations[END_REF]].

4. 1

 1 .1.1 Atom-detailed simulations Expression (1.4) was fitted to the numerically lateral center-of-mass mean square displacements computed for both the shorter (15 ns) simulation and and the longer (150 ns) simulation. The results are presented in Fig. 4.1. In both cases the MSDs exhibit clearly a sublinear behavior. This strongly contrasts with normal di↵usion in case of simulated SPC/E water, for which the W(t) evolves linearly with time (inset in Fig. 4.1a), as one expects for normal di↵usion. The fitted values of the parameter ↵ for 15 ns (↵ = 0.668)

Figure 4 . 1 :

 41 Figure 4.1: The MSD, W(t), calculated for the centres of the masses of the POPC lipids simulated for 15ns (a) and 150 ns (b). The inset in (a) shows the MSD for normal di↵usion (SPC/E water).

  12 m 2 /s) 2.1 ± 0.7 4.0 ± 0.8 7.0 ± 1.0 temperature than the all-atom calculations. According to the temperature dependence of the POPC di↵usion constant shown in Table

Figure 4 .

 4 Figure 4.2: The MSD, W(t), calculated for the centres of the masses of the lipids simulated for 600 ns.

4. 1

 1 .1.3 Comparison between OPLS-AA and MARTINI force-fields In Figure 4.3 the lateral mean squared displacements for the simulations with the 15 ns OPLS-AA and the 600 ns MARTINI force field (NV T and N Ap z T ensembles in the latter case) are presented. The MSDs are compared for a maximal time-lag of 1.5 ns. For all simulations the calculated MSDs have a sublinear form. The curve of the all-atom MSD

Figure 4 . 3 :

 43 Figure 4.3: Comparison of the centre of mass MSDs calculated for the 15 ns all-atom and the 600 ns coarse-grained simulations (the maximal time-lag is 1.5 ns).

  [2009]] on a time scale of 10 ns and normal di↵usion for POPC with the MARTINI ↵ by Niemelä et al.[2010], who analyzed a 500 ns simulation. In this work the subdi↵usive nature of lipids motions is preserved for 15 ns and up to 50 ns

Figure 4 . 4

 44 shows the mean square displacement calculated for the reference SPC/E bulk water, which was simulated for 100 ps in canonical thermodynamic conditions at the temperature of 310 K (section 3.3). The di↵usion constant calculated from the fit of Relation (1.3), gives the value of 3.528 • 10 9 m 2 s (3.528 nm 2 ns ). In comparison to experimental findings, where the di↵usion constant is D = 2.919 • 10 9 m 2 s[START_REF] Mills | Self-di↵usion in normal and heavy water in the range 1-45 C[END_REF]] for temperature of 307.15 K, the obtained value is reasonably close. According to studies of[START_REF] Mark | Structure and Dynamics of the TIP3P, SPC, and SPC/E Water Models at 298 K. The journal of physical chemistry[END_REF], simulations using the SPC/E water model yield best di↵usion constants (compared to other models, e.g. TIP3P and SPC), but still the obtained value is about 20 % too large for the given temperature.A very interesting subject is the behavior of water around lipid membranes. Experimental findings show that water molecules close to the surface of biological membranes interact with lipid molecules, forming a particular hydration layer[START_REF] Arnold | E↵ect of poly(ethylene glycol) on phospholipid hydration and polarity of the external phase[END_REF]; Fitter & Lechner [1999]; Gawrisch et al. [2007]; Pearson & Pascher [1979]]. They create a network of intermolecular bonds between lipids, which are linked to each other via hy-drogen bonds of water molecules, so-called water-bridges, and thus stabilize the structure of membrane [Binder et al. [1999]; Mombelli et al. [2003]; Pasenkiewicz-Gierula et al. [1997]; Pearson & Pascher [1979]; Róg et al. [2009]; Steinbauer et al. [2003]].

Figure 4 . 4 :

 44 Figure 4.4: W(t), calculated for the centres of the masses of the bulk SPC/E water simulated for 100 ps.

Figure 4 . 5 :

 45 Figure 4.5: The MSD, W(t), calculated for the centres of the masses of the water around POPC bilayer in case of 15 ns all-atom simulation (a) and 600 ns MARTINI (b).

Figure 4 . 6 :

 46 Figure 4.6: Comparison of the centre of mass Velocity Autocorrelation Functions calculated for the 15 ns all-atom and 600 ns coarse-grained simulations.

Figure 4 . 7 :

 47 Figure 4.7: Comparison of the integrals over the centre of mass VACFs calculated for the 15 ns all-atom and 600 ns coarse-grained simulations. 57

. 27 )

 27 Here the C vv (n) is defined in range of n 2 ( N + 1, N 1) and additionally C vv (n) = C vv (n + l(2N 1)) for l 2 Z, as the DFT implicitly assumes a periodic signal. Here N is the length of the trajectory. Then the approximation of the Fourier integral (4.26) is given by Cvv (!) ⇡ t Ĉvv (k),(4.28)

Figure 4 . 8 :

 48 Figure 4.8: Calculated densities of states from the VACFs Fourier transforms for the 600 ns MARTINI simulations (CG) and 15 ns OPLS-AA simulation (AA -the inset).

Figure 4 . 9 :

 49 Figure 4.9: Figures show the fits of Expression for the DOS (4.24) for the MARTINI (a) and the OPLS-AA (b) simulations. The insets present values of the function s(!), Eq. (4.25), for the ↵ parameters obtained from the fits (the red line is the value of D ↵ from the MSD calculations). 63

D

  Fig. 4.9). The functions were computed according to Eq. (4.24). The discrepancy to the red line (D ↵ from MSD) is attributed to the fact that ↵ parameters obtained by the DOS and MSD methods, respectively, do not have the same values.

C

  vv (|j k|)a k (j = 1, 2, ..., M) ( 4 .4 5 )

Fig. 4 .

 4 Fig. 4.12 are presented the fits of Relation. 4.24 to the estimated g(!) from the MARTINI (Fig. 4.12a) and OPLS-AA simulations Stronger fluctuations of the g(!) can be seen for both cases, but the characteristic curves of the DOS for subdi↵usion are preserved. The fitted fractional di↵usion constants and the ↵ parameters are very close to the computed values from the Fourier transformed VACFs (Table4.7). This is in particular true in case of the coarse-grained calculations, for which the D ↵ have exactly the same values of 0.059 nm 2 ns ↵ . The ↵ parameters are practically the same and the ↵-parameters for the all-atom simulations obtained from the AR model and by classical Fourier transform are also very close. On the other hand, D ↵ = 0.023 nm 2 ns ↵ overestimates the lateral di↵usion of the POPC lipids even somewhat stronger compared to the analysis with the Fourier

Figure 4 . 10 :Figure 4 . 11 :

 410411 Figure 4.10: The figures show the densities of states for MARTINI (a) and OPLS-AA (b) estimated by AR model. In the insets are presented the results of the FT method for comparison.

.051 @nm 2 êns a D

  model for gHwL Ê Ê ÊÊ ÊÊÊ ÊÊ Ê ÊÊ ÊÊÊÊ ÊÊÊ ÊÊÊÊÊ ÊÊÊÊÊ ÊÊÊÊÊ ÊÊ ÊÊÊ ÊÊÊÊÊÊÊ ÊÊÊ ÊÊÊÊ Ê D a =0.018 @nm 2 êns a D

Figure 4 . 12 :

 412 Figure 4.12: Figures show the fits of Eq. (4.24) to the autoregressive DOS for the MARTINI (a) and OPLS-AA (b) simulations. The insets compare the evolution of s(!) with the value of D ↵ obtained from the MSD fits.

Figure 4 . 13 :

 413 Figure 4.13: The schematic visualising the definition of the pair distribution function 11 .

.55) 11 Source:Figure 4 . 14 :

 11414 Figure 4.14: The pair distribution function, g(r), for the bulk SPC/E water.

Figure 4 . 15 :

 415 Figure 4.15: Three dimensional plots of the distant time-dependant pair correlation functions, G D (r, t), for the all-atom POPC 15 ns (a) and 150 ns (b) simulations.

Figure 4 .

 4 Figure 4.16: 3D plots of the G D (r, t), for the 600 ns NV T (a) and N Ap z T (b) MARTINI simulations.

Figure 4 .

 4 Figure 4.17: 3D plot of the computed G D (r, t) for the simulated bulk SPC/E water.

[

  1996]]. Moreover the first peak at r = 0.2 nm agrees with experimentally measured pair

Figure 4 . 18 :

 418 Figure 4.18: Plots show the calculated G D (r, t) of the OPLS-AA 15 ns (a) and 150 ns (b) simulations for three time-lags. Visible slow time decay of the first shell of neighbours in long time-tail analysis, the log-log plots in the insets, indicate power-law behavior.

Figure 4 .

 4 Figure 4.19: Pair correlation function for the centre of mass of SPC/E water molecules for three time slices : t = 0 (solid line), t = 300 fs (dashed line) and t = 1.5 ps (dotted line). The first shell of neighbours, in log plot shown on the inset, decays rapidly and exponentially.

  4.18 show log-log plots of the G D (r, t) values for the first peak in function of time. The long-time part, chosen from 400 ps for the 15 ns and from 3ns for the 150 ns simulations, shows a power-law decay dependency. A corresponding function was fitted to the data points,G D (t) = At . ÊÊ Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê ÊÊ ÊÊ ÊÊ Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê

Figure 4 .

 4 Figure 4.20: The G D (r, t) of the MARTINI NV T (a) and N Ap z T (b) simulations for three time lags : t = 0 (solid line), t = 500 ps (dashed line) and t = 1.5 ns (dotted line).Similarly to the all-atom case, the long-time tails (the insets) decay as a power-law.

Figure 4 .

 4 Figure 4.21: Figures show the parameter dependence on the time-lag, defining the beginning of the G D (r, t) long-time tail for the simulations : (a) 15ns OPLS-AA, (b) 150 ns OPLS-AA, (c) 600 ns NV T MARTINI and (d) 600 ns N Ap z T MARTINI. The maximal time-lag for which the G D (r, t) functions were computed is denoted by t MAX .

  Avec le modèle MARTINI, cette décroissance est beaucoup plus rapide, ce qui peut indiquer les limites de ce type de champ de force pour reproduire entièrement la dynamique réelle. Membrane simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.1 OPLS-AA force-field . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.1.1 System equilibration . . . . . . . . . . . . . . . . . . . . . 3.2.1.2 Short production run -15 ns . . . . . . . . . . . . . . . . 3.2.1.3 Extended simulation -150 ns . . . . . . . . . . . . . . . . Pair Correlation Function . . . . . . . . . . . . . . . . . . . . . . . 4.2.2 Visualization of the cage e↵ect . . . . . . . . . . . . . . . . . . . . . 4.2.3 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

	Contents 5 Conclusions List of Figures
	6 Publications
	References
	Abstract	i
	Résumé	iii
	Résumé Substantiel	v
	Contents	xiv
	List of Figures	xvii
	1 General Introduction	1
	xiii xvi

1.1 Biological membranes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Transport of molecules in cells and membranes . . . . . . . . . . . . . . . . 7 1.3 Modeling anomalous di↵usion . . . . . . . . . . . . . . . . . . . . . . . . . 1.4 Recent simulation work . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2 Concept of Molecular Dynamics Simulations 17 2.1 Classical Molecular Dynamics Simulations . . . . . . . . . . . . . . . . . . 2.1.1 Concept and potential function . . . . . . . . . . . . . . . . . . . . 2.1.2 Periodic boundary conditions and treatment of electrostatic interactions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Treatment of thermodynamic conditions . . . . . . . . . . . . . . . . . . . xiv 3.1 System setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.1 All-atom OPLS force-field . . . . . . . . . . . . . . . . . . . . . . . 3.1.2 Coarse-grained MARTINI force-field . . . . . . . . . . . . . . . . . 3.2 3.2.2 MARTINI force-field . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.2.1 System equilibration . . . . . . . . . . . . . . . . . . . . . 3.2.2.2 600 ns production run in the canonical ensemble . . . . . 3.2.2.3 600 ns production run in the N Ap z T ensemble . . . . . . 3.3 Bulk water reference simulation . . . . . . . . . . . . . . . . . . . . . . . . 4 Results 4.1 Single molecule dynamics in a lipid bilayer . . . . . . . . . . . . . . . . . . 4.1.1 Mean Square Displacements . . . . . . . . . . . . . . . . . . . . . . 4.1.1.1 Atom-detailed simulations . . . . . . . . . . . . . . . . . . 4.1.1.2 Coarse-grained force-field results . . . . . . . . . . . . . . 4.1.1.3 Comparison between OPLS-AA and MARTINI force-fields 4.1.1.4 Bulk water comparison with membrane hydration water . 4.1.2 Velocity autocorrelation function analysis . . . . . . . . . . . . . . . 4.1.3 Density of States (DOS) . . . . . . . . . . . . . . . . . . . . . . . . 4.1.3.1 DOS as Fourier transform of VACF . . . . . . . . . . . . . 4.1.3.2 Density of the States from the Autoregressive Model . . . 4.1.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Collective motions of lipids . . . . . . . . . . . . . . . . . . . . . . . . . . . xv CONTENTS 4.2.1 1 La complexité d'une membrane biologique 1 . . . . . . . . . . . . . . . . . . vi 2 Comparaison entre les modèles tout-atome OPLS (AA) et gros-grain MAR-TINI (CG) de membrane de POPC. . . . . . . . . . . . . . . . . . . . . . . viii 3 Les fits de l' éq. 1 de déplacement carré moyen pour les simulations: a)

Table 1

 1

	.1: Lipid composition in the heart, liver, erythrocytes and plasma (weight per-
	centage of the total) [Christie [1985]].
	LIPID CLASS		TISSUE
	Heart	Liver	Erythrocytes	Plasma
	Cholesteryl esters Triacylglycerols Cholesterol Diacylglycerols Phosphatidylethanolamine 33.44 ± 0.21 19.91 ± 0.08 20.75 ± 0.21 0.22 ± 0.04 1.47 ± 0.10 -3.77 ± 0.23 6.66 ± 0.15 -4.06 ± 0.16 5.40 ± 0.13 30.23 ± 0.34 5.65 ± 0.10 16.04 ± 0.30 48.65 ± 0.71 0.65 ± 0.07 -0.38 ± 0.09 0.54 ± 0.06 -Phosphatidylinositol 3.69 ± 0.06 4.43 ± 0.08 -3.42 ± 0.14 Phosphatidylserine ---3.11 ± 0.46 Phosphatidylcholine 38.57 ± 0.20 55.18 ± 0.13 32.03 ± 0.14 24.13 ± 0.59 Sphingomyelin 1.76 ± 0.50 2.09 ± 0.06 8.22 ± 0.17 1.99 ± 0.07 Lysophosphatidylcholine --0.89 ± 0.15 1.35 ± 0.06
	Phospholipids have three main moieties (Fig. 1.2):
	1. A hydrophilic headgroup, which consists of the phosphate group and additional
	groups, e.g. a N (CH) 3 group. The latter is called a choline group and is a part of
	phosphatidylcholines. Another example for a headgroup is the ethanolamine rest in
	phosphatidylethanolamine.

Table 1 .

 1

	LIPID CLASS	WHOLE		MEMBRANE
		TISSUE Nuclei Mitochondria Microsomes Plasma
		[A]	[B]	[C]	[C]	[C]
	Phosphatidylethanolamine 25.2 ± 0.5 26.1 Phosphatidylinositol 7.2 ± 0.2 3.9 Phosphatidylserine 3.2 ± 0.1 5.5 Phosphatidylcholine 50.8 ± 0.4 57.3 Sphingomyelin 4.2 ± 0.2 6.3 Lysophosphatidylcholine 1	33.6 ± 1.3 6.6 ± 1.5 0.9 ± 1.3 40.5 ± 2.9 2.4 ± 0.9	21.6 ± 2.8 8.2 ± 1.5 3.9 ± 2.8 58.7 ± 4.0 4.0 ± 0.5	19.8 6.5 3.7 43.1 23.1

2: Lipid abundance in di↵erent organelles in rat liver (weight percentage of the total) [A -Wuthier [1966], B -Gurr et al. [1965], C -Colbeau et al. [1971]].

Table 1 .

 1

	16 : 0	23	26	21	29	21	10
	16 : 1	5	8	6	3	11	16
	18 : 0	6	5	6	18	-	-
	18 : 1	35	31	49	41	23	26
	18 : 2	19	9	16	8	1	2
	18 : 3	2	18	2	-	-	-
	20 : 1	-	-	-	-	10	14
	20 : 5	-	-	-	-	9	7
	22 : 1	-	-	-	-	10	7
	22 : 6

3: Fatty acids composition in di↵erent animal species (weight percentage of the total) [A -Brockerho↵ et al. [1966], B -Christie & Moore [1970], C -Brockerho↵ & Hoyle [1963]]. FATTY ACID SPECIES Rat [A] Horse [A] Duck [A] Pig [B] Herring [C] Seal [C]

Table 4 .

 4 1: Di↵usion fractional constants, D ↵ and ↵ parameters obtained from the fit of the Eq. (1.4) to the CM MSDs of POPC lipids for 15 ns and 150 ns OPLS-AA simulations.[START_REF] Lewis | Thermotropic Phase Behavior of Model Membranes Composed[END_REF];[START_REF] Litman | Packing characteristics of highly unsaturated bilayer lipids: Raman spectroscopic studies of multilamellar phosphatidylcholine dispersions[END_REF]]. In comparison, the experimental value of D ↵ = 0.088 ± 0.007 nm 2 ns ↵ for ↵ = 0.701 has ⌧ vacf = 0.117 ps, which is twice larger than for POPC and three times smaller than for DOPC. In this context it is important to underline that the experimental results in Ref.Schwille et al. [1999a] concern studies of a di↵erent type of lipid -DLPC (dilauroyl-sn-glycero-3-phosphocholine), but computed⌧vacf for the simulations of POPC are of the right order of magnitude. Due to the lack of any experimental works on fractional di↵usion constants for POPC lipids, one can only compare results for the normal di↵usion. NMR studies on the temperature dependence of normal di↵usion constant (↵ = 1) for POPC give values of D = 0.012 ± 0.004 nm 2 ns for the temperature of 308.15 K [Lindblom et al. [1981]] and D = 0.014 ± 0.002 nm 2 ns for 323.15 K [Köchy & Bayerl [1993]], which correspond to ⌧

		OPLS-AA
		15 ns 150 ns
	↵	0.668 0.701
	D ↵ (nm 2 /ns ↵ ) 0.018 0.016

vacf = 0.0017 ps and ⌧ vacf = 0.0021 ps, respectively. These values are much smaller than ⌧ vacf obtained for the simulations, which means that di↵usion of POPC lipids observed in the NMR experiments is much slower.

Table 4 . 3

 43 , for which both MSDs are practically identical. The parameters from the fits of Eq. (1.4), which are presented in Table4.3, give very similar values. The fractional di↵usion constants D ↵ from the fits up to 50 ns are the same for the NVT and N Ap z T simulations: 0.051 nm 2 ns ↵ . The corresponding ↵ parameters of ↵ = 0.558 and ↵ = 0.571, respectively, have also very similar values, and lie in the range of subdi↵usion.

	↵	0.515	0.508	0.571	0.558
	D ↵ (nm 2 /ns ↵ )	0.057	0.058	0.051	0.051
	(Fig. 4.2). The MSDs are very similar, especially for the shorter time-lags of 1.5 ns (inset

: Fractional di↵usion constants, D ↵ and ↵ parameters obtained from the fit of Eq. (1.4) to the CM MSDs of coarse-grained POPC lipids for 600 ns N Ap z T and NVT simulations. MARTINI N Ap z T 1.5 ns NVT 1.5 ns N Ap z T 50 ns NVT 50 ns

Table 4 .

 4 4: Fractional di↵usion constants, D ↵ and ↵ parameters obtained from the fit of Eq. (1.4) to the CM MSDs of coarse-grained POPC lipids for 600 ns N Ap z T and NV T simulations.

		OPLS-AA	N Ap z T CG	NV T CG
		15 ns 150 ns 1.5 ns ⇤ 50 ns ⇤ 1.5 ns ⇤ 50 ns ⇤
	↵	0.668 0.701	0.515	0.571	0.508	0.558
	D ↵ (nm 2 /ns ↵ ) 0.018 0.016	0.057	0.051	0.058	0.051
	of MARTINI lipids (13 beads in comparison to 134 atoms building real POPC molecule)

gHwL Data points -FT Fit of the model for gHwL

	0.20
	0.15
	0.10
	0.05
	0.00	1	0.2	0.3 w @THzD	0.4	0.5	0.6

Table 4 .

 4 6: Comparison between D ↵ and ↵ parameters for OPLS-AA and MARTINI POPC lipids, obtained from the fit of Eq. (1.4) to the MSDs and of the Eq. (4.24) to the DOS for small !.

@THzD gHwL Data points -AR Fit of the model for gHwL

	0.20
	0.15
	0.10
	0.05
	0.00	1	0.2	0.3	0.4	0.5	0.6
			w

Table 4 .

 4 7: Comparison between D ↵ and ↵ parameters for OPLS-AA and MARTINI simulations obtained from the fit of Eq. (1.4) to the MSDs and of Eq. (4.24) to the DOS for small !. The DOS-FT means the g(!) calculated by the Fourier transform and DOS-AR for g(!) computed within the AR process.

			OPLS-AA			MARTINI
		MSD DOS-FT DOS-AR MSD DOS-FT DOS-AR
	↵	0.668	0.437	0.429	0.571	0.467	0.466
	D ↵ (nm 2 /ns ↵ ) 0.018	0.020	0.023	0.051	0.059	0.059

Source : http://www.nature.com/embor/journal/v5/n1/fig_tab/7400056_F3.html

Source : http://openi.nlm.nih.gov/detailedresult.php?img=2741594_pbio.1000203.g012&req=4

Acknowledgements

I would also like to acknowledge financial support by the Agence Nationale de la Recherche.

all-atom calculations [START_REF] Stachura | Anomalous lateral di↵usion in lipid bilayers observed by molecular dynamics simulations with atomistic and coarse-grained force fields[END_REF]]. The comparison between all-atom and coarse-grained models is presented in Fig. 3.3.

Membrane simulations

All simulations in this thesis were conducted with the GROMACS package [START_REF] Hess | GROMACS 4: Algorithms for Highly E cient, Load-Balanced, and Scalable Molecular Simulation[END_REF]] and the resulting trajectories were analysed with the nMoldyn program [START_REF] Hinsen | nMoldyn 3: Using task farming for a parallel spectroscopy-oriented analysis of molecular dynamics simulations[END_REF]; [START_REF] Kneller | nMOLDYN -a program package for a neutron scattering oriented analysis of molecular dynamics simulations[END_REF]]. All simulated systems (MARTINI and OPLS-AA ↵s)

were energetically equilibrated, before the production runs. Equilibration does not only imply constant total energy but also stable box fluctuations and surface areas per POPC lipid, which should be close to experimental data on average. All the simulations were performed above the main phase transition temperature of POPC lipids (270 K), which means that bilayers were in liquid phase [START_REF] Lewis | Thermotropic Phase Behavior of Model Membranes Composed[END_REF]; [START_REF] Litman | Packing characteristics of highly unsaturated bilayer lipids: Raman spectroscopic studies of multilamellar phosphatidylcholine dispersions[END_REF]].

OPLS-AA force-field

The integration time-step in the all-atom simulations was generally set to 1 fs. An exception was the very first stage of the system equilibration, where the time-step was 2 fs.

All bonds involving hydrogen atoms were constrained using the LINCS algorithm [START_REF] Hess | LINCS: a linear constraint solver for molecular simulations[END_REF]]. Periodic boundary conditions were applied in all three directions with the minimum image convention. The Coulomb interactions were treated using the Particle-Mesh Ewald method with a real cut-o↵ radius, r cutof f , of 1.2 nm. The same r cutof f was chosen for the van der Waals interactions. During simulations a constant temperature of As r tends to infinity, the spatial correlations between the position of a tagged central particle ↵ and the partner particles decay and therefore their density in the shell between r and r + r approaches the mean density, ⇢, such that lim r!1 G D (r, t) = 1.

(4.56)

The static pair correlation function is obtained through

Visualization of the cage e↵ect

The lateral time-dependant pair correlation functions, G D (r, t) were calculated according to Eq. (4.55) choosing a maximal time lag of 10% of the respective trajectory length.

For comparison, a similar analysis, but in three dimensions, was performed for the bulk SPC/E water system simulated for 1 ns in NVT ensemble. The biggest radius value was set to be smaller than half the box length in either x or y direction for the all-atom POPC lipids, i.e. 4.4 nm, with the radial step, r, of 0.1 nm. The same parameters were used for the MARTINI simulations. For the water the maximal radius was 2.0 nm, which was less than half the edge of cubic simulation box, with r = 0.01 nm.

The distinct time-dependant pair correlation functions, G D (r, t), for the all-atom and the