DPD Exercise 2: Lipid bilayers, micelles and vesicles

Introduction

Following on from DPD Exercise 1, we can extend DPD by joining beads together into molecules using e.g. harmonic springs. Since we can vary conservative force parameters between different bead species, we can therefore construct molecules that have different regions of behaviour, e.g. amphiphiles with hydrophilic (water-loving) head groups and hydrophobic (water-hating) tails. In solution, these molecules can readily assemble into large-scale structures such as micelles, vesicles, bilayers and membranes.

Among others, Shillcock and Lipowsky [Shillcock2002] have proposed DPD as a viable modelling method for biological molecules such as lipids and proteins. In particular, DPD is able to adequately model the essential topology of lipid molecules in terms of local hydrophobicity and bond interactions, which can also include controlling the angles formed between pairs of bonds.

Using atomistic or coarse-grained molecular dynamics to model how lipid structures form is possible but time-consuming due to the large number of atoms required, mostly for the solvent. This is less of a problem for DPD, since the softer potentials available allow for the use of larger timesteps to get to equilibrium structures more quickly.

Aim

In this exercise you will start with a simple DPD-based model for an amphiphilic lipid molecule. This consists of a single hydrophilic bead as a head group and 6 hydrophobic beads as a tail, all joined together with bonds. The solvent (water) is represented by unbonded beads.

Our model for the lipid molecule includes harmonic bonds between the beads:

(1)\[U_{ij} = \frac{\kappa}{2} \left(r_{ij} - r_0 \right)^2\]

with a spring force constant \(\kappa\) and an equilibrium bond length \(r_0\). We are also using a cosine potential between pairs of bonds to control the angle between them:

(2)\[U_{ijk} = A \left[ 1 + \cos{\left( m\theta_{ijk} - \theta_{0} \right)} \right]\]

where \(\theta_{ijk}\) is the angle between particles \(i\), \(j\) and \(k\), \(A\) is the energy level, \(m\) is a multiplier and \(\theta_0\) is an equilibrium angle.

The straightness of bonds and the concentration of lipid molecules will determine the mesoscopic structures that they form. The two most likely kinds of structures for these molecules are micelles – generally spherical structures with hydrophobic tail groups pointing inwards towards the centre – and bilayers consisting of two planar layers of lipids with the tail groups pointing towards each other. Vesicles (or liposomes, if they are artificially created) are spherical bilayers that include some of the solvent, which is accommodated by the existence of hydrophilic head groups inside the structures.

Any initial configuration for the lipids and solvent can be used in a DPD simulation, as they will eventually arrange themselves into the lowest potential energy configuration, i.e. an equilibrium state.

Instructions

For this exercise, you will need the main DL_MESO_DPD executable dpd.exe - ideally the parallel version with MPI - and either traject_vtf.exe or traject_xml.exe to visualise the simulations along with VMD or OVITO respectively. Optionally, you can also make use of the dlmresultviewer.py Python script to look at calculated system properties, although you can also use graph plotting software that can read tabulated text files (e.g. Gnuplot, Excel).

Download CONTROL and FIELD into your working directory: these are input files for a smaller version of one of the test cases supplied with DL_MESO, which you can find in dl_meso/DEMO/DPD/LipidBilayer [1]. Both versions include three species of bead - solvent beads W, head beads H and tail beads C - and defines molecules (HC6), each consisting of one H bead and 6 C beads. The harmonic bonds between the beads - described by (1) - use a spring force constant \(\kappa = 128 k_B T / r_c\) and equilibrium bond length \(r_0 = 0.5 r_c\), while the bond angles described by the cosine potential in (2) use an energy level \(A = 20 k_B T\), multiplier \(m = 1\) and equilibrium angle \(\theta_{0} = 0{^\circ}\).

Since we want to run some variants of this simulation, create four directories inside your working directory, one for each simulation you are about to run:

mkdir LIPID1
mkdir LIPID2
mkdir LIPID3
mkdir LIPID4

and copy the CONTROL and FIELD files into each of these directories. The CONTROL file can be used for all four simulations without modification, while we are going to modify the FIELD file for three of the simulations to modify the angle interactions and lipid molecule concentrations.

