Well, that's a tricky question without any information on your biological system. If both sites are symmetrical and filled with the same ligand, I would first think of an artifact and check the available biochemistry data.

Well, that's a tricky question without any information on your biological system. If both sites are symmetrical and filled with the same ligand, I would first think of an artifact and check the available biochemistry data.

Information about the MD package, force-field and the total simulation time would also be helpful.

Based on your current description, one could suggest that you run a longer simulation to follow the development of your observation in time. It might occur in the other active site too, but slower.

Additionally, you might want to look at the ions and the solvent molecules.

As for the charges, you can compare your ligand with other ligands. As long as the parameters are meant for the same simulation package, you should be fine.

Thanks Karine and Balder

Sorry for the concise intro with the question.

This is a homodimer (same sequence), not symmetric but the structural aspects are almost same in both the monomers. Same ligands are bound in each of the monomer with similar sorts of interactions with the neighboring residues.

I'm runnig the simulations in Gromacs-4.5.5 with amber94 force field. The RESP charges were derived from RED tools and/or server. Parameters for the ligands were assigned based on the chemical analogy with similar kind of molecules present in the amber force field.

Up to now I ran the simulation for 5 ns and continuing to 50 ns.

To add a bit to what both Karine and Balder have said, it would help to know how you derived your original model (i.e., what crystal structure was it from and what was the associated space group), and also a bit more detail about how you prepared the structure (what equilibration protocols did you effect).

I think what we are all saying is that there are conditions under which the two components of a homodimer can really (biologically) behave in a somewhat distinct non-equivalent fashion, but also that there are lots of ways to artificially (or accidentally) bias the system toward unrealistic inequivalence.

If you started with a high symmetry system and were careful in your model preparation, the resulting simulation should not stray far from equivalent monomer behavior and deviations should revert to the mean.

Depend on the size of your dimer as well as the force field package. Take into account the computational cost of your simulation. As G. Lushington said, you should be careful with your model preparation, once, it is a very important step to get a correct information about your system.

Are there water molecules involved in the ligand binding , probably they have a role in your observation

I agree with Karine Loth's answer - it is probably an artifact. This being said, it is definitely relevant to include water in your simulation - it is not clear from your question whether you have done this. Is your molecule a dimer in the crystal or is it a dimer of dimers - if the latter is true, you will get artifacts when simulating only the dimer.

You need to carefully prepare the structure. Make sure that both have identical sequence, sometimes certain residues may be missing from one monomer and if this residue is critical for stability it will be a big problem. Second NEVER use AMBER94 forcefield, it is very outdated, you have the much better AMBERff99SB and its ILDN version which i think will do the job properly. You need also to consider to run two simulation one with and one without the ligand.

You can also consider running simulation for one monomer. And you have to do careful minimization and long equilibration , etc..

Hi Tarak , You should first check your starting structure. Are you starting with Crystal structure ? Check ligand in both monomer , sometimes bond is not correctly assigned. Also check that you are using same topology file for both ligand Also check size of solvation box and position restrain setting. If System is stable in NPT equilibration then you can go to Production simulation.

It depends on preparation of input file, may be initial structure is far from the final optimum structure. Some times all phase space for the structure are not accessible. For the time scale problem increasing time is not enough toll for the simulation. It is better to use the annealing technique due to find global optimization structure in your molecular dynamics simulation.

What program and/or are you using? Is there a substrate on the active sites? Are there solvent orany other small molecule(sugar, phosphate,etc) added orleft out?

Hi,

To put a simple answer to this question is

"There is something wrong with the way you are simulating."

Key point to focus on is proteins can be simulated for really small time scales to obtain energetic fingerprints. I would term small times scales to be around 2-3 ns. To study changes on the conformational level you need to simulate for very longer times scales.

Without knowing much about the system -- since you say they are exactly alike - then the binding sites are exactly alike as well. In this case unless you simulate for way longer time periods you should not be seeing what you see. Moreover changes should be symmetrical.

Therefore the simulation process is flawed.

Now at this point we can all guess wildly at what the problem could be.

So here is my guess.

