Exploration of the Substrate Preference of Lysine Methyltransferase SMYD3 by Molecular Dynamics Simulations.

SMYD3, a SET and MYND domain containing lysine methyltransferase, catalyzes the transfer of the methyl group from a methyl donor onto the Nε group of a lysine residue in the substrate protein. Methylation of MAP3 kinase kinase (MAP3K2) by SMYD3 has been implicated in Ras-driven tumorigenesis. The crystal structure of SMYD3 in complex with MAP3K2 peptide reveals a shallow hydrophobic pocket (P-2), which accommodates the binding of a phenylalanine residue at the -2 position of the substrate (F258) is a crucial determinant of substrate specificity of SMYD3. To better understand the substrate preference of SMYD3 at the -2 position, molecular dynamics (MD) simulations and the MM/GBSA method were performed on the crystal structure of SMYD3-MAP3K2 complex (PDB: 5EX0) after substitution of F258 residue of MAP3K2 to each of the other 19 natural residues, respectively. Binding free energy calculations reveal that the P-2 pocket prefers an aromatic hydrophobic group and none of the substitutions behave better than the wild-type phenylalanine residue does. Furthermore, we investigated the structure-activity relationships (SAR) of a series of non-natural phenylalanine derivative substitutions at the -2 position and found that quite a few modifications on the sidechain of F258 residue could strengthen its binding to the P-2 pocket of SMYD3. These explorations provide insights into developing novel SMYD3 inhibitors with high potency and high selectivity against MAP3K2 and cancer.

Introduction

The lysine methylation is a principal regulatory mechanism that influences the protein activity, stability, and function.[1] SMYDs (SET and MYND domain containing proteins) catalyze the transfer of the methyl group from S-adenosyl-l-methionine (SAM) onto the Nε group of a lysine residue in the target protein substrate. Both histones and nonhistone proteins have been shown to be the substrates of SMYD family proteins. Their reported histone targets include H3K4 (SMYD1 and SMYD3), H3K36 (SMYD2), and H4K5 (SMYD3),[2−5] functioning in the regulation of chromatin structure and gene expression.[6−11] Nonhistone substrates can also be methylated by SMYD proteins (Figure A). SMYD2 has a broad spectrum of substrates including tumor suppressor protein p53 (K370),[12] retinoblastoma (Rb, K860),[13] estrogen receptor α (ERα, K266),[14] heat shock protein 90 (HSP90, K615),[15] and period circadian protein homolog 2 (PER2, K798).[16] In the case of SMYD3, vascular endothelial growth factor receptor 1 (VEGFR1, K831) and MAP3 kinase kinase (MAP3K2, K260) have been reported to be the substrates of SMYD3.[17,18]

Figure 1

Comparison between the substrates of SMYD3 and SMYD2. (A) Amino acid sequence alignment of the reported substrates of SMYD3 and SMYD2. The −2- and 0-position residues in substrates of SMYD3 and the −1- and 0-position residues in substrates of SMYD2 are highlighted in the red boxes, respectively. (B) Electrostatic potential and substrate-binding cleft on the surface of SMYD3. The SET, MYND, post-SET, and C-terminal domains of SMYD3 are shown in green, blue, yellow, and pink, respectively. SAH and the MAP3K2 peptide are shown in a stick model (carbon, gray, and cyan, respectively; nitrogen, blue; oxygen, red). Zinc ions are shown as gray spheres. The −2-position F258 residue in MAP3K2 peptide is highlighted in the red circle. (C) Structural alignment of the shallow hydrophobic pocket of SMYD3 and SMYD2 with their substrates MAP3K2 and p53 peptide, respectively. The residues at the −2 position in SMYD3, at the −1 position in SMYD2, in MAP3K2 peptide, and in p53 peptide are shown in a stick model with carbon in green, orange, cyan, and yellow, respectively.

SMYD3 has been implicated in a variety of cancers including liver, colon, and breast cancer related to its methyltransferase activity on nonhistone targets.[19−22] Methylation of VEGFR1 K831 by SMYD3 enhances its autophosphorylation and kinase activity in the cell and promotes carcinogenesis.[23−25] Methylation of MAP3K2 at K260 in the cytoplasm activates the MAP kinase signaling module and promotes RAS-driven tumorigenesis.[18,26] These results suggest that targeting the methyltransferase activity of SMYD3 would be a useful strategy in anticancer therapy.

