Structure and dynamics of liposomes designed for drug delivery: coarse-grained molecular dynamics simulations to reveal the role of lipopolymer incorporation

In this work, coarse-grained molecular dynamics simulations are carried out in NPTH and NVTE statistical ensembles in order to study the structure and dynamics properties of liposomes coated with polyethylene glycol (PEG). The considered liposomes are made by membrane bilayer DPPC with DPPC-PEG incorporated lipopolymers, in an aqueous environment. We have described the two essential PEG conformation regimes, mushroom and brush, and their properties which depend on the grafting density. The effects of grafting density on the structure and dynamics of the membrane were also studied. Our simulations were then discussed by comparing with the available experimental results and by referring to the suitable theoretical models. The results from the NPTH simulations agree with the experimental data of X-ray diffraction and with scale and mean-field theories in terms of thickness of the PEG layer and thickness of the DPPC bilayer membrane. The results from NVTE simulations are found in good agreement with the experimental results from studying the diffusion of the DPPC bilayer membrane and the PEG. The analysis of the mean square displacement revealed that the dynamics of the membranes in the plane show a subdiffusion due to the cage effect and that the grafted PEG dynamics is better described by the Rouse diffusion-mode. Thus, from a macroscopic viewpoint, the incorporation of DPPC-PEG plays an important role in the protection and lubrication of the liposome.


Introduction
Supported lipid bilayer membranes have a growing interest in studying the structure and function of membrane proteins and receptors. [1][2][3] In order to increase the spacing between the hydrophilic support and the bilayers, a random polymer is introduced into this space in a controlled manner. This is supposed to increase the dynamic exibility of the supported bilayer and allow the incorporation of large membrane proteins. This new concept has been used to incorporate integral membrane protein functions. 4,5 Polymer-graed lipid bilayers have several roles in drug delivery, 5,6 such as the stabilization of liposomes, 6 synthesis of supported bilayers for biomaterial applications, 7 modication of the surface of implanted medical devices and protection of cytotoxic drugs. 8 It has also been found that gra lipopolymers stabilize liposomes sterically and provide protection against immune system attacks and antibiotic resistance. 9,10 Polymers play an important role in the synthesis of supported membranes where the bilayer rests on a polymer cushion tethered to a solid substrate. These polymer-supported membranes have shown remarkable stability in air and have the potential to be used for in vitro biosensor applications. 11,12 A polymer commonly used in graing studies is polyethylene glycol (PEG), which is hydrophilic and biocompatible. Cell adhesion is a phenomenon that is essential for most physiological cellular phenomena such as, survival, differentiation, migration, and activation. [13][14][15] It is also essential for pathological situations such as the formation of cancerous metastases, tissue invasion by a pathogenic agent, inammation, and reaction of the organism with respect to biomaterials. [16][17][18] Liposomes are made out of the same material as a cell membrane, so they are supposed to the adhesion phenomena. The incorporation of lipopolymers can be an efficient solution to prevent the adhesion between liposome-liposome and between liposome-human cells. The PEG coating acts as a shield against the hydrophobic properties of the target organ cells that tend to stick the liposomes on their membranes. Everything happens as if it added some kind of lubricant. The quantitative study of the physical properties of the polymer serves to nd the most effective parameters for acting on the target cells. Thus, PEG covered liposomes would be an alternative to the use of nanoparticles to deliver treatments of different pathologies, especially anticancer treatments.
For an experimental viewpoint, several microscopy techniques have been used to examine supported lipid bilayers, such as quartz crystal microbalance (QCM-D), 19 surface plasmon resonance, 20 neutron reectivity, 21 microscopy to atomic force, 22,23 spectroscopy and ellipsometry, 24 and the uorescence interference contrast microscopy (FLIC). 25 Quartz crystal microbalance (QCMD) is a microscopy technique that tracks adsorption, fusion, and rupture of vesicles in different types of surfaces in real time. The neutron reectivity with the labeled deuterium components is a relatively well-established technique for measuring different distances of the diffusion density layer in the supported membranes. However, this technique requires relatively large amounts of hardware, expensive equipment, and data analysis depending on the model. Neutron reectivity measurements also have a poor lateral resolution. Surface plasmon resonance and quartz microbalance techniques simply detect the refractive index and surface mass changes and are therefore insufficient to solve the structural details of the supported bilayers. The uorescence interference contrast microscopy (FLIC) was developed to measure cell/ substrate distances but don't give information about the nanostructure, dynamics and thermodynamic information. Volker Kiessling and Lukas K. Tamm measured the distance of polyethylene glycol between the supported lipid bilayers and the surface of the oxidized silicon chips by uorescence interference contrast microscopy (FLIC). They found a distance ranging from 1.7 to 3.9 nm, depending on the molecular weight of graed polymers, 25 Lambacher and Fromherz measured the distance of poly (2-methyl-2-oxazoline) between the substrate and bilayers supported by uorescence contrast microscopy (FLIC). They found that, for a polymer fraction of 0.5, the increase in the degree of polymerization from n p ¼ 33 to n p ¼ 104 resulted in an increase in the distance from d ¼ 2.3 AE 0.7 nm to d ¼ 4.8 AE 0.6 nm. 26 Note that there are several prerequisites, especially for making and controlling adhesive surfaces dened on a subnanometric scale, and for understanding the Brownian motion of a particle in the vicinity of a wall. To this end, we propose to associate the molecular dynamics simulation (MD) method with the traditional tools of cell biology, to measure and exploit the phenomenon of liposomes adhesion, this method is promising in biophysics to explore the details of the most important physical properties with a good agreement with the available experimental techniques. MD also opens a wide range of studies by modifying experimental parameters and conditions, particularly for complex biological systems. For a dynamic point of view, the subdiffusion behavior, predicted by several theoretical approaches and by experimental methods designed to probe motion at picosecond timescales, can be evaluated using allatom MD simulation and as well as by coarse-grained MD simulations. For the time being there is no simple and satisfactory theoretical model to describe the subdiffusion phenomena caused by the cage effect, as the later can depend on several physicochemical parameters. For this consideration, MD method, which is based on empirical interaction potentials, is a preferred tool for evaluating the anomalous and slow dynamical behavior, as it separately determines the motion of each particle in the simulated system. The document is organized as follows: in Section 2, we present the details of the performed simulations. In Section 3, the distance between two membranes with PEG-graed polymers is estimated by using molecular dynamics simulations in the NPHT statistical ensemble and is compared with that obtained by experiment, then a study of the dynamic and structural properties of the DPPC-PEG lipopolymers and the membrane lipids is carried out in the NVTE statistical ensemble.
2 System model and simulation method

