SLCO2A1 simulation analysis

This repo contains files and codes relevant for analyzing MD simulations of SLCO2A1 systems presented in the mansucript Structural basis for prostaglandin and drug transport via SLCO2A1. Note that processed simulation trajectories are available at [Open Science Framework (OSF)].

Prerequisites

Directory structure

Here is the directory structure of the repo.

Analysis protocol

Index group preparation

• Identification of binding pocket residues
  • Convert the cyro-EM PDB into GRO format

    gmx editconf -f [cryo_EM_PDB].pdb -o cryo_EM.gro
    

  • Use identify_pocket.py to identify the binding pocket residues.

    python identify_pocket.py -i cryo_EM.gro -l [LIGAND_RESNAME] --c 6.0
    

    The command will also print convenient PyMOL and NDX selection commands. Here is an example STDOUT:

    Identified 33 residues in the binding pocket around the ligand 'ZLK' within 6.0 Å:
    GLN42, GLN45, LEU46, GLN188, PHE212, ALA215, VAL216, PRO219, ALA220, TYR223, PHE337, SER338, ILE341, LEU344, SER345, THR346, LEU348, ASN349, ASN364, ILE367, GLY368, ASN371, LEU372, ALA375, ALA376, MET379, ALA526, CYS530, HIS533, PHE557, MET560, ARG561, TRP565
    
    PyMOL command binding pocket residues:
    select pocket, resi 42+45+46+188+212+215+216+219+220+223+337+338+341+344+345+346+348+349+364+367+368+371+372+375+376+379+526+530+533+557+560+561+565
    
    NDX command to select the pocket backbone:
    a N CA C O & r 42 45 46 188 212 215 216 219 220 223 337 338 341 344 345 346 348 349 364 367 368 371 372 375 376 379 526 530 533 557 560 561 565
    

• Index file preparation We run the following command to create an index file index.ndx.

  gmx make_ndx -f ../production/rep_1/sys.gro -o index.ndx
  

  Default index groups will be shown in the interactive prompt. Assuming that there are 18 default groups (indices 0-17), here we define the following additional index groups useful for downstream analyses and visualizatoins. Here is where the command generated by identify_pocket.py becomes handy.

  13 & ! a H*  # Assuming the ligand has and index of 13
  name 18 LIG_heavy
  
  r DPPC | 1  # Protein and DPPC membrane, assuming the protein has an index of 1
  name 19 MemProt
  
  19 | 13
  name 20 MemProtLIG
  
  a N CA C O & r 42 45 46 188 212 215 216 219 220 223 337 338 341 344 345 346 348 349 364 367 368 371 372 375 376 379 526 530 533 557 560 561 565
  name 21 pocket_bb
  
  21 | 18
  name 22 pocket+LIG
  
  4 | 18
  name 23 bb_LIG_heavy
  
  q
  

MD trajectory preprocessing

• Concatentation of simulation trajectories To concatenate the trajectories from the three replicates, run the following command and select c (continue) in the interactive prompt:

  gmx trjcat -f ../production/rep_1/md.xtc ../production/rep_2/md.xtc ../production/rep_3/md.xtc -o md_all.xtc -settime
  

• Removal of jumps across periodic boundaries Using the following command, we select Group 4 (Backbone) for centering and Group 20 (MemProtLIG) for output.

  gmx trjconv -s ../production/rep_1/md.tpr -f md_all.xtc -o md_all_nojump.xtc -center -pbc nojump -n index.ndx -dt 200
  

• Recentering the trajectory Using the following command, we again select Group 4 (Backbone) for centering and Group 20 (MemProtLIG) for output.

  gmx trjconv -s ../production/rep_1/md.tpr -f md_all_nojump.xtc -o md_all_center.xtc -center -pbc whole -ur compact -n index.ndx -dt 200
  

• Convert the TPR file We need to create a new TPR file that matches the system of the processed trajectory. We select Group 20 (MemProtLIG) for output here.

  gmx convert-tpr -s ../production/rep_1/md.tpr -n index.ndx -o ../production/rep_1/md_system.tpr
  

Binding pose clustering

Given the concatenated and preprocessed trajectory md_all_center.xtc, we then cluster the binding poses of the ligand using the following command. We select Group 22 pocket+LIG for leqst squres fit and RMSD calculation, and Group 23 (bb_LIG_heavy) for output.

gmx cluster -f ../md_all_center.xtc -s ../../production/rep_1/md_system.tpr -n ../index.ndx -method gromos -cutoff 0.13 -g -cl

Interaction analysis

• General analysis We ran the following command to perform general interaciton analysis. This includes generating interaction fingerprints, plotting interaciton frequencies, plotting interaction time series, and analyzing Arg561-Glu78 salt bridge interactions.

  python analyze_interactions.py 
  

• Interaction maps The interaction maps generated by ProLIF cannot be exported as an image, so we used the jupyter notebook interaction_maps.ipynb to generate them, then reproduce them using PowerPoint for better quality.

RMSD analysis

To calculate backbone aligned ligand RMSD, we run the following command:

gmx rms -s ../../production/rep_1/md_system.tpr -f ../md_all_center.xtc -o bb_aligned_ligand_rmsd.xvg -n ../index.ndx

In the interactive prompt, we selected Group 4 (Backbone) for least squared fit and Group 18 (LIG_heavy) for RMSD calculation.

SASA analysis

To calculate the solvent accessible surface area (SASA) of the ligand, we run the following command:

gmx sasa -s ../../production/rep_1/md_system.tpr -f ../md_all_center.xtc -o sasa_ligand.xvg -n ../index.ndx
``
In the interactive prompt, we selected Group 18 (`LIG_heavy`) for SASA calculation.


## Reproducing results from the manuscript