A Coarse-Grained Molecular Dynamics Simulation Using NAMD Package to Reveal Aggregation Profile of Phospholipids Self-Assembly in Water

The energy profile of self-assembly process of DLPE, DLPS, DOPE, DOPS, DLiPE, and DLiPS in water was investigated by a coarse-grained molecular dynamics simulation using NAMD package. The self-assembly process was initiated from random configurations. The simulation was carried out for 160 ns. This study presented proof that there were three major self-assembled arrangements which became visible for a certain duration when the simulation took place, that is, liposome, deformed liposome, and planar bilayer. The energy profile that shows plateau at the time of these structures emerge confirmed their stability therein. Our findings have highlighted the idea that liposomes and deformed liposomes are metastable phases which eventually will turn into planar bilayer, the stable one.


Introduction
Solution of phospholipid molecules can demonstrate more than one micellar structures, namely, spherical micelles, rodlike structures, liposomes, bilayers, and others due to their surfactant-like features [1][2][3][4].These structures play important role in drug delivery systems as well as in biological systems [5,6].Micellar structures rely on the molecular species, composition, and also on the self-assembly pathways affected by the initial configuration [7][8][9][10].A great deal of experimental research has been done to study self-assembly of phospholipid molecules.Nevertheless the dynamics information about the liposome formation is still hard to achieve experimentally.
Molecular dynamics computer simulation has the ability to deliver more detailed information.It is an impressive device to investigate the mechanism of self-assembly [11][12][13][14].Conventional molecular dynamics (MD) uncovers maximum features however they are limited to small time scales.
It demands a long time to equilibrate a real physical system.Hence a coarse-grained molecular dynamics (CGMD) method was built as a simplified model to carry out molecular dynamics.The CGMD models have been used to explore a variety of structural and dynamic properties in large molecular systems.The CGMD method has offered significant outcome when exploring time and length scales further than what is viable with conventional MD.While CGMD has brought important findings for understanding the phospholipid self-assembly [12,14,15], there is still limited information accessible for a theoretical perceptive of phospholipid self-assembly pathway.
An all atomic simulation on phospholipid aggregation by Marrink and coworkers [13] has shown a typical pathway for bilayer formation.However, it did not state the formation of liposomes during the course of aggregation.Applying CGMD method on DLPE, DLPS, DOPE, DOPS, DLiPE, and DLiPS we demonstrate for the first time the aggregation profile of  phospholipid molecules which clearly show the formation of liposome as the metastable phase.These phospholipids have been reported as the main phospholipid component of coconut, sesame, and candlenut endosperm [16] which produced liposomes and planar bilayer during aggregation [2].

Methodology
The structure of all phospholipids used in the simulation is presented in Figure 1.The model molecule was prepared using the Open Babel package [17].In this study, 256 phospholipid molecules were placed randomly in a cube-shaped box with a size of 8 nm using Packmol package [18].Residuebased coarse-graining was applied on the system based on Martini Force Field ver.2.0 [15,19] using VMD package [20].The force field was parameterized to reproduce accurate thermodynamic properties [21].Each phospholipid molecule was represented as 10-14 beads.Water molecules were modeled by hydrophilic beads; each one represented four real water molecules.Figure 2 presents coarse-grained structure of all phospholipids used in the simulation.In Martini Force Field each bead interacts with the pair wise Lennard-Jones potential (LJ).Screen Coulomb interaction is used to model the electrostatic interaction between the zwitter ionic head groups of phospholipids.
Molecular dynamics simulations of phospholipid system began with energy minimization using NAMD package [22].Energy minimization was done to adjust the structure to the force field, the distribution of solvents, and especially to reduce the steric clashes that might occur in the system.This phase provided the system with the lowest energy to do the simulation.It had been marked by the achievement of energy convergency at the end of minimization, 0.6 ns.
After the minimization simulations were performed with 40 fs time step integration during the effective time of 160 ns.Simulations were conducted on periodic boundary conditions (PBC).The duration of liposome formation and the total energy systems were analyzed from the simulation results.Visualization during the simulation process was also done by VMD.The simulation was also undertaken for larger systems, that is, 1500 phospholipid molecules.
To evaluate our system, before running the simulation on the phospholipids we conducted the simulation on DLiPC which has been recognized to form liposome [23].