Our previous work solved the crystal structure of SMYD3 in complex with MAP3K2 peptide.[27] The structure reveals that there is an amphiphilic cleft for substrate binding on the surface of SMYD3. A shallow hydrophobic pocket (P-2), which accommodates the binding of a phenylalanine residue at the −2 position of the substrate (F258) is a crucial determinant of the substrate specificity of SMYD3 (Figure B). This P-2 pocket is composed of several hydrophobic residues including L104, V178, I179, and V195, with two serine residues, S101 and S182, standing at the perimeters of the pocket (Figure C). SMYD2 contains a similar shallow hydrophobic pocket to hold the leucine residue at −1 position of the substrate p53, with L108, V179, N180, S196, T105, and G183 taking the place of L104, V179, I179, V195, S101, and S182 of SMYD3, respectively (PDB IDs: 3TG5 and 5EX0, Figure C). It is interesting to find that a majority of the reported substrates of SMYD2 include a leucine residue at the −1 position[12−16] (Figure A). And the importance of this hydrophobic pocket in substrate binding preference of SMYD2 has been proved by several experimental and computational studies.[16,28−32]

In the past decade, several computational protein design protocols have been reported to predict the specificity of protein–protein or protein–ligand interactions.[33−36] Coupled protein docking and conformational ensemble methods[37−39] were further developed to enhance the accuracy of prediction. For example, a multistate computational procedure combining backbone ensemble, energy minimization, amino acid substitution, and fitness calculation was developed to probe the substrate specificity of SMYD2.[32] The result reveals that SMYD2 exhibits narrow specificity for substrates at the −1 and +2 positions. Besides these approaches, molecular dynamics (MD) simulations are valuable adjuncts to reveal substrate recognitions by conformational dynamics. The motions of substrate and protein induced by their mutual interactions could be directly observed through trajectories.[40−42]

To better understand the substrate preference of SMYD3 at the −2 position, different kinds of natural amino acid substitutions are first used as probes to explore the binding preference of the P-2 pocket in this work. Binding free energy calculations of MD simulations from the MM/GBSA method reveal that the P-2 pocket plays a crucial role in substrate binding with narrow specificity. It prefers an aromatic hydrophobic group and none of the substitutions behave better than the wild-type phenylalanine does. Furthermore, a series of non-natural phenylalanine derivative substitutions at the −2 position was investigated by similar MD simulations procedures. Results show that quite a few modifications on the sidechain of F258 residue could further stabilize the sidechain conformation, and thus strengthen its binding to the P-2 pocket of SMYD3. These results provide insights into developing novel SMYD3 inhibitors with high potency and high selectivity against MAP3K2 and cancer.

Results and Discussion

First, we performed a 100-ns MD simulation to observe the stability of SMYD3-MAP3K2 complex. SMYD3 is composed of three domains, an N-terminal catalytic domain (SET and post-SET), a MYND domain inserts into the SET domain, and a C-terminal domain (CTD) (Figure B). The crystal structure of SMYD3-MAP3K2 complex reveals that the SET and post-SET domains of SMYD3 contribute most of the binding force to substrate binding other than the MYND and CTD domains. root mean square fluctuations (RMSF) analyses of each residue of the complex after 100-ns MD simulation show that the MYND and CTD domains of SMYD3 are more flexible than the SET and post-SET domains (Figure A). In the case of MAP3K2 peptide, residues from −2 to +2 positions (F258 to G262) are more stable than the remaining N- and C-terminal residues (Figure B). Figure C exhibits time evolutions of root mean square deviations (RMSD) of the stable regions of SMYD3 and MAP3K2. The RMSD values of both SET and post-SET domains of SMYD3 and F258 to G262 residues of MAP3K2 reach a stable plateau after ∼10 ns of simulations and keep a position around an average value of 2.3 and 1.0 Å during the remaining MD simulations, respectively. These results indicate that a 50-ns MD simulation unit performed on the stable regions of SMYD3 and MAP3K2 in the following study is sufficient for detecting reliable binding free energies and conformational dynamics of the complex.

Figure 2

Structural deviation and flexibility in the 100-ns MD simulation. (A) RMSF as a function of the residue index of SMYD3. The SET, MYND, post-SET, and CTD domains are shaded in green, blue, yellow, and pink, respectively. (B) RMSF as a function of the residue index of MAP3K2 peptide. (C) Time evolutions of RMSD of the SET and post-SET domains in SMYD3 and the residues F258–G262 in MAP3K2 peptide.

