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

Abstract

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.

---

1 Introduction

Supported lipid bilayer membranes have a growing interest in studying the structure and function of membrane proteins and receptors.^1–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 flexibility 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-grafted lipid bilayers have several roles in drug delivery,^5,6 such as the stabilization of liposomes,⁶ synthesis of supported bilayers for biomaterial applications,⁷ modification of the surface of implanted medical devices and protection of cytotoxic drugs.⁸ It has also been found that graft 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 grafting 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–15 It is also essential for pathological situations such as the formation of cancerous metastases, tissue invasion by a pathogenic agent, inflammation, and reaction of the organism with respect to biomaterials.^16–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 find 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),¹⁹ surface plasmon resonance,²⁰ neutron reflectivity,²¹ microscopy to atomic force,^22,23 spectroscopy and ellipsometry,²⁴ and the fluorescence interference contrast microscopy (FLIC).²⁵ 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 reflectivity 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 reflectivity 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 fluorescence 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 fluorescence interference contrast microscopy (FLIC). They found a distance ranging from 1.7 to 3.9 nm, depending on the molecular weight of grafted polymers,²⁵ Lambacher and Fromherz measured the distance of poly (2-methyl-2-oxazoline) between the substrate and bilayers supported by fluorescence 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 ± 0.7 nm to d = 4.8 ± 0.6 nm.²⁶ Note that there are several prerequisites, especially for making and controlling adhesive surfaces defined 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 all-atom 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-grafted 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

2.1 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 finite ensemble of CG beads, by integrating numerically the classical equations of motion F(x_i) = −∇U(x_i) = M_iV̇_i(t). 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 field, 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 field describing the interaction potentials between the simulated system particles is of undeniable interest.

2.2 Coarse-grained model

We simulate two adjacent DPPC membrane bilayer with PEG grafted polymers in an aqueous environment using the MARTINI force-field. In fact, several coarse-grained models have been proposed to study lipids and polymers. The first 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²⁹ or the antimicrobial action of a polymer chain.³⁰ 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–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 field appears to be the best choice because the membranes are under the action of amphipathic assembly forces. The MARTINI force field is, therefore, the most suitable field 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 first 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 angle-interaction-potential that governs the rigidity and orientation of the lipid-tails. The MARTINI force field 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 defined, 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 field containing terms typical for other classical force fields. 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.

Fig. 1 Schematic representation of the studied system, (a) initial distribution of the simulated model used in the molecular dynamics simulation, the lengths of the simulation box are lx = 480 Å, ly = 480 Å and lz = 465 Å. (b) DPPC lipid model. (c) Lipopolymer (DPPC-PEG) model.

2.3 Interaction potentials in the MARTINI force field

The proposed model for membranes supported via graft 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 field 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 field is one of the widely used Coarse-Grained (CG) models in MD simulations. More details about the MARTINI CG force field for lipids and polymers can be found in the literature.^34–40 Here, we will present the model briefly. In the MARTINI CG model, all particle pairs (in the Martini force field) i and j at distance r interact via a shifted 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):

V_b = K_b(r − r₀)² (4)

With the interaction parameters between PEG monomers K_b(EG–EG) = 20.3 kcal mol⁻¹ Å⁻², and r_0(EG–EG) = 3.3 Å. The bonded interaction parameters for bond DPPC bead are set to K_b(ij) = 24 kcal mol⁻¹ Å⁻² and r_0(ij) = 4.7 Å. The angle interaction between triplets of DPPC beads use a cosine squared harmonic potential V_a1(θ)

V_a1 = K_a1(cos(θ)−cos(θ₀))² (5)

With K_a1 = 3 kcal mol⁻¹ for all angles, θ₀ = 120° for the Q_a–Na–Na equilibrium angle and θ₀ = 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.

With K_a2 = 20.3 kcal mol⁻¹, and θ₀ = 130°.

The parameters of the dihedral interaction are given in Table 2.

2.4 Simulation details

