GROMACS version: 2023.3
GROMACS modification: No
I simulated a small metalloprotein with Cu and Cd ions in Charmm36 force field. I tried different lengths of equilibration and md with following parametrs:
NVT (tcoupl= V-rescale, define = -DPOSRES) 0,1-1 ns
NPT (pcoupl= C-rescale, tcoupl= V-rescale, define = -DPOSRES) 0,1-1 ns
MD (pcoupl= C-rescale, tcoupl= V-rescale) 50-200 ns
Then PBC effects were removed and protein was centered
But visualization of trajectory and calculation of the RMSD showed several problems:
protein unfolds completely after some time (<50 ns)
all metal ions are moving from the binding site
I add my mdp files and some images of protein before simulation and after unfolding/
I’ll be grateful for advice if this problems can be fixed or this behavior of the protein may be expected? nvt.mdp (2.4 KB) npt.mdp (2.6 KB) md.mdp (2.5 KB)
That’s a lot of highly charged species in tight spaces. Have you properly protonated the coordinating residues (e.g. anionic Cys, appropriate His, etc)? Normally one would apply covalent bonds to these kinds of coordinations, which you would have to add yourself via specbond.dat entries.
Yes, this metalloprotein has ~20 cystein amino acids that coordinate 7 divalent metal cations (zinc or cadmium).
At the first time I used general command to create a topology (and no bonds were created)
gmx pdb2gmx -f protein.pdb -o protein_processed.gro -water tip3p -ter
Then according to your advice, I tried to create bonds with metal ions (in this case cadmiun). For this purpose I added two files to my working directory:
residuetypes.dat with additional line
CD2 Protein
specbond.dat with additional line
CYS SG 1 CD2 CD 4 0.3 CYS2 CD2
I measured distance between sulfur and cadmium, which is 3.0 A.
Then I tried to create my topology with a modified command:
gmx pdb2gmx -f protein.pdb -o protein_processed.gro -water tip3p -ignh -ter -merge all
But I still don’t see any changes in topology file or or produced protein_processed.gro file.
I add my files to this message. May be I am missing something, so no additinal bonds are generated in topology. specbond.dat (374 Bytes) residuetypes.dat (2.7 KB) topol.top (214.6 KB)
I also add the end of my protein.pdb file, where ligands are specified as the .pdb extension is not allowed on this site:
ATOM 384 N CYS A 59 11.829 -8.504 6.671 1.00132.16 N
ATOM 385 CA CYS A 59 10.678 -9.138 7.248 1.00132.16 C
ATOM 386 CB CYS A 59 9.583 -8.114 7.588 1.00132.16 C
ATOM 387 SG CYS A 59 10.200 -6.820 8.706 1.00132.16 S
ATOM 388 C CYS A 59 10.096 -10.150 6.317 1.00132.16 C
ATOM 389 O CYS A 59 9.516 -11.137 6.765 1.00132.16 O
ATOM 390 N CYS A 60 10.174 -9.928 4.995 1.00148.29 N
ATOM 391 CA CYS A 60 9.565 -10.889 4.126 1.00148.29 C
ATOM 392 CB CYS A 60 8.341 -10.289 3.458 1.00148.29 C
ATOM 393 SG CYS A 60 7.277 -9.709 4.799 1.00148.29 S
ATOM 394 C CYS A 60 10.562 -11.306 3.098 1.00148.29 C
ATOM 395 O CYS A 60 10.247 -11.490 1.923 1.00148.29 O
ATOM 396 N ALA A 61 11.811 -11.505 3.542 1.00 48.50 N
ATOM 397 CA ALA A 61 12.825 -11.922 2.628 1.00 48.50 C
ATOM 398 CB ALA A 61 14.091 -11.052 2.699 1.00 48.50 C
ATOM 399 C ALA A 61 13.191 -13.337 3.042 1.00 48.50 C
ATOM 400 O ALA A 61 12.248 -14.152 3.215 1.00 48.50 O
ATOM 401 OXT ALA A 61 14.411 -13.624 3.189 1.00 48.50 O
TER 402 ALA A 61
ATOM 403 CD CD2 B 62 -3.443 2.400 -8.436 1.00 16.05 CD
ATOM 404 CD CD2 C 63 -5.195 5.082 -10.843 1.00 16.50 CD
ATOM 405 CD CD2 D 64 -5.175 5.985 -6.744 1.00 12.10 CD
ATOM 406 CD CD2 E 65 8.984 -7.285 4.859 1.00 11.34 CD
ATOM 407 CD CD2 F 66 3.168 -4.634 3.144 1.00 16.09 CD
ATOM 408 CD CD2 G 67 5.408 -8.528 2.824 1.00 20.57 CD
ATOM 409 CD CD2 H 68 7.228 -4.839 1.900 1.00 14.94 CD
END
Please provide the PDB coordinates with a .txt extension. 3.0 Å sounds like a very long bond, even for something involving these species, but I guess it could be possible.
You also probably want CYM, not CYS2, as the resultant ligand. Divalent metals are Lewis acids so they lead to deprotonation (ionization), which is not the same as the oxidation occurring upon disulfide formation.
A distance of 0.3 nm in specbond.dat should work here. Everything is certainly within 0.3 ± 0.03 nm as far as I can tell. Check the terminal output from pdb2gmx carefully to see the distance matrix it prints and whether it indicates linkages. Your Cd2+ ions are in separate chains from the protein, so unless you use the -merge option with pdb2gmx, the bonds won’t be created. Or put everything in chain A.
Thank you very much. It finally worked for me with command line:
gmx pdb2gmx -f protein.pdb -o protein_processed.gro -water tip3p -ignh -ter -merge all
I defined terminus for protein NH3+ and COO-, and for Cadmium None.
But now i get errors at the grompp stage like:
No default Bond types
No default U-B types
No default Proper Dih. types.
They corresponds to the lines where CD is interacting.
Changing some input data (like ATOM to HETATM for CD or terminus options hasn’t worked).
I can guess (but I am not sure) that there is no specific data for CYS-CD bonds in the ffbonded.itp file. But I have no idea how to fix it correctly.
For Zn ion there is no error (may be because of the default data available in ffbonded.itp file).