Conformational Dynamics of F258 and Substitutions at the −2 Position

Binding preference at the −2 position was detected by the conformational dynamics of the F258 residue and each of its substitutions. A 50-ns MD simulation was performed three times. Averaged binding free energy of SMYD3 to F258 residue or its substitution from the last 2-ns was calculated, as shown in Figure A and Table , and time evolutions of its binding free energy in each repetition is depicted in Figure S1. The binding free energy tends to convergence, which means the complex has reached a relatively stable state during MD simulations. The results show that wild-type F258 of MAP3K2 has a most negative binding free energy of −9.64 kcal/mol toward SMYD3, indicating the most favorable substrate–protein recognition. Replacement of F258 to a leucine or a methionine residue results in an increased binding free energy of −5.46 and −5.31 kcal/mol, respectively. Other hydrophobic residue substitutions like I, V, C, and Y show tolerable but much weaker binding energy toward SMYD3. F258Y introduces an additional hydroxyl group to the phenyl ring, which lowers its binding affinity to the P-2 hydrophobic pocket. While in the case of the remaining F258 substitutions, the calculated binding free energy may be too high for efficient recognition. These results are in consistence with our previous experimental data[27] and similar to the substrate specificity of SMYD2 at −1 position proved by computational and experimental approaches.[41]

Figure 3

Conformational dynamics of F258 and substitutions at the −2 position. (A) The binding free energies between SMYD3 and F258 or its substitutions. The averaged value and standard deviation of three repetitions in each system are depicted in the histogram, respectively. The same scheme is used in the following figures unless explicitly specified. (B) The initial conformation in the crystal structure. F258 in MAP3K2 peptide and the residues around F258 in SMYD3 are shown in a stick model. The carbon atoms are shown in magenta and green, respectively. The main chain of the residue is hidden for brevity if not necessary. The same scheme is used in the following figures unless explicitly specified. The distances between F258 and the marked residues in SMYD3 are shown in the color map. (C) Conformational alignment of F258 after MD simulation with carbon in magenta and F258 in the crystal structure with carbon in cyan. (D–H) The conformations of other representative substitutions after MD simulations.

Table 1

Averaged Binding Free Energies and Standard Deviations between SMYD3 and F258 or Its Substitutions from Three Times-MD Simulations, Respectivelya

			repetition
residue type	average	std. dev.	1	2	3
PHE	–9.64	0.13	–9.64	–9.52	–9.77
LEU	–5.46	0.95	–5.47	–6.41	–4.50
MET	–5.31	1.19	–4.44	–6.66	–4.82
ILE	–4.32	0.98	–5.42	–3.52	–4.02
VAL	–3.99	0.32	–3.95	–3.69	–4.33
CYS	–3.89	1.50	–3.49	–2.62	–5.55
TYR	–3.88	0.74	–3.64	–3.28	–4.71
PRO	–2.46	0.80	–2.40	–1.70	–3.29
THR	–2.14	1.10	–2.13	–3.25	–1.06
HIS	–1.92	0.75	–2.74	–1.26	–1.77
ALA	–1.74	0.87	–2.46	–1.98	–0.77
TRP	–0.85	1.13	0.44	–1.42	–1.58
ASN	–0.04	1.41	0.70	–1.66	0.85
GLY	–0.03	0.42	–0.01	–0.46	0.38
GLN	0.03	1.18	1.05	0.32	–1.27
SER	0.22	0.79	0.16	1.04	–0.54
ARG	0.65	0.84	–0.24	1.42	0.78
GLU	1.03	0.91	1.99	0.89	0.19
LYS	1.98	0.56	2.62	1.68	1.64
ASP	2.37	0.84	1.91	3.35	1.86

All binding free energies are in kcal/mol.