The principle of MD simulation consists of generating the trajectories of a finite 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: where A designs any physical property and Ā_t represents its temporal mean-expectation-value, and 〈A〉 accounts for an equilibrium-average calculated in the adopted statistical ensemble. To study the physical properties of lipid membranes with grafted polymers using the molecular dynamics simulation method based on the MARTINI coarse-grained model, we have simulated a lipid membranes with grafted 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 software package,⁴¹ In order to eliminate the contributions of surfaces that affect the physical properties of the system, periodic boundary conditions are imposed.⁴² 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 first 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⁴³ to the Langevin thermostat⁴⁴ for one million time-step, with a time-step of 10 fs, we remark that, after a 2 × 10⁵ 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². 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,⁴⁵ 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.

3 Results and discussions

3.1 Polymer layer and membrane thickness

Polymer physics theories, starting with Flory⁴⁶ for the free chains and passing to de Gennes⁴⁷ and Alexander,⁴⁸ for grafted polymers, can be used to characterize statistical configurations 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” configurations of the graft polymer chain. The “mushroom” and “brush” regimes are qualified 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 qualified 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 graft 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.,⁴⁹ 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 findings. 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 configuration 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 grafted onto a surface with little density, they do not interact and behave almost like isolated grafted chains as shown in Fig. 2a and b. It is customary to call this regime the mushroom regime. As a first 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 ν is the excluded volume by monomer. In general ν = a³(1 − 2χ). Where χ is the Flory–Huggins interaction parameter with an athermic solvent, for which non-steric intramolecular interactions may be neglected ν ≈ a³. The first term of eqn (8) represents the excluded volume interactions, which are represented in the mean field 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.⁴⁸ Free energy in eqn (8) with respect to R, gives the following expression for the Flory radius (three-dimensional):

Hristova and Needham⁵⁰ find that for the PEG monomer unit of oxyethylene: a ≈ 0.35 nm. From the adhesion measurements, Evans et al.⁴⁹ found that the data for PEG lipids is better adjusted by a value of a ≈ 0.43 nm. In this work, for a reason to compare the coarse-grained MD simulation results with the experimental findings, we simulate PEG polymers with a degree of polymerization n_p = 45 and different molar fractions, similar to the works of Evans et al.⁴⁹

3.1.2 Transition from the mushroom regime to the brush regime. The transition between mushroom and brush regimes occurs at the graft lipid concentration for which the polymer chains associated with the surface begin to overlap. This condition is satisfied approximately to polymer molar fractions such as:where A_l is the area of the membrane by lipid (A_l = 0.6–0.7 nm²) for lipid in the fluid phase and (A_l = 0.432 − 0.48 nm²) in the gel phase. For a graft polymer with molecular weights of 2000, i.e. (n_p = 45) for the PEG polymer, the transition between the mushroom and brush configurations must take place at for fluid phase membranes.

3.1.3 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 graft polymer increases, the polymer head groups begin to interact and present a more stretched configuration (like a brush) in which the polymer chains extend from the surface of the membrane (Fig. 3 and 4).

Fig. 2 Captions of the NPHT simulation at T = 300 K and P = 1 atm of the mushroom configuration regime. The model membranes composed of DPPC lipid bilayer are decorated with a molar fractions X_p = 0.005 of lipopolymer DPPC-PEG with the same degree of polymerization n_p = 45 and the whole is hydrated by water molecules. After one nanosecond, the simulated system is equilibrated and the membrane–membrane distance remains generally constant (for clarity, the water molecules are not showing). (a) MD Simulation in the mushroom for X_p = 0.005. (b) Lipopolymer conformation in the mushroom regime.

Fig. 3 Captions of the NPHT simulation at T = 300 K and P = 1 atm of the critic configuration regime. The model membranes composed of DPPC lipid bilayer are decorated with a molar fractions X_p = 0.014 of lipopolymer DPPC-PEG with the same degree of polymerization n_p = 45. (a) MD Simulation in the critic regime for X_p = 0.014. (b) Lipopolymer conformation in the critic regime.

