How to avoid an 'isopeptide bond' to break in GROMACS?

GROMACS version: 2020.4
GROMACS modification: No

Hi,
I am trying to simulate a diubiquitin chain, where the NZ of a lysine in one ubiquitin (atom 2464) is covalently linked to the C-ter of the other (atom 1962). In the PDB file this ‘isopeptide bond’ is defined at the end by:
CONECT 1962 2464
CONECT 2464 1962

The CONECT lines are eventually ignored when pdb2gmx writes the topology, and apparently when H atoms are added in both molecules somehow the double bonding of C=O at the C-ter breaks down into C-O (x2). So now the “extra oxygen” clashes with NZ, breaking the isopeptide bond I am trying to follow here. Carrying out the simulation this way both of the involved side chains displace both molecules considerably which of course is not what I am looking.

Trying to overcome this I’ve followed a few strategies. At first, and ignoring the appearance of the extra oxygen mentioned above, I applied position restraints at the C-ter and NZ atoms in an attempt to mimic the isopeptide bond as much as possible. After failure I tried to restraint the whole residues (GLY_C-ter and LYS_NZ) but given that they are in separate chains, numbering issues regarding the atoms result in an error message similar to ‘Atom index is out of bounds’ (even when the position restraints were included in the topology file at the right position). Finally, I managed both molecules as a single chain and again tried it with position restraints and even as ‘frozen groups’; although in both cases the atoms remained well fixed, now an artifact is formed between C-ter atom of the 1st ubiquitin and both N-ter and NZ atoms of the second (giving place to this weird ‘pentavalent carbon’).

Is there any other strategy I could try? Or a manner in which I may force this covalent bond to stay? I’m frankly lost.

Thank you very much,
Fidel.

Branched bonds such as these require the use of specbond.dat to define the bond. You will have to create a new lysine .rtp entry with a different name and only one HZ proton, and use that residue name in specbond.dat as the renamed residue, and then define the bond between NZ and the carbonyl C of the C-terminus. You will also have to choose ‘None’ as the C-terminal patch type, though I don’t know if pdb2gmx will complain about that case (missing atoms). But that’s the only way to attempt making such a bond.

Thank you very much for your reply @jalemkul. I tried your suggestion in accordance with the manual. So I created a new lysine (LYS1) which was defined within the aminoacids.rtp file (I copied those data from LYS and simply erased the extra H atoms):

[ LYS1 ]
[ atoms ]
N NH1 -0.47 0
HN H 0.31 1
CA CT1 0.07 2
HA HB 0.09 3
CB CT2 -0.18 4
HB1 HA 0.09 5
HB2 HA 0.09 6
CG CT2 -0.18 7
HG1 HA 0.09 8
HG2 HA 0.09 9
CD CT2 -0.18 10
HD1 HA 0.09 11
HD2 HA 0.09 12
CE CT2 0.21 13
HE1 HA 0.05 14
HE2 HA 0.05 15
NZ NH3 -0.30 16
HZ1 HC 0.33 17
C C 0.51 20
O O -0.51 21
[ bonds ]
CB CA
CG CB
CD CG
CE CD
NZ CE
N HN
N CA
C CA
C +N
CA HA
CB HB1
CB HB2
CG HG1
CG HG2
CD HD1
CD HD2
CE HE1
CE HE2
O C
NZ HZ1
[ impropers ]
N -C CA HN
C CA +N O
[ cmap ]
-C N CA C +N

Further I declared LYS1 as Protein in the residuetypes.dat and updated specbond.dat by adding the line:

GLY C 1 LYS1 NZ 1 0.2 ISO ISO

My guess here is that LYS1 would now appear as an option when I run pdb2bmx with -inter to define the protonation state of each lysine, but this is not the case. I am using CHARMM27 as forcefield. Is there any other file that I miss to modify? Since I am not adding any new atom or bond types I did not edit the .itp files regarding bonded or non-bonded interactions, but perhaps I am supposing it wrong.

Besides, by choosing ‘None’ for the C-terminal patch type the error message regarding ‘dangling bond at the terminal end’ actually appeared. So I tried to modify the c.tdb file in a similar way that LYS1, tagging this new option as ‘CO-’. But again, no option appeared available in the list after -inter was used.

I would thank you for any other advice,
Fidel.