Besides binding free energies, binding modes at the P-2 pocket could be traced after MD simulations. The initial conformation of F258 in the crystal structure is shown in the left panel of Figure B. Distances of sidechain-centroids between F258 and each of the six residues forming the P-2 pocket of SMYD3 (numbered 1–6) are shown in a graduated color map at the right panel, from blue (closest) to red (farthest). Moreover, the averaged distance in the color map would help to estimate the occupancy rate of F258 in the P-2 pocket. Average conformations of F258 and its 19 replacements in the P-2 pocket during the last 2-ns MD simulations are shown in Figures C–H and S2. From the results, we can see that the conformation of F258 is rather stable, and its averaged distance to the P-2 pocket residues changed from 6.37 Å to a little bit more favorable 6.24 Å after simulations (Figure C). Among the 19 substitutions of F258, only leucine and methionine residues insert into and fully occupy the P-2 pocket like the phenylalanine residue does, with a slightly larger distances of 6.57 and 6.92 Å, respectively (Figure D,E). In the case of other hydrophobic substitutions like I, V, C, A, and W, their sidechains are not symmetrically presented (I), too short (V, C, A), or too large (W) to fully insert into the P-2 pocket (Figures F,G and S2). Other residues with hydrophilic or charged sidechains like arginine could not be stably located inside the pocket and the unfavorable bindings lead to their sidechains flipping out of the pocket during MD simulations (Figures H and S2).

The P-2 pocket of SMYD3 is a shallow hydrophobic pit mainly formed by L104, V178, I179, and V195; its binding preference is associated with the hydrophobic property and the sidechain volume of the interacting residue. From our calculation results, −2 position residue with aromatic hydrophobic sidechain shows a much negative binding free energy to the P-2 pocket than residues with small or hydrophilic sidechains. Meanwhile, nonpolar sidechains occupy the −2 pocket better and demonstrate a much favorable binding mode than polar or charged ones. These are consistent with the size and hydrophobic property of the P-2 pocket.

Our previous work showed that the catalytic activity of SMYD3 toward MAP3K2 peptide decreased ∼6-fold when F258 was mutated to a leucine residue, and no methylation activity was detected when an F258R mutation was introduced to the substrate.[27] Combined with binding free energy calculations and structural analyses, this work provides good explanations of the experimental results. Compared with the wild-type phenylalanine residue, a leucine residue has a 4.2 kcal/mol higher binding free energy to the P-2 pocket, indicating a lower binding capacity (−5.46 to −9.64 kcal/mol). MD simulations demonstrate that a leucine residue can occupy the P-2 pocket well but has a lower occupancy rate than phenylalanine (6.57–6.24 Å). These may be the major reasons causing a lower activity of SMYD3 on F258L mutant of the substrate. This can be further proved by the substitution of the −2 residue of VEGFR1, another reported substrate of SMYD3. VEGFR1 happens to contain a leucine residue at the −2 position (L829, Figure A), and the activity of SMYD3 on VEGFR1 is lower than that on MAP3K2.[27] When L829 is substituted by a phenylalanine residue, the binding free energy of F829 to the P-2 pocket decreases from −4.49 to −8.52 kcal/mol (Table S1), which confirms the importance of −2 phenylalanine to the activity of SMYD3.

To further prove that the binding affinity is related to the catalytic activity directly, in other words, is not affected by the change of the catalytic property of SMYD3 when catalyzing on mutated substrates, we first measured Kcat value of SMYD3 on wild-type MAP3K2 peptide and F258L mutant by in vitro methyltransferase assay. Results show similar Kcat values of SMYD3 on wild-type MAP3K2 (5.51 ± 0.85 h–1) compared with that on F258L mutant (5.13 ± 0.58 h–1), indicating a constant catalytic property of SMYD3 (Table S2). In addition, we demonstrate computationally that the active site structures remain the same for different mutants during MD simulations by calculating the averaged RMSD value of the key residues both in the active site and in the P-2 pocket of SMYD3 (Table S3). Results show small RMSDs (from 0.58 to 1.21 Å) with reasonable standard deviations (0.08–0.36 Å), indicating that the structures of the active site of SMYD3 remain the same for all of the 19 mutated substrates during MD simulations.

In the case of F258R mutation, the positively charged large sidechain of an arginine residue can hardly fit into the P-2 pocket after the MD simulations, and binding free energies rose to an even positive charged value with big deviations among paralleled simulations (Table ). The failure of recognition at the −2 position may disturb substrate binding, hindering the insertion of the target lysine residue into the catalytic pocket. The improper distance or uncomfortable circumstance between the Nε group of lysine and the methyl group on SAM will prevent an efficient transfer of the methyl group.

SARs of Phenylalanine Derivatives at the −2 Position