Fig. 4 Captions of the NPHT simulation at T = 300 K and P = 1 atm of brush mushroom configuration regime. The model membranes composed of DPPC lipid bilayer are decorated with a molar fractions X_p = 0.1 of lipopolymer DPPC-PEG with the same degree of polymerization n_p = 45 and the whole is hydrated by water molecules. After one nanosecond, the simulated system is equilibrated and the membrane–membrane distance remains generally constant (for clarity, the water molecules are not showing). (a) MD Simulation in the brush regime for X_p = 0.1. (b) Lipopolymer conformation in the brush regime.

3.1.3.1 Mean field theory for the brush regime. Two different treatments have been adopted to study the polymer chain bond in the brush regime. The first is based on the average field 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 field theory is presented in this section and the theory of scaling laws in the next section. In the brush regime, the polymer chains are extended and the problem is essentially unidimensional. The chains are confined in the normal direction to the grafting surface (Fig. 2a and b). For a grafted polymer chain, the one-dimensional version of the average field expression for free energy is:⁵¹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 confined 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,52where D is the diameter of the confinement, i.e. the distance between the grafting 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 grafting points is given by Eqn (10) then becomes:

This is exactly the same result as obtained from the mean field 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.⁵³ and DPPC bilayer thickness measured by N. Kucerka et al.⁵⁴ 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 profile of the simulated system using the Ovito visualization tool.⁵⁵ We find that the thickness of the polymer increases with the increase of the molar fraction of the graft 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 liposome–liposome even with a small molar fraction of grafted polymers, and they form a layer which protects the liposome from supposed system immune attack (Fig. 5 and Table 3).

Table 1 Non-bonded interaction potential parameters

Pair type	σ (Å)	ε (kcal mol^−1)
Q0–Q0, C1–C1	4.7	0.836521
Qa–Qa, Qa–EG, MW–MW	4.7	1.195030
Na–Na, Q0–Na, Qa–Na, Q0–EG	4.7	0.956024
EG–EG	4.3	0.806644
Q0–Qa	4.7	1.075527
Q0–C1, Qa–C1	6.2	0.478012
MW–C1	4.7	0.478012
Na–C1, EG–C1	4.7	0.645316
Q0–MW, Qa–MW	4.7	1.338434
Na–MW	4.7	0.956024
EG–Na	4.7	1.075527
EG–MW	4.7	1.075527

---

Fig. 5 The particle-count in the z-direction of the simulated system for the two essential regimes of polymers configuration, where h is the membrane thickness and d the distance between the two adhesive membranes. The polymer thickness is calculated as L = d/2. (a) The particle-count in the z-direction in the mushroom regime. (b) The particle-count in the z-direction in the brush regime.

Table 2 Dihedral interaction potential parameters

ϕ0(°)	Ki (kcal mol^−1)	ni
180	0.46845	1
0	0.046845	2
0	0.078872	3
0	0.028681	4

---

Table 3 Thickness of the polymer layer obtained by MD simulations for DPPC-PEG_np lipopolymer in DPPC bilayer membranes

Xp	Regime	LExp (nm)	LDM (nm)	hDM (nm)
0.005	Mushroom	3.75	3.6 ± 0.2	4.1 ± 0.05
0.014	Critic	—	4.3 ± 0.2	4.05 ± 0.05
0.1	Brush	6.6	6.4 ± 0.2	4.0 ± 0.05

---

3.2 Radial distribution function and mean-field interaction potential

The investigation of the three-dimensional distribution of the beads of DPPC lipids and PEG polymer, for the two grafting regimes, is carried out by analyzing the radial-distribution-function (RDF) from MD simulations in the NVTE statistical ensemble. RDF can be defined as the probability of finding 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 defined by the following thermal-average:where ρ is the number-particle-density. With ρ = 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 finding 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 confirms 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 graft 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-field interaction potential computed as U_eff/k_BT = −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 grafted 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 configuration. 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.

