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
You can run AMOEBA simulations using Tinker9 on GPUs directly. Same Tinker (CPU) input files. It has the (best) Langevin piston pressure control.
ReplyDeletehttps://github.com/TinkerTools/tinker9
https://biomol.bme.utexas.edu/tinkergpu/index.php?title=Tinkergpu:Tinker-tut#How_to_run_Tinker9_using_GPU_(2021)
Thank you for the contribution.
DeletePlease note that we would strongly encourage that "all communications would be written such that they are assignable to a person and affiliation", see the post On credits.
Thanks for the comment.
DeleteIndeed, we have been able to run AMOEBA using Tinker9 on GPUs directly with a decent performance. What is missing thought is the semi-isotropic pressure coupling which would be optimal for lipid bilayers. Nevertheless, I think that the fully anisotropic is ok as long as the system does not explode.
The problem in our simulations with Tinker9 is that they go to gel-like phase, while AMOEBA lipid bilayer simulations with OpenMM in correct phase has been reported in the literature in the publications cited in the original post. Question is now if this is due to the parameters we have or some algorithms used. Personally I would guess the parameters, but it would be really good to exclude the possibility of algorithm.
We are not associated with the lipid parameters published by Li. It is in AMOEBA style but we do not distribute them with AMOEBA parameter files in anyway. I won't be surprised they need to be worked on/improved (all parameters straightly came out of QM of small fragments). We only did some limited tests with them so far. There are issues but haven't seen gel. We plan to work on an AMOEBA lipid library with Blake (see his comments below). -Pengyu Ren
DeleteBatuhan,
ReplyDeleteI've been working on reproducing the results from the Li parameters as part of my sabbatical project for the past 18 months. I've been running DMPC simulations in both Tinker8-OpenMM and Tinker9, which has the langevin barostat included. This may do a better job of treating the system, as I don't see any transition to the gel phase upon visual inspection (am calculating the area per lipid this week). I would be interested in participating in the Zoom discussion when it happens, as I have a vested interest in seeing a polarizable lipid force field be successful (we want to use it for the protein systems we are investigating). You can get in touch with me at blake.mertz@mail.wvu.edu. Thanks for taking on this project, it is an ambitious one!
Blake
Thanks for the information!
DeleteDid you use anisotropic pressure scaling? I think that I had problems to use it with other than Berendsen barostat in Tinker9.
Hello,
DeleteThank you for the information, Blake.
I tried the POPE system with the Langevin barostat. I have tried different time steps from 0.1 to 2 fs but all crashed with induced dipoles not converging. Did you have similar issues before?
You're more than welcome. I'll have to reach out the Tinker9 developers about the ability to run anisotropic or semi-isotropic (the documentation is something I'm trying to contribute to as well, as it's pretty sparse). I need to go back and look at my data -- I've been benchmarking the systems with both Nose-Hoover and Langevin barostats (NH requires a 1 fs timestep, which is unfortunate). The DMPC system I mentioned that is stable was run with the NH, so I hope I didn't speak too soon. I do question the published results -- the DMPC system I've simulated severely underestimates the 2H NMR order parameters for most of the acyl tail, from C2 to C11, compared to experimental results.
ReplyDeleteRegarding the semi-isotropic pressure coupling, I already contacted the developers
Deletehttps://github.com/TinkerTools/tinker9/issues/181
For the POPE lipids, the C2 order parameter for the SN-1 chain has an error in the original publication. However, the other tail group order parameters, as I remember, were reproducible.
Yes it can run anisotropic (ANISO-PRESSURE will allow anisotropic relaxation). I suggest asking questions on https://github.com/TinkerTools/tinker9 (issues). We will all see it.
ReplyDeleteIf dipole convergence issue occurs right away, the structure or prm files are probably problematic. If it happens after a while, check the biomol utexas page for options. The latest Tinker9 should be very stable regarding the dipole scg.
ReplyDeleteHello,
ReplyDeleteI think at the moment it would be useful for us to try to get our simulations running on OpenMM as well. While looking into this, I have found the Li parameters already available in OpenMM format (https://github.com/Inniag/openmm-scripts-amoeba/). I unfortunately do not have a working experience with OpenMM. Would anyone be interested in creating the simulation input files with OpenMM?