SMD Pulling of DNA Toward Asymmetric Bilayer

User discussions
1

GROMACS version: 2024.4
GROMACS modification: No

Hi everyone,

I’m simulating the interaction between a DNA duplex and an asymmetric lipid bilayer, and I’m running into some unexpected behavior during steered molecular dynamics (SMD). The goal is to pull the DNA through the bilayer along the Z-axis, but so far, it only reaches the interface — and the pulling reaction coordinate starts from a negative value. I’d appreciate any feedback on the setup and results.

System Overview

  • A DNA helix was placed ~5 nm above the bilayer COM, roughly parallel to the membrane surface.
  • Built and visualized in VMD.
  • Lipid bilayer is asymmetric and simulated with CHARMM36 force field.

Equilibration Protocol

I used a multi-step equilibration with gradually released restraints:

Step Duration Restraints (DNA/Lipid/Dihedral) (kJ/mol/nm²)
EM – 1000 / 1000 / 1000
NVT 1 ns 1000 / 1000 / 1000
NPT1 1 ns 1000 / 1000 / 1000
NPT2 1 ns 500 / 400 / 400
NPT3 1 ns 200 / 200 / 100
NPT4 5 ns No restraints

  • Thermostat groups:

    tc-grps = DNA Membrane Water_and_ions

Pulling Setup (SMD)

DNA was pulled toward the bilayer along Z (membrane normal):

pull                    = yes
pull_ncoords            = 1
pull_ngroups            = 2
pull_group1_name        = Membrane
pull_group2_name        = DNA
pull_coord1_geometry    = direction-periodic
pull_coord1_dim         = N N Y
pull_coord1_vec         = 0 0 -1     ; pulling down
pull_coord1_rate        = 0.001      ; 1 nm/ns
pull_coord1_k           = 2000
pull_coord1_start       = yes
pull-group1-pbcatom        = 11114
pull-group2-pbcatom        = 43856
pull-pbc-ref-prev-step-com = yes
pull-nstxout               = 50
pull-nstfout               = 50

  • Pulling time: 10 ns

  • Used semiisotropic pressure coupling, with Z compressibility = 0:

    pcoupltype              = semiisotropic
ref_p                   = 1.0 1.0
compressibility         = 4.5e-5 0.0
refcoord_scaling        = com

Key Observations

  1. DNA Does Not Penetrate Bilayer
  • After 10 ns of pulling, DNA reaches the headgroup interface, but never crosses the bilayer.
  • It flattens against the surface instead.
Frame View Description
First Side DNA placed clearly above bilayer (~5 nm)
Final Side DNA stops at bilayer interface
Final Top DNA is no longer centered in XY
  1. Pulling Force (pullf.xvg)
  • Force reaches ~400–500 kJ/mol/nm, indicating strong resistance.
  • Yet, DNA does not insert — likely pointing to a significant energy barrier.
  1. Unexpected Reaction Coordinate (pullx.xvg)

  • Pulling distance starts from a negative value, even though DNA was above the membrane.

  • Could this be:

    • A misinterpretation of the pull group COMs?
    • A wrapping or periodic image issue?
    • Axis flipped?

My Questions

  1. Why does the pulling start from a negative value if DNA starts above the membrane?
  2. Is direction-periodic the best choice for this SMD geometry?
  3. Should I be concerned about XY drift? If so, how should I fix it?
  4. Could the DNA’s failure to insert be physical (electrostatics, packing, etc.), or due to setup artifacts?
  5. Why do some lipids move so oddly? Should I apply restraints on lipids even during production MD, or might this be an artifact of removing restraints quickly?

Thanks so much in advance for any insights or suggestions! I’m happy to provide more details or files if needed.

Grace: pullx.xvg_001
Grace: pullx.xvg_0011328×736 10.6 KB

Grace: pullf.xvg_001
Grace: pullf.xvg_0011328×736 12.8 KB

First_frame
First_frame669×736 138 KB

Final_frame_top
Final_frame_top669×736 136 KB

Final_frame
Final_frame669×736 110 KB

2

I think that you are confusing direction of the reaction coordinate with the pulling direction. Use 0 0 1 as direction and set a negative pull rate if you want to pull down. Furthermore, you should use geometry=cylinder where a membrane is involved.

But if your molecule flattens, I think there that the reaction coordinate is simply not good enough. If you are lucky pulling slower might help.

3

Thank you very much for the clear explanation and the crucial corrections!

I now understand the confusion between the direction vector and the pulling rate, and I’ll immediately adjust to using pull_coord1_vec = 0 0 1 with a negative pull rate. I will also try the geometry = cylinder option, as you suggested.

May I ask two quick follow-up questions?

  1. When using geometry = cylinder, should I ensure that the Z-axis box length is at least twice the length of the reaction coordinate? I’m concerned about box size and water molecules.
  2. For the next step, which is umbrella sampling, can I use a shorter equilibration phase (like a single NPT step) instead of the multi-step protocol I used for SMD? Would this still be reliable for the umbrella windows?
4

I’ve also noticed that water molecules reappear in the bilayer core during equilibration and production steps, even after I manually removed them following solvation. Could you suggest reliable methods to prevent water from re-entering the membrane core?

Unfortunately, the membrane rotates during simulation. It starts in the xy plane (z-axis normal) and ends in the xz plane (y-axis normal). Should I be concerned about the membrane rotating from z-normal to y-normal? Would this affect the accuracy of PMF or SMD analysis?

5

When using geometry = cylinder, should I ensure that the Z-axis box length is at least twice the length of the reaction coordinate? I’m concerned about box size and water molecules.

No, there is no such requirement.