Fig. 6 Radial distribution function (RDF): (a) RDF of the group of lipopolymer molecules, (b) RDF of the group of lipid molecules.

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.

3.3 Dynamic aspect

3.3.1 Lipids and polymer dynamic. Membrane bilayer systems containing lipopolymers are idealized model systems designed to study fundamental physical properties of liposome. From an experimental point of view, these systems can indeed be synthesized and analyzed by different tools, such as neutron scattering or fluorescence correlation and pulsed-gradient spin-echo nuclear magnetic resonance. The MD simulation is capable of determining the type of diffusion in which the model membranes are decorated with a different molar fraction of membrane lipopolymers.
3.3.1.1 Lipids dynamic. The lipids dynamic can be described by analyzing the mean square displacement defined as:

The 〈Δ²r_CM〉 is found proportional to the diffusion coefficient D and the diffusion exponent α for a system of dimension d by:⁵⁶

〈Δ²r_CM〉 = 2dD_αt^α (17)

which can be expressed as:

log(〈Δ²r_CM〉) = log(2dD_α) + α × log(t) (18)

where the log function refers to the common logarithm. From the log–log plot of MSD versus time we determine the diffusion coefficients D_α and the exponent α < 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 significant times, that is, beyond a characteristic time depending on the specific details of the broadcast process and the host support structure. This long-term behavior deviates from the linear dependence on time found for Brownian motion.⁵⁷ Here, D_α represents the generalized diffusion coefficient or (the fractional diffusion coefficient). The latter is expressed in length²/time^α 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 significant times, that is, beyond a characteristic time depending on the specific details of the broadcast process and the host support structure. In general, a particle is said to be subdiffusive if the condition 〈Δ²r_CM(t)〉(t)/t→0, for t → +∞, is satisfied (very slow diffusion). This explains why the exponent α must be in the range 0 < α < 1. From an experimental point of view, subdiffusion has been observed in many scientific fields.⁵⁷ Subdiffusion originates from the cage effect,⁵⁷ where the tracer is trapped in a cage formed by its neighbors. It has been found that the length of stay, Δt_c, of a tracer in the cage, composed of N neighbors, is noted: Δt_c ∼ N^1/α, where α is the abnormal exponent.

3.3.1.2 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.⁵⁸ The dynamics of the PEG is modeled by the Brownian motion of EG-monomers. The polymer dynamics model, in the coarse-grained scale, was first 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₁, R₂,…, R_N ≡ R_n) is the position of the beads, k = 3k_BT/a² is the entropic spring force and f_n a random force. We assume that each bead feels the same friction ξ, that its motion is overdamped, and that the diffusion coefficient D = k_BT/ξ is independent of the position R_n of the bead.

The random force f_n verifies the conditions:⁵⁹

〈f_n(t)〉 = 0 (20)

〈f_n(t)f_m(t′)〉 = 6ξk_BTδ(t−t)δ_nm (21)

The solution of the Langevin equation gives the following behavior of the mean-square-displacement as a function of time.

with τ 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 center-of-mass is not normal. In this case, the mean-square-displacement of a monomer is given by:

〈Δr_CM²〉 = 〈[R_G(t) − R_G(0)]²〉 ≃ 〈[R_n(t) − R_n(0)]²〉 ∝ √t, t < τ (23)

3.3.2 Analysis of MD results. To study the influence of the molar fraction of the grafted polymer on their dynamics, and then on the membrane dynamics, we keep their degree of polymerization fixed 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 fixed at 300 K. The diffusion coefficients D_α, 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 first remark that after a short initial regime of subdiffusive motion (t = 100 ps), with a diffusion exponent α = 0.5, the polymer dynamics reach the normal diffusion, after this time MSDs are straight lines of the same slope α = 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 Experimental values of the diffusion coefficient of the PEG polymer were found by Waggoner et al.⁶⁰ For a low molecular weight, the diffusion coefficient is of the order of 10⁻⁶ cm² s⁻¹ and for a large molecular weight, it is of the order of 10⁻⁷ cm² s⁻¹. Our coarse-grained simulation for branched PEGs is in good agreement with the experimental findings. 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 (α < 1), α 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⁻⁷ cm² s⁻¹,⁶¹ MD results show that the polymer coating to the membranes has no important influence 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 ≳ 10) ns, as observed by Flenner and coworkers in an all-atom MD simulation of the DMPC bilayer membrane.⁶² 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 after 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 viewed as dispersed colloidal microspheres, the polymer plays an important role in preventing the liposomes aggregation and in facilitating their transport (Table 5).⁶³

