Mdrun Instability with deprotonated ARG topology with CHARMM36

GROMACS version: 2022.4
GROMACS modification: No

When altering the protonation state of a positively charged ARG to a neutral ARGN using topologies from pdb2gmx with the -arg flag enabled with the charmm36 FF, there seems to be massive mdrun instability the simulation ending with the dreaded LINCS or XTC warning as the system explodes. grompp has no complaints and the instability is centered around the ARGN atoms HH11 and HH21.

The system I am testing on is a coordinate file from PDB:1UBQ with no alterations and waters removed. I am running in a 999nm^3 box to simulate vacuum conditions (see Konermann, L. et al., Methods 144, 104–112 (2018) if curious), but have also run with smaller box sizes to no avail. I have also listed all the files used in my test below. This failure happens both on native windows and linux (HPC) builds with/without GPU acceleration.

Through testing as far as I can tell this is not specific to the coordinate file I am using (PDB:1UBQ), it is not specific to any one ARG residue (I have tested each ARG residue in 1ubq.pdb), and it is not specific to the charmm36 version I am using (I have tested charmm36-jul2017.ff, charmm36-mar2019.ff, and charmm36-jul2022.ff). All other protonable amino acids including LYS, GLU, ASP, HIS, and C/N-termini can handle altered protonation states fine so it seems to be specific to ARG. Altering ARG → ARGN works fine with other FF’s like OPLS or GROMOS, so it also seem to be specific to CHARMM36. I am guessing it is some improper parameters within the charmm36.ff files, which I am currently looking into, but hesitant to alter given the complexity and effect it may have on proper parametrization.

Altering protonation states with pdb2gmx manually is quite tedious so I built a python script, protonation_test.py copied below to automate the process as well as automate mdrun for em/NVT/NPT equilibration by using the subprocess module. The .py file runs on windows, but might not work properly on mac/linux without some alterations or you can just run the gmx calls normally without using the .py file. If you do run the .py file, start it in a dir with 1ubq.pdb and the respective mdp files which are also listed below.

When run, the NVT equilibriation always fails around step 600-800. The trajectory file is quite interesting as the deprotonated ARG behaves normally until HH11 and HH12 pull together seemingly randomly at step ~600 before oscillating to explosion. If you run the prodrun.mdp file first you get similar results. I have tried editing nearly every parameter of these two files to no help.

protonation_test.py file with deprotonated 42ARG (see comments in code below)

#!/usr/bin/env python
import subprocess

#expands coordinate file to 10.0 nm^3 box and removes waters
subprocess.run(['gmx', 'editconf', '-f', '1ubq.pdb', '-o', 'prot.gro', '-box', '10.0'], shell=True)
subprocess.run(['echo', 'q', '|', 'gmx', 'make_ndx', '-f', 'prot.gro', '-o', '1ubq.ndx'], shell=True)
subprocess.run(['echo', '1', '|', 'gmx', 'editconf', '-f', 'prot.gro', '-o', 'prot.gro', '-n', '1ubq.ndx'], shell=True)

#list of ubiq amino acid residues from PDB:1UBQ
proton_res = ['6LYSN', '11LYSN', '27LYS', '29LYS', '33LYS', '48LYS', \
    '63LYS', '42ARG', '54ARG', '72ARG', '74ARG', '21ASP', '32ASP', '39ASP', \
        '52ASP', '58ASP', '16GLU', '18GLU', '24GLU', '34GLU', '51GLU', '64GLU', \
            '68HISE', 'Nterm', 'Cterm']

#dict of protonation states to be fed into pdb2gmx 
#flip any of the ARG numbers from 1 (protonated ARG) -> 0 (deprotonated ARGN)
#currently 42ARG is set to be deprotonated as it is set to 0
#if 42ARG set to 1, the simulation will run with no errors to completion
proton_map = {'6LYSN': 1, '11LYSN': 1, '27LYS': 1, '29LYS': 1, '33LYS': 1, \
    '48LYS': 1, '63LYS': 1, '42ARG': 0, '54ARG': 1, '72ARG': 1, '74ARG': 1, '21ASP': 0, \
        '32ASP': 0, '39ASP': 0, '52ASP': 0, '58ASP': 0, '16GLU': 0, '18GLU': 0, '24GLU': 0, \
            '34GLU': 0, '51GLU': 0, '64GLU': 0, '68HISE': 1, 'Nterm': 1, 'Cterm': 0}

#calling pdb2gmx and feeding values from dict
proton_map_vals = []
for i in range(len(proton_res)):
    proton_map_vals.append(str(proton_map[proton_res[i]]))