For the next step, which is umbrella sampling, can I use a shorter equilibration phase (like a single NPT step) instead of the multi-step protocol I used for SMD? Would this still be reliable for the umbrella windows?

If you can get a reasonable trajectory with SMD, which is doubtful, you can extract starting frames from that and you will need no equilibration at all.

6

Thanks a lot for your response.

7

Dear @obZehn and @alevilla

I truly appreciate the valuable insights shared in this discussion.

I’ve noticed from your previous contributions that you both have deep expertise in membrane simulations and have helped clarify many nuanced issues related to bilayer systems, pulling setups, and system preparation.

Given your experience, I would be very grateful to hear your thoughts on the challenges I described above.

Thank you in advance for your time and consideration!

8

Hi @mostafagh98

I do not have much to add with respect to what @hess already pointed out about the SMD dragging of the DNA strand, in terms of system set-up. I might add a few points here and there that might help you

NPT4 5 ns No restraints

In general, 5ns of equilibration or a lipid bilayer seems extremely short. This won’t change much of your outcome, but still a point to raise.

I’ve also noticed that water molecules reappear in the bilayer core during equilibration and production steps, even after I manually removed them following solvation. Could you suggest reliable methods to prevent water from re-entering the membrane core?
Unfortunately, the membrane rotates during simulation. It starts in the xy plane (z-axis normal) and ends in the xz plane (y-axis normal). Should I be concerned about the membrane rotating from z-normal to y-normal? Would this affect the accuracy of PMF or SMD analysis?

These are a bit more worrying. Do you expect water molecules to be there? I do not know the nature of the membrane. By the images you posted, it seems like it doesn’t span the whole box but it seems more a nano disk of some sort. In that case then the curvature at the borders might let some water in, as the chains are much less packed. If this is a lipid disk then there are no lipids spanning the whole box and the slab can rotate freely, as all free structures diffusing in water can do, and this can explain why you see the slab reorienting.

Regarding the DNA flattening on the bilayer, I guess that makes perfectly sense. I am not competent in DNA simulations or biology, but can a DNA strand physiologically partition through a bilayer? I think DNA is (heavily?) charged and will prefer by far to stick on the membrane surface, interacting with the lipid heads and water, rather than go through a lipophilic slab. On top of that, you are trying to pull a huge strand that has a ton of degrees of freedom. What is the expected behavior that you want to simulate? Is is that of the strand along its long axis penetrating the bilayer? In this set up you are applying just a force to the center of mass of a large molecule, you have zero control on the rest, so the first DNA atoms are going to stick to the surface of the bilayer and then the applied force will push the rest down and eventually stick everything to the surface, making you lose whatever configuration you were aiming for. Moreover, the DNA strand is as large as the lipid patch, so I would expect very large changes also in the lipid slab. I don’t know if with such a molecule partitioned inside the slab will be even able to hold together or just fully reorient to minimize the interaction of the strand with the lipid tails. As said in the latter, maybe if you pull very very slowly you might end up in a quasi-equilibrium regime where you gradually pull through the lipids and from there you can take the windows you need for sampling, but I honestly have a hard time seeing it given the nature of the molecules. I think you are very much limited by you collective variable (a simple distance) and might need much more refined CVs that actually also accelerate slow degrees of freedom in the membrane, as well as probably enforcing with some kind of restraints the DNA and/or the slab. But again, is this expected to happen? Is this a physiological process/an experimental set up of some kind?

9

Thank you very much for your thoughtful and detailed feedback.

I definitely do not expect water molecules to enter the bilayer core. I tried removing them manually after solvation, but they reappear during equilibration. I’ve read about possible strategies like:

Flat-bottomed position restraints on water molecules,
Or using SETTLE constraints in the topology file (e.g., for water hydrogens),

But I’m unsure whether either method would introduce artifacts or negatively affect the membrane dynamics. Do you have experience with such techniques or any recommendations for keeping water out of the bilayer?
Do you think I may have contributed to this by setting compressibility = 0.0 in the z-direction? I found it based on a suggestion here.

My membrane is not a nanodisc. It’s composed of:

Upper leaflet: 128 POPC, 27 Cholesterol, 15 GM3
Lower leaflet: 110 POPC, 27 Cholesterol, 33 POPS
(170 lipids per leaflet, total 340)

I made sure the area per lipid (APL) is balanced between the two leaflets before and after equilibration and that the bilayer plane spans at least 1.5 nm more than the DNA helical axis on each side. Do you think I should still increase the membrane size?

You’re absolutely right to ask — no, DNA does not physiologically partition through a bilayer on its own. In experiments, we use naked DNA as a negative control.

My actual goal is to later simulate DNA + peptide complexes, where different numbers of peptides are bound to the DNA. The current naked DNA simulation serves as a baseline for PMF comparison, to show that the DNA alone has a higher free energy barrier compared to the complexes.

I’ve seen similar PMF studies done using symmetric bilayers. But in those cases, the neutral charge of the membrane tends to favor DNA over complexes. That’s why I’m using an asymmetric bilayer to reflect more realistic surface charge conditions.

I’m intending to pull the DNA along its long (helical) axis.

Regarding your point about applying force only to the center of mass (COM) of the DNA, I understand that this provides little control over the molecule’s orientation. I’ve seen COM-based pulling used in several studies, but I’m open to better alternatives.
Would you recommend a more localized pulling group or a more refined collective variable?

The membrane rotation (from z-normal to y-normal) is concerning.

I’m considering adding the following lines in the .mdp to suppress membrane rotation:

nstcomm       = 100
comm_mode     = angular
comm_grps     = MEMB

Would this be appropriate for a bilayer system, or would it introduce unintended side effects?

Thanks again for all your valuable insights.