In the first part of the NMRlipids VI project on polarisable force fields, we have focused on the CHARMM-Drude polarisable force field and analyzed the POPC bilayers at various NaCl and CaCl2 concentrations as well as POPE bilayers without salt. Our results have been summarized in a previous blog posts and the manuscript in GitHub. As a quick recap, we found that CHARMM-Drude predicts a forking for the head group beta and alpha order parameters, which is not correct and not present in CHARMM36. We confirmed that this is due to the fitting of order parameters to the wrong target values. We further found that sodium binding to POPC bilayers is too strong. As the fitting-target data for the CHARMM-Drude force field is not correct, we could not get useful information on the effect of polarizability. During the latest NMRlipids meeting in September 2021 in Prague, we decided to include POPE and POPC simulations with the AMOEBA force field into the project if possible. This was mostly due to the public interest in the AMOEBA force field, and to get a more complete picture of the available polarisable models; furthermore, we wanted to see if ion binding is good in this polarizable model — as hinted by the recent developments in the electronic-continuum-correction (ECC) force fields. Here we summarize our current experience in attempting to run AMOEBA simulations.
Since the NMRlipids meeting in Prague, we have worked on obtaining the AMOEBA force field parameters, cross-checking them for errors, and setting up the same systems we had for CHARMM-Drude. After the NMRlipids Databank meeting online in December 2021, we finalized the list of simulations to run and distributed them among the NMRlipids contributors interested in running them. Our plan was to simulate pure POPE and POPC bilayers, as well as POPC bilayers at 350 mM, 450 mM, and 1000 mM NaCl and CaCl2 concentrations. We obtained the POPC force field parameters from the authors ahead of publication and created the initial structures for the aforementioned systems.
Our initial understanding was that the AMOEBA simulations would need to be run using either Tinker or the Tinker/OpenMM interface. To this end, we prepared all the files compatible with the Tinker software. As Tinker does not support semi-isotropic pressure coupling, we used the fully anisotropic coupling. However, after running the simulations for a few nanoseconds, we observed a drastic decrease in the area per lipid values, and visual inspection of the trajectories showed that the membranes transitioned into the solid gel-like phase instead of being at the liquid-disordered phase.
To test if the ordering of the membrane originated from the parameters or from the algorithms used, we started to search alternative options to run the AMOEBA membrane simulations. It turned out, contrary to our initial understanding, that AMOEBA simulations can indeed be run with OpenMM. For that, the Tinker force field files need to be converted into OpenMM format, which is achievable via publicly available scripts. The conversion is only possible, unfortunately, for the Allinger-type out-of-plane-bending angles and energies. Allinger type has been used in the force field for the POPE, but POPC bilayer uses Wilson–Decius–Cross (WDC) formulation of the out-of-plane-bend angle, which is not supported by the current version of OpenMM. Therefore, the POPE bilayer simulations using AMOEBA seem possible to run with the OpenMM engine — but the POPC not. We are in contact with the OpenMM developers to see if the WDC formulation can be implemented, but we believe it is not realistic to expect a very rapid fix to this problem.
Interestingly, according to the original publication for DOPC lipids, OpenMM was used for the simulations, although WDC has never been part of OpenMM. It is not yet clear if the authors used an in-house implementation to add WDC functionality to OpenMM, or if there is simply a typo in the publication.
In conclusion, we should be able to run the POPE simulations with the correct pressure coupling and barostat — but for POPC this does not seem possible. We have not yet been able to run AMOEBA simulations of lipid bilayers that would remain in the liquid disordered phase. We are not yet sure if this is due to the parameters we have or the simulation algorithms we can use. We are still planning to use some time to understand this, but we are also open to alternative options to focus on ion-binding affinity. For example, one option would be to run simulations in the NVT ensemble with area per lipid fixed to a reasonable value — yet such a membrane would be under tension, and the area per lipid could not equilibrate upon ion binding, which could potentially affect the results. However, such simulations could still give valuable information on how polarizability affects ion binding (while we wait for fully functional methods for polarizable membrane simulations to become available). In addition, the NVT ensemble might be able to give us useful results on the head group order parameters. In any case, until the WDC angle gets implemented in OpenMM it, unfortunately, seems that we have no way to create the exact equivalents of our CHARMM-Drude simulation setups.
In order to discuss these findings and possible ways to continue this project, we will have a Zoom discussion on 10th of March 2022 at 15:00 CET. If you are interested to join, send email (b.kav at fz-juelich.de) for the Zoom link.
Batuhan Kav
