Probing the interaction mechanism of small molecule inhibitors with matriptase based on molecular dynamics simulation and free energy calculations

Matriptase is a serine protease associated with a wide variety of human tumors and carcinoma progression. Up to now, many promising anti-cancer drugs have been developed. However, the detailed structure–function relationship between inhibitors and matriptase remains elusive. In this work, molecular dynamics simulation and binding free energy studies were performed to investigate the biochemistry behaviors of two class inhibitors binding to matriptase. The binding free energies predicted by MM/GBSA methods are in good agreement with the experimental bioactivities, and the analysis of the individual energy terms suggests that the van der Waals interaction is the major driving force for ligand binding. The key residues responsible for achieving strong binding have been identified by the MM/GBSA free energy decomposition analysis. Especially, Trp215 and Phe99 had an important impact on active site architecture and ligand binding. The results clearly identify the two class inhibitors exist different binding modes. Through summarizing the two different modes, we have mastered some important and favorable interaction patterns between matriptase and inhibitors. Our findings would be helpful for understanding the interaction mechanism between the inhibitor and matriptase and afford important guidance for the rational design of potent matriptase inhibitors.


Introduction
Tumor metastasis is the main event leading to death in individual with cancer (Coghlin & Murray, 2010). In this process, proteases play an important role. They are involved in mediating degradation of extra-cellular matrix and inter-cellular cohesive structures allowing migration of cells into the extra-cellular surroundings, and stimulate various angiogenic factors and cell growth (Deryugina & Quigley, 2006;Duffy, 1987). Matriptase, as a Type II transmembrane serine protease, is overexpressed in many human cancers and proved to contribute to the tumor development (Szabo et al., 2003;Uhland, 2006). The oncogenic and metastatic activities of matriptase mainly come from its ability that can degrade extracellular matrix and activate various pro-oncogenic and pro-metastatic factors such as precursor of hepatocyte growth factor/scatter factor (pro-HGF/SF) (Goswami et al., 2014;Uhland, 2006). Thus, a possible method to prevent tumor development and metastasis is to use specific matriptase inhibitors.
Matriptase has been recognized for the participation in numerous physiological and pathological processes such as inflammation, tumor invasion, and metastasis (Antalis, Buzza, Hodge, Hooper, & Netzel-Arnett, 2010;Heutinck, ten Berge, Hack, Hamann, & Rowshani, 2010;Hooper, Clements, Quigley, & Antalis, 2001). Matriptase is overexpressed in many epithelial tumors and often correlated with clinical outcome. A modest overexpression of matriptase in epidermis cells would cause spontaneous squamous cell carcinoma and increase carcinogen susceptibility (List et al., 2005). Matriptase has played an important role in maintaining the epithelial integrity. The matriptase deletion disturbs epithelial development by rapid dehydration and impairs hair follicle formation, and a mutation in the coding region of matriptase would cause autosomal recessive ichthyosis with hypotrichosis (Basel-Vanagaite et al., 2007;Desilets et al., 2008). Moreover, matriptase often coexpressed with hepatocyte growth factor activator inhibitor-1 (HAI-1), which has been proposed to exert essential function in placentation and other embryonic processes (Lee, Kiyomiya, Benaud, Dickson, & Lin, 2005;Lin, Anders, Johnson, & Dickson, 1999;Szabo, Molinolo, List, & Bugge, 2007). The unbalanced ratio between matriptase and HAI-1 result in many cancers such as breast cancer (Kang et al., 2003) and prostate cancer (Saleem et al., 2006). In summary, it is helpful to consider that matriptase may be a central regulator of cell immigration and cancer invasion. Consequently, this protease may appear as a potential target to develop anti-cancer drugs.
In the last years, many promising anti-cancer drugs have been developed, such as peptide inhibitors and small molecule inhibitors (Colombo et al., 2012;Enyedy et al., 2001;Galkin et al., 2004;Goswami et al., 2013Goswami et al., , 2014Gray et al., 2014;Jiang et al., 2007;Li et al., 2007;Long et al., 2001;Quimbar et al., 2013;Steinmetzer et al., 2006). In 2006, Steinmetzer et al. identified a series of amidinophenylalanines as potent and selective small molecule matriptase inhibitors (Steinmetzer et al., 2006). In 2013, Goswami et al. also designed a series of sulfonamide and amide derivatives of pyridyl bis(oxy)benzamidine as matriptase inhibitors (Goswami et al., 2013). However, the analysis of binding modes for inhibitors cannot fully provide insight into the structural features in the active site of matriptase based on available experimental data. In the present work, four inhibitors which are from two class inhibitors with different scaffolds were selected to investigate atomistic details binding to matriptase. Our study is aimed to compare the different binding modes of the two class inhibitors binding to the same protein, as well as to find favorable interaction patterns using molecular dynamic (MD) method and molecular mechanics generalized born surface area (MM-GBSA) theory. This work provides detailed atomistic insight into the understood structurefunction relationships of matriptase-inhibitor complexes and provides meaningful information for the further design of novel, potent, and selective matriptase inhibitors.

