Tuesday, April 7, 2020

NMRlipids VI: Polarizable force fields

The primary goal of the NMRlipids Project is to find atomic resolution MD simulation force fields that correctly capture the lipid headgroup structures of biologically relevant lipids, and their interactions with ions and other biomolecules.

Current results from the NMRlipids Project indicate that none of the existing force fields correctly captures the lipid headgroup structures (NMRlipids I and IVa). However, the differences between PC, PE, PG and PS headgroup structures are roughly reproduced in CHARMM36, and the description of ion binding to PC and PS headgroups is substantially improved when electronic polarizability is implicitly included using the electronic continuum correction (ECC).

So far, the NMRlipids Project has focused on force fields that lack electronic polarizability. However, the number of available polarizable lipid force fields is increasing, and our PC and PS simulations with ECC suggest that the electronic polarizability may be an essential player in lipid–ion interactions. For these reasons, Batuhan Kav, an active NMRlipids contributor, has suggested the NMRlipids community to make a systematic review and benchmark study of the available polarizable lipid force fields. To this end, we hereby launch the NMRlipids VI project. It will follow the normal NMRlipids rules, with the exception that Batuhan Kav will mainly push the project and thus be the corresponding author.

Note that besides its primary goal, the NMRlipids Project has produced the largest publicly available collection of lipid bilayer MD simulations (indexed also at www.nmrlipids.fi) and evaluations of lipid headgroup force field quality against NMR experiments. To strengthen this side of NMRlipids, we will in NMRlipids VI test a new data contribution, indexing and analysis protocol that paves the way toward the planned leap to the NMRlipids Databank.

As in all NMRlipids projects, the contributions to NMRlipids VI can be made by commenting blog posts related to this topic or by contributing to the related GitHub repository. The GitHub repository already contains a draft review by Batuhan Kav on the published simulations on polarizable lipid force fields. In addition, a python script and stepwise instructions on how to contribute data are available.