To start with, go into the LIPID1 directory and take a look at the CONTROL and FIELD files with your preferred text editor. The current form of the FIELD file with the specified interactions and numbers of particles - both the solvent (W) and the amphiphilic molecules (HC6) - should eventually produce a lipid bilayer.

Now change into the LIPID2 directory and open the FIELD file using a text editor. We first want to see the effect of removing the angle interactions, so either delete the line beginning with angles and the five lines that follow, or change the energy levels from 20.0 to 0.0. (The latter will allow you to monitor the bond angles but will not apply the potential.) After saving this file and exiting the text editor, copy this FIELD file into the LIPID4 directory to overwrite the file already in there.

Change into the LIPID3 directory and open the FIELD file. We now want to observe the effect of lipid molecule concentration on the structure that forms, so change the number in the line beginning with nummols to a lower value. Note that we still want the same total number of beads in the simulation and the same particle density, so add some more water beads in the appropriate line under SPECIES.

Change into the LIPID4 directory, open the FIELD file and make the same changes to the numbers of molecules and water beads as you did for the third simulation.

Now run each of the simulations in turn by entering each directory and using one of the following commands:

../dpd.exe
mpirun -np 4 ../dpd.exe

to launch the executable found in the parent directory [2].

When each calculation has finished, run the required trajectory conversion utility to produce your file(s) to launch in VMD or OVITO from the trajectories in the HISTORY files. Again you can run the utility (located in your working directory) in each sub-directory:

../traject_vtf.exe
../traject_xml.exe

Analysis

  1. Take a look at the first simulation: visualise it with VMD or OVITO, and also take a look at the total and/or potential energies stated in OUTPUT or CORREL, using either dlmresultviewer.py or graphing software for the latter.
    • Does a lipid bilayer eventually form? If so, does its formation coincide with the lowest energy levels for the system?
    • The final column in CORREL gives averaged bond angles in degrees: how close are they to the desired 180°?
  2. Now take a look at the second simulation.
    • What effect does removing the angle interactions have on the structure that forms, compared with the first simulation? Can you think of a reason why it makes a difference?
    • What are the averaged bond angles in this case?
  3. Take a look at the third and fourth simulations with reduced molecular concentrations.
    • What structures do the lipid molecules form? Does either structure include water inside, i.e. does a vesicle form?
  4. You might be interested to see how each of these structures forms: do you observe any intermediate (metastable) structures before the system reaches an equilibrated state?
  5. Optionally, take a look at the original input files supplied with DL_MESO in DEMO/DPD/LipidBilayer and compare them with those you have been using.
    • If you have time, try a simulation with an intermediate volume and total number of particles, i.e. try doubling the system you have been running so far, and try different numbers of lipid molecules (but keep the total number of beads the same). How many molecules do you need to form a complete bilayer (without any holes)?
    • Given the differences in cubic volumes between these three simulations, all of which produce complete bilayers, can you devise some kind of relationship between the simulation box size and the number of amphiphilic molecules needed to form a bilayer?

Footnotes

[1]You could use the original input files from the LipidBilayer test case for this exercise if you have access to a more powerful computer, e.g. a high-performance computing (HPC) platform. Exactly how to log onto the machine, transfer the DL_MESO source code onto it, compile DL_MESO_DPD and launch calculations will depend on how it is set up: you may need to consult the available documentation and/or the system administrator to find out this information.
[2]If you are using a computer with multiple processor cores, we would strongly recommend running the simulations on at least four cores: each calculation should take about 10 minutes to complete on a four-core laptop, based on tests carried out on an available machine. Even if your computer only has two physical processor cores, it might be capable of hyperthreading, i.e. using each of your computer’s threads as a virtual (logical) core, and you may therefore still get some extra performance by specifying 4 cores in the mpirun command. If you do need to hyperthread, depending on your MPI implementation you may need an additional command line option with mpirun to use hyperthreading or to ‘oversubscribe’ tasks.

Reference

[Shillcock2002]JC Shillcock and R Lipowsky, Equilibrium structure and lateral stress distribution of amphiphilic bilayers from dissipative particle dynamics simulations, Journal of Chemical Physics, 117, 5048-5061 (2002), doi: 10.1063/1.1498463.