Insights into the key interactions between human protein phosphatase 5 and cantharidin using molecular dynamics and site-directed mutagenesis bioassays

Serine/threonine protein phosphatase 5 (PP5) is a promising novel target for anticancer therapies. This work aims to uncover the key interactions at the atomic level between PP5 and three inhibitors (cantharidin, norcantharidin and endothall). We found that, unlike previous report, Arg 100 contributes less to PP5-inhibitor binding, and the residues His 69, Asn 128, His 129, Arg 225, His 252 and Arg 250 are of importance to PP5-inhibitor binding. The hydrophobic interactions established between the residues Val 254, Phe 271 and Tyr 276, especially Glu 253, are very important to enhance the inhibitive interaction. We suggested that, to increase the inhibitory activity, the interactions of inhibitor with three negatively charged unfavorable interaction residues, Asp 99, Glu 130 and Asp 213, should be avoided. However, the interactions of inhibitor with favorable interaction residue Arg 250 could enhance the inhibitory activity. The Manganese ion 2 (MN2) unfavorably contribute to the total interaction free energies. The coordination between MN2 and chemical group of inhibitor should be eliminated. This work provides insight into how cantharidin and its analogs bind to PP5c at the atomic level and will facilitate modification of cantharidin-like chemicals to rationally develop more specific and less cytotoxic anti-cancer drugs.


Supplementary Tables
Table S2. Calculated average distances between manganese and its ligated atoms of three complexes (in Å) compared with X-ray experimental values a a The labels of residues and the names of the atoms adopted from PDB convention are given in Figure S4.

