MM-PBSA setup for peptide/DNA–membrane interaction

[Mostafa Gheydi]

Hello everyone,

I’m working on a system involving a peptide/DNA complex interacting with a multi-component lipid bilayer, and I’m interested in estimating the free energy contributions from each membrane leaflet separately. Specifically, I’d like to analyze how the outer vs. inner leaflet contributes to the free energy barrier associated with penetration, mainly to help predict the peptide’s cell penetration efficiency.

I know that umbrella sampling is the most common and appropriate method for such analysis, but I’d like to explore whether MM-PBSA using gmx_MMPBSA can give insight in this case.

A few questions I hope you can help with:

1. Is it valid to define a whole leaflet (e.g., all outer-leaflet lipids) as the “receptor” group in the index file, and the peptide/DNA complex as the “ligand” (using -cg)? I’d like to treat the leaflet–complex interaction as a type of binding.

2. My system contains a full explicit bilayer (all atoms included). In this case, is it recommended to set memopt = 2 to enable the heterogeneous implicit membrane model in the PB calculation? Would that interfere with the real membrane already present in the system?

3. I saw in the membrane-protein example (here) that the ligand must be provided as an Antechamber .mol2 file. However, since I am not using a small molecule as a ligand (but a peptide–DNA complex), can I ignore this requirement and just define everything using the -cg flag and an index file?

My machine has enough resources to handle larger systems, so I’m open to heavier calculations if needed. 

Any guidance or suggestions from the developers or users would be greatly appreciated!

Thanks in advance,
Mostafa Gheydi
University of Tehran

[mariosergi...@gmail.com]

Sorry for the late response
1. Theoretically, yes. Other users have already done so. But keep in mind that the results are uncertain because these methods can considerably overestimate the energy in these contexts.
2. No. This consideration is an approximation for estimating the dielectric effect of the membrane on an embedded protein without including it in the calculation. In your case, the protein/DNA is not embedded, so it is not necessary.
3. The example is somewhat old and outdated. The option to use mol2 is obsolete. Simply use a calculation like any other, where you define the groups according to your requirements and define the topology, as shown in the gmx_MMPBSA.log file you sent me.

You must take into account the chains in your structure. Since there are several residues without chain IDs and the same number, it can result in overlapping, resulting in errors.

[Mostafa Gheydi]

Dear Mario,
Thank you for your helpful support and for developing such a valuable tool.

I wanted to briefly share that I ran gmx_MMPBSA to evaluate DNA interaction with an anionic membrane leaflet, and obtained a binding free energy (ΔTOTAL) of about 51.55 kcal/mol (101 frames). Given the strong electrostatic character of both the DNA and membrane surface, I believe this result is reasonable; I'd appreciate any thoughts you might have.

During my run, gmx_MMPBSA issued a warning saying no reference file was found, and that the provided structure file (.pdb from -cs) also lacked chain IDs, so it proceeded to assign them automatically. My question is, does gmx_MMPBSA check both the reference file (-cr) and the structure file (-cs) when assigning chain IDs? For example, if I don’t use -cr but supply a .pdb with correct chain IDs in -cs, will those be preserved?

Thanks again for your time and support.

kind regards,

Mostafa Gheydi
University of Tehran

[mariosergi...@gmail.com]

I wanted to briefly share that I ran gmx_MMPBSA to evaluate DNA interaction with an anionic membrane leaflet, and obtained a binding free energy (ΔTOTAL) of about 51.55 kcal/mol (101 frames). Given the strong electrostatic character of both the DNA and membrane surface, I believe this result is reasonable; I'd appreciate any thoughts you might have.

This is tricky. These methods tend to overestimate the energy and have some problems with the electrostatic component in highly charged systems.
Here, you can vary parameters like the internal dielectric constant and use nonlinear PB.
But remember that this is largely experimental; I don't think many case studies have these characteristics.

 

During my run, gmx_MMPBSA issued a warning saying no reference file was found, and that the provided structure file (.pdb from -cs) also lacked chain IDs, so it proceeded to assign them automatically. My question is, does gmx_MMPBSA check both the reference file (-cr) and the structure file (-cs) when assigning chain IDs? For example, if I don’t use -cr but supply a .pdb with correct chain IDs in -cs, will those be preserved?

That's correct. gmx_MMPBSA checks if there are chain IDs in the structure file (-cs); if not, it checks the reference file, if it was defined (-cr). We had to set it this way because, in many cases, without chain IDs, parmed can fail to identify different residues (recently we saw an example of this, the structure included carbohydrates without the ID, and since they had the same name and the same number, parmed overwrote them).
Let's say gmx_MMPBSA is reasonably good at identifying chains, but sometimes the structures we receive are extremely inconsistent, so there's no way to resolve the issues.

[Mostafa Gheydi]

Thanks, Mario, that clarifies a lot.

