NMRlipids VI: First results from polarizable Charmm-Drude force field

After analyzing the simulation data for the POPC lipids at various NaCl and CaCl₂ concentrations using the Drude polarizable model, I have compiled all the results in the manuscript. You can access all the figures, scripts, and the source files used in the analysis in the GitHub repository. Here I’d like to summarize the main results we have obtained so far.

1) The headgroup and glycerol backbone order parameters predicted by Drude force field are very different than the ones with Charmm36 force field.

The Fig. 1 in the manuscript show that Drude model predicts a forking at the beta and alpha order parameters, which are not predicted by the Charmm36 force field. 

Figure 1: Headgroup and glycerol backbone order parameters of POPC without ions from Charmm-Drude and Charmm36 simulations

We further see that the response of these order parameters to changing Na⁺/Ca²⁺ concentrations also deviate significantly from the Charmm36 results, as shown in Fig. 2.

Figure 2: Changes in the headgroup order parameters upon addition of ions in Charmm-Drude and Charmm36 simulations.

 

2) The dihedral distributions around the lipid headgroup and glycerol backbone have completely different profiles compared to the Charmm36 force field.

In Figs. 3 and 4, we show the dihedral angle distributions around the lipid headgroup and glycerol backbone at different ion concentrations. We observe the largest differences, as a function of the ion concentration, at the O_g3-P-O_α-C_α and the g₃-O_g3-P-O_α pairs. For the remaining of the distributions, we don’t observe large changes.

When compared to the Charmm36 results [Fig.2 of NMRlipids IVb manuscript], we see that the Drude model yields completely different dihedral angle distributions.

3) The density profiles for the ions suggest that the sodium binding is too strong.

Following Samuli’s comment: "Based on rough comparison of calcium ion density profiles to the results in Fig.5 in http://dx.doi.org/10.1021/acs.jpcp.7b12510, it seems 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 line with the experiments."

Overall, our calculations show that the Charmm Drude polarizable model does not improve on the results of the Charmm36 force field. In fact, we see that the Drude polarizable model yields worse results than Charmm36, as evident from the forking of the beta and alpha order parameters in Drude model. This suggests that inclusion of the polarizability does not automatically improve the headgroup conformational ensembles and ion binding affinities to the membranes. Based on the current evidence, we are not planning to run any more simulations using the Drude model (for the POPC lipids) and switch our focus to the other polarizable force fields.

From now on, we are planning to concentrate on obtaining simulation results for the AMOEBA and other polarizable force fields. Still, we think it could be very helpful to have data from Charmm-Drude simulations (also older versions) for cross-checking our results. To be consistent with the Drude simulations, our first focus will be the pure POPC membranes at different CaCl₂ and NaCl concentrations. It seems like we will be obtaining some trajectories soon. However we kindly invite anyone who might have some trajectories with the AMOEBA and other force fields to share them with us.

Any comment on the current results or data contributions, as always, will be greatly appreciated.

 

Batuhan Kav

NMRlipids IVb: Toward submission of the manuscript with PE and PG results

After the previous post on the NMRlipids IVb manuscript, I have made further analysis on the structural differences between PC, PE, PG and PS headgroups using the development version of the NMRlipids databank and CHARMM36 simulations, which give the headgroup and glycerol backbone order parameters closest to experiments for all lipids. Based on these analyses, I have reformulated the NMRlipids IVb manuscript to focus on the headgroup conformations of different lipids.

Main points of the manuscript are now:

1) Differences in NMR order parameters between different headgroups can be explained by the changes in dihedral angle distributions, suggesting that similar conformations are accessed by all headgroups, but with different probabilities.

CHARMM36 simulations approximately reproduce the main differences of headgroup order parameters between PC, PE, PG and PS lipids (forking of the ⍺-carbon in PS and positive β-carbon value in PG in Figs. 1 A and B). These differences can be explained by the differences in the heavy atom dihedral distributions (Fig. 1 C).

Figure 1: Headgroup and glycerol backbone order parameters of POPC, POPE, POPG and POPS from (A) experiments and (B) CHARMM36 simulations. (C) Heavy atom dihedral angle distributions from CHARMM36 simulations.  