In the published literature on polarizable force fields (see Batuhan's draft review), the acyl chains have already been evaluated against NMR order parameter data, but the quality of headgroup structures and ion binding remains largely untested. As a first step to do this test, we need trajectories from polarizable lipid bilayer simulations. During initial attempts to run lipid bilayer simulations with polarizable force fields, it turned out to be significantly more complicated than for non-polarizable force fields. Therefore, we specifically ask contributions from people, who already have data from lipid bilayer simulations with polarizable force fields, or know how to run these in practise. The polarizable lipid force fields that we are aware of are:
If you have access to lipid bilayer simulations using these or other polarizable lipid force fields, and are willing to contribute to the project, please contact us and/or contribute the data according to the instructions. Currently our analysis script works for data from Gromacs and NAMD, but our goal is to accept data from any MD engine. Do not hesitate to give us feedback on the data contribution and analysis, or to develop it further by yourself on GitHub.

Batuhan Kav
Markus Miettinen
Samuli Ollila


  1. Hello,

    I have generated a trajectory of pure POPC using the Drude polarizable model with OpenMM simulation package. The trajectory files are at


    and the analysis results are at


    The initial data shows that the Drude and Charmm36 force fields behave very similarly when it comes to the order parameters without any ions.

    One thing we can try to check is how these force fields differ with respect to the ion binding. From NMRLipids Lipid-Ion interactions project, we already know how Charmm36 behaves. My suggestion would be to run POPC simulations with the Drude model at 0, 350, 450, 650, and 1000 mM NaCl or CaCl2 concentrations. If no one is interested in running these simulations, I will run them.

    1. Where is the data for CHARMM36 POPC taken from? It does not look like a typical results we have got from this force field. Alpha and beta carbon in CHARMM36 are typically quite close to experiments and do not exhibit any forking.

    2. Hi Samuli,

      There was a mistake in the plot for the Charmm36 results. I have updated the plots on GitHub.

      Prior to a simulation with Drude force field, I need to run 100-200 ns of equilibration with Charmm36. I calculated the order parameters from this simulation, using the last 100 ns of a total 200 ns trajectory.

      The order parameters for the Charmm36 now seem close to what we expect.

    3. Thanks. It is interesting that the CHARMM Drude model has this large forking in the order parameters. In the parametrization publication they have tuned the headgroup dihedrals to capture these order parameters correctly. However, they report only absolute values and I cannot see any forking in their results.

      Anyway, I agree that these simulations with ions would be highly interesting. Are there any publications where simulations are ran with ions using this model?

    4. I'm not sure how they calculate the order parameters for alpha and beta hydrogens. I am also not aware of any publication for Drude lipids + ions.

      Initially I will look at the NaCl binding behavior of the Drude POPC model.

  2. Hello,

    I finished running the simulations of pure POPC (72 x 72) membrane at 350mM, 450mM, 650mM, and 1000mM NaCl (or CaCl2 -- 1000mM simulations for CaCl2 not available yet) concentrations using the Drude polarizable force field. The trajectories are available under Zenodo (links follow below). My plan is to finish analyzing them by early October.


    1. And the order parameters are uploaded to https://github.com/NMRLipids/NMRlipidsVIpolarizableFFs/tree/master/Figures

    2. Thanks, this is interesting. Here are some comments:

      1) The forking of alpha and beta order parameters seen already in figure without any ions that you linked in earlier comment is quite extreme compared to any other models we have seen (I am comparing to figure 2 in NMRlipids I publication). This indicates that headgroup may have very peculiar conformations in these simulations. Calculation of the headgroup region dihedral angles as we have recently done in NMRlipids IV (https://github.com/NMRLipids/NMRlipidsIVPEandPG/blob/master/Figs/DIHEDRALS-eps-converted-to.pdf) would be useful to check this. Because the data is in the databank format, this could be doable with the codes that calculate all the dihedrals through the databank:
      However, the first one has some performance issues (especially memory) and the latter works only for gromacs trajectories.

      2) I would be happy to see the changes of headgroup order parameters plotted as a function of sodium and calcium ion concentration against experimental data, as in figure 2 in NMRlipids II. Also ion density profiles along membrane normal would be highly interesting. Because we have nowadays the ion density profiles that reproduce the correct order parameter responses (Figure 5 in http://dx.doi.org/10.1021/acs.jpcb.7b12510), we could also compare to those.

      3) Data from POPC C36 450 mM CaCl2 seems an outlier. Maybe it would be worth of double check this dataset?

    3. Samuli wrote: Because we have nowadays the ion density profiles that reproduce the correct order parameter responses (Figure 5 in http://dx.doi.org/10.1021/acs.jpcb.7b12510), we could also compare to those.

      Are these ion density profiles available somewhere as text files?

    4. They should be in this repository, but I cannot find them now: https://github.com/ohsOllila/NMRlipids_VI-NewIonModel

      Maybe Josef Melcr could specify the location if they are there?

    5. The profiles are there, but it takes some search (sorry!) → you can reach them (including other profiles and analysis outputs) in the following folder:


      under the folder-names starting with "sim22a_*".

      �� For your convenience, I have collected the files (with sensible names) also into the ECC-lipids repository at:



  3. I have uploaded dihedral distributions, ion number density profiles, and order parameter changes results to


    For each Drude system I have only one trajectory. Therefore, currently the error bars in the ion density plots correspond to the standard deviation of the ion number per bin over the trajectory.

    1. Thanks again.

      It seems clear to me now that the CHARMM-Drude model gives very different headgroup and glycerol backbone conformational ensembles than the non-polarizable CHARMM, and that the order parameters from polarizable simulations are far from experimental values.

      I believe that the inaccurate headgroup conformational ensembles lead to the weird responses of order parameters to the bound ions shown in your figures.

      Based on rough comparison of calcium ion density profiles to the results in Figure 5 in http://dx.doi.org/10.1021/acs.jpcb.7b12510, it seems to me that the binding affinity of calcium could be similar to the model which gives the correct order parameter responses.

      However, based on both order parameter changes and density profiles, the sodium binding to membranes seems to be clearly stronger than calcium, which is not in line with experiments.

      In conclusion, it seems to me that the inclusion of polarizability do not automatically improve the headgroup conformational ensembles and ion binding affinitied to membranes.

      I think that next it would be highly interesting to have data also from other polarizable force fields than CHARMM.

  4. I also agree, based on the current results, that the addition of polarizability does not readily improve the head group conformations.

    Regarding the ion distributions, I need to say that I have used only one trajectory per system. The errors shown in the figure (https://github.com/NMRLipids/NMRlipidsVIpolarizableFFs/blob/master/Figures/ca_number_density.png) are essentially the standard deviation of the number density for a given bin throughout one trajectory. Do you think that having maybe one or two more simulations will help us to get a better error estimate?

    For the dihedral distributions, I realized that my results are differing from the data shown in (https://github.com/NMRLipids/NMRlipidsIVPEandPG/blob/master/Figs/DIHEDRALS-eps-converted-to.pdf) by 180 degrees in the x-axis. I have used my own Charmm36 POPC trajectories to confirm that the dihedral distributions I calculate for the Drude model are correct, however I always get a 180 degree shift in my results. This could be from a very obvious reason but my mind seems to be blocked here to understand why this is happening.

    One other model that we can consider for the polarizable models is the AMOEBA force field. That'd be really helpful if anyone has some data for that force field and share with us.

    I'll be updating the manuscript (https://github.com/NMRLipids/NMRlipidsVIpolarizableFFs/blob/master/Manuscript/manuscript.pdf) to include the latest results.

    1. In simulations with calcium, also the slightly different densities in different leaflets indicate that the statistics could be better. More simulations should help, but I am not sure how accurate results we need from this for this project.

      In your dihedral plot the x-axis goes from -180 to 180, and in NMRlipids the axis is from 0 to 360. Here, -180 to 0 means angles from 180 to 360. I am not sure if your dihedral are still shifted by 180 after doing this. I am not sure how you calculate dihedrals, but another issue may the order of atoms given in the code. Many atom selection tools, including gmx make_ndx, takes atoms in the order they appear in coordinate/topology files. This may not be the correct order for the dihedral calculation. For this reason, I make the atom selection in a slightly complicated way in the dihedral calculation codes: https://github.com/NMRLipids/NMRlipidsIVPEandPG/blob/master/scripts/calcDIHEDRALS.py, https://github.com/NMRLipids/NMRlipidsIVPEandPG/blob/master/scripts/calcDIHEDRALSgromacs.py


  5. Hello,
    I will prepare the input files for pure POPC (64 x 64) membrane at 350mM, 450mM, 650mM, and 1000mM NaCl ( or CaCl2 ) and neutral salt concentration, using AMOEBA FF.
    I will submit the setup files to Batuhan.

    1. Hello,

      Thanks, Suman, for the contribution.

      Just to clarify any possible confusion: me and Suman Samantray (IBI-7, FZ-Juelich) will be working on the AMOEBA simulations of the POPC lipids together.

      Again, thank you!

    2. Hi Batuhan, I have been wanting to test out OpenMM myself, also in combination with AMOEBA as previously discussed. Being new to both, this combination however has a learning curve to it.. I found some useful examples with lipids (POPC), but I need to update my OpenMM for the AMOEBA support (version>7.5) on my uni's cluster, which must be done by our support staff. I can let you know how I progress.


  6. Hi Jesper, thanks again for the contribution! I'm looking forward to seeing the data.


Please sign in before writing your comment.