With the rise in computing power and cluster accessibility, computer simulations experience an increase in popularity amongst all kinds of sciences. They are used to explore problems from the macroscopic to the microscopic world with the help of many different methods, including Molecular Dynamics (MD) simulations. Simulations enable research of new topics like movements and interactions between atoms, which can be specifically useful for large biomolecules like lipids. While these simulations are limited to a very short period of time, they deliver valuable information. Interactions between molecules, and thus mechanisms-of-action concerning different biological processes, can be examined in a more detailed fashion compared to standard practice. Running simulations under different initial conditions with a precise control over all parameters helps reduce noise and promotes a clearer analysis of the systems behavior. This is only possible in silico as some conditions are tough to create and sustain in a lab environment (Hollingsworth & Dror, 2018; Karplus & Petsko, 1990).
The procedure of a MD simulation revolves around calculations of forces that act upon the individual atoms in the system. Forces of attraction and repulsion result in a net force that can be calculated and then used in order to gain a function describing the movement of an atom between two points in time. The movement follows the Newtonian Law of Motion, where the force applied to the atoms results in them being accelerated according to their mass:
A common approach for calculating the movement of atoms in MD simulations is Verlet integration. It is a numeric method used to solve the kinetic equation, without the need to know the precise velocities of the simulated atoms, as it can be described by only the acceleration as well as current and previous position of the atom. As we calculate all forces acting on said atoms, we can use the Newtonian Law of Motion to gain insight into the acceleration of atoms, given their mass. This gives us a relatively fast approach for calculating movements, without the error becoming overwhelmingly big (Hollingsworth & Dror, 2018; Karplus & Petsko, 1990).
In order to calculate the aforementioned forces, a so-called forcefield is needed. Force fields are collections of equations which can be used to model the bonds and geometries of the simulated molecules. We used the Martini Force-Field, which is a general purpose force-field specialized in coarse grain simulations. In a Martini model the molecules simulated are coarse-grained in a 4:1 ratio, meaning that 4 atoms which are heavier than hydrogen are combined to one node. This transformation requires specialized force-fields with specific parameters (Hollingsworth & Dror, 2018; Karplus & Petsko, 1990).
In our case, studying the hydrodynamic behavior and formation mechanisms of liposomes in silico was helpful to gather insights which were valuable for our wet lab experiments, including production of liposomes and developing a microfluidic device. By simulating the interactions between lipids and their liposome formation mechanism, we were able to evaluate the lipid composition used by the wet lab in a proof-of-concept approach. The time needed for the lipids to form liposome-like structures is a great first indicator for the stability and usability of these liposomes for our use-case. Additionally, further analysis like the evolution of total energy or root-mean-square deviation were conducted to examine the simulation’s output. We opted to simulate liposome formation using an explicit solvent. Using this approach, the positions of all water molecules from the solvent in the water box created around the lipids are also calculated. While this requires more computing power and time, it delivers more accurate results. Especially during our research, we were interested in the interaction between the water and the newly formed vesicles, which led to us using the explicit solvent. In contrast, an implicit water model would replace the individual water molecules by a homogeneously polarizable medium. While this reduces the runtimes and the needed computing power, it also results in a less accurate representation of the interactions between water molecules and the simulated molecules. In some cases this approximation can be acceptable and using this approach one could then profit from faster simulations (Hollingsworth & Dror, 2018; Karplus & Petsko, 1990).
For the purpose of simulating lipids, we used an approach previously mentioned by Parchekani et al. (2022). As we focused on the macroscopic formation of vesicles, the added complexity of an all-atom simulation was not needed. The best choice for this kind of simulation is a coarse-grain approach. The lipid chains are represented by a reduced number of nodes, making the simulation more effective due to less computing needed. The coarse-grained data of lipids was obtained from the web-service CHARMM-GUI (Jo et al., 2008; Brooks et al., 2009; Lee et al., 2016). The molecules were then placed into a cubic water-box of 10 nm in length, together with 0.15 M NaCl.
Also relying on CHARMM-GUI, we used the Martini Maker (Qi et al., 2015; Hsu et al., 2017) based on the Martini force 2.2 field to construct all necessary files for the implementation of the simulation, including processes like charge neutralization, solvatation and equilibration. The Martini Random Maker distributes coarse-grained lipids randomly in a defined cubic box. Afterwards, we performed the simulation using the GROMACS software package (version 2022.2-cuda-11.6), an open-source simulation suite which can be used for a variety of molecular dynamic simulations. All necessary files for running the simulation can be found on our gitlab page. The lipids simulated by GROMACS were chosen according to the lipids we planned to work with in the wet lab in order to evaluate the choice. There, DOTAP, DOPE, mPEG2000 and cholesterol are used. For our simulation, lipids with similar physical properties and structures were chosen for a proof-of-concept approach. Therefore, we used DOPC, DOPE and cholesterol as shown in table 1 to run a simulation for 1000 ns at 303,15 K. Next to the lipids, 6339 water molecules and 73 NaCl particles were included.
| Lipid | Concentration [mM] | # Molecules in MD simulation |
| DOPC | 46 | 28 |
| DOPE | 42 | 25 |
| Cholesterol | 100 | 60 |
Following, the simulation output was analyzed with the Python package MDanalysis (version 2.2.0) (Michaud-Agrawal et al., 2011). The main usage of this package was to extract the trajectory data from the compressed files created during the MD simulation and to create coordinate files (.pdb) from the trajectory file (.xtc). This enables us to comprehend the movement of the lipids during the simulation and evaluate the formation of vesicles. A python notebook including these steps can be found on our gitlab page as well.
Next to those analysis steps, different simulation parameters like total energy, potential and root-mean-square deviation of atomic coordinates, or just root-mean-square deviation (RMSD), were performed with the purpose of evaluating the simulation’s output.
The simulation can be visualized at different times to examine the movement of lipids. For this purpose, the VMD software (http://www.ks.uiuc.edu/Research/vmd/) was used. These visualizations can be seen in figure 1. As explained above, every node represents one coarse-grained group. In figure 1D, for example, three nodes can be seen at the top of the image. Due to the periodic boundaries, lipids that leave the box at the bottom appear at the top, which explains this observation, also appearing in other images.
The initially randomly distributed lipids (A) seem to form three distinct liposome-like structures after 10 ns (B), which continue to merge at around 30 ns (C). After about 50 ns all lipids are included in one round structure (D). In the following 950 ns, the structure moves around the box without any major structural changes. The lipid distribution after 1000 ns can be seen in figure 2.
Also relying on VMD, a simulation of the first 100 ns was produced and can be seen below. Again, bonds ranging through the whole box can be explained by periodic boundaries, which let parts of lipids leave the box on one side and re-enter it on the opposite side. Due to file-size and file-type, we were only able to upload the animation you can see below. Another, more elaborate animation can be found here: md_animation_final.
To evaluate the simulation’s output, the total energy, enthalpy and temperature as well as the pressure was plotted against the simulation time. This is shown in figures 4 and 5. For the temperature and pressure plots, the median value is included.
The formation of liposomes is driven by non-covalent interactions, which are attractive Van-der-Waals and repulsive electrostatic interactions. While every single one of these interactions has a low energy content itself, the addition of many of them is having a great impact on the formation of liposomes. During the formation process, the total energy and the enthalpy of the system decrease first, and reach a constant value when the system is stable and equilibrated. Therefore, the formation of liposomes and their stability can be observed via the total energy and enthalpy changes in the course of the simulation (Parchekani et al., 2022).
According to the diagrams in figure 4, the total energy as well as the enthalpie is decreasing drastically in the beginning while the liposome is formed, indicating a suitable formation of the lipids. Following, the energy diagrams show no more drastic changes which shows the stable state of the formed structure and the reached equilibrium. However, there is an intense fluctuation visible in both diagrams after reaching the equilibrium.
In figure 5, the temperature and pressure diagrams can be seen. The temperature shows a median value of 302.7 K, which is slightly less then the set temperature of 303.15 K. The pressure’s median value is 0.47 bar.
Another important parameter of the simulation is the root-mean-square-deviation (RMSD). The RMSD of two structures which are superimposed and have the same amount of atoms (or in our case, the same amount of coarse-grained nodes) is computed with the following equation:
where di is the distance between the i-th atom in the two structures and N is the total number of equivalent atoms (Sargsyan et al., 2017). Therefore, the RMDS shows how great the deviation of atoms or atom groups is compared to the reference point. When choosing the initial state of the system as the reference point, the RMSD slope over time indicates the structural stability of the liposome during the simulation. When the slope is close to zero after reaching an equilibrium, the formed structure is stable. However, if the slope still changes drastically, the simulated structure will be more unstable (Parchekani et al., 2022). Figure 6 shows RMSD diagrams of cholesterol (A), DOPC (B), DOPE (C) and of the whole system (D). In all cases, the equilibrium is reached quickly in less than 100 ns according to the results mentioned above. The slopes of diagrams 5A, 5B and 5C do not change drastically after reaching the equilibrium, however, there are still fluctuations visible. When viewing the diagram 5D (whole system), no such fluctuations are visible after the liposome formation.
MD simulations offer the possibility to visualize movements and formation mechanisms of biomolecules like proteins or lipids in order to predict their in vivo or in vitro behavior. Here, lipids similar to the composition used for the production of liposomes in the wet lab were simulated for the purpose of evaluating the chosen lipids. The simulation output showed the formation of a liposome-like structure during the first 100 ns as well as the persistence of this structure over the remaining 900 ns. The analysis of different simulation parameters like the total energy or the RMSD made it possible to take a closer look at the stability of the simulated model.
The simulation was performed in a cubic box of 10 nm in length. Altogether, 113 lipids were put into the box together with 6339 water molecules, all in coarse-grained state. Due to the small size, the results cannot be used to exactly picture the formation of liposomes as conducted in the wet lab. However, the quick formation of a round structure and its persistence indicate that the chosen composition could be suitable in larger dimensions for liposome production as well. The analysis of energies and RMSDs also showed the quick setting of an equilibrium, which seems to be quite stable. However, fluctuations were clearly visible regarding the total energy, enthalpy and RMSD after the formation of the round structure indicating that the lipid composition could be improved further before being used for production. As stated further above, the lipid composition simulated is not the same as the one proposed for liposome production in the wet lab. For example, this proposed composition also includes PEG-lipids, which are reported as a major component of LNP-siRNA delivery systems (Kulkarni et al., 2019; Evers et al., 2018).
To conclude, the simulation of lipids similar to our wet lab experiment proved the concept of the approach in general. Moreover, especially the animation and the visualizations of different time frames helped us gain new insights on how the liposomes’ formation process works. Moreover, even though our work only includes a simulation in small dimensions, it still can be a manual for simulating lipids for the iGEM community. All steps taken during the simulation and analysis process are documented and can be found on our gitlab page including input files for the simulation, commands used and instructions to install all necessary packages.