System setups
The crystal structures of matriptase with four small molecule inhibitors were obtained from Protein Data Bank (PDB). The PDB entries used in this study are 2GV7 (matriptase-CJ-672) (Steinmetzer et al., 2006), 2GV6 (matriptase-CJ-730) (Steinmetzer et al., 2006), 4JYT (matriptase-N4A) (Goswami et al., 2013), and 4JZ1 (matriptase-F4D) (Goswami et al., 2013). The amidinophenylalanine inhibitors include CJ-672 and CJ-730, and the pyridyl bis (oxy)-benzamidine inhibitors are N4A and F4D, respectively. The structures of four inhibitors are shown in Figure 1. These four structures are used as the initial models for MD simulation. All of the missing atoms were added to complexes with the t-Leap procedure in Amber 11 software package (Case et al., 2010). H++ server (Gordon et al., 2005) has been used to determine the protonation states at pH 7.0, which can predict the pK a values of protein residues at a given pH. By calculating the pK a value of residues, His57 was assigned to be fully protonated at both nitrogen atoms, and the other ionizable residues to be their default protonation states at a neutral pH for four systems. The ff99SB force field was applied to produce the force field parameters for the protein and crystal water molecules (Hornak et al., 2006). The general Amber force field (Wang, Wolf, Caldwell, Kollman, & Case, 2004) was used to establish the potentials for the four inhibitors. Appropriate sodium ions were added to keep the whole systems neutral. The atom charges and atom types were assigned by the Antechamber tools in the Amber software package. Finally, TIP3P explicit water boxes with a 10.0 Å distance around the solute were added to the four systems.

MD simulations
Energy minimizations and MD simulations were carried out using the SANDER module of AMBER 11 package (Pearlman et al., 1995). Energy minimizations were carried out for two steps. Initially, the solvent and ions were subjected to 2500 steps of steepest decent minimization followed by 2500 steps of conjugate gradient minimization with the protein and ligands fixed with a 500 kcal/ (mol Å 2 ) constraint. Then, each system was totally minimized for another 5000 steps with no restraint (2500 steps of steepest decent minimization and 2500 steps of conjugate gradient minimization). After minimization, the system restrained with a force of 10 kcal/(mol Å 2 ) in the Cα atoms was gradually heated from 0 to 300 K during 150 ps in the canonical (NVT) ensemble. Then, each system went through 500 ps equilibrium at an NPT ensemble conditions (300 K, 1 atm). Finally, 50 ns MD simulations were carried out for each system with the isothermal isobaric (NPT) ensemble. Particle Mesh Ewald (PME) was employed to treat the long-range electrostatics interactions (Darden, York, & Pedersen, 1993). The short-range interactions were set to 10.0 Å. The time step was set to 2 fs for the MD simulations and PME method. SHAKE algorithm was used for constraining bonds involving hydrogen atoms (Kräutler, van Gunsteren, & Hünenberger, 2001). A hydrogen bond is defined by the distance between the heavy atoms using a cut-off of 3.5 Å and the angle between the donor and acceptor atoms using a cut-off 120 0 . The distance cut-off for ionic bond is set to 4.5 Å. MD simulations of 50 ns were performed for the four systems, respectively. The cluster analysis of MD trajectory was carried out to group together coordinate snapshots from the trajectory into distinct sets using average linkage as the clustering algorithm and using the root-mean-square deviation (RMSD) of backbone atom as the distance metric (Shao, Tanner, Thompson, & Cheatham, 2007). In this contribution, module ptraj of AMBER 11 (Case et al., 2010) was used for the cluster analysis of the four trajectories. The representative structures of the clusters with highest occurrences for matriptase-CJ-672, matriptase-CJ-730, matriptase-N4A, and matriptase-F4D systems were employed to depict structural information.

