Tutorial for the calculation of per-residue contributions to the energy barrier in CP2K

Mutagenesis Scan

This protocol is designed to perform a mutagenesis scan with residues such as ALA or GLY (other residues are also allowed). The mutated residues are optimized while keeping the rest of the system fixed, and then single-point calculations are performed on the provided reactant and transition state structures.

The protocol requires the following inputs:
● File listing the residues to be mutated
● Topology file in *prmtop format
● Reactant and transition state structures in *pdb format
● LEaP and CP2K input templates
● Selection of the QM region

Additionally, the following software packages are needed:
● CP2K
● PyMOL
● VMD
● CPPTRAJ (from AmberTools)
● LEaP (from AmberTools)
● ParmED (from AmberTools)

I - Input Preparation

To properly construct the QM system, a selection in the format of qm_selection.dat is required. Note that atom numbers will change upon residue deletion, while residue IDs (resid's) will remain the same in the PDB files. The preparation of the selection is explained in the I - Input Preparation section of the residue deletion protocol.

The mut_scan_qmmm_cp2k.sh script has the following usage:

mut_scan_qmmm_cp2k.sh <residue_list> <scan_type> <topology> <reactant_structure> <ts_structure> <selection> <leap_template> <cp2k_template> <qm_selection>

user@machine:~$ mut_scan_qmmm_cp2k.sh residue_list.dat GLY hpla2_ee.prmtop R.pdb TS.pdb :1-124 leap_template.in cp2k_template.inp qm_selection.dat

It prepares a directory for each residue in the list where the input files for CP2K will be output. The <scan_type> argument represents the three-letter code of the residue used in the scan. ALA and GLY are suitable for assessing sidechain contributions, whereas residues with greater volume and more degrees of freedom may complicate the analysis. Note: when mutating a CYX residue in a disulfide bridge, the other CYX is converted to CYS. The <selection> argument defines the residue range of the enzyme (its parameters are updated, while the rest of the system remains unchanged). The LEaP input should be consistent with the original parameterization. The cp2k_template.inp file is used to produce geometry optimization and single-point point energy input files.

To generate the structures and topology of the mutated enzyme, the sp_mutation.sh script is called by the mut_scan_qmmm_cp2k.sh script and has the following syntax:

sp_mutation.sh <number> <residue> <topology> <reactant_structure> <ts_structure> <selection> <leap_template>

Since mutating residues changes the atom numbering, the QM/MM settings must be updated for each mutation. The mut_qm_sel.sh script is called by the mut_scan_qmmm_cp2k.sh script to check how the mutated residue should be included in the QM layer and modify the qm_selection.dat file accordingly. It has the following usage:

mut_qm_sel.sh <number> <residue> <topology> <qm_selection>

The vmd_forceeval.tcl script is called within the mut_scan_qmmm_cp2k.sh to produce a file with the configuration of the QM layer, defined by the selection in the qm_selection.dat file.

The calculations can then be run using a for loop:

user@machine:~$ for i in GLY_*; do cd "$i" ; cp2k.popt -i opt_res_R.inp -o opt_res_R.out ; cp2k.popt -i sp_res_R.inp -o sp_res_R.out ; cp2k.popt -i opt_res_TS.inp -o opt_res_TS.out ; cp2k.popt -i sp_res_TS.inp -o sp_res_TS.out ; cd .. ; done

II - Output Processing

After running the optimization and single-point calculations, the following command allows us to extract the absolute energies and calculate the R->TS energy barrier for each residue mutation:

user@machine:~$ paste <(for i in GLY_*; do echo "$i" | sed 's/GLY_//g'; done) <(for i in GLY_*; do echo $(grep "Total FORCE" "$i"/sp_res_TS.out | tail -n -1) ; done | awk '{print $9}') <(for i in GLY_*; do echo $(grep "Total FORCE" "$i"/sp_res_R.out | tail -n -1) ; done | awk '{print $9}') | awk '{print $1,($2-$3)*627.509}' | sort -n -k1,1 > energy_differences_mut.dat

The energy barriers can be compared to the ones calculated with the residue deletion protocol using the E_diff_comparison_bar_plot.py script:

user@machine:~$ python E_diff_comparison_bar_plot.py energy_differences_del.dat energy_differences_mut_gly.dat

The energy barriers exhibit a similar trend in both cases, however, backbone deletion results in higher barriers compared to mutation by GLY. The difference is more pronounced with CYX43, as deletion results in a deprotonated cysteine residue, whereas the mutation yields a protonated cysteine. As expected, either deletion or mutation of the catalytic ASP91 results in a significantly higher energy barrier. Residues that significantly lower the energy barrier become key hotspots for targeted mutagenesis studies. It is important to note that different conformations in the R and TS states might break the adiabatic connection and lead to jumps in the PES.