Principle of molecular dynamics
Coarse-grained Molecular Dynamics (CG-MD) simulation is a dynamical simulation technique, in which the equations of motion for a system of particles are numerically integrated over time to obtain particle trajectories in phase space. In contrast to all-atom MD simulation, which usually describes microscopic particles (atoms, or at most small united atom clusters), CG-MD aims can describe particle systems at the coarse-grained mesoscopic level. Accordingly, the point particles in a model describe center-of-mass positions of atoms, clusters, molecules or subsections of polymeric molecules. The principle of CG-MD simulation is particularly simple and consists of generating the trajectories of a nite ensemble of CG beads, by integrating numerically the classical equations of motion F( . The potential energy function U(x i ) of the bead (i) is a function of the particle coordinates x i , it is the applied interaction by the other beads in the chosen force eld, and M i is the mass of the CG bead (i). Each particle trajectory is resolved separately. Thus, the determined trajectories are used to evaluate the static and dynamic properties, passing by some statistical physics equations. It should be noted that, the temporal averages, at very large-times, coincide with the statistical-averages, for ergodic systems. Indeed, the evolution of the positions and velocities of the N beads is obtained by solving a system of motion equations. Such an equation can be integrated numerically using the Verlet algorithm, to give the trajectories of the beads in phase space in each time-step.
At each time-step, as well as for the particle positions, any physical quantity takes instantaneous value. When the system is at thermodynamic equilibrium, the static properties becomes stationary, and its average value can be determined within a time interval. Thus, MD makes the study of the behavior of the considered system possible, on average, by calculating its temporal evolution numerically and averaging the quantity of interest over a long enough time. In the framework of MD simulation, it is necessary to compute the forces on all the beads and to update all the positions. The principle is to discretize, with respect to time, the Newton equations of motion.
For the MD simulation requirement, the development of a force eld describing the interaction potentials between the simulated system particles is of undeniable interest.

Coarse-grained model
We simulate two adjacent DPPC membrane bilayer with PEG graed polymers in an aqueous environment using the MARTINI force-eld. In fact, several coarse-grained models have been proposed to study lipids and polymers. The rst GG model used to represent lipid membranes and polymers has been proposed, in 2001, by Klein et al. 27,28 This model was used to study the self-assembly of lipids/peptides, the insertion of hydrophobic nanopores in a lipid bilayer 29 or the antimicrobial action of a polymer chain. 30 This set of works combines real statistical mechanics with empirical modeling and included considerable insight. These approaches might not be considered fully rigorous, but they have nevertheless a considerable impact on the results they have produced. Then, Marrink et al. proposed another CG model of lipids, which makes it possible to reproduce several lipid phases. [31][32][33] This model is an essentially empirical approach to CG modeling. It uses simple functional forms to describe interaction potentials. The parameters of the interaction potentials are determined by adapting to the experimental data, such as the sharing of free energies between oil and water. To simulate membranes and membrane-related phenomena, the MARTINI CG force eld appears to be the best choice because the membranes are under the action of amphipathic assembly forces. The MARTINI force eld is, therefore, the most suitable eld for modeling such systems. Within the framework of this model, the molecules are not represented by individual atoms, but rather by "pseudo-atoms" approaching groups of atoms, such as whole amino acid residues. By decreasing the freedom degrees, much longer simulation times can be studied at the expense of the molecular details. Fig. 1(b) shows the coarse-grain representation of a DPPC-lipid. In this representation, the PC-head-group consists of two hydrophilic groups: choline (blue) and a phosphate-group (pink). The rst carries a positive charge (+e), while the second, carries a negative one (Àe). The Na particles (yellow) refer to the intermediate polarity glycerol moiety, and the tail double-bonds can be efficiently modeled using slightly less hydrophobic beads, as well as a change in the angleinteraction-potential that governs the rigidity and orientation of the lipid-tails. The MARTINI force eld is based on a mapping of one to four, which means that, on average, four heavy atoms, including the associated hydrogens, can be represented by a single coarse-grained bead. Consequently, one coarse-grained water bead corresponds to four water molecules. To properly reproduce the chemical nature of the modeled systems, four main types of the coarse-grained particles are dened, namely, metallic, polar, nonpolar and charged. The four main types of the coarse-grained particles are divided into sub-types, based on hydrogen-bonding capabilities (donor, acceptor, both or none) and polarity (ranging from 1-low polarity to 5-high polarity), giving a total of 18 unique "building blocks". The described mapping scheme provides a relatively straightforward and effective way of switching from all-atom to coarse-grained representation for a wide range of the biological systems. Interactions between the coarse-grained particles are described by a force eld containing terms typical for other classical force elds. In the MARTINI coarse-grained model, each monomer in polyethyleneglycol PEG is represented by a single bead (EG). Fig. 1(c) show schematic representation of a lipopolymer DPPC-PEG, where the DPPC-PEG bond is assumed to be rigid.

Interaction potentials in the MARTINI force eld
The proposed model for membranes supported via gra polymers consists of lipid bilayers supported by PEG polymer chains, the bilayers are consisting of dipalmitoylphosphatidylcholine lipids (DPPC), a schematic representation is presented in Fig. 1. The simulation used the Martini coarse-grained force eld developed by J. Marrink for use with DPPC lipids, water (MW), and PEG polymer. 34,35 The importance of using the coarse-grained model is to allow the extend of the spatial and temporal scales of the simulations compared to the all-atom model. The MARTINI force eld is one of the widely used Coarse-Grained (CG) models in MD simulations. More details about the MARTINI CG force eld for lipids and polymers can be found in the literature. [34][35][36][37][38][39][40] Here, we will present the model briey. In the MARTINI CG model, all particle pairs (in the Martini force eld) i and j at distance r interact via a shied Lennard-Jones V LJ and coulombic V C interactions: For the DPPC lipids and PEG polymer the bond interactions between adjacent pairs of beads are described by a harmonic potential V b (r): With the interaction parameters between PEG monomers K b(EG-EG) ¼ 20.3 kcal mol À1ÅÀ2 , and r 0(EG-EG) ¼ 3.3Å. The bonded interaction parameters for bond DPPC bead are set to K b(ij) ¼ 24 kcal mol À1ÅÀ2 and r 0(ij) ¼ 4.7Å. The angle interaction between triplets of DPPC beads use a cosine squared harmonic potential V a1 (q) With K a1 ¼ 3 kcal mol À1 for all angles, q 0 ¼ 120 for the Q a -Na-Na equilibrium angle and q 0 ¼ 180 for the others angles. The PEG polymer, the a harmonic angle interaction V a2 is used, and Fourier dihedral interaction potential V d is added to represent the PEG chains torsion.
The parameters of the dihedral interaction are given in Table 2.

Simulation details
The principle of MD simulation consists of generating the trajectories of a nite set of particles, by integrating numerically the classical equations of motion. The trajectories thus determined are used to evaluate the static and dynamic properties by temporal averages, which, at very large-times, coincide with the statistical-averages, for ergodic systems. The motion equation of each bead can be integrated numerically using MD simulation, to give the trajectories of each bead in phase space. At each point of the trajectory, any physical quantity, A, takes instantaneous value, A(t). When the system is at thermodynamic equilibrium, A becomes stationary, and its average value can be determined within a time interval. Thus, DM makes possible the study of the behavior of the considered system, on average, by calculating its temporal evolution numerically and averaging the quantity of interest over a long enough time. Then, for larger times, physics respects the ergodicity principle stipulating that: lim where A designs any physical property andĀ t represents its temporal mean-expectation-value, and hAi accounts for an equilibrium-average calculated in the adopted statistical ensemble. To study the physical properties of lipid membranes with graed polymers using the molecular dynamics simulation method based on the MARTINI coarsegrained model, we have simulated a lipid membranes with graed polymers model consisting of 15 376 DPPC molecules and different PEG lipopolymer molar fraction (0.005, 0.014, and 0.1), fully hydrated by 112 659 water beads. Molecular dynamics simulations were performed using the LAMMPS soware package, 41 In order to eliminate the contributions of surfaces that affect the physical properties of the system, periodic boundary conditions are imposed. 42 We performed molecular dynamics simulations in the NPHT and NVTE ensembles. N is the number of coarse-grained beads, V is the volume of the simulated system, T is the temperature, H is the enthalpy, E is the internal energy and P is the external pressure. In the rst simulation we plan to calculate the equilibrium distance between the two membranes for different molar fractions of lipopolymer at ambient conditions of temperature and pressure. The preformed initial distribution (Fig. 1) is adjusted in local potential energy minimum by iteratively adjusting atom coordinates. Then, the system is equilibrated under the NPHT conditions, by coupling the NPH barostat 43 to the Langevin thermostat 44 for one million time-step, with a time-step of 10 fs, we remark that, aer a 2 Â 10 5 time-step, the volume remains practically constant and then the polymer layer in the z-direction can be calculated. The NPHT simulation updates the position and speed for each time step and changes the volume freely to achieve the imposed equilibrium conditions of external pressure and temperature (T ¼ 300 K, P ¼ 1 bar). For these conditions, the membrane bilayer is in a gel phase giving an area per lipid A l around 0.7 nm 2 . The second simulation is a Brownian dynamics (BD) performed in the NVTE ensemble, it is done by coupling the Langevin thermostat which models an interaction with an implicit background solvent with T ¼ 300 K, to the NVE integration scheme which performs a constant NVE simulation to update the position and the speed of atoms, 45 in this simulation we investigate the structure and dynamics of both DPPC lipids and PEG polymers. The NVTE simulations were run for two million time-step, with a time-step of 10 fs, giving a simulation time of 20 ns.

Polymer layer and membrane thickness
Polymer physics theories, starting with Flory 46 for the free chains and passing to de Gennes 47 and Alexander, 48 for graed polymers, can be used to characterize statistical congurations and the thermodynamics of polymer chains. The study of polymers treated on solid surfaces shows the existence of two regions of concentration of polymers. These are characterized by the so-called "mushroom" and "brush" congurations of the gra polymer chain. The "mushroom" and "brush" regimes are qualied from the visualization of the simulated system and by comparing the polymer layers of the MD simulation with those of the theory.
The mushroom and brush regimes are qualied by qualitatively from the visualization of the simulated system and quantitatively by comparing the polymer layer to that established theoretically and experimentally.
The mushroom regime is due to low concentrations of gra polymer and the brush regime refers to higher concentrations. The reason for using the same parameters to present the PEG used in the experimental work of Evans et al., 49 in term of the degree of polymerization n p ¼ 45 and the monomer radius a ¼ 4.3Å is to validate the MD simulations by comparing its results to the experimental ndings. We note that, the investigation of other degrees of polymerization and molar fractions for PEG or for other polymers, polycations or polyanions is can be conducted using the present MD simulations protocol. For a comparison to the theoretical predictions, the conguration of the chain in the two regimes can be treated by the methods of polymer physics, as developed below: 3.1.1 Mushroom regime. When the chains are graed onto a surface with little density, they do not interact and behave almost like isolated graed chains as shown in Fig. 2a and b. It is customary to call this regime the mushroom regime. As a rst approximation, they are not affected by the presence of the surface and their average height is the order of R F which is determined by the degree of polymerization n p and the size a of the monomer unit. The dimensions of the free polymer are determined by steric interactions (the excluded volume effects) that tend to extend the polymer. These are balanced by the entropic effects unfavorable to the stretching of the coiled polymer chain. In Flory's theory, the free energy of a real chain is given by: where R is the end-to-end distance of the chain and n is the excluded volume by monomer. In general n ¼ a 3 (1 À 2c). Where c is the Flory-Huggins interaction parameter with an athermic solvent, for which non-steric intramolecular interactions may be neglected n z a 3 . The rst term of eqn (8) represents the excluded volume interactions, which are represented in the mean eld approximation and are proportional to the square of the local concentration. The second term of the eqn (8) represents the free energy of stretching, which is derived from the entropy of an ideal polymer chain whose extremities undergo a random (Gaussian) walk. 48 Free energy in eqn (8) with respect to R, gives the following expression for the Flory radius (threedimensional): Hristova and Needham 50 nd that for the PEG monomer unit of oxyethylene: a z 0.35 nm. From the adhesion measurements, Evans et al. 49 found that the data for PEG lipids is better adjusted by a value of a z 0.43 nm. In this work, for a reason to compare the coarse-grained MD simulation results with the experimental ndings, we simulate PEG polymers with a degree of polymerization n p ¼ 45 and different molar fractions, similar to the works of Evans et al. 49 3.1.2 Transition from the mushroom regime to the brush regime. The transition between mushroom and brush regimes occurs at the gra lipid concentration for which the polymer chains associated with the surface begin to overlap. This condition is satised approximately to polymer molar fractions X c/b p such as: where A l is the area of the membrane by lipid (A l ¼ 0.6-0.7 nm 2 ) for lipid in the uid phase and (A l ¼ 0.432 À 0.48 nm 2 ) in the gel phase. For a gra polymer with molecular weights of 2000, i.e. (n p ¼ 45) for the PEG polymer, the transition between the mushroom and brush congurations must take place at X c/b p ¼ 0:014 for uid phase membranes.

Brush regime.
A polymer brush consists of polymer chains, one and only one end of which is covalently bonded to a solid substrate, as schematically shown in Fig. 4a and b. As the concentration of the gra polymer increases, the polymer head groups begin to interact and present a more stretched conguration (like a brush) in which the polymer chains extend from the surface of the membrane (Fig. 3 and 4).
3.1.3.1 Mean eld theory for the brush regime. Two different treatments have been adopted to study the polymer chain bond in the brush regime. The rst is based on the average eld theory which, in essence, is similar to that presented in eqn (8). The second is the theory of scale laws that has been proposed by de Gennes. The average eld theory is presented in this section and the theory of scaling laws in the next section. In the brush  This journal is © The Royal Society of Chemistry 2020 RSC Adv., 2020, 10, 3745-3755 | 3749 regime, the polymer chains are extended and the problem is essentially unidimensional. The chains are conned in the normal direction to the graing surface ( Fig. 2a and b). For a graed polymer chain, the one-dimensional version of the average eld expression for free energy is: 51 where X p is the molar fraction of the lipid polymer and A l is the area per lipid. Minimizing the free energy of eqn (11) with respect to R (end-to-end distance of the chain) gives the following expression for the equilibrium length of the polymer brush: 3.1.3.2 Scale theory for the brush regime. The essential determinant in scale law theory is the one-dimensional nature of the problem for conned polymer chains in the brush regime ( Fig. 2a and b). Under these conditions, the length L of the polymer chains is linear in degree of polymerization n p . For long polymer chains, the scale law for the equilibrium length is given by: 47,52 where D is the diameter of the connement, i.e. the distance between the graing points of the polymer. The scale of length in the problem is the Flory radius R F3 , given by eqn (6). The requirement that L TE be linear in n p gives for the exponent the value m L ¼ 2/3. In terms of the molar fraction of the polymer lipid, the distance D between the graing points is given by : Eqn (10) then becomes: This is exactly the same result as obtained from the mean eld theory given in eqn (9). In Table 1 the thickness of the polymer layers available experimentally, for the PEG polymer with a degree of polymerization n p ¼ 45 and different molar fractions, by A. K. Kenworthy et al. 53 and DPPC bilayer thickness measured by N. Kucerka et al. 54 is compared to that obtained by the present MD simulation at the NPHT statistical ensemble under the ambient condition of pressure and temperature, for different molar fractions of lipopolymer (0.005, 0.014 and 0.1) using a lipopolymer composed of a DPPC molecule bound to a PEG polymer chain of degree of polymerization n p ¼ 45 and monomer radius a ¼ 0.43 nm, for this parameters the transition mushroom-brush lipopolymer molar fraction is X p ¼ 0.014, from Fig. 2, it is clear that for X p ¼ 0.005 the lipopolymer shows a mushroom regime, for X p ¼ 0.014 shows a critical regime and for X p ¼ 0.1 shows a brush regime. From molecular dynamics simulation, as the membranes are more congested for comparison to the diluted polymers, we can deduce the polymer and membrane thickness by analyzing the density prole of the simulated system using the Ovito visualization tool. 55 We nd that the thickness of the polymer increases with the increase of the molar fraction of the gra polymers and that the thickness of the membrane decreases a little from the mushroom regime to the brush regime. Thus, the polymer layer, which covers the liposomes, plays an important role in its protection. The incorporation of lipopolymers prevent the adhesion liposomeliposome even with a small molar fraction of graed polymers, and they form a layer which protects the liposome from supposed system immune attack (Fig. 5 and Table 3).

Radial distribution function and mean-eld interaction potential
The investigation of the three-dimensional distribution of the beads of DPPC lipids and PEG polymer, for the two graing regimes, is carried out by analyzing the radial-distribution-  function (RDF) from MD simulations in the NVTE statistical ensemble. RDF can be dened as the probability of nding a coarse-grained bead distributed around a given (central) bead, in three dimensions. In fact, RDF provides information on the local distribution of the coarse-grained beads in a group among the various beads of the system. Analytically, RDF is dened by the following thermal-average: where r is the number-particle-density. With r ¼ N/V, where N is the number of constituents (DPPC-molecules or PEG polymer beads) and V is the volume of the simulation box. The above sum is performed over all pairs of constituents. Fig. 6a and b shows the radial distribution function (RDF), which gives the probability of nding a coarse-grained bead at a distance r from another bead of the same considered group of molecules (lipids or polymers). Inspection of this radial distribution function conrms that, for brush and mushroom regimes, the membrane bilayer is in the gel phase, the main peaks of the RDF between lipids increase proportionally to the gra polymer molar fraction, i.e. the membrane lipids are more congested for the brush regime, for the reason that the interaction between lipids is more attractive in the brush regime than the mushroom regime. For polymers, in the brush regime, the polymer RDF has narrow peaks corresponding to an ordered monomers distribution, and in the mushroom regime, the polymer RDF has large peaks, which means that it has a messy structure. Fig. 7a and b shows the effective mean-eld interaction potential computed as U eff /k B T ¼ Àlog(g(r)) for polymers and lipids in the brush and mushroom regimes. In the mushroom regime, the interaction between the monomers is an attractive type. However, when the concentration of the graed polymer increases from the mushroom regime to the brush regime, the interaction between the monomers becomes less attractive and the polymer head groups present a stretched conguration. The interaction between lipids is more attractive in the brush regime. Thus, the stability of DPPC lipids bilayer increases with the increase of lipopolymer molar fraction.    Table 2 Dihedral interaction potential parameters The hD 2 r CM i is found proportional to the diffusion coefficient D and the diffusion exponent a for a system of dimension d by: 56 which can be expressed as: where the log function refers to the common logarithm. From the log-log plot of MSD versus time we determine the diffusion coefficients D a and the exponent a < 1 related to the different diffusion regimes. We underline that subdiffusion is a feature of overcrowded systems, where the trajectories of their mobile constituents are strongly correlated. Note that the scale relationship above is valid for signicant times, that is, beyond a characteristic time depending on the specic details of the broadcast process and the host support structure. This longterm behavior deviates from the linear dependence on time found for Brownian motion. 57 Here, D a represents the generalized diffusion coefficient or (the fractional diffusion coefficient). The latter is expressed in length 2 /time a unit. We underline that subdiffusion is a characteristic of saturated systems, where the trajectories of their mobile constituents are strongly correlated. Note that the scale relationship above is valid for signicant times, that is, beyond a characteristic time depending on the specic details of the broadcast process and the host support structure. In general, a particle is said to be subdiffusive if the condition hD 2 r CM (t)i(t)/t/0, for t / +N, is satised (very slow diffusion). This explains why the exponent a must be in the range 0 < a < 1. From an experimental point of view, subdiffusion has been observed in many scientic elds. 57 Subdiffusion originates from the cage effect, 57 where the tracer is trapped in a cage formed by its neighbors. It has been found that the length of stay, Dt c , of a tracer in the cage, composed of N neighbors, is noted: Dt c $ N 1/a , where a is the abnormal exponent.

Polymer dynamic.
In the coarse-grained representation, PEG chains can be understood as a sequence of N beads, EG-monomers, joined by a stretching potential. 58 The dynamics of the PEG is modeled by the Brownian motion of EGmonomers. The polymer dynamics model, in the coarsegrained scale, was rst proposed by Rouse. In the framework of the Rouse model, the equation of motion of the EG-beads is described by the Langevin equations for segmental motion: where, (R 1 , R 2 ,., R N h R n ) is the position of the beads, k ¼ 3k B T/ a 2 is the entropic spring force and f n a random force. We assume that each bead feels the same friction x, that its motion is overdamped, and that the diffusion coefficient D ¼ k B T/x is independent of the position R n of the bead.
The random force f n veries the conditions: 59 The solution of the Langevin equation gives the following behavior of the mean-square-displacement as a function of time.
with s is a characteristic time for wich the motion of the center of mass fallow a normal diffusion law, in a space of dimension d, with a diffusion constant D.
On the other hand, the polymers can show a subdiffusion due to the cage effect. For this reason, the diffusion of the Fig. 7 The mean-field interaction potential (U eff (r)), (a) U eff (r) of the group of lipopolymer molecules, (b) U eff (r) of the group of lipid molecules.
center-of-mass is not normal. In this case, the mean-squaredisplacement of a monomer is given by:

Analysis of MD results.
To study the inuence of the molar fraction of the graed polymer on their dynamics, and then on the membrane dynamics, we keep their degree of polymerization xed to the values: n p ¼ 45, and vary the molar fraction. Fig. 8 and 9 show the results of the molecular dynamics simulations in the NVTE statistic ensemble, for both brush and mushroom regimes, the temperature is xed at 300 K. The diffusion coefficients D a , that are determined by extrapolating the log-log plot of the MSD as a function of time, are given in Table 4. In Fig. 8, we depict the log-log plot of MSD against time for polymers. We rst remark that aer a short initial regime of subdiffusive motion (t ¼ 100 ps), with a diffusion exponent a ¼ 0.5, the polymer dynamics reach the normal diffusion, aer this time MSDs are straight lines of the same slope a ¼ 1. The PEG dynamics, for the mushroom regime, is in perfect agreement with the analytical prediction in the framework of the Rouse model. However, for the brush regime, PEG dynamics, at short time, show a deviation from the Rouse model. In this regime, the polymer density is high, and can cause a cage effect that blocks the monomers movement and therefore, the exponent a\ 1 2 : Experimental values of the diffusion coefficient of the PEG polymer were found by Waggoner et al. 60 For a low molecular weight, the diffusion coefficient is of the order of 10 À6 cm 2 s À1 and for a large molecular weight, it is of the order of 10 À7 cm 2 s À1 . Our coarsegrained simulation for branched PEGs is in good agreement with the experimental ndings. In Fig. 9 we depict the log-log plot of MSD against time for lipids, It is well known that the membrane bilayer at ambient temperature appears to be in a gel phase characterized by a slow lateral diffusion of lipids, the diffusion of lipids is anomalous. This is demonstrated in MSD curves depicting the time-dependent anomalous diffusion exponent (a < 1), a remains <1 even for long time. We remark a slight decrease of the membranes dynamics passing from the mushroom to the brush regime. The values obtained experimentally of the diffusion coefficient in the plane of a gel membrane are of the order of 10 À7 cm 2 s À1 , 61 MD results show that the polymer coating to the membranes has no important inuence on the lipids dynamics properties. The dynamics behavior of lipids can be described, as well as polymers, by the Rouse model. But, the normal diffusion expected by Rouse appears at very long time (t T 10) ns, as observed by Flenner and coworkers in an all-atom MD simulation of the DMPC bilayer membrane. 62 In this work, we are interested in the subdiffusion phenomenon. The simulation time is of 20 ns and it is sufficient to explore the anomalous dynamics, which takes place for a short time t < 10 ns. At long simulation time, this regime disappears and the normal diffusion takes place. We recall that lipids show a long time subdiffusion. In contrast, the polymers have a normal diffusion aer a short regime of subdiffusion, even in the brush regime, where the molar fraction of coating polymers is important. Furthermore, it should be noted that, for a macroscale viewpoint where the liposomes can be    This journal is © The Royal Society of Chemistry 2020 RSC Adv., 2020, 10, 3745-3755 | 3753 viewed as dispersed colloidal microspheres, the polymer plays an important role in preventing the liposomes aggregation and in facilitating their transport (Table 5). 63

Conclusion
We studied the polymer layer, the structural and dynamic properties of a liposome model designed for drug delivery using coarse-grained molecular dynamics simulations at two different statistical ensembles. From the NPTH MD simulations, it has been observed that the distance of the polymer layer increases with the increase of the molar fraction of the gra polymers and shows a transition from the mushroom regime to the brush regime. When polymer chains are graed onto the surface of the membrane with a small molar fraction, they do not interact and behave almost like isolated chains. It is customary to call this regime the mushroom regime. As a rst approximation, the gra polymers are not affected by the presence of the surfaces, for ambient conditions of pressure and temperature, and their average height is of the order of R F which is determined by the degree of polymerization n p and the size of the monomer unit.
As the concentration of the gra polymer increases, the polymer beads of different chains interact and exhibit a more stretched conguration (such as a brush) in which the polymer chains extend from the surface of the membrane. The results obtained using the MD simulation are in good agreement with the experimental ndings from X-ray diffraction. From the NVTE MD simulations, the analysis of the radial distribution function has shown that the modication of the molar fraction of lipid gra polymers changes the interaction of the hydrophobic chain. For the brush regime, the interaction of the hydrophobic chain is more attractive within the bilayer, resulting in a more close-packed and rigid membrane. The analysis of the meansquare-displacement versus time revealed that the dynamics of the membrane slowing during the transition from the mushroom regime to the brush regime. The polymers move freely following a normal diffusion and then, from a macroscopic scale viewpoint, increase the dynamic exibility of the liposomes by preventing the macroscopic aggregation of liposomes. The diffusion coefficients calculated using NVTE MD simulation are of the order of those measured experimentally using pulsed-gradient spin-echo nuclear magnetic resonance method and single-particle tracking experiments for PEG polymers and DPPC lipids, respectively. We recall that the present nanoscale results can be extended to the industrial scale and proposing techniques of encapsulation of fragile pharmaceutical compounds in nanocapsules. Because of their properties, liposomes covered with biocompatible polymers are ideal for protecting and transporting sensitive pharmaceutical compounds such as those used in cancer therapy.

Conflicts of interest
There are no conicts to declare.