MM-GB/SA calculations
MM-GB/SA methods, implementing in AMBER11 (Case et al., 2010), were applied to calculate the binding free energies of the four matriptase-inhibitor systems. A total of 1000 snapshots, which were selected at an interval of 10 ps from the last 10 ns stable MD trajectories, were used to calculate the binding free energy. The free energy was calculated by the following equations: G complex , G protein , and G inhibitor are the free energy of complex, protein, and small molecule inhibitors, respectively. ΔG bind is written as the sum of molecule mechanics energy. In Equation (2), where E MM is the molecule mechanics component in gas phase, G sol stands for the salvation free energy, and TS is a vibrational entropy term. E MM is given as a sum of E int , E ele , and E vdw , which represent the internal, coulomb, and Van der Waals interaction term, respectively. Similarly, the solvation free energy (G sol ) can be divided into electrostatic solvation free energy (G GB ) and nonpolar solvation free energy (G SA ). G GB can be computed using Generalized Born (GB) Method (Hawkins, Cramer, & Truhlar, 1996), while G SA is considered to be proportional to the molecular solvent accessible surface area (Hou, Li, Li, & Wang, 2012;Weiser, Shenkin, & Still, 1999).

Calculation of the inhibitor-residue interaction
Per-residue binding free-energy decomposition was performed using the MM-GB/SA method in Amber 11 (Case et al., 2010). It can be given as the sum of Van der Waals contribution (ΔE vdW ), electrostatic contribution (ΔE ele ), and solvation contribution (ΔG ele, sol ). The solvation contribution can also be divided into the polar (ΔG ele, sol ) and the nonpolar (ΔG nonpol,sol ) parts. 3. Results and discussion 3.1. Stability during MD simulations MD simulations of 50 ns were performed for each system in the explicit solvent. In order to evaluate the stability of the MD trajectories, the RMSD of Cα atoms with respect to the initial frame has been monitored. As shown in Figure 2, each system gradually reached the equilibrium and remained quite stable after 25 ns, where RMSD value of each system converged to an average value around 1.0 Å. The result suggests that four complexes have stable protein structures and it is realizable for the further analysis.
To measure the fluctuations of a certain residue over the simulation time, the root-mean-square-fluctuation (RMSF) of Cα atoms in matriptase was calculated. The result is shown in Figure 3. It is easy to notice that the RMSF values of most residues are less than 1 Å. The residues near the active site (including His57, Phe99, Asp189, Ser190, Trp215, Gly216, and Gly219) show low RMSF values in the four complexes. It indicates conformation of matriptase active site hardly ever change upon ligand binding. In the meantime, there are four residues showing high values, including Arg60c, Val129, Gln165, and Arg179. These four residues show certain flexibility, which may be due to the existing of the long side chain. On the above structural analysis, the global protein structure, especially the active site, is reserved.