Fig. 8 Mean square displacement (MSD): log–log of MSD as a function of time for DPPC molecules.

Fig. 9 Mean square displacement (MSD): log–log of MSD as a function of time for lipopolymers molecules.

Table 4 Diffusion coefficient D_α of DPPC lipids obtained by MD simulations in the NVTE statistical ensemble

Regime	t (fs)	α	(Dα) (cm^2 s^−α)
Mushroom	[0, 10^5.3]	0.50	7.4 × 10^−4
	[10^5.3, 10^6.4]	0.63	1.2 × 10^−4
	[10^6.4, 10^7]	1.00	1.4 × 10^−6
Brush	[0, 10^5.6]	0.42	6.9 × 10^−4
	[10^5.6, 10^7]	0.63	9.6 × 10^−5
	[10^6.6, 10^7]	1.00	8.9 × 10^−7

---

Table 5 Diffusion coefficient D_α of polymer obtained by MD simulations in the NVTE statistical ensemble

Regime	t (fs)	α	(Dα) (cm^2 s^−α)
Mushroom	[0, 10^5.2]	0.50	2.8 × 10^−4
	[10^5.2, 10^7]	1	1.4 × 10^−6
Brush	[0, 10^5.4]	0.34	2.2 × 10^−4
	[10^5.4, 10^7]	1	8.5 × 10^−7

---