Our research reveals that the P-2 pocket has a high binding preference for phenylalanine residue. The structure–activity relationship (SAR) studies of modifications on the −2 phenylalanine would be important for future inhibitor design targeting the P-2 pocket. For this propose, a series of substitutions of −2 residue by non-natural phenylalanine derivatives was investigated by similar MD simulation procedures. A total of 51 phenylalanine derivatives or mimics was chosen from the SwissSidechain database,[43] which provides a unified web resource for hundreds of non-natural sidechains with biochemical parameters. Time evolutions of binding free energy (Figure S3) show that each complex has reached a relatively stable state during MD simulations. Binding free energy calculations show that quite a few modifications on the phenyl ring and two cycloalkyl mimics have more negative values than phenylalanine, indicating a more favorable binding to the P-2 pocket (Figure A, Table S4).

Figure 4

Conformational dynamics of F258 modifications at the −2 position. (A) The binding free energies between SMYD3 and F258 modifications. (B–P) The conformations of the representative F258 modifications. The residue labels of SMYD3 are shown in (B) and omitted in the other figures for brevity. Hydrogen bond is depicted as dash line if it exists.

Among these beneficial modifications, halogen substitutions prove to be good candidates. And the less hydrophilic substitution results in a more favorable recognition as exampled by monosubstitution at the para position of the phenyl ring by different halogens. Binding free energy results show that substitutions by I (PHI, −10.53 kcal/mol), Br (4BF, −10.49 kcal/mol), and Cl (F20, −10.04 kcal/mol) are more suitable to the pocket than wild-type F258 (−9.64 kcal/mol), but F (PFF, −9.51 kcal/mol) is not a good candidate of substitution (Table , Figure B–E). These results correspond with the hydrophobic property of the P-2 pocket. Electronegativity of iodine is very close to that of carbon, and thus C–I bond is less polarized than other carbon–halogen bonds, which makes iodine substitution the most favorable. On the contrary, the electronegativity of fluorine is much higher than that of carbon, so the C–F bond is strongly polarized, which results in the unfavorable binding of fluorine substitution. Furthermore, monosubstitution by halogen at the para position of the phenyl ring seems to be more effective than that at the ortho or meta position (Table ), which may due to the spatial locations of residues forming the P-2 pocket. As shown in Figure B, the distance between −2-position F258 and L104 residue of SMYD3 leaves adequate space to allow suitable modifications at the para position on the phenyl ring. In contrast, the sidechain hydroxyl group of S182 at the periphery of the P-2 pocket is rather close (∼4 Å) to the ortho and meta sites on one side of the phenyl ring, and thus it may hinder further modifications at those positions, especially for big hydrophobic halogen atoms. However, monosubstitution by the hydroxyl group at the ortho position is better than that at the meta or para position for recognition because its hydrophilic property is more suitable for interacting with S182 residue (Table ). Double substitutions by halogens (except fluorine) at both para and meta sites result in a more favorable binding to the P-2 pocket as both of the halogen atoms insert into the pocket in good conformations (I34, B34, and CP3 in Table , Figure F–H). However, the para–ortho-site double modifications (CP2) is not as good as the para–meta one (CP3) because of the rather outer position and more hydrophilic property of the ortho site compared to the meta site in the pocket, though it is still better than a monosubstitution at the para site (F20) (Table , Figure I).

Table 2

Averaged Binding Free Energies between SMYD3 and the Representative Monosubstituted F258 Modificationsa

	binding free energy
substituent group	ortho	meta	para
–I	–9.69	–9.82	–10.53
–Br	–9.52	–9.66	–10.49
–Cl	–8.96	–9.46	–10.04
–F	–7.63	–8.65	–9.51
–CH3	–9.45	–9.88	–10.17
–OH	–8.91	–6.70	–3.88

All binding free energies are in kcal/mol.

Table 3

Averaged Binding Free Energies between SMYD3 and the Representative Multisubstituted F258 Modifications or Cycloalkyl Derivativesa

	substituent group/position
name	2	3	4	5	6	binding free energy
I34		–I	–I			–12.79
B34		–Br	–Br			–12.29
CP3		–Cl	–Cl			–12.13
CP2	–Cl		–Cl			–9.97
BOH	–OH	–Br	–Br			–12.14
BO2	–OH		–Br	–Br		–11.49
BMO	–Br	–OH	–Br			–10.40
DBY		–Br	–OH	–Br		–9.67
ALC	cyclohexyl sidechain					–12.64
CPA	cyclopentyl sidechain					–10.82

All binding free energies are in kcal/mol.