Clustering analysis
Clustering analyses for the four systems were carried out by the ptraj module of AMBER11. Each trajectory was divided into five clusters (C0, C1, C2, C3, and C4). In matriptase-CJ-672 and matriptase-CJ-730 systems, the most popular clusters are C0 and C1, which contain 30.5 and 54.5% of the total population, respectively. In matriptase-N4A and matriptase-F4D systems, the most popular clusters are C3 (70.5%), and C2 (54.2%), respectively. Thus, the most representative structure of cluster for each complex is chosen to investigate the structural characteristics. Partial residues interacting with four ligands are shown in Figure 4. It is clear that both CJ-672 and CJ-730 adopt stable chair conformations (Figure 4 A and B). The proximal piperazine (CJ-672) or piperidinyl (CJ-730) ring can fold an integral moiety with the sulfonylated aryl ring. The moiety is found to face to a hydrophobic surface cleft in space formed by Phe97, Asp98, Phe99, Gln175, and Trp215. Moreover, the distal positively charged amidino group (CJ-672) and guanidinoethyl group (CJ-730) are surrounded by a negatively charged "cation cleft" (Figure 5). Due to the complementarity of charge, both CJ-672 and CJ-730 can be stabilized well in the "cation cleft." These structural features mentioned above are consistent with the previous study of matriptase complexes (Steinmetzer et al., 2006).
As shown in Figure 4(C) and (D), the benzamidines of N4A and F4D extend into the S2/S4 hydrophobic cleft. At another ends of N4A and F4D, both the naphthyl of N4A and the fluoro-phenyl moiety of F4D orient toward 60-loop, which do not make contact with the "cation cleft" due to the two groups with electrically neutral property. Obviously, the distal moieties of the two class inhibitors extend into the different site of S1' region. It is clear that the two class inhibitors exist different binding modes interacting with the same protein. Furthermore, the benzyl of Phe99 goes through an obvious rotation in different degrees for each complex. It may be  the reason that the extended S2/S4 cleft is not a rigid structure. The Phe99 benzyl group drifts away from its original site to enlarge the hydrophobic cleft to avoid the steric hindrance, when the groups of small molecules approach the S2/S4 cleft.

Hydrogen bond analysis
Hydrogen bonds play an essential role in stabilizing protein-ligand complexes. The key residues participating in the hydrogen bond interactions are listed in Supplementary Table S1 and displayed in Figure 6. It is clear that Gly219 forms hydrogen bonds with all ligands. For matriptase-CJ-672 and matriptase-CJ-730 systems, four residues (Gly219, Gly216, Ser190, and Asp189) have important contributions to inhibitors through the formation of hydrogen bonds in high frames. Thereinto, Gly219, Gly216, and Ser190 can form the hydrogen bonds with the two inhibitors, which has been proven by experiment (Steinmetzer et al., 2006). CJ-672 and CJ-730 can be stabilized well by the formation of stable hydrogen bonds between inhibitors and two residues (Gly219 and Gly216) in high hydrogen bond occupancies in the process of MD simulations. The side chain hydroxyl of Ser190 forms hydrogen bonds with the terminal benzamidines of the two inhibitors in approximately 55.12 and 59.38%, respectively. In addition, we have observed the terminal benzamidines of CJ-672 and CJ-730 can form hydrogen bonds with the carboxylate oxygens of Asp189 at the bottom of S1 pocket. It can be considered that the two inhibitors can be firmly fixed at the S1 pocket by the hydrogen bonds formed by the inhibitors and four residues. For matriptase-N4A system, His57 and Gly193 form hydrogen bonds with N4A in high hydrogen bond  occupancies over 99% in our MD simulation. It is obvious that the two hydrogen bonds play a significant role in stabilizing N4A. Moreover, Gly219, Ser195, Ser190, and Asp189 also can form hydrogen bonds with N4A in approximately 55.04, 60.66, 31.58, and 39.82% frames. For matriptase-F4D system, there is a little difference from N4A. Because of the inward deflection of fluorophenyl moiety of F4D, it is found to form a hydrogen bond with Gly219 in 96.08% occupancy, which makes major contributions to stabilize F4D. In addition, Gly193 and His57 can form hydrogen bonds with F4D in 61.82 and 35.10% frames, respectively. It can be considered that the N4A and F4D may be fixed by hydrogen bonds at both ends of N4A and F4D.
The above analyses elucidate that the two class inhibitors exist different binding modes, the hydrogen bonds have significant effects on the binding modes of matriptase.