I work with CHARMM. I have never faced a problem like this - but if i had - i would have looked into my pre-processing. In my protocol i fix aminoacids around the binding sites that might have variable charges at the ph relevant to where i am working. if neutral then you might need to fix Histidine. Apart from that some amino acids in the PDB have alternate conformations - you sure that if there were any alternate conformers in the PDB of the crystal structures you took you picked the same same for both proteins (the dimer.)

I have seen some people compromising on water molecules. If the binding site shows issues like yours is doing maybe waters are relevant .. one guess can be that you stripped the waters from the PDB - and in your solvation that is if you are doing it - added your own shell. Maybe that shell somehow kept essential waters for part A of dimer in tact and didn't for part B and maybe that's why part B is showing instability.

another guess would be that you parent structure was missing atoms which your simulation protocol guessed. In which case you should have run some form of minimization prior to your simulation.

I cant think of a reason why i would need to simulate a dimer. I mean there is no reason to do that if there are two ligand that are bound to two proteins presenting two distinct binding sights. Are you sure in nature it works as a dimer and just not in the PDB. The author of the PDB might have published 2 same structures and some parts of one might be missing. You perhaps need to review your need for simulating a dimer?

Going on further - if there are two structures in the pdb and you are using one simulation to simulate the entire thing - then why do you expect them to generate the same results. Maybe its just a relaxation event that will stabilize in a while. Maybe the structure showing the difference is right and the one that's not moving is stuck in a minimum.

I think i should stop like i said - it could be any ones guess.

Try looking into the alternate conformers, charged/uncharged amino acids and the water guesses. Maybe 1 of them might work for you.

Good Luck.

Best,

/A

"This is a homodimer (same sequence), not symetric but the structural aspects are almost same in both the monomers.".... Hum, very strange, it could be helpfull to have the crystallographic space group, resolution, Rfac and Rfree, it's a homodimer or you have two monomer in asymetric unit, in this case the cristal packing could made some differences in the two structures.

after that, if it's really an asymetric homodimer, it's probably because the two monomers do not play the same role.

Indeed, homodimers are predominantly symmetric. Asymmetry can arise from binding events specially with DNA or to promote 2:1 ratio. If I have well understood, it is not what you are doing, so I agree with all suggestions : check the pdb, check your structure preparation/equilibration, change the force field...

It seems to be artifacts. You could search for the availability of wet laboratory data for molecular interaction between your protein(?) and ligand. If both are known probably you can have the help of one chemistry person for there possible biochemical interactions. After that you can compare your simulated results with all the possible biochemical interactions. If both your compounds are known then you wont get much problems. However, it is always anticipated to not ask tricky questions in general if you need very particular answers from the experts. If you want privacy of your compounds then you can submit your questions personally to the experts of your choice rather than putting questions in open forums without any details.

If you could demonstrate sufficient sampling of the corresponding trajectory, then, yes, what you observe is a 'real' consequence of the given system+forcefield+simulation protocol. In the absence of demonstrated sufficient sampling, any interpretation is necessarily tentative (especially when you attempt to imply biological relevance). My twocents.

And remember when running MD simulations: if you see something happening once, it basically never happened... I.e. statistical value zero...

Ideally you should run a few simulations with different starting conditions (e.g. different seed for random velocities).

Just to add to all the other excellent comments.

Try multiple (8-10) short simulations (2-3ns) with proper heating protocol and check RMSD, RMSF, protein volume, pre/post equilibration contact maps of protein-ligand. If you still see the same behavior, take couple of simulations and extend them as long as your computers let you (ideally tens of nanoseconds). If the behavior still persists then you can be reasonably sure that it is not an artifact. Without knowing more about your system this is the best suggestion I have. Good luck!

As Pavan said check the intermediate snapshots and see what happened to the whole system and RMSD too.

What program are you using? without any information on your biological system--since you say they are exactly alike - then the binding sites are exactly alike as well.

Is ligand present in both monomers when you carried out teh simulation? If only one monomer has ligand occupancy and the other one does not, then the active site of the unoccupied monomer can fluctuate a lot depending on the protein.

Thanks to all for your nice comments and suggestions.

I'll cross check the pdb structure (space group, symmetry and all) before running the simulations.

As suggested by so many experts the best way to deal with this kind of problem is to run few more independent trajectories and then get the statistical interpretation.