Substitution by the methyl group seems to be an alternative strategy that the mono-meta site (APD, −9.88 kcal/mol), mono-para site (4PH, −10.17 kcal/mol), and double-meta–para-site substitution (MP3, −9.80 kcal/mol) show a slightly lower binding free energy than F258 (−9.64 kcal/mol), respectively. A mono-ortho-site methyl modification is not as good as the para or meta substitution somehow (MPH, −9.45 kcal/mol) because the CH3 group at the ortho position may have space clashes with the aforementioned S182 residue of the P-2 pocket (Figures J and S4). Larger methyl group derivatives like the C(CH3)3 group behaves better than the methyl group at the para position (−10.46 to −10.17 kcal/mol) because of its greater hydrophobic property (Table S4, Figure K).

As one can predict, none of the hydroxyl group substitutions are suitable for the P-2 pocket, though the ortho substitution seems to be better than meta or para ones (Table ). However, triple substitution by the hydroxyl group and halogen seems to be a good option as exampled by hydroxyl and halogen combinations. Substitutions by two iodine or bromine atoms and one hydroxyl group exhibit good binding free energies to the P-2 pocket and in the meantime attenuate the strong hydrophobic property introduced by the two halogen atoms (Table ). Among these, BOH and BO2 with para-, meta-double bromine and ortho-hydroxyl group substitutions show pretty good binding free energy to the P-2 pocket as two bromine atoms insert into the pocket in good conformations, though BOH with a hydroxyl group positioned in the same side of the phenyl ring to bromine atoms behaves better than the opposite-side BO2 does (Figure L,M). If the hydroxyl group is positioned at the meta (BMO) or para (DBY) position, the binding free energy is higher and the hydrogen bond is formed with D100 instead of S182 residue, requiring a big conformational change in D100 (Table , Figures N and S4). Iodine may not be a good choice in these mixed triple substitutions because of its big size, though the binding free energy is rather outstanding (Table S4 and Figure S4). These findings throw light on future inhibitor designs by using substrate mimics at the −2 position, since the solubility of the small inhibitor can be increased by introducing a hydrophilic group but not at the expense of losing too much binding affinity.

Furthermore, substitutions of the whole phenyl ring by cycloalkyl mimics also show rather good results. This may be connected with the flexibility of naphthene that it could fit into the P-2 pocket better than a more rigid phenyl ring. Cyclohexyl ring is better than the cyclopentyl ring as it resembles the phenyl ring more in shape (Table , Figure O,P).

Conclusions

In this work, we use MD simulations and the MM/GBSA method to study the substrate preference of methyltransferase SMYD3. The P-2 pocket of SMYD3 is a well-defined shallow hydrophobic subpocket adjacent to the catalytic pocket. The −2-position F258 residue of MAP3K2 was substituted to each of the other 19 natural residues followed by a 50-ns MD simulation repeated three times, respectively. The energetic and structural results show that the P-2 pocket prefers an average-sized hydrophobic group. None of the substitutions behave better than the wild-type phenylalanine residue.

As the P-2 pocket plays a significant role in substrate binding, we investigated the SARs of a series of non-natural phenylalanine derivative substitutions at the −2 position. Results show that the P-2 pocket has adequate space to allow hydrophobic modifications at the para positions on the phenyl ring, especially for halogen atoms like bromine and chlorine. However, the hydrophilic S182 residue in the P-2 pocket confines the size and hydrophobicity of the substitution at the ortho or meta position of the ring. Multiple substitutions by hydroxyl group and halogen combinations seem to be a good guidance for future inhibitor design to achieve both solubility and affinity.

Recently, EPZ030456 was developed as an SMYD3 inhibitor with a nanomolar inhibition concentration.[44] However, the crystal structure shows that except for the headgroup that inserts into the lysine channel, the compound makes few specific interactions with SMYD3. To our knowledge, the P-2 pocket has not been shown as a target pocket in developing SMYD3 inhibitors. Investigations of binding preference of the P-2 pocket, as well as the SARs of modifications on the −2-position phenyl ring in this study, will provide insights into the rational drug design to develop novel SMYD3 inhibitors with high potency and selectivity.

Materials and Methods

Systems for MD Simulations