Ionic bonds
The ionic bond analysis was carried out to investigate the interaction characteristics between protein and ligands ( Table 1). The carboxylate oxygens of Asp189 can form two stable ionic bonds with the two amidine nitrogens of four inhibitors with high occupancy (>75%), respectively. The result is in accordance with the experimental data (Goswami et al., 2013;Steinmetzer et al., 2006). Therefore, Asp189 is a crucial residue which is responsible for constructing hydrogen bond and ionic bond interactions between the protein and inhibitors. Moreover, we have observed that CJ-672 can form two additional ionic bonds with Asp60b in the 60-loop region, which is not observed in crystal structure. The two ionic bonds can play an important role in stabilizing the matriptase-CJ-672 system. From another side of view, we consider that the length of the chain which extended into the S1' region may have influence on the contact with 60-loop region. These nonbonding interactions, including hydrophobic, hydrogen bond, ionic bond interactions, are crucial for the stability of the four complexes.

Binding free energy calculations
To investigate the binding free energies between the protein and ligands, the MM-GB/SA calculations were performed. Using the competitive inhibition equation, K i = IC 50 /(1 + [S]/K m ), the K i value of the four inhibitors to matriptase were obtained by Steinmetzer (Steinmetzer et al., 2006) and Goswami (Goswami et al., 2013). Based on K i value from experiment, the ΔG exp was determined. All of energy terms and the entropy changes are listed in Table 2. The binding free energies of CJ-672, CJ-730, N4A, and F4D systems are −37.21, Figure 6. Percentage occupancy of hydrogen bonds in different systems. All these bonds are monitored during the last 10 ns of each system. −34.46, −34.62, and −28.63 kcal/mol, respectively. It is encouraging to find that the ranking of predicted binding free energies are in good agreement with that of the experimental data (Goswami et al., 2013;Steinmetzer et al., 2006). The Van der Waals (ΔE vdW ), electrostatic (ΔE ele ), and nonpolar solvation terms (ΔG SA ) are favorable for four ligands binding to the matriptase. Thereinto, the Van der Waals energy is found to be the main driving force of binding affinity. Overall, the results indicate amidinophenylalanine inhibitors have slightly higher ability to bind to the matriptase than pyridyl bis (oxy)benzamidine inhibitors. It has shown the highest value of Van der Waals in matriptase-CJ-672 system and the highest value of electrostatic interactions in matriptase-CJ-730 system. In the case of matriptase-F4D system, the value of electrostatic binding free energy is almost half as that of matriptase-CJ-730 system. Moreover, it can be noticed that matriptase-N4A and matriptase-F4D systems lose more entropy during two inhibitors binding to matriptase than the other two systems. It may be the reason that CJ-672 and CJ-730 can adopt stable chair conformations and have obvious favoring structures in the binding process of protein-ligand.