MN1
Mutation sites are indicated by lowercase. All results are presented as the mean ± S.D. for three independent experiments.         Derivation of manganese-related force field parameters. The Mn (II)-Mn (II) center that derived from our three systems was bound with three His residues, two Asp residues, and one Asn residue on the protein side. For convenience, the residues bound to the manganese ion were approximated and we used three ethyl imidazole molecules to mimic the side chains of the three His residues. Two propionic acid molecules were used to mimic the side chains of two Asp residues, and a propionamide molecule was used to mimic the side chain of the Asn residue. On the ligand side, the integral structures of three small molecules (Cantharidin, Norcantharidin, and Endothall) were retained in our three systems. Cantharidin and Norcantharidin were present in the hydrolyzed form. We established three models (Model I) summarized in Figure S3 to perform QM computations. All QM computations were performed using the Gaussian 09 program package. 2 The center of PP5c contains two metal Mn (II), five electrons occupied the 3d orbitals of each Mn (II).Due to the total energy for the antiferromagnetic coupling being higher than for the ferromagnetic coupling by calculating the stable geometries of energy minimization, only the ferromagnetic coupling was considered for the Mn (II)-Mn (II) center. The parallel spins on each Mn 2+ were ferromagnetically coupled. 3 The charge of cantharidin, norcantharidin, and endothall is -2 in three models. The total charge in the three models is 0, and the total spin S is 11/2 for the ferromagnetically coupled state. Three models were optimized with the long-range corrected hybrid density-functional theory (DFT) method WB97XD. 4 The basic set 6-31G** was used for the carbon (C), hydrogen (H), oxygen (O), and nitrogen atoms (N). The basic set LANL2DZ was employed for manganese ions (Mn). After the geometry optimization, frequency analysis was performed at the same level to confirm the existence of true local minima without any imaginary frequency. The frequency analysis also produced the Cartesian Hessian matrices that were required for the following step.
The bond and angle force constants (K r and K θ ) were obtained using Seminario's method. K r and K θ are derived from the Cartesian Hessian matrices that obtained from the frequency analysis.
The derivations of manganese-related force field parameters are compatible with the AMBER force field. The energy function is given by equation (1 The four terms in equation (1) relate to the energies of bond stretching, angle bending, dihedrals and nonbonded van der Waals and electrostatic interactions, respectively. 5 In our study, we derived bond-stretching and angle-bending parameters for manganese ion. According to the Seminario's method, the bond-stretching force constants of bond A-B can be extracted the 3×3 submatrix from the Hessian matrix in Cartesian coordinates.  Calculations of three models as shown in Figure 1 were performed as described above using the MTK++ program. 6 Neglecting the dihedral parameters is a common procedure for this symmetric geometry of the metal coordination center. 7 Atomic single charge calculations. Large manganese cluster models containing all atoms of a bound residue were constructed to keep the native crystallographic geometry. These models (Model II) were capped with acetyl (ACE) and N-methylamino (NME) residues ( Figure S4).
Sampling was conducted by the Merz-Singh-Kollman (MK) 8 while the restrained electrostatic potential (RESP) method was used to derive atom-centered partial charges. Electrostatic potential (ESP) charge fitting was performed by QM computations using Gaussian 09, with keywords WB97XD/6-31G** for C, H, O and N, but LANL2DZ for Mn. A Van der Waals radius of 1.69 Å was assigned to the manganese center. 9 This task was conducted based on the RESP fitting protocol implemented in the MTK++ program. This protocol restrains the backbone heavy atoms (CA, N, C, O) to those values found in the AMBER parm94 force field and has given the most impressive performance in previous reports.

Calculated vibrational frequencies.
To check the quality of the force field parameters, we were analyzed with the PTRAJ module of the AMBER package. Hydrogen bonds (H-bond) were also assigned using the PTRAJ on the basis of the following criteria: the distance between the proton donor and acceptor atoms was ≤3.5 Å, and the angle formed by the donor, hydrogen, and acceptor was ≥120º . 13 The force field parameters obtained in our study were applied to the molecular dynamics (MD) simulations of three systems in explicit solvents ( Figure S2). The RESP fitting charges were only applied to the Mn (II)-Mn (II) center of three systems indicated in Figure   2. The partial charges from the force field ff99SB set were used to assign the rest of parts on three systems. A bond was created to mimic the coordination interaction between the manganese ion and manganese-related atom both on the protein side and ligand side. These systems were soaked in a rectangular box of TIP3P water 14 with the smallest distance between the protein surface and cell boundary set to 10 Å. An appropriate number of counterions were added to neutralize the global charge of the entire system. For each system, we performed three individual 50ns MD at different seeds for production phase without any restraint was performed. This was preceded by 2000 steps energy minimization with a weak positional restraint to eliminate unfavorable contacts, 500 ps slowly heating in the canonical ensemble (NVT) from 0 K to 300 K and 500 ps of density equilibration to adjust the solvent density under 1 atm pressure in the isothermal-isobaric ensemble (NPT), and followed by a 5 ns constant pressure equilibration which was performed unrestrained at 300 K. All MD simulations were carried out using a Langevin thermostat 15 with a collision frequency of 2.0 ps -1 dynamics for temperature control. Constant pressure was controlled with an average pressure of 1 atm using a Berendsen barastat. 16 The time step was set to 2 fs. The periodic boundary condition was enabled during MD simulation. The particle mesh ewald (PME) method 17 was used to handle the long-range electrostatic interactions. The distance cutoff for the real-space nonbond interactions was set to 12 Å. All bonds with hydrogen atoms were constrained using the SHAKE algorithm. 18 The three individual 50ns MD trajectories for each system was recorded every 10 ps for subsequent analysis.
Binding free energy calculation and spectrum of free energy decomposition. The binding free energy for both systems was estimated by the Molecular Mechanics-Poisson-Boltzmann Surface Area (MM-PBSA) approach as implemented in AMBER12 using the same force field parameters as described above. MM-PBSA calculations were performed on 1000 snapshots exacted from 40~50 ns production trajectories with a time interval of 10 ps. The binding free energy ( bind G  ) of each system was evaluated as follows: The molecular mechanics energy, TS is the change of the conformational entropy upon ligand binding, which was not considered here because our aim is to compare the binding affinity of three ligands in our three systems.
In order to estimate the contribution of key residues on the ligand binding, the protein-ligand interaction spectrum of each complex was decomposed based on a per-residue method 21

SASA.
The decomposition energies for each residue in the complex are further broken down into backbone, sidechain, and total energy contributions. All energy components in equation 5 were calculated using the same snapshots as the binding free energy calculation.
Computational alanine scanning mutagenesis. The computational alanine scanning mutagenesis (ASM) protocol has been widely used in structure-based drug design and protein engineering for evaluating the contributions of individual residue side chains to protein−protein or protein−ligand binding free energy and understanding the structural and energetic characteristics of the hot-spots.
It has been shown to be an effective and reliable method, and can now be applied with an accuracy of 1 kcal/mol. ASM is capable of anticipating the experimental results of site-directed mutagenesis and achieves an overall success rate of 80% and an 82% success rate in residues for which alanine mutation causes an increase in the binding free energy > 2.0 kcal/mol (warm-and hot-spots). 22 To further evaluate the impact of mutational effects on the ligands binding in our three systems, ASM was applied to estimate the relative binding free energy of different PP5c mutants to the ligands.
The binding free energies for the complex and for the alanine mutants were calculated using the MM-PBSA method. We used the same 1000 snapshots taken from the production phase MD trajectory for performing the alanine scanning. The alanine mutant trajectory was initially generated from the wild type molecular dynamics trajectory by truncating the side chains of the mutated residue at Cγ replacing them with a hydrogen atom and setting the Cβ-H direction to that of the former Cβ-Cγ. The relative binding free energy is the free energy difference between the wild-type and alanine mutants and is defined as: Site-directed mutagenesis. The sequence encoding the ORF of Human PP5 (HuPP5) was isolated from the liver cancer cell line SMMC-7721 by PCR using primers that introduced a Nde I site and 6xHis tag at the forward one, and a Xho I site at the reverse one. The PCR products were cloned into pCR2.1 (Invitrogen, Carlsbad, CA, USA) and transformed into DH5α (TaKaRa, Dalian, China) and then sequenced (AuGCT, Inc., Beijing, China). The mutants were generated using the QuikChange® site-directed mutagenesis kit (Stratagene, La Jolla, CA, USA) following the instruction manual. The mutagenic primers were designed based on the sequence of wild-type HuPP5 and are listed in Table S5. The mutants were sequenced to confirm the presence of the desired mutations and the lack of undesired mutations. Inhibition assays were carried out by adding compounds to the assay mixture 5 min before starting the reaction with the addition of p-NPP as described above. Compounds were dissolved in DMSO (dimethyl sulphoxide) to obtain 1 M stock solutions and diluted to desired concentrations with assay buffer before use. A non-enzyme reaction was used as the background control, while a non-compound reaction was taken as the full-activity control. The inhibition ratio was calculated as the percentage of A 410 values of the inhibition assay reaction divided by the full-activity control, having subtracted the background control A 410 value for both. The IC 50 values were obtained using GraphPad Prism 5.0.