There are many possible reasons for this:

• Your starting structure may contain (an) error(s). What resolution was it solved at? Anything above 2.5 Angstroms begins to steeply increase the likelihood of problems. What is the B-factor of residues surrounding the phosphate? I'd highly recommend learning how to generate the electron density map (it's quite straightforward in Phenix), and have a look at your problem region overlaid with the map. A great way to check for issues.
• Many phosphate binding sites involve metal ions (particularly magnesium). Is this the case for yours - and if so, is the magnesium present?
• At the beginning of your simulation, if you don't give sufficient time for the water to soak in (with your protein/ligand complex) fixed in space, the initial jolt can cause some substantial structural problems.
• Is your ligand a standard molecule, or did you have to generate a topology in PRODRG or ParamChem or similar? If so, is the generated topology reasonable?

Hi Dr. Croll,

Thanks for looking into my problem.

My X-ray structure is with 1.5 ang resolution.

Beta-factors are within 17 to 22 for the residues surrounding the phosphate group.

Few terminal residues are missing, which are not very near to the active site. Do they make any difference in the dynamics of the active site pocket?

No, in my case there is no metal in  the phosphate binding pocket.

I kept few water molecules, deeply buried in the pocket,  in the starting structure.

My ligand is a standard nucleotide, force field parameters are well-established within charmm or amber category.

I'll try to learn generating the electron density map and correlate the results with the event I'm presently seeing.

Thanks again for your suggestions.

Tarak

Hmm. At 1.5 angstroms there's very little chance of any serious error in the structure. Does the environment around the phosphate group "make sense" (i.e. in terms of charge compatibility with the phosphate groups)? Was this structure crystallised at a particularly unusual pH (i.e. where you might need to consider protonation states)?

Would you be willing to share the PDB accession ID?

Of course, there's always the possibility that the terminal phosphate residues is naturally mobile...

Yes, there are several hydrogen bonding partners in that pocket. You are right, there may be some effects due to the pH, few residues in that pocket may have different protonation states, need to look at it carefully.

I'm extremely sorry that I could not share this protein details as it is a collaborative work with some experimentalists. :(

Thanks for you support.

I'll check the things as you've suggested.

Building up on what has already been said: 

You said you are in equilibration mode,

how far into the simulation does the ligand start to drift away?

if your system is THAT unstable you should see this behaviour during minimization. 

Did you perform minimization?

If you did - you can get an idea of the interaction energy between the ligand and the 5A binding site residues. Is the value favorable ? It should be. Because if you got to the point of the crystal structure in a set of known conditions- it should remain that way in the simulation. 

Binding is also dependant on the type of ensemble you use. Look into that as well. 

Also during the course of the simulation you can do simple correlation analysis to see what residues are moving the most. Phosphate carries a charge - ensure the site its sitting in and in its vicinity is oppositely charged.

Another thing which I just remembered - how did you make the psf/pdb which you used with NAMD. If you used VMD, it currently uses topology/parameters from CHARMM27 and if you give it CHARMM36 it complains when you try to build psf for nucleotides - did you work through this problem? 

On a final note for the electrostatics to work out correctly you need to have a neutral system at the start. What ions are you using to neutralize? Do any of them fall close to your system?

There are many more things that you should double check apart from some that I have pointed to, however start with these.

Hope this helps.

Good Luck.

/A

Thanks Ashar.

Yes, I minimized the system prior to the NPT equilibration. After few ps it is coming out.

I used charmm36, vmd did not show any warning or error.

To neutralize the system I use 0.15 mol of NaCl. None of the ions comes close to to the binding pocket.

I'm currently checking few things related to the charge density and electrostatic plots.

In the minimization - does the ligand stay in place ? What do the energy values look like if you do a selected analysis of the binding pocket only?

Is the binding pocket oppositely charged ?

Few ps means almost at the start. This would mean that either the minimization pushed it out of its energy minimum or the charge complementarity is not working out for you.

How many steps do you minimize for? One way of being sure is not to do any minimization at all. Try this.

If you used VMD, did you use psfgen (script based) or autopsf (gui based). If you used automatic psf builder, it be defaults read charmm27. Look into that as well. 

Good Luck.

Hi Ashar,

During the minimization ligand was in the pocket, I didn't do any selected analysis yet, I'll do. The binding pocket has some residues which may in principle form hydrogen bonds, but there are few residues which are with neutral.

Minimization has been carried out for around 20000 steps.

In VMD, autopsfgen has an option to make it read user specific force field. I have deleted the default one and provide the charmm36 and the ligand.rtf file to generate the psf file for the complex.

Thanks for your suggestions. I'll surely do some analyses accordingly.

Thanks,

Tarak

Hey, 

20,000?? That big a number is really unnecessary. Usually its between 50 to 500 steps. You can try reducing this number. You can also justify this reduction by seeing that if you plot the energy during minimization with the number of minimization steps - your plot would start to converge around this mark. So start of with this.

Sorry - I read your original post as phosphate from a DNA part. I realize now that its not related to the DNA, perhaps thats why you didn't get the error. The answer still applies though. Phosphate ions have a -ve charge and the primary reason for it to stay there is because of charge complementarity.   You are looking for arginines and lysines in your binding site and look for those that move away from phosphate during the simulation. If that is happening that could explain why the phosphate sees a drift. 

Good Luck. 

Best,

/A

I hope that my answer is not too late, since it's been 2.5 months.

IMO the most plausible explanation is the absence of Mg2+ ion inside the binding cavity of phosphate. It is a common feature of ATP, ADP, GTP binding protein/enzyme to have Mg2+ ion inside the binding cavity, which provide the bridge between negative residue inside the binding cavity (aspartate or glutamate) with phosphate ion. My best guess is you're removing every non-protein atom from the crystal structure, including the Mg2+ ion. That way the negative residue close to the phosphate ion is rejecting the negative phosphate ion due to the charge similarity.

Hope that help.

Radif