The previously solved crystal structure of SMYD3 and MAP3K2 complex (PDB: 5EX0) was used as the starting structure. F258 and G262 at the −2 and +2 positions of the MAP3K2 peptide were substituted to each of the other 19 natural residues, respectively. F258 was also substituted, respectively, to several non-natural residues chosen from the SwissSidechain database,[43] which provides a unified web resource for hundreds of non-natural sidechains with biochemical parameters. The selection of non-natural residue types follows the rules: (1) Sidechains were mainly chosen from the set of phenylalanine derivatives. (2) Sidechains with large substituted groups on the phenyl ring such as phenyl group, phosphate group, nitro group, cyan group, and guanidine group were not considered in case of steric conflicts with the P-2 pocket. (3) Sidechains with modifications on the CB atom such as hydroxyl- or carbonyl substitutions were not considered to simplify the system of comparison. (4) There is no special requirement for the position (ortho, meta, or para) or the number (mono or multi) of the substituted group. (5) Two alanine derivatives were selected with cycloalkyl sidechains to resemble the phenyl ring. Names of these non-natural residues were shown by three-letter abbreviations. (6) Eleven phenylalanine derivatives with halogen-substitutions, which are not included in the database, were also employed in the SARs study. They are I34, B34, BOH, IMO, IO2, BO2, BMO, MIF, MBF, OIF, and OBF. Their structures were modified from 4-bromo-phenylalanine (4BF) in the database by changing the substituted groups manually. Finally, 51 non-natural sidechains were investigated in MD simulations.

Parameters in the MD Simulations

The LEaP module of the AMBER18 package[45] was used for residue substitution and adding hydrogen atoms for each complex. The cofactor S-adenosyl-l-homocysteine (SAH) and three zinc ions in the crystal structure were reserved in the MD simulations. The chelating cysteine and histidine residues were deprotonated. The amber FF14SB force field[46] was used for protein and natural residues. The force field parameters of non-natural residues and SAH were derived from the SwissSidechain database[43] and the general Amber force field (GAFF),[47] respectively. TIP3P water model[48] was used to solvate the complex in a hexagonal explicit water box under the periodic boundary condition. The distance between the edges of the box and the closet atoms of the complex was 12 Å. Na+ was added as a counterion to neutralize each system.

For each solvated system, 5000-step energy minimization was performed, followed by a combined equilibration process with a 500-ps constant volume ensemble to heat the system from 0 to 300 K and a 500-ps constant pressure ensemble at a constant pressure of 1 bar. The Langevin thermostat[49] and Berendsen barostat[50] were used for temperature and pressure control, respectively. During equilibration, a force constant of 10 kcal mol–1 Å–2 as a harmonic constraint was applied. Then, 50-ns MD simulation of each system was performed in constant pressure ensembles at 300 K with the constraint released. The time step was set to 2 fs. The SHAKE algorithm[51] was used to restrain all of the bond lengths involving hydrogen atoms. The particle mesh Ewald (PME) method[52] was used to calculate the long-range electrostatic contributions. The cut-off value of the van der Waals interactions was set to 10 Å. Each system was performed three times under the same conditions to confirm the results. Besides, an extra 100-ns MD simulation of SMYD3 in complex with wild-type MAP3K2 peptide was performed to investigate their stability. The AMBER18 software package[45] was employed to perform MD simulations. Coordinates were saved every 10 ps.

Analysis of the Simulations

The interactions between SMYD3 and the focused residues were calculated by the MM/GBSA method, which combines molecular mechanics (MM), the Generalized Born (GB) equations, and surface accessible (SA) calculations to calculate the binding free energy.[53−55] The calculated binding free energy is given bywhere ΔGbind, vac is the free energy of the complex in the gas phase, including nonbonded and bonded potentials calculated by the MM method, and ΔGcom, sol, ΔGrec, sol, and ΔGlig, sol are the free energies of the complex, receptor, and ligand in the aqueous phase, respectively. The solvated free energies include the polar and nonpolar free energies, which are computed by the GB and SA methods, respectively. In this work, the binding free energy between SMYD3 and the residue in the peptide substrate was calculated by a total of 200 snapshots extracted from the final 2 ns trajectory of each system. All of the parameters were set as default values in the calculations. In each system, the residue in the MAP3K2 peptide to be calculated and the whole protein with SAH and zinc ions were considered as ligand and receptor, respectively. Since we are mainly interested in the differences of the binding free energies, entropy was ignored in the calculations.

The cpptraj module[56] of AmberTools18 was used to calculate the averaged structure, the averaged distance, root mean square fluctuations (RMSF), and root mean square deviations (RMSD). The illustrated representative conformation was averaged from the last 2-ns MD simulation and derived from the first repetition of each system. The distance was measured by sidechain-centroids between residues unless explicitly specified.