Again thank you all.

Ofcourse you need to rerun simulations with different starting velocities and get statistics. But assymetry has been suggested to occur in homodimers before, have a look at this paper and also at papers citing it:

http://www.ncbi.nlm.nih.gov/pubmed/17089205

"Dynamic models of G-protein coupled receptor dimers: indications of asymmetry in the rhodopsin dimer from molecular dynamics simulations in a POPC bilayer."

In case you are to model a homodimer and you have an apprpriate template, I suggest to use the Modeller. In case, you do not have the appropriate template, you may want to try a protein-protein docking program. To this end, you should have information about the particular amino acids, that participate in coupling the homodimer. At least a few of them....

Dear Tarek,

A simulation is just a simulation!

You are not "touching a reality".

Run simulations with different initial velocities and check for trends in what is coming out. Put global angular momentum and check again.

After a while you will be more confident about the type of results the model delivers.

Keep us informed if you have time for that.

Thank u Sir @Prof. Orlando

It's my pleasure to have all experts near me. Comments, advices from more experienced people in this field always encourage new comers like me to do better job and deal with problems. I'll update this post with the observations I'll get after running few independent trajectories.

Thanks to each of the experts for guiding me.

In addition to my previous suggestions, you may also want to run fast dynamics perturbation analysis (BMC Struct Biol. 2008 Jan 30;8:5. doi: 10.1186/1472-6807-8-5.)

In order to perform coarse grain simulation we use Gromacs program with MARTINI force field.

Tarak from your question itself a number of other questions arise. For eg.

1) The asymmetry may have its origin in your starting file so is your structure a modeled structure or is it a crystal structure. In both cases asymmetries might result from missing loops, residues etc. I just hope you have thoroughly looked for any inherent discrepancies in your input file and filled the missing parts. Also I hope that you have run at least a few steps of energy minimization before simulating your dimer.

2) The asymmetry might indicate how the thing unfolds in real time. I also work with a homodimer and both monomers may not at one time react equally to an activation event. I just think you will have to analyze the trajectory of the monomers individually and try and make sense out of it by coupling it to experimental data.

3) I would also do one other thing. I would generate two seperate PDB files one with each monomer(keeping the PDB co-ordinates same i.e deleting one monomer at a time and saving the rest as an individual PDB file). Then I would simulate both files separately and analyze the trajectories. If you still see the same issue in one of the files then the problem/asymmetry

lies with your input file itself.

Thanks for the elaborative answer

I'm runnig the simulations and as you suggested I'll analyse each of the monomer individually and then try to interpret the observations.

After running few trajectories I got some information about my protein dimer. The overall rmsd of the dimer is ~1.5 A which is very good. However, one important residue present near the ligand binding region is flexible. This residue is expected to interact with the ligand and should stabilize it inside the active pocket. Sometimes it is moving away from the ligand molecules.

Is there any charge problem for the ligand?

Note: The charges of all the atoms present in ligand molecules are assigned based on the RESP model in REDtools or Server.

Hi,

Isnt that how its suppose to behave biologically.

All ligands have binding affinities and disassociation constants.

Disassociation is a measure of the unbinding of a ligand.

Chemically speaking that is saying that the ligand undocks --- which happens because the attaction between the site and the ligand weakens allowing the ligand to diffuse away.

If the residue would n't move away - the energy would never go to a point where the ligand can disassociate.

You need to understand - what makes the molecule move away (assuming the "Away" state is significant meaning its keep showing up - rather than being visited only a few times)

Once you have ascertained this - consider an area (its actually volume) of 5~7 Ang around your residue of interest in the binding site. Get all the residues in that region and plot individual energies between each pair.

When you over-lay this spctrum plot with your simulation - timestep RMSD plot (of th residue of interest) you will see what the spike in RMSD is augmented with in the Spectrum. If there is an association between the plots - you can perhaps call it a significant state. Which is what (i think happens in real - allowing the ligand to disassociate).

The above assumes that there are no fundamental flaws with your data.

If not, with this you can safely categorize the motion of residue of interest as a significant or non-significant state.

Hope this was relevant and helpful.

Best,

/A

Thank u Ashar for nice suggestions.