Results and Discussion
Molecular dynamics simulation of phospholipid molecules was able to provide an overview of the mechanism of the aggregation process and the relationship between aggregate structure and the total energy of the system.Simulation was performed on a system with 256 molecules of phospholipids in aqueous medium with density of 0.00609 atom/A 3 using NAMD package.For all molecules used in this report the simulation began with a random position.In general the molecules then started to form various clusters of phospholipids with hydrocarbon tails directed to their interior.After that they began to form liposomes or half-liposomes which were followed by deformed liposomes or planar bilayer formation at the end.We think the formation of half-liposome is due to shortage of phospholipid molecules supplies in the system.
To examine the influence of the number of molecules the simulation was also performed on the system with 1500 molecules with the same density.These simulations have shown the aggregation process better.It was preceded by the formation of small clusters of phospholipids which then was followed by mergers into larger aggregates, in the form of worm-like, cup-like, tube-like, and other structures, leading to formation of liposome or planar bilayer.For some phospholipids liposomes remained stable until the simulation ends.For others the process was then followed by deformation to produce deformed liposomes that lasted until they all became a planar bilayer.Examples of various aggregate structures observed in the aggregation process are presented as snapshots on Figure 3.
Observation of the total energy changes during the aggregation process showed that the process was accompanied by a decrease in the total energy of the system.The energy decrease occurred in stages before reaching a minimum when the simulation was terminated.This means that metastable structure is formed which is subsequently followed by a stable structure with minimum energy state.The simulations show that liposome and deformed liposome are metastable structures while a planar bilayer is a stable structure.This finding supports the views expressed by Luisi, Lasic, and Laughlin [24][25][26] that liposomes are metastable structures.For that liposomes also have a finite lifetime.The change of aggregation state which was accompanied by decreases of total energy system throughout the 160 ns simulation time  is represented by the aggregation of 256 DOPE molecules.The DOPE molecules showed an energy decrease when the system formed liposome, deformed liposome, and finally planar bilayer.These self-assembly structures appeared at the first, the second, and the last energy plateau on the energy profile of the simulation (Figure 4).It also showed clearly that DOPE liposome is a metastable system with lifetime for 24 ns.Simulations with a larger number of molecules showed that liposome formation occurs at a lower energy.
The transformation and the lifetime of the aggregate structures as well as the changes of system total energy throughout the period of simulation were varied for each molecular species.The data for simulation of 256 molecules of each phospholipid molecular species are presented in Table 1.
Based on the results in Table 1 we propose the aggregation mechanism of phospholipid molecules is through the stages as in Figure 5. Starting from a random configuration, there is a rapid decline in total energy and it is accompanied by formation of irregular phospholipid clusters during the decline.
After that the liposome is formed at a relatively stable total energy system, liposomes then change shape (deformed, nonspherical structure) followed by a decrease in the total energy of the system.Deformation of liposomal structures occurs for a certain time and at a relatively stable total energy.In the final stage, deformed liposome releases its total energy to form a planar bilayer which is stable until the end of the simulation.
The simulation suggests that liposomal structure is a metastable structure and the aggregation eventually will produce a planar bilayer at its final stage.It is in line with the notion suggested by Zana [27] that liposome formed from a single species of phospholipid cannot be thermodynamically stable since the bilayer will have high bending energy in the liposomal structure.On the contrary a planar bilayer can be thermodynamically stable due to its zero curvature.
The simulations also reveal how the hydrophilic head group, the number of carbon atoms in each carbon chain, and the number of double bond affect the aggregation stages  experienced by each phospholipid species and the lifespan of liposome produced.Observation on the total system energy of liposomes with ethanolamine (PE) head group shows that they have total system energy higher than the serine (PS) and choline (PC).These results are consistent with Martens and McMahon [28] findings on mechanism of membrane fusion which suggest that phospholipids with ethanolamine head group prefer negative curvature, while serine and choline are positive.

Conclusions
This work simulated the vesiculation pathway of several phospholipid species, namely, DLPE, DLPS, DOPE, DOPS,  DLiPE, and DLiPS from random initial configuration.We employ Marrink's coarse-grain model and NAMD package for this study.We have provided further evidence that there are three main self-assembled structures which materialize for a certain period of time during the simulation, that is, liposome, deformed liposome, and planar bilayer.The energy profile associated with the appearance of these structures indicated their stability.Our results have lent power to the hypothesis that liposomes and deformed liposomes are metastable phases which will then become the stable one, planar bilayer.

Figure 1 :
Figure 1: The molecular structure of phospholipids in the simulation.

1, 2 -Figure 2 :
Figure 2: Coarse-grained structures of phospholipids in the simulation.Bead of atoms are presented by colored beads.Ethanolamine head group are blue bead, brown for phosphate, glycerol backbone pink, green for hydrocarbon tail groups, and purple for the double bonds.
(a) Clusters of 1500 molecules DLiPE at  = 14.4 ns.Water molecules are not presented for clarity purpose only (b) Cup-like structure of 1500 molecules DLPE at  = 124.8ns.Water molecules are not presented for clarity purpose only (c) Tube-like structure of 1500 molecules DLPE at  = 96.96ns.Water molecules are not presented for clarity purpose only (d) Liposome of 1500 molecules DOPE at  = 21.6-45.6ns.This picture clearly shows some water molecules inside the liposome (e) Deformed Liposome of 1500 DLiPS molecules at  > 101.6 ns.This picture clearly shows some water molecules that still reside in the deformed liposome (f) Planar bilayer of 1500 molecules DLiPE at  = 133.6 ns

Figure 3 :
Figure 3: Various aggregate structures of phospholipid molecules during 160 ns simulation time.Bead of atoms in phospholipid molecules are represented by colored lines.Phospholipid head groups are green or blue lines, brown for phosphate, glycerol backbone pink and green hydrocarbon tail groups, and purple for the double bonds in the tail groups.Water molecules are presented by light blue beads.

Figure 4 :
Figure 4: The total energy profile of 256 DOPE molecules during 160 ns of simulation.

Figure 5 :
Figure 5: Stages of phospholipids aggregation with its accompanied changes of total energy system.

Table 1 :
Formation of aggregate structures of 256 phospholipid molecules during 160 ns simulation.