AMOEBA POPE simulation about openmm-scripts-amoeba

[Inniag]

Hi Batukav, the link you posted requires a Google Doc account, so I can't look at it right now. However, my hunch would be that your coordinate files contain more than 10000 atoms in the POPE residues.

I had a similar problem in the past: bonds are specified as pairs of indices in PDB files, with only 4 digits reserved for each index. So with more than 10000 atoms, you overflow the index field and OpenMM can no longer understand which atoms are bonded.

For standard residues like water and amino acids, OpenMM has (heuristic?) rules for inferring bonds, but for "non-standard" lipids this is not implemented afaik.

[Inniag]

Take a look at my simulation script: https://github.com/Inniag/openmm-scripts-amoeba/blob/master/scripts/simulation/simulate_amoeba.py

In line 32 I define a function to workaround this problem. Perhaps it helps you as well.

[batukav]

Hi @Inniag , thank you for the quick reply.

Sorry for the file location, I attach a zipped version here. I only have 72 lipids, which makes only 9000 atoms for the system.
pope.zip

I am already using the simulate_amoeba script to go through the PDB file and I am already using the pdb_file_nonstandard_bonds function.

[Inniag]

I am on mobile at the moment, can't really check the zip ATM.

My next guess would be that your system is not properly equilibrated yet. The guess_bonds makes heuristic assumptions on atom proximity I think and CHARMM GUI generated bilayers might not place all lipids in energetically sensible configurations (e.g. overlapping vdW spheres).

I solved that by doing a short (10ns) atomistic simulation using CHARMM36m. This is likely need anyway, as AMOEBA is very sensitive to non-physical starting conditions and the electronic degrees of freedom will likely diverge of you don't equilibrated atomistically first.

[batukav]

Thank you for the suggestion! I will try running the scripts after equilibrating the bilayer. I'll update here on how it goes.

Update:

Indeed, it was the minimisation. After running the initial structure using the CHARMM36 for a few ns, simulate_amoeba worked. Thank you!

[batukav]

Hello,

I have managed to run some simulations of pure POPE bilayers using the parameters and scripts provided here. Thank you again for providing these tools.

In my simulations, I tried three different options: 1- default settings, 2-using Langevin multistep integrator, and 3- using Langevin multistep integrator with PME cutoff at 1.2 nm (as opposed to the default 0.8 nm). The results from all these simulations are very similar to each other.

One thing I realised is that the area per lipid values I obtain (66.3 +- 1.8 A^2) are quite larger than the originally published values (61.1 A^2) [https://doi.org/10.1080/00268976.2018.1436201] .

I would like to ask if you have calculated the area per lipid values (without the protein) from your simulations? I could not find this information in the published JACS paper.

[Inniag]

Glad you got your simulations to work.

Your observations regarding integrator and PME are in line with my experiences. I would generally use the multi-step integrator for performance reasons alone and set the PME cutoff to whatever value was used in force field parameterization (if possible - sometimes lipids and proteins were parameterized with different settings).

Regarding area per lipid: I never measured that because I was not looking into lipid-protein interactions or the like. If you have already ruled out electrostatics as a cause for the difference here, my next step would be to check thermostat and - in particular - barostat parameters. Are they consistent between the source paper and your simulations? Are they using the same algorithm (I don't recall anymore and its paywalled to me now). Barostats work by adjusting box size, so it's quite likely they impact something like area per lipid.

The other thing to check: Is your system properly equilibrated? Does area per lipid no longer change over time (for tens of ns at least would be good) and is area per lipid independent of the initialisation conditions? I never checked, but I could imagine that polarisable force fields have a longer relaxation time due to the increased "stickyness" caused by induced dipole interactions.

[batukav]

Thank you for the explanations.

Regarding the area per lipid: I already ruled out the thermostat as I have tried both Andersen and Langevin and both gave the similar results. The barostat is the same barostat in the original publication (membrane Monte Carlo barostat with the pressure along xy-plane is independent from the z-component). Therefore, I am expecting to be able to get the same area per lipid values published in the original paper.

I would expect the same in terms of the relaxation time. This is something we have also seen with the CHARMM-Drude force field that it takes quite longer for the system to relax. At the moment, within the first 100 ns of the simulations the area per lipid fluctuates around its average value without any trend of going down. This could be due to the initial structure and/or slow relaxation dynamics.

I will now first run a long equilibration with CHARMM36, then I will start the simulations with the correct area per lipid as the initial condition. Maybe it will make a difference; I will edit this comment with the results.

Thanks for the suggestions.