Sim of multiple small molecules with protein simultaneously

[agoliaei]

Hi,
Thank you for sharing small molecule parameters with the community.
I was wondering about your experience in simulating many small molecules with one protein at the same time. Have you done such simulations? If so how did you solve self assembly of small molecules?
Is there a way to turn interaction between small molecules off? If so, how is it done in Gromacs?
I am curious to hear your thoughts on this.
Thanks again.

[agoliaei]

I found a way for doing this. I'm going to close it as it is no longer required.
Thank you again for sharing these parameters.

[ricalessandri]

@agoliaei I'm very sorry for missing the notification on this discussion! Thanks for reaching out here. You're welcome, happy to hear the parameters are useful.

For what my answer would have been: these are my first thoughts:

1. Do you have any reason (best if experimental evidence) that the aggregation shouldn't happen? If so, then this may indicate that the model is over-aggregating and there might be something missing. If not, consider that aggregation may be "normal". Of course, you might be using artificially high concentrations of ligands to increase binding chances so that would be OK. Just wanted to spell this out.
2. One way to achieve non-interacting small molecules, at least the easiest that comes to mind right now, would be to assign all the beads of the small molecule with a new bead name (TC5-SC3-TN5-P4 become TZ5-SZ3-TY6-H4, in which I replace C by Z, N by Y and P by H and kept the bead sizes and degree of polarity label the same), and then set to 0 all the interactions between these beads. You can set them to zero by modifying the martini_v3.0.0.itp file so that every line with the above beads will be:

   TZ5    TZ5  1 0.000000e+00    0.000000e+00
   

   This has the downsides that it won't describe intra-molecular interactions within the small molecule. This is often not an issue. It may be an issue if if the small molecule is large-ish, but you could model such interactions via bonded terms.

Not sure if this is the way you used, feel free to post how you did, if you want.

[agoliaei]

@ricalessandri thanks so much for getting back to me. I appreciate your suggestions. I don't have any experimental evidence for this. I just wanted to be able to simulate large number of small fragments to save simulation time. I ended up doing this:

1. use gromacs editconf to move the gro file to the center and right at the center add one extra atom (virtual). for example for caffeine this is the final caff.gro:
One caffeine molecule
8
1CAFF R1 1 0.116 -0.164 0.233
1CAFF O2 2 -0.123 -0.104 0.203
1CAFF R3 3 -0.298 0.005 0.075
1CAFF O4 4 -0.145 0.092 -0.108
1CAFF R5 5 0.067 0.008 -0.033
1CAFF N6 6 0.277 0.010 -0.095
1CAFF R7 7 0.108 0.150 -0.275
1CAFF V1 8 0.000 0.000 0.000
10.00000 10.00000 10.00000

2. in the martini_v3.0.0.itp file under section [ atomtypes ] add one extra bead like this:
VS 0.0 0.000 V 0.0 0.0
and under section [ nonbond_params ] add one extra interaction that is only repulsive
VS VS 1 0.000000e+00 4.900000e+00

3. and finally in the caff.itp file, change according to this:
[moleculetype]
; molname nrexcl
CAFF 1

[atoms]
; nr type resnr residue atom cgnr charge mass
1 TN1 0 CAFF R1 1 0 63
2 TP1a 0 CAFF O2 2 0 0
3 TN1 0 CAFF R3 3 0 63
4 TP1a 0 CAFF O4 4 0 0
5 TC5 0 CAFF R5 5 0 0
6 TN5a 0 CAFF N6 6 0 63
7 TN1 0 CAFF R7 7 0 63
8 VS 1 CAFF V1 8 0 0

[bonds]
#ifndef FLEXIBLE
[constraints]
#endif
; i j funct length
1 3 1 0.474 1000000 ; cog
1 6 1 0.405 1000000 ; cog
3 6 1 0.599 1000000 ; cog
3 7 1 0.555 1000000 ; cog
6 7 1 0.283 1000000 ; cog

[dihedrals]
; improper
; i j k l funct ref.angle force_k
1 3 6 7 2 180 200

[virtual_sites2]
; sites positioned at 1-a
; site from funct a
5 3 6 1 0.635 ; cog

[virtual_sites3]
; peripheral virtual sites
; site positioned as a linear combination of 3 atoms
; the site is in the same plane
; site from funct a b
2 5 1 3 1 0.640 0.605 ; cog
4 5 3 7 1 0.650 0.600 ; cog

[virtual_sitesn]
; site funct constructing atom indices
8 1 1 2 3 4 5 6 7

[exclusions]
1 2 3 4 5 6 7
2 3 4 5 6 7
3 4 5 6 7
4 5 6 7
5 6 7
6 7
8 1 2 3 4 5 6 7

the virtual site is defined as one that represents all atoms (1 2 3 4 5 6 7), and the exclusions removes intra contributions. this VS is only there to push two molecules apart when they get closer to one another.

I don't know if there is a right way of doing this. I think all of these methods are artificial, but I think the possible issue I see with your approach is that multiple ligands could end up occupying the same pocket in the protein. am I wrong? my approach is also not great. the ligands push each other out of the active site if they get close to each other in the pocket. so I think at least it's nice we have option. we can choose accordingly.

I'm curious to hear your thoughts on my approach.
thanks again. this is very useful and interesting.

[ricalessandri]

@agoliaei thanks for sharing your approach. Yes, you're right about the shortcoming of mine. I think your approach is better than mine, even though I do see the problem of the ligands pushing each other out of the pocket. Perhaps tuning the repulsion down might help with this issue? I'm thinking there might be a sweet spot that prevents (excessive) aggregation while not affecting the residence in the binding pocket. But I'm afraid overall there's no "clean" solution, at least that I can think of.