SMD Pulling of DNA Toward Asymmetric Bilayer

[mostafagh98]

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

2. 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.

3. 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.

[hess]

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.

[mostafagh98]

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?

[mostafagh98]

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?

[hess]

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.

[mostafagh98]

Thanks a lot for your response.

[mostafagh98]

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!

[obZehn]

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?

[mostafagh98]

Dear @obZehn, thank you very much for your thoughtful and detailed feedback.

@obZehn

Do you expect water molecules to be there?

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,
But I’m unsure whether this 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.

@obZehn

I do not know the nature of the membrane.

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?

@obZehn

can a DNA strand physiologically partition through a bilayer?

@obZehn

is this expected to happen? Is this a physiological process/an experimental set up of some kind?

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.

@obZehn

What is the expected behavior that you want to simulate?

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.

[obZehn]

But I’m unsure whether this 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?

I would avoid restraining water in general. A flat bottomed potential with a nanodisk-like system will not work as the potential will keep the water away from a full slab of space, also where there are no lipids.

Do you think I may have contributed to this by setting compressibility = 0.0 in the z-direction?

Honestly I do not think so. Maybe if you had a normal “infinite” bilayer and you decoupled the xy plane, starting with a bad lipid configuration you might end up with a ruptured/too much compressed bilayer. However, along z I can’t see problems, especially considering that you have a nanodisk and not a full bilayer; the system should be able to relax properly. Setting off in the z-direction is useful to decouple the pressure if you need to use SMD to pull groups farther away than half the box size, otherwise you will have ill-defined PBCs distances and the system will crash. You can read it here. If you do not need to pull groups that far away then you do not need the decoupling (but given how big is the DNA strand you most likely will need it).

The water molecules going in the bilayer core are to me or i) something expected because in some way this is a poorly hydrophilic/somewhat solvatable bilayer or ii) symptoms or a bad starting structure. I honestly would tend for the second one, which is also the easier to verify (re-build the box and try anew).

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?

If you want to simulate the partition of a DNA in a bilayer - not on a nanodisk - then honestly I don’t see a reason for not having a bilayer! :) A nanodisk will just over complicate your life.

Would you recommend a more localized pulling group or a more refined collective variable?

Depends. What do you want to simulate? If it is the DNA going through along its long axis then I feel like this is never going to happen without further refining the collective variables you are using. Already the COM itself is super degenerate; there are infinite directions in which the strand can be oriented which share the same COM position. So yes, I would 100% thing about something more refined, but here I can’t help as this goes beyond my expertise.

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

I am not completely sure, but my take on the nstcomm and following key words is that this a collective motion removal to avoid artifact like the flying ice cube effect or similar (which I guess are not really present anymore anyway), and in general avoid the transfer of energy from some degrees of freedom to those of the box. As such, I don’t think this will fix your problem, but actually make it worse, as you will subtract angular and linear velocities to a subgroup of the system, but not to the rest. Here I might be mistaken, but I’d rather fall back on more common position restraints on the lipids that you able to control effectively much more in depth.

Anyway I think that what you are trying to do is very very hard, as the DNA strand has a ton degrees of freedom. So my best suggestion would be: has anyone done this type of PMF estimates before in literature? Can you follow roughly their approach/methods?

[mostafagh98]

Hi @obZehn,
Thanks a lot for your reply and all the helpful suggestions.

@obZehn

So my best suggestion would be: has anyone done this type of PMF estimates before in literature? Can you follow roughly their approach/methods?

Regarding your question about similar studies: yes, I came across a few papers that inspired my approach, including:

• https://doi.org/10.1021/acs.jpcb.9b05110
• https://doi.org/10.1021/acs.jpcb.5b01256
• https://doi.org/10.1039/C6RA05607E

Originally, I meant to work with a full planar bilayer, but I accidentally made the box larger than the membrane’s XY size, which caused the system to behave more like a nanodisk. I’ve fixed that now.

I also updated my pulling setup based on @hess’s recommendations. Here’s what I’m currently using:

; Pull options: constant-velocity SMD along Z-axis
pull                    = yes
pull_ncoords            = 1
pull_ngroups            = 2
pull_group1_name        = Membrane
pull_group2_name        = DNA
pull_coord1_type        = umbrella
pull_coord1_geometry    = direction-periodic
pull_coord1_dim         = N N Y
pull_coord1_groups      = 1 2
pull_coord1_start       = yes
pull_coord1_vec         = 0 0 1
pull_coord1_rate        = -0.0001
pull_coord1_k           = 2000