One follow-up: I understand that the non-linear PB solver isn’t compatible with solvopt 1-4 and bcopt = 10 (periodic).
But I also saw in this Amber list post (http://archive.ambermd.org/202006/0099.html) that periodic boundaries are recommended for membrane systems.

So I’m not sure how to reconcile the two; if I switch to non-linear PB, what boundary-condition and solver settings should I use for a membrane?
Any practical advice would be greatly appreciated.

Many thanks again!

Best regards,
Mostafa

[mariosergi...@gmail.com]

It's certainly a problem because the membrane typically extends from side to side of the box, so not considering the BCs can give inconsistent results. Since in your case, the membrane is your receptor, you would then be left with the option of testing whether NonLinear PB gives you acceptable results compared to the alternative, which would be to play with the internal dielectric constant using Linear PB.
In either case, you have a large implicit uncertainty because these methods have not been parameterized to consider the membrane as a receptor.

[Mostafa Gheydi]

Dear Mario,
Thank you again for your thoughtful and clear explanation, much appreciated!

To better understand how parameters affect binding energy in membrane systems, I ran several tests by varying the non-polar solvation model (inp) and the internal dielectric constant (indi). I used the same general input across all runs:

&general sys_name="DNA-Memb", startframe=1, endframe=2001, interval=20, PBRadii=7, temperature=300, /
&pb memopt=2, indi=..., mctrdz=80.325, mthick=42.85, radiopt=0, istrng=0.150, fillratio=1.25, inp=..., sasopt=0, solvopt=2, ipb=1, bcopt=10, nfocus=1, linit=10000, eneopt=1, cutfd=7.0, cutnb=10.0, maxarcdot=5000, npbverb=1, /

Here are the results I obtained:

1. inp = 2, indi = 4
   ΔTOTAL = +51.55 kcal/mol
   ΔGGAS = -39.43
   ΔGSOLV = +90.98
   ΔEDISPER = +184.83

2. inp = 1, indi = 20
   ΔTOTAL = -109.11 kcal/mol
   ΔGGAS = -88.90
   ΔGSOLV = -20.22
   ΔEDISPER = 0.00

3. inp = 1, indi = 8
   ΔTOTAL = -85.87 kcal/mol
   ΔGGAS = -65.65
   ΔGSOLV = -20.22
   ΔEDISPER = 0.00

My system involves a DNA near the membrane surface (outer leaflet).

I’d really appreciate your input on a couple of things:

1. How should I choose inp in this kind of system? It seems inp = 1 is used for shallow or interfacial pockets (e.g., water/lipid headgroup interface), while inp = 2 might be more suited for deeply buried pockets. Since the DNA interacts with the membrane surface, inp = 1 feels more appropriate.
   However, I came across this thread, where it was mentioned that the choice of inp can be force-field dependent. I’m using CHARMM36, but I wasn’t sure how this should influence my decision; any clarification would help.

2. How should I properly choose indi (internal dielectric)? I’ve seen a wide range of values from 1 to 20 used in the literature, but it's still unclear how to select a physically reasonable one. For a highly charged system like DNA interacting with an anionic membrane, what would be a sensible value (or range)?

Thanks again for your time and guidance.

Best regards,

Mostafa Gheydi
University of Tehran

[Mostafa Gheydi]

Dear Mario,

I understand you may be busy, but I would greatly appreciate your insight whenever you have a moment, as I really need it to move forward with my work.

Thank you very much for your time and consideration.

[marioe911116]

Please check this paper of ours where we discussed the effect of the internal dielectric on the calculations (there are also other works on this): 

https://pubs.acs.org/doi/10.1021/acs.jpcb.2c07079

given that you are working with a highly charged system like DNA, indi of ~10 might be appropriate

[Mostafa Gheydi]

Thank you again for all your help. I genuinely appreciate the time you took to carefully address my questions.
> I’m still unsure about one point from my earlier message, so I’d like to clarify it briefly.

   How should I choose inp in this kind of system? It seems inp = 1 is used for shallow or interfacial pockets (e.g., water/lipid headgroup interface), while         inp = 2 might be more suited for deeply buried pockets. Since the DNA interacts with the membrane surface, inp = 1 feels more appropriate.
   However, I came across this thread, where it was mentioned that the choice of inp can be force-field dependent. I’m using CHARMM36, but I wasn’t sure     how this should influence my decision; any clarification would help.
 

My main concern is how to select the inp.
In my tests (DNA on the outer leaflet, CHARMM36 FF) inp = 1 gave
ΔTOTAL ≈ –86 kcal/mol (indi = 8) while inp = 2 gave +52 kcal/mol (indi = 4).

Could you share your opinion on how I should choose inp in this kind of system?

Thank you once again for your generous and helpful replies.

Kind regards,

Mostafa Gheydi
University of Tehran