proc = subprocess.Popen(['gmx','pdb2gmx' ,'-f', 'prot.gro', '-o', 'prot.gro', '-p', 'topol.top', '-ff', 'charmm36', \
    '-water', 'none', '-ter', '-lys', '-arg', '-asp', '-glu', '-his', '-ignh'], stdin=subprocess.PIPE, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
for input in proton_map_vals:
    proc.stdin.write(f"{input}\n".encode())
proc.stdin.flush()
output, errors = proc.communicate()
print(output)

#em equilibriation
subprocess.run(['echo','q','|','gmx','make_ndx','-f','prot.gro','-o','prot.ndx'], shell=True)
subprocess.run(['gmx','grompp','-f','em.mdp','-c','prot.gro','-r','prot.gro','-p','topol.top','-n','prot.ndx','-o','em.tpr','-maxwarn','10'], shell=True)
subprocess.run(['gmx','mdrun','-v','-deffnm','em'], shell=True)

#running protein in vacuum so need to expand to max box size
#see Konermann, L. et al., Methods 144, 104–112 (2018) if curious
subprocess.run(['gmx', 'editconf', '-f', 'em.gro', '-o', 'expanded.gro', '-box', '999'], shell=True)

#nvt equilibriation
subprocess.run(['echo','q','|','gmx','make_ndx','-f','expanded.gro','-o','prot.ndx'], shell=True)
subprocess.run(['gmx','grompp','-f','nvt.mdp','-c','expanded.gro','-r','expanded.gro','-p','topol.top','-n','prot.ndx','-o','nvt.tpr','-maxwarn','10'], shell=True)
subprocess.run(['gmx','mdrun','-v','-deffnm','nvt'], shell=True)

#npt equilibriation
subprocess.run(['echo','q','|','gmx','make_ndx','-f','nvt.gro','-o','prot.ndx'], shell=True)
subprocess.run(['gmx','grompp','-f','npt.mdp','-c','nvt.gro','-r','nvt.gro','-p','topol.top','-n','prot.ndx','-o','npt.tpr','-maxwarn','10'], shell=True)
subprocess.run(['gmx','mdrun','-v','-deffnm','npt'], shell=True)

#prodrun
subprocess.run(['echo','q','|','gmx','make_ndx','-f','npt.gro','-o','prot.ndx'], shell=True)
subprocess.run(['gmx','grompp','-f','prodrun.mdp','-c','npt.gro','-r','npt.gro','-p','topol.top','-n','prot.ndx','-o','prodrun.tpr','-maxwarn','10'], shell=True)
subprocess.run(['gmx','mdrun','-v','-deffnm','prodrun'], shell=True)

#deletes tpr/log files so mdruns don't rerun old files if grompp fails
#this is with windows syntax, change del to rm if on mac/linux
subprocess.run(['del', '*log'], shell=True)
subprocess.run(['del', '*tpr'], shell=True)

em.mdp

integrator	        = steep		
emtol		        = 0.01 
emstep              = 0.01 
nstxout-compressed  = 100     
nsteps		        = 10000	  	
nstlist		        = 10		    
cutoff-scheme       = Verlet  
rlist		        = 1.2		
coulombtype	        = cut-off		   
rcoulomb	        = 1.2		 
vdwtype             = cut-off
vdw-modifier        = force-switch
rvdw-switch         = 1.0
rvdw		        = 1.2		
pbc                 = xyz 		 
DispCorr            = no

nvt.mdp file

integrator		        = md		
nsteps			        = 10000		
dt			            = 0.001		
nstxout			        = 0		
nstvout			        = 0		
nstenergy		        = 0		
nstlog			        = 0	
nstxout-compressed      = 1
continuation	        = no		
constraint_algorithm    = lincs	   
constraints	            = h-bonds	
lincs_iter	            = 1		    
lincs_order		        = 4		 
comm-mode 		        = linear   
cutoff-scheme   	    = Verlet	
nstlist			        = 20
rlist 			        = 333	    
vdwtype                 = cut-off
rvdw                    = 333    
DispCorr		        = no  
coulombtype             = cut-off     
rcoulomb                = 333  
tcoupl			        = nose-Hoover
nh-chain-length         = 1	          
tc-grps			        = system	
tau_t			        = 2.0		        
ref_t			        = 370		     
pcoupl		            = no 	   	
gen_vel			        = yes		
pbc                     = xyz 
gen_temp        	    = 370  
gen_seed         	    = -1

npt.mdp file

integrator              = md      
nsteps                  = 10000    
dt                      = 0.001    
nstxout                 = 0       
nstvout                 = 0       
nstenergy               = 0       
nstlog                  = 0      
nstxout-compressed      = 0  
continuation            = yes     
constraint_algorithm    = lincs     
constraints             = h-bonds   
lincs_iter              = 1         
lincs_order             = 4       
comm-mode 		        = linear  
cutoff-scheme           = Verlet     
nstlist                 = 20  
rlist 			        = 333
vdwtype  		        = cut-off   
rvdw                    = 333      
DispCorr                = no
coulombtype             = cut-off       
rcoulomb                = 333 
tcoupl                  = Nose-Hoover  
nh-chain-length         = 1  
tc-grps                 = System  
tau_t                   = 2.0          
ref_t                   = 370            
pcoupl                  = Parrinello-Rahman   
pcoupltype              = isotropic            
tau_p                   = 4.0                
ref_p                   = 0.0        
compressibility         = 4.5e-5              
refcoord_scaling        = com
pbc                     = xyz      
gen_vel                 = no

prodrun.mdp

integrator              = md        
nsteps                  = 10000 
dt                      = 0.001   
nstenergy               = 0 
nstcalcenergy           = 0
nstlog                  = 0 
nstxout-compressed      = 0   
continuation            = yes       
constraint_algorithm    = lincs    
constraints             = h-bonds 
lincs_iter              = 1      
lincs_order             = 4          
comm-mode 		        = linear
cutoff-scheme           = Verlet   
nstlist                 = 20        
rlist                   = 333
vdwtype                 = cut-off
rvdw                    = 333 
DispCorr		        = no     
coulombtype             = cut-off     
rcoulomb                = 333      
tcoupl                  = Nose-Hoover  
nh-chain-length         = 1	                  
tc-grps                 = System    
tau_t                   = 2.0                       
ref_t                   = 370
pcoupl                  = no                    
pbc                     = xyz    
gen_vel                 = no

Any help would be greatly appreciated.
Thanks,

Michael Cordes
Baylor University