Interestingly, with this setup the DNA actually inserted into the membrane and even crossed it, but I noticed that the pulling force went negative, which I didn’t expect and found a bit confusing.

A few questions I was hoping to get your thoughts on:

1. Do negative pulling forces impact the PMF in umbrella sampling? Should I be worried about that?
2. I was thinking of using ~0.2 nm spacing between windows along the Z-axis. In your experience, would 0.1 nm be safer to ensure better overlap for WHAM?
3. When calculating distances using gmx distance or extracting configurations for windows, should I use the raw pull.xtc or a PBC-corrected trajectory?
4. If I notice poor overlap during analysis, would you recommend using tighter window spacing from the beginning, or is it okay to just go back and add windows where needed?
5. And finally, Is it possible to switch geometries (e.g., from direction-periodic to direction) in umbrella sampling (USMD) without introducing artifacts?

Thanks again for all your help.

[hess]

The pulling force has to be negative to pull the molecule down against the upward force resisting the pulling.

I would still suggest using geometry=cylinder, with a sufficiently large radius.

My guess is that the main issue is going to be properly sampling rearrangements in the bilayer and the DNA. This might or might not work with umbrella sampling. In your simulation above the work is about 2000 kJ/mol, so you are extremely far from equilibrium.

[mostafagh98]

Thanks again for taking the time to help.
About the cylinder geometry, I totally understand your recommendation. I actually avoided using it because, according to the GROMACS manual, the cylinder radius needs to be smaller than half the box size. Since my DNA is about 7 nm long, I’d need at least a 15–16 nm box in X and Y, and to maintain a truly infinite membrane without edge effects, I’d have to double the membrane area, which seemed too demanding in terms of computational cost. Am I thinking about this the right way? Or do you think there’s a smarter way to manage this without dramatically increasing the system size?
Also, just to clarify, I tried to pull the DNA into the membrane by one of its termini, but it starts out lying parallel to the membrane.

[hess]

But the pore that is created has a diameter much smaller than half the box size, or?

[obZehn]

I think @hess addressed the most important of the points, which is that it is going to be extremely challenging to converge any kind of free energy given how slow many of these degrees of freedom will be.

Given your points:
ii) No idea, as this is system dependent. 0.1nm is a common choice, if the papers you link go for one I would start from there;
iii) I would use the pullx.xvg file that is generated during the pulling, as that is the distance between the COM of the strand and that of the bilayer;
iv) I think I do not understand the question. Formally there is no difference, in the sense that if I want more sampling in a region you go back to the SMD trajectory, pull out new configurations that cover the CV space you want to be covered, and you add those to the sampling/analysis;
v) I wouldn’t expect many differences if you only re-couple along z, but this is a bit up to you; I would probably re-couple the pressure and just run equilibration of the windows.

[mostafagh98]

Hi @hess,

Sorry for the late reply, and thank you again for your valuable suggestions. I tried switching to geometry = cylinder as you recommended (with pull-cylinder-r = 5), and it really improved the behavior of the simulation. I’ve attached a few snapshots to show the results.

I also increased pull_coord1_k from 2000 to 5000 kJ/mol/nm². I noticed some related studies use high force constants. Is there a specific reason for choosing such large values?

Additionally, is there any relationship between pull_coord1_k and pull_coord1_rate? For example, now that I’m using a stiffer spring (k = 5000), should I also consider adjusting the pulling rate? Right now I’m using -0.0001 nm/ps (0.1 nm/ns).

[mostafagh98]

Dear @obZehn
Thank you again for your generous and helpful replies. I truly appreciate your time and the clarity you’ve provided on all my questions; it’s been incredibly valuable.

[hess]

You are generating very large forces, so you will need a high force constant. The force constant and rate are not directly related. The rate needs to be sufficiently small to stay close to equilibrium. For you case I think this will be several order of magnitude slower than the rate you are using now and likely practically unreachable.

If you want a real estimate: repeat the same SMD procedure once or twice. Then integrate the force over the distance to get the work. To get reasonable results the variation in this work should be on the order of a few kT. I expect this to be tens to hundreds of kT in your case.

[mostafagh98]

Appreciate your continued support. I’m learning a lot from your advice.