2) The changes in headgroup order parameters upon addition of charged molecules can be explained by the changes in dihedral angle distributions close to the phosphate region.

CHARMM36 simulations reproduce the experimental changes in headgroup order parameters upon addition of cationic surfactants which are related to the substantial tilt of P-N dipole (Fig. 2 A) and changes in dihedral distributions close to the phosphate region (Fig. 2 B). 

Figure 2: Response of POPC headgroup A) order parameters, P-N vector and B) dihedral angles to the addition of cationic surfactant from CHARMM36 simulations compared with experimental data.

3) Wide range of headgroup conformations are observed also in lipids that are directly bound to proteins, suggesting that the specific binding of lipids to proteins is dominated by the intermolecular lipid-protein interactions rather than the differences in conformational restrictions of lipids.

Dihedral angles calculated from protein bound lipid structures from the protein data bank (PDB) have wide distributions independently on the headgroup chemistry (Fig. 3).

Figure 3: Dihedral distributions of protein bound lipids analysed from the protein data bank (PDB).

Comparison of headgroup and glycerol backbone order parameters between different force fields and experiments, and results from lipid mixtures are now in the supplementary information because conclusions from these are essentially the same as in NMRlipids I and IV publications. None of the force fields correctly capture the conformational ensembles of lipid headgroups, but CHARMM36 is closest, and simulations containing charged lipids are complicated by the overestimated binding affinity of counterions. 

I believe that the manuscript is approaching the submission stage. Because MD simulations from the NMRlipids project are able to explain the differences in experimental order parameters between different lipids, thereby giving an interpretation of the lipid conformational ensembles, my current plan is to submit the manuscript to the ACS Central Science.

Most critical things to do before submission are related to the methods, see todo points in the manuscript and especially in the supplementary information. There are also some issues still open in the GitHub. All kind of comments on the manuscript are welcomed.

Online meeting about the NMRlipids databank

The first NMRlipids online meeting, focused in the development of the NMRlipids databank, was held at 16.00 CET on Monday 29th of June 2020. The meeting started with a presentation by Samuli Ollila about the current status of the databank (slides available in here). The schematic structure of the new databank is shown in Fig. 1.

Figure 1: Schematic structure of the new databank. Beta versions of the Databank, Databank builder and Databank analyzer codes are available.

The presentation was followed by a highly useful discussion, thanks to more than 20 participants. The discussion was mainly focused on urgent issues brought up in the presentation that were complemented by additional points raised by the participants. The outcomes of the meeting and some decisions based on the discussions are listed here.

1. What information will be stored into the dictionary files composing the databank? Current plan is to include information requested from contributor that cannot be read afterwards (force field information, trajectory length, etc.), and information necessary for using (file names and sources) and searching (number of molecules and temperature) the data. Note, however, that the tpr (or corresponding) and trajectory files are accessible through the databank. Thereby all the information of each simulation is available even thought everything is not written directly into the dictionary. For detailed discussion, see the GitHub issue.
2. How molecules will be named? When writing and searching the data from the databank, we need unique machine readable names for molecules. There will be a list of molecule names (for example, POPC, POT, TIP3P, etc.) that will be used by default. If the uploaded simulation has different names, user has to tell those. For detailed discussion, see the GitHub issue.
   
3. Unique convention for the atoms within the molecules. For now, we will use the idea of mapping files updated with a third column that tells the residue name for each atom. This should be useful in situations where parts of one lipid are named with different residue names, such as in the current Amber force field convention. For detailed discussion see the GitHub issue.
4. File format for the dictionary. If practically feasible, we will consider saving dictionary in yaml format instead of json. For detailed discussion see the GitHub issue.

As a first step, we will build a prototype databank containing simulations from NMRlipids IVb manuscript and use this to analyze, for example, P-N vector and dihedrals angles required for the manuscript. Therefore all the related issues are now in the GitHub repository of NMRlipids IVb.

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:

• CHARMM-Drude
• AMOEBA for zwitterionic and anionic lipids
• CHARMM-Fluctuating Charge (FQ)

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