III - Multipoint mutagenesis
This protocol is designed to evaluate activation energies of multipoint mutant variants. The mutated residues are submitted to molecular dynamics and geometry optimization to produce a more natural conformation while keeping the rest of the system fixed.
The protocol requires the following inputs:
file listing the mutations in the variants;
topology file in *prmtop format;
reactant and transition state structures in *pdb format;
LEaP and CP2K input templates;
selection of the QM region;
file listing cpptraj selections of atoms to keep free during GEO_OPT and MD.
Additionally, the following software packages are needed:
CP2K;
CPPTRAJ (from AmberTools);
LEaP (from AmberTools);
ParmEd (from AmberTools);
Python 3 with NumPy and Matplotlib.
Input preparation
A selection in the format presented in qm_selection.dat is required to adequately build the QM system (the atom numbers will change upon mutation, resid's remain the same). The preparation of the selection is explained in the I - Input Preparation section of the residue deletion protocol.
The mut_mp_qmmm_cp2k.sh script has the following usage:
mut_mp_qmmm_cp2k.sh [options] <mutant_list> <topology> <reactant_structure> <ts_structure> <selection> <leap_template> <cp2k_template> <qm_selection> <selection_free> user@machine:~$ mut_mp_qmmm_cp2k.sh --qmmm-reference forceeval_qmmm.inc variant_list.dat hpla2_ee.prmtop R.pdb TS.pdb :1-124 leap_template.in cp2k_template.inp qm_selection.dat free_selection.dat
It prepares a directory for each mutant variant in the list where the input files for CP2K will be output, and a WT directory with the single-points of the unmodified system. A given mutant variant can be specified as <name> <number><residue> in the file mutant variant list (e.g. MUT_001 41LYS 53ASP 74GLU 92ASP). 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, which has to be the first block of atoms of the topology (its parameters are updated, while the rest of the system remains unchanged). The LEaP input has to load the same force field that was used for the original topology, as explained in the mutagenesis scan protocol. The cp2k_template.inp file is used to produce molecular dynamics, geometry optimization and single-point energy input files. The <selection_free> argument requires a file containing cpptraj selections for each mutant variant, specifying the atoms that remain unfrozen (e.g. MUT_001 :41,53,74,92|((((:41,53,74,92)&!@N,H,CA,HA,C,O)<:2.5)&!(:1-124&@N,H,CA,HA,C,O)), i.e. the mutated residues and every residue within 2.5 Å of their sidechains, without the protein backbone). Backbone atom names should only be excluded within the protein residue range, since water oxygens are also named O; the script reports solvent molecules that end up only partly free. The selection is evaluated on the transition state structure and the same atoms are free in both states. The mp_summary.dat file reports what was done to each variant.
The first optimization and the molecular dynamics are classical, and the force field cannot describe the bonds that are being formed and broken in the transition state. For this reason, the atoms of the QM layer that are selected by <selection_free> and do not belong to a mutated residue are kept fixed (they are listed in mp_summary.dat). The --allow-free-qm option turns this off, and should only be used if none of those atoms takes part in the reaction. The length, time step and temperature of the molecular dynamics can be changed with the --md-steps, --timestep and --temperature options.
To generate the structures and topology of the mutated enzyme, the build_mutant.sh script is called by the mut_mp_qmmm_cp2k.sh script and has the following syntax:
build_mutant.sh [options] <mut_name> <topology> <reactant_structure> <ts_structure> <selection> <leap_template> <mutation> [<mutation> ...]
The mutated residues keep their backbone atoms (and CB, unless GLY is involved), and LEaP rebuilds the rest of each residue from its template. The place-sidechains subcommand of qmmm_mut_tools.py then samples the sidechain torsions of the new residues (staggered values for sp3 torsions, 30° steps for planar groups, followed by a local refinement) and keeps the conformation with the fewest steric overlaps with the rest of the system, including the solvent and the substrate. This is done separately for the reactant and the transition state, and the result of each residue is written to placement.log in the mutant directory. The --no-placement option keeps the geometry built by LEaP.
Since mutating residues changes the atom numbering, the QM/MM settings must be updated for each mutation. The qmmm_mut_tools.py script checks how the mutated residues should be included in the QM layer and modifies the qm_selection.dat file accordingly. It has the following usage:
qmmm_mut_tools.py mutate-selection -w <topology> -m <mutant_topology> -q <qm_selection> -o <new_qm_selection> <mutation> [<mutation> ...]
The qmmm-section subcommand of qmmm_mut_tools.py then writes a file with the configuration of the QM layer (QM cell, QM atoms and link atoms), defined by the selection in the qm_selection.dat file. The resulting QM layer is then checked against the topology (link atoms, charge, number of electrons), and no inputs are written for a variant that fails these checks.
By default, the QM cell is built around the QM atoms of each system (5 Å on each side) and the link atom parameters (ALPHA_IMOMM and QMMM_SCALE_FACTOR) are estimated from the charges of the MM atoms. If the QM/MM section used to obtain the original reaction path is available, it should be given with --qmmm-reference forceeval_qmmm.inc (the file used for this example is in the hpla2_example folder). The QM cell is then copied from that file, and every link atom that also exists in the original system keeps its original parameters, so that the WT single points reproduce the original setup and all variants share the same QM cell. Only links that did not exist before (for instance, around a mutated residue of the QM layer) get estimated parameters, and they are listed in the qmmm_section.log file. A system whose QM atoms do not fit in the reference cell (less than 2 Å from its edges, see --min-margin) is skipped.
The calculations can then be run using a for loop:
user@machine:~$ cd WT ; cp2k.popt -i sp_res_R.inp -o sp_res_R.out ; cp2k.popt -i sp_res_TS.inp -o sp_res_TS.out ; cd .. user@machine:~$ for i in MUT_*; do cd "$i" ; cp2k.popt -i opt_md_res_R.inp -o opt_md_res_R.out ; cp2k.popt -i md_res_R.inp -o md_res_R.out ; 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_md_res_TS.inp -o opt_md_res_TS.out ; cp2k.popt -i md_res_TS.inp -o md_res_TS.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
Output processing
After running the optimization and single-point calculations, the extract_barriers.py script extracts the absolute energies and calculates the R->TS energy barrier for each mutant variant and for the unmodified system:
user@machine:~$ extract_barriers.py --wt-dir WT -o energy_differences_multi.dat MUT_*
Entries whose optimization stopped at the maximum number of steps without converging are reported by the script, and left out with the --strict-opt option.
The energy barriers can be plotted with the E_diff_multi_bar_plot.py script:
user@machine:~$ E_diff_multi_bar_plot.py energy_differences_multi.dat
The calculated energy barriers of the mutant variants can be compared with the original energy barrier to see if the introduced mutations stabilize or destabilize the transition state of the reaction step. Since the reactant and the transition state of each variant are relaxed independently, the mutated residues and their surroundings may end up in different conformations in the two states, which adds a contribution to the energy barrier that is not related to the reaction. In this example the effects are close to additive: MUT_001 lowers the barrier by 5.8 kcal⋅mol-1, and the GLY scan values of the same four residues add up to 6.0 kcal⋅mol-1; for MUT_003 the two values are 7.6 and 7.7 kcal⋅mol-1. It is advisable to compare the optimized reactant and transition state structures of each variant, and to repeat the procedure from different molecular dynamics lengths or velocities before drawing conclusions from small differences.