MD run explodes

[Lee Xiong An]

Hi Sam,

I was doing an MD run, but it has exploded. For speed purposes, i rigidified most of the RNA, but kept residues 51-69 free and included default MD parameters for these residues and a few others more. However, the MD run was exploding. Is there a way to reduce the time step ? I started from a crystal structure, so i suppose it is a reasonably equilibrated structure. ]

Xiong An

Protocol
firstStage 125
lastStage 125


reportingInterval .005
numReportingIntervals 100
RNA A 1 GUUCUUAUCAAGAGAAGCAGAGGGACUGGCCCGACGAAGCUUCAGCAACCGGUGUAAUGGCGAAAGCCAUGACCAAGGUGCUAAAUCCAGCAAGCUCGAACAGCUUGGAAGAUAAGAAC

mobilizer Rigid A 1 50
mobilizer Rigid A 70 119
mobilizer FREE A 51 69

baseInteractionScaleFactor 200
baseInteraction A 24 WatsonCrick A 30 WatsonCrick Cis
baseInteraction A 23 WatsonCrick A 31 WatsonCrick Cis
baseInteraction A 22 WatsonCrick A 32 WatsonCrick Cis
...
... plus the other 50 base pairs
...
baseInteraction A 4 WatsonCrick A 116 WatsonCrick Cis
baseInteraction A 3 WatsonCrick A 117 WatsonCrick Cis
baseInteraction A 2 WatsonCrick A 118 WatsonCrick Cis

setDefaultMDParameters
includeResidues A 41 79

[Samuel Flores]

I would use "Default" rather than "Free". Free confers 6 degrees of freedom (DOFs) for each bond. Default makes most bonds either Torsion (1 DOF) or Rigid (0 DOF). Rigid is for groups like phenols, or the peptide bond, that MMB prefers to treat as though they had no internal DOFs. Anyway, Free costs a lot more.

I would prefer you use the capitalization in the tutorial (e.g. Free rather than FREE). Some of the parameters are case-insensitive, but I have not gotten around to making ALL of them case insensitive, so alternate capitalizations are not supported.

I would not refer to this as an MD run in the sense of GROMACS or NAMD, because so much of the molecule is rigid, and for other reasons. I think of it as a prototype MD, you might still want to do formal MD with explicit solvent using GROMACS to confirm your results. I think with MMB you can get a fast idea of what your system is doing, and generate lots of cool hypotheses in a short time.

Try turning Free into Default and see if it is still exploding. I can't check your work because you did not provide me with a PDB ID.

One more suggestion: instead of instantiating the RNA explicitly you can just do:

loadSequencesFromPdb

The advantage of this is that it can handle insertion codes, deletions, etc. There are still a few bizarre atom name and residue number practices that some PDB files have which will make MMB puke, but the latter has gotten more robust in the trunk version (not yet available as binaries). Let me know if you have any further trouble. I love to hear what my users are up to.

Sam

PS I see you are at A-Star. Are you a PhD student? What is your thesis topic?

[Lee Xiong An]

Hi sam,

I am not running MD, but i have been trying to run some homology modelling, but the insertions and base pairs are not properly preserved nor is their structure 'reasonable' after quick visual inspection. I ran this 'MD' step just as a debugging step for MMB and i found the forcefield or timestep problematic, causing the run to explode.

It seems to explode when i rigidfy two portions of template, 1-50, 70-119, leaving 51-69 as default with 1 DOF. Kept MD parameters for residues 41-79.
Also, rigidfying 10-110 and leaving both ends as default and keeping MD pararmeters also had the same explosions.

The pdb id is 4kqy

Xiong An

I am a research asst, but i believe i am on an equal standing as a PHD student at least currently . My current project is developing RNA homology modelling protocol for insertions and deletions, and effectively the 'thesis' is finding antiobitics for a TPP riboswitch with lots of insertions (Pseudomas A) .

[Samuel Flores]

Hi Xiong An,

I made a few modifications to your script. Assuming you are using the latest build, it will be easier to work this way.

First, you can just issue "mobilizer Rigid" to rigidify everything, and then set 51-69 back to being flexible. You can use nucleicAcidDuplex to make base pairs in a double helix, for any stretch that is all watson-crick base pairs with no bulges. You should constrain the two rigid segments to each other to keep them from drifting apart.

You can use "physicsRadius" to just include all residues within some distance of the flexible pieces.

I append my script. it does not explode at all. However it is not really at the point where I would say it is informative. I would add a few more base pair forces. the flexibility zone is not quite right. You may wish to add solvent and ions. Depending on what you want to do , it could be a one semester or one year project. What is it you want to do, exactly?

Sam


firstStage 2
lastStage 2
loadSequencesFromPdb 4KQY.pdb

mobilizer Rigid
mobilizer Default A 51 69
constraint A FirstResidue Weld A LastResidue


baseInteractionScaleFactor 200
nucleicAcidDuplex A 56 60 A 69 65



setDefaultMDParameters
physicsRadius 1.0




reportingInterval .05
numReportingIntervals 100

[Lee Xiong An]

Hi sam,

This MD run was not meant to be informative, but was just a debugging step just to ensure that any molecule fed into MMB retains proper sterochemistry. For all my scripts, i found this was not the case. I suspect i am missing out on a few crucial commands or failed to understand how certain things work.

Areas where my script breaks down
1. Many areas seem to unphysically relax if it is not held rigid, and it looks unreasonable based on my eyes

less serious problem
2. . base pairs command seems to fail and residues fail to basepair , though sometimes it works

What i am really trying to do is benchmark your homology approach against mine. However, i have failed to model and deletions or insertions using MMB properly, mainly due to problem 1, the insetion or deletions unphysically relax, even in the current script.

Xiong An

[Samuel Flores]

Sorry I seem to be missing these messages.

1. You can expect unphysical results if parts of your molecule are flexible and physics is not turned on, or is not turned on for enough of the nearby residues. For homology modeling you have to be very careful to keep the MD forces turned OFF for the template and ON for the target, or at least the part of the target that matters. You can also turn the MD off for the template, so it only feels the threading forces. But you have to be aware that the atoms thus become invisible to each other which may cause the issue you are seeing.

2. The base pairs might not be strong enough, so try increasing the strength. Sometimes a factor of 1000 is good. It could also be that they do not form because of a steric barrier or other constraint. Maybe you are applying forces that are not consistent with each other, or which are impossible to satisfy simultaneously.

Let me know if this doesn't work. you might email me directly if I don't answer on here .. samuelfloresc@gmail.com

Sam