4 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 graft polymers and shows a transition from the mushroom regime to the brush regime. When polymer chains are grafted 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 first approximation, the graft 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 graft polymer increases, the polymer beads of different chains interact and exhibit a more stretched configuration (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 findings from X-ray diffraction. From the NVTE MD simulations, the analysis of the radial distribution function has shown that the modification of the molar fraction of lipid graft 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 mean-square-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 flexibility 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 conflicts to declare.

Notes and references

1. M. Eeman and M. Deleu, Biotechnol., Agron., Soc. Environ., 2010, 14, 719–736.
2. J. T. Groves and R. Parthasarathy, J. Phys. Chem. B, 2006, 110, 8513–8516.
3. M. L. Wagner and L. K. Tamm, Biophys. J., 2000, 79, 1400–1414.
4. S. C. Hayden, A. Junghans, J. Majewski and M. A. Firestone, Biomacromolecules, 2017, 18, 1097–1107.
5. M. Hu, F. Stanzione, A. K. Sum, R. Faller and M. Deserno, ACS Nano, 2015, 9, 9942–9954.
6. R. N. Majzoub, K. K. Ewert and C. R. Safinya, Philos. Trans. R. Soc., A, 2016, 374, 20150129.
7. M. Tanaka and E. Sackmann, Nature, 2005, 437, 656.
8. J. L. Paris, M. V. Cabanas, M. Manzano and M. Vallet-Regi, ACS Nano, 2015, 9, 11023–11033.
9. K. P. Miller, L. Wang, B. C. Benicewicz and A. W. Decho, Chem. Soc. Rev., 2015, 44, 7787–7807.
10. Y. Kakimoto, Y. Tachihara, Y. Okamoto, K. Miyazawa, T. Fukuma and R. Tero, Langmuir, 2018, 34, 7201–7209.
11. M. Krishnamoorthy, S. Hakobyan, M. Ramstedt and J. E. Gautrot, Chem. Rev., 2014, 114(21), 10976–11026.
12. H. Liu, T. L. Doane, Y. Cheng, F. Lu, S. Srinivasan, J. J. Zhu and C. Burda, Part. Part. Syst. Charact., 2015, 32, 197–204.
13. D. P. Julien, A. W. Chan, J. Barrios, J. Mathiaparanam, A. Douglass, M. A. Wolman and A. Sagasti, J. Neurogenet., 2018, 32(4), 336–352.
14. Y. Sun, W. Z. Liu, T. Liu, X. Feng, N. Yang and H. F. Zhou, J. Recept. Signal Transduction, 2015, 35, 600–604.
15. W. A. Comrie, S. Li, S. Boyle and J. K. Burkhardt, J. Cell Biol., 2015, 208, 457–473.
16. L. Seguin, J. S. Desgrosellier, S. M. Weis and D. A. Cheresh, Trends Cell Biol., 2015, 25, 234–240.
17. T. T. Lee, J. R. Garcia, J. I. Paez, A. Singh, E. A. Phelps, S. Weis, Z. Shafiq, A. Shekaran, A. del Campo and A. J. Garcia, Nat. Mater., 2015, 14, 352–360.
18. J. M. Nam, Y. K. Lee, S. Kim, and K. Kim, US Patent application, 10,509,030, 17 Dec 2019.
19. P. Parkkila, M. Elderdfi, A. Bunker and T. Viitala, Langmuir, 2018, 34, 8081–8091.
20. M. Kim, M. Vala, C. Ertsgaard, S. H. Oh, T. P. Lodge, F. S. Bates and B. J. Hackel, Langmuir, 2018, 34, 6703–6712.
21. A. Luchini, Y. Gerelli, G. Fragneto, T. Nylander, G. K. Palsson, M.-S. Appavou and L. Paduano, Colloids Surf., B, 2017, 151, 76–87.
22. J. D. Unsay, K. Cosentino and A. J. Garcia-Saez, J. Visualized Exp., 2015, e52867.
23. B. Gumi-Audenis, L. Costa, F. Carla, F. Comin, F. Sanz and I. M. Giannotti, Membranes, 2016, 6, 58.
24. P. A. De Beule and A. Miranda, Anisotropy Imaging of Supported Lipid Bilayers via Spectroscopic Imaging Ellipsometry, in Optical Trapping Applications, Optical Society of America, 2015, p. JT3A-42.
25. V. Kiessling and L. K. Tamm, Biophys. J., 2003, 84, 408–418.
26. A. Lambacher and P. Fromherz, J. Opt. Soc. Am. B, 2002, 19, 1435–1453.
27. J. C. Shelley, M. Y. Shelley, R. C. Reeder, S. Bandyopadhyay, P. B. Moore and M. L. Klein, J. Phys. Chem. B, 2001, 105, 9785–9792.
28. J. C. Shelley, M. Y. Shelley, R. C. Reeder, S. Bandyopadhyay and M. L. Klein, J. Phys. Chem. B, 2001, 105, 4464–4470.
29. C. F. Lopez, S. O. Nielsen, P. B. Moore and M. L. Klein, Proc. Natl. Acad. Sci. U. S. A., 2004, 101, 4431–4434.
30. C. F. Lopez, S. O. Nielsen, G. Srinivas, W. F. DeGrado and M. L. Klein, J. Chem. Theory Comput., 2006, 2, 649–655.
31. S. J. Marrink, A. H. Vries and A. E. Mark, J. Phys. Chem. B, 2004, 108(2), 750–760.
32. S. J. Marrink, H. J. Risselada, S. Yefimov, D. P. Tieleman and A. H. De Vries, J. Phys. Chem. B, 2007, 111(27), 7812–7824.
33. S. J. Marrink and A. E. Mark, Biophys. J., 2004, 87, 3894–3900.
34. S. J. Marrink and J. Siewert, et al., J. Phys. Chem. B, 2007, 111, 7812–7824.
35. H. Lee, A. H. de Vries, S. J. Marrink and R. W. Pastor, J. Phys. Chem. B, 2009, 113(l), 13186–13194.
36. S. J. Marrink, A. H. De Vries and A. E. Mark, J. Phys. Chem. B, 2004, 108, 750–760.
37. L. Monticelli, S. K. Kandasamy, X. Periole, R. G. Larson, D. P. Tieleman and S.-J. Marrink, J. Chem. Theory Comput., 2008, 4, 819–834.
38. R. Baron, A. H. De Vries, P. H. Hunenberger and W. F. van Gunsteren, J. Phys. Chem. B, 2006, 110, 15602–15614.
39. G. Rossi, P. F. J. Fuchs, J. Barnoud and L. Monticelli, J. Phys. Chem. B, 2012, 14353–14362.
40. M. Pannuzzo, D. H. de Jong, A. Raudino and S. J. Marrink, J. Chem. Phys., 2014, 140, 124905.
41. S. Plimpton, P. Crozier and A. Thompson, Sandia Natl. Lab., 2007, 18, 43.
42. M. P. Allen, Computational soft matter: from synthetic polymers to proteins, 2004, vol. 23, pp. 1–28.
43. W. Shinoda, M. Shiga and M. Mikami, Phys. Rev. B: Condens. Matter Mater. Phys., 2004, 69, 134103.
44. R. L. Davidchack, T. E. Ouldridge and M. V. Tretyakov, J. Chem. Phys., 2015, 142, 144114.
45. T. F. Miller, M. Eleftheriou, P. Pattnaik, A. Ndirango, D. Newns and G. J. Martyna, J. Chem. Phys., 2002, 116, 8649–8659.
46. P. J. Flory, Principles of polymer chemistry, Cornell University Press, 1953.
47. P. G. De Gennes and P.G. Gennes, Scaling concepts in polymer physics, Cornell University Press, 1979.
48. S. Alexander, J. Phys., 1977, 38, 983–987.
49. E. Evans, D. J. Klingenberg, W. Rawicz and F. Szoka, Langmuir, 1996, 12, 3031–3037.
50. K. Hristova and D. Needham, Macromolecules, 1995, 28, 991–1002.
51. S. T. Milner, T. A. Witten and M. E. Cates, Europhys. Lett., 1988, 5, 413.
52. M. Daoud and P. G. De Gennes, J. Phys., 1977, 38, 85–93.
53. A. K. Kenworthy, K. Hristova, D. Needham and T. J. McIntosh, Biophys. J., 1995, 68, 1921–1936.
54. N. Kucerka, J. Pencer, M. P. Nieh and J. Katsaras, Biochim. Biophys. Acta, Biomembr., 2011, 1808, 2761–2771.
55. A. Stukowski, Modell. Simul. Mater. Sci. Eng., 2009, 18, 015012.
56. R. Kubo, M. Toda and N. Hashitsume, Statistical Physics II: nonequilibrium statistical mechanics, Springer Science and Business Media, 2012, vol. 31, p. 31.
57. M. Badia, S. El-Moudny, M. Benhamou and M. El Ossmani, J. Mol. Liq., 2017, 240, 1–30.
58. E. Rolls, Y. Togashi and R. Erban, Multiscale Model. Simul., 2017, 15, 1672–1693.
59. D. Osmanović and Y. Rabin, Soft Matter, 2017, 13, 963–968.
60. R. A. Waggoner, F. D. Blum and J. C. Lang, Macromolecules, 1995, 28, 2658–2664.
61. R. Metzler, J. H. Jeon and A. G. Cherstvy, Biochim. Biophys. Acta, Biomembr., 2016, 1858, 2451–2467.
62. E. Flenner, J. Das, M. C. Rheinstadter and I. Kosztin, Phys. Rev. E, 2009, 79, 011907.
63. A. Iborra, G. Diaz, D. Lopez, J. M. Giussi and O. Azzaroni, Eur. Polym. J., 2017, 87, 308–317.

---