ALK_files_2021_PNAS

README: supporting data for 2021 manuscript PNAS Patil et. al. There are three folders F_unbiased that includes the files to setup and perform Unbiased Molecular Dynamics, F_metad contains files to setup and perform metadynamics, F_post_processing contains files to perform post processing

Folders

F_unbiased: This folder has the input gromacs .mdp scripts required to perform Unbiased Molecular Dynamics and the active and the inactive structures of ALK
1. The topology files .top are generated using Gromacs function : "pdb2gmx" and we implement charmm27 forcefield
2. The equilibration and production procedure is as mentioned in the SI
3. The outputs (.gro, .cpt, .top, .tpr) of the unbiased MD simualtion are inputs to Metadynamics. The active and inactive ALK structures at 101 ns of MD simulation were chosen as reference in the Metadynamics run. These structures are included in F_metad as active_wt_alk_CA_meta_min.pdb and inactive_wt_alk_CA_meta_min.pdb

F_metad: This folder has the input scripts required to perform Metadynamics and also the output free energy file that is obtained by summing over HILLS files

1. plumed.dat - a plumed script that needs to be run for performing Metadynamcs
   
2. active_wt_alk_CA_meta_min.pdb - the reference active structure, one of the two collective variables is RMSD from this structure
   
3. inactive_wt_alk_CA_meta_min.pdb - the reference inactive structure, one of the two collective variables is RMSD from this structure
   

PS note: Additionally, you would need the usual gromacs file: .mdp file, .top file, .cpt file, .tpr file which are the outputs of unbiased MD

4. fes_650.dat - The free energy file which has free energy at every coordinate on CV space obtained by summing over HILLS files. The reason this is named fes_650.dat, since each walker runs 650 ns, so the aggregate time that we have is 2.6us (4 walkers)
   

F_post_processing: This folder contain the python scripts used for post-processing trajectories data to calculate hydrogen bond occupancies, plot free energy landscapes, and check convergence of metadynamics free energy zones

ALK_final_frames: This folder contains the .pdb structures of active and inactive ALK wildtype and 42 variants obtained from last frames of the Molecular Dynamics simulations.

Forcefield implemented

Our MD and metadynamics simulations use charmm27.ff from GROMACS 5.0.7

Procedure

1. MD simulation in GROMACS

The MD simulations are performed in GROMACS. The paper mentions Biophyscode which actually is a wrapper on Gromacs. Please feel free to try out Biophyscode, a creation from Radhakrishnan lab. https://biophyscode.github.io. Please note: Following is the walk-through on how to run a simulation using Biophyscode:

"https://biophyscode.github.io/molecular_dynamics_lab/".

All mutations are introduced using BioPhysCode Automacs routine based on MODELLER. Please see MODELLER here:

https://salilab.org/modeller/

All the requisite files needed to setup and run the simulations are included in the F_unbiased. We followed Bevan Lab Tutorials: "Lysozyme in water" example for equilibration and production
http://www.mdtutorials.com/gmx/lysozyme/index.html

2. Enhanced sampling simulation of Metadynamics using PLUMED

We patch GROMACS 5.0.7 with PLUMED 2.3.5 and use multiple walker (num_of_walkers=4) metadynamics. For installation procedure of PLUMED and the subsequent patching with GROMACS, see here: https://www.plumed.org/doc-v2.6/user-doc/html/_installation.html

In the folder F_metad we have the scripts required to run the metadynamics. Please note: Metadynamics has to be run using the output .gro and .cpt files of the unbiased MD simulations. The folder contains the .mdp file, plumed.dat and the reference active and inactive ALK alpha Carbons structure files

plumed.dat metadynamics script is run in GROMACS 5.0.7 patched with PLUMED 2.3.5 using multiple walker (num_of_walkers=4) metadynamics. We set in the PLUMED script, the energy in kcal/mol and length in Å. The parameters used in this study to perform Well-Tempered multiple walker metadynamics (WTMD) are:bias factor γ = T +∆T/T = 20, height = 0.6 and pace = 500.

Aggregate of 2.6us metadynamics run in GROMACS 5.0.7 patched with PLUMED 2.3.5 using multiple walker (num_of_walkers=4) metadynamics until the convergence criterion is met for the zones of interest (See section C1 in the SI of Patil et. al. PNAS 2021

The output of the metadynamics run are the HILLS files. We will have four of them since we used four walkers. To get the free energy file from the HILLS use plumed sum_hills action of PLUMED (https://www.plumed.org/doc-v2.5/user-doc/html/sum_hills.html)

3. Description of the post-processing codes

Analysis of the trajectories was done using python and mostly MDanalysis package in python: https://www.mdanalysis.org

The folder F_post_processing contains three codes:

1. Free_energy_convergence.ipynb is the code to plot the free energy landscapes using fes.dat file and also to give the free energy estimate of the zone of interest using the equation 4 under section C1 in SI.
2. hbonds.py is the code that takes in trajectory files from the zones as input and outputs things required for the next code - plot_zone.py (to reproduce plots in Fig 3A and 3B in main text)
3. plot_zone.py takes in the output from hbonds.py and reproduces the plots in Fig 3A and 3B in main text)

Citations

If you have any suggestions or queries please feel free to reach out at : patilk@seas.upenn.edu
If you found the above scripts and/or codes helpful in your work, please cite:

1. Jordan, E. Joseph, et al. "Computational algorithms for in silico profiling of activating mutations in cancer." Cellular and Molecular Life Sciences 76.14 (2019): 2663-2679.
2. Patil, Keshav, et al. "Computational studies of anaplastic lymphoma kinase mutations reveal common mechanisms of oncogenic activation." Proceedings of the National Academy of Sciences 118.10 (2021).