LigParGen

OPLS/CM1A Parameter Generator for Organic Ligands

NAMD Protein Ligand Complex Simulations

1. Preparing Protein-Ligand System

LigParGen server provides OPLS-AAM templates with CM1A/CM1A-LBCC charges for small organic molecules. In general, molecular dynamics simulations are focused on protein/NA-ligand interactions rather than just small molecules. For this reason, in this tutorial, a robust protocol to prepare Gromacs protein/NA-ligand systems using LigParGen server will be explained in detail.

NAMD, Chimera and python must be installed in your computer to perform this tutorial.

1.1 Obtain a complex

The tutorial uses an example structure of T4 Lysozyme L99A with a Benzene molecule bound (PDB ID: 4W52.pdb, see figure below) which contains just one chain A.

T4 Lysozyme L99A with Benzene

T4 Lysozyme L99A with Benzene

1.2 Cleaning & Preparing Protein for MD simulations

To prepare your system, this tutorial provides the python code below which uses Chimera dockprep functionality to clean the protein and to add missing residues using the Dunbrack rotamer Library. It is important to notice that hydrogens are not added to the protein in this step because it will be done automatically by a VMD script to avoid issues with atom names differences. On the other hand, for ligands, hydrogens are added using custom forcefield parameters explicitly.

## PDB_FILE SHOULD THE COMPLETE PATH OF THE FILE
## REPLACE BNZ with LIGAND resname 
## USAGE: Chimera --nogui --script "prep_prot_lig.py 4w52.pdb BNZ" 
import chimera
from DockPrep import prep
import Midas
import sys
import os
PDB_file = sys.argv[1] 
lig_name = sys.argv[2]
os.system('grep ATOM %s > %s_clean.pdb'%(PDB_file,PDB_file[:-4]))
os.system('grep %s  %s > %s.pdb'%(lig_name,PDB_file,lig_name))
protein=chimera.openModels.open('%s_clean.pdb'%PDB_file[:-4])
ligand=chimera.openModels.open('%s.pdb'%lig_name)
prep(protein,addHFunc=None,addCharges=False)
prep(ligand)
Midas.write(protein,None,"protein_clean.pdb")
Midas.write(ligand,None,"ligand_wH.pdb")

Run this code by using Chimera. Make sure you have Chimera installed and can be called from Command Prompt.

Chimera --nogui --script "prep_prot_lig.py 4w52.pdb BNZ"

This produces protein_clean.pdb and ligand_wH.pdb that will be used further. Make sure that the right number of hydrogens are added and the ligand charge matches the expected value.

1.3 Preparing Ligand for MD simulations

Upload the ligand_wH.pdb to LigParGen server to get the parameter files. For protein-ligand simulation, before uploading, make sure that ligand residue number is changed to 1 (BOSS, the core of LigParGen server, only works for a limited number of residues). Before to submit the structure to the server, user must specify two options: optimization and charge model (1.14*CM1A or 1.14*CM1A-LBCC). Optimization will perform a quick minimization of the submitted structure using the OPLS-AA and the charge model parameters generated by the server. If the ligand submitted is not neutral, just CM1A model can be applied (the charge must be specified by the user) and the scale factor 1.14 will not be used. Then, user can submit the structure and download the NAMD files, i.e BNZ.prm (coordinate file) and BNZ.rtp (topology) files.

2. NAMD simulations

2.1 Combining Protein-Ligand coordinates

Combine the PDB coordinates of the apo-protein and ligand using: cat protein_clean.pdb BNZ.pdb > complex.pdb, remove any line with TER and the connectivity information. Moreover, change the residue number of Benzene ligand for continuity according to the last residue number of the protein.

2.2 Setting up the system with VMD

System preparation for NAMD simulations is described in the following ten steps:

  1. Load complex.pdb
  2. Open Extension -> Modeling -> Automatic PSF generator
  3. Remove default Topology file.
  4. Download topology file for OPLS-AA/M from Jorgensen group and add the path to topology files.
  5. Add the path of BNZ.rtf file too.
  6. Click Load input files
  7. Click Guess and Split Chains it will ask you for location of complex.pdb again, so add the location
  8. Click Create Chains
  9. Solvate the protein using Extension -> Modeling -> Add Solvation box and use PSF and PDB files generated from AutoPSF generator.
  10. Neutralize the system using Extension -> Modeling -> Add Ions

TODO: Add the NAMD scripts for equilibration and production

Optional:

Users can set up the system using the VMD command line with the code lines below:

package require psfgen   
topology top_opls_aam.inp
topology BNZ.rtf 

pdbalias HIS HSD
pdbalias atom SER HG HG1
pdbalias residue HIS HSE     
pdbalias atom ILE CD1 CD

segment A {pdb complex.pdb}  
coordpdb complex.pdb A   
guesscoord   
writepdb complex_autopsf.pdb     
writepsf complex_autopsf.psf 

package require solvate  
solvate complex_autopsf.psf complex_autopsf.pdb -t 5 -o complex_wb 
package require autoionize
autoionize -psf complex_wb.psf -pdb complex_wb.pdb -neutralize -o ionized
set ubq [atomselect top all]
measure minmax $ubq   ## USE THE DIFFERENCE BETWEEN MIN AND MAX FOR CELL BASIS VECTOR DIMENSION
measure center $ubq   ## USE THE CENTER FOR CELLBASIS ORIGIN
exit

Run it using $VMD -dispdev text -e UNK.pgn and it will generate ionized.psf and ionized.pdb files that will be used for performing MD simulations.