Decomposition analysis of the binding free energy
The per-residue binding free energy decomposition analysis was explored to gain more-detailed insight into the effects of individual residues on the inhibitor bindings. Lower than −1 kcal/mol interaction energies are regarded as essential for ligand bindings. The results are shown in Supplementary Table S2 and displayed in Figure 7. For matriptase-CJ-672 and matriptase-CJ-730 complexes, the analyses of interaction energies highlight nine residues within the active site: His57, Asp60b, Asp96, Phe99, Asp189, Ser190, Trp215, Gly216, and Gly219. Thereinto, Asp189, Ser190, Gly216, and Gly219 can firmly fix the inhibitors at the S1 pocket by hydrogen bond interactions. Especially, Gly216 and Gly219 make outstanding contributions to stabilize the CJ-672 (−2.64 and −2.53 kcal/mol) and CJ-730 (−3.21 and −1.86 kcal/mol), respectively. Moreover, the distal positively charged amidino group (CJ-672) and guanidinoethyl group  can be stabilized by a negatively charged "cation cleft," which is consistent with the previous data (Steinmetzer et al., 2006). According to the binding free energy decomposition, the electrostatic energies of Asp60b and Asp96 for CJ-672 and CJ-730 are −5.57 and −1.49 kcal/mol, −2.52 and −2.27 kcal/mol, respectively. Though the total binding free energies of both Asp60b and Asp96 make no contributions to the binding process due to the solvation effects, they still play significant role in stabilizing the inhibitors in favorable positions. The most significant energetic contribution is found for the residue Trp215, whose total energetic contribution to CJ-672 and CJ-730 are −4.50 and −4.29 kcal/mol, respectively. As can be observed, both CJ-672 and CJ-730 are surrounded by the hydrophobic residues, particularly Phe97, Phe99, Trp215, and His57, collectively form a hydrophobic cleft and provide hydrophobic interaction to stabilize the inhibitors. Phe99 is observed to adopt a suitable position to accommodate the inhibitors by a rotation of its side chain and subsequently has indirect impact on ligand binding. For matriptase-N4A and matriptase-F4D complexes, the key residues that contribute to the binding of the ligands: His57, Phe99, Asp189, Ser190, Gly193, Gly195, Trp215, Gly216, and Gly219. The hydrophobic residues including Trp215 and Phe99 have important contributions to binding between matriptase and N4A and F4D because of the low interaction energies. The benzamidines of N4A and F4D can be fixed in the hydrophobic cleft. Gly219 and Asp189 have played significant roles in stabilizing the inhibitors in a favorable position by hydrogen bond and ionic bond interactions at the S1 pocket. Due to contribution of electrostatic energy made by His57, Gly193, and Gly195, the naphthyl of N4A and the fluoro-phenyl moiety of F4D, which orients toward 60-loop, are fixed to be packing against the Cys42-Cys58 disulfide bridge and Ile41 side chain. This structural feature that we have observed is in concert with the previous study (Steinmetzer et al., 2006). The Van der Waals energies of Ile41 are −1.64 kcal/mol for N4A, and −1.06 kcal/mol for F4D, respectively. It can be considered that the distal moieties of N4A and F4D can be stabilized well in the S1' region. In summary, the two class inhibitors exist different binding modes due to the different structural features. For the four inhibitors, the crucial residues with the hydrophobic side chains, especially Trp215 and Phe99, may have great impact on the inhibitors binding.

Conclusion
The application of MD simulations has become indispensable in computational area of particular emphasis with regard to the protein-ligand interactions Wang et al., 2013). In this work, an integrated computational approach by combining MD simulations, MM/GBSA binding free energy calculations, and MM/ GBSA binding energy decomposition analysis was used to find out the binding modes and the interaction mechanisms of matriptase with small molecule inhibitors. In the results of MM-GB/SA, the contribution of Van der Waals interactions is the main force to give impetus for the binding. Our results elucidate the two class inhibitors exist different binding modes due to the different structural features. Through summarizing the two different modes, we have mastered some important and favorable interaction patterns between matriptase and inhibitors: Firstly, the electrostatic interactions between the distal group of inhibitor and the negatively charged "cation cleft," such as amidino, guanidino, and so on. Secondly, the hydrophobic interactions between inhibitor and Trp215 and Phe99. Thirdly, we have paid close attention to the position that should be responsible for the formation of hydrogen bonds, one position is at the bottom of S1 pocket, the other is in the site of the sulfonyl group. Accordingly, for the rational design of more potent inhibitors of matriptase, the sulfonyl group is essential. In addition, several key residues responsible for achieving strong binding have been identified. Especially, Trp215 and Phe99 have an important impact on active site architecture and ligand binding. It is actually significant to consider these favorable protein-ligand interaction patterns for designing more potent inhibitors of matriptase. On basis of a series of analysis, we can obtain the important binding sites of protein and the potent functional groups of inhibitors at the atomic level. This work will provide a better structural understanding of small molecules targeting the matriptase, the information will be valuable not only for obtaining the binding mode and inhibition mechanism but also for designing more novel and potent matriptase inhibitors.

Supplementary material
The supplementary material for this paper is available online at http://dx.doi.org/10.1080/07391102.2016. 1160259.

Disclosure statement
No potential conflict of interest was reported by the authors.