Effect of mutations on drug resistance of smoothened receptor toward inhibitors probed by molecular modeling

Effect of mutations on drug resistance of smoothened receptor toward inhibitors probed by molecular modeling

Chemical Physics Letters 741 (2020) 137126 Contents lists available at ScienceDirect Chemical Physics Letters journal homepage: www.elsevier.com/loc...

2MB Sizes 0 Downloads 19 Views

Chemical Physics Letters 741 (2020) 137126

Contents lists available at ScienceDirect

Chemical Physics Letters journal homepage: www.elsevier.com/locate/cplett

Research paper

Effect of mutations on drug resistance of smoothened receptor toward inhibitors probed by molecular modeling ⁎

Zhenmei Gaoa, Jingxiao Baoa, Shuhua Shib, Xiongwen Zhanga, Ya Gaoa,c, , Tong Zhua,d,

T



a

Shanghai Engineering Research Center of Molecular Therapeutics & New Drug Development, School of Chemistry and Molecular Engineering, East China Normal University, Shanghai, China b School of Science, Shandong Jianzhu University, Jinan, China c School of Mathematics, Physics and Statistics, Shanghai University of Engineering Science, Shanghai, China d NYU-ECNU Center for Computational Chemistry at NYU Shanghai, Shanghai, China

H I GH L IG H T S

drug resistance mechanisms of Vismodegib, LY2940680, and L4 were investigated. • The D473H mutation can weaken the binding of Vismodegib. • The • For LY2940680 and L4, hydrogen bond networks in binding pocket are preserved.

A B S T R A C T

The smoothened receptor (SMO) plays a critical role in embryonic development. Vismodegib was approved by the FDA to treat metastatic basal-cell carcinoma. However, The D473H mutation can cause obvious drug resistance toward Vismodegib. In this work, the drug-resistance mechanisms of three molecules were explored by theoretical calculations. The results indicate that the D473H mutation can disrupt the hydrogen-bonding network in the binding pocket, which can weaken the binding of Vismodegib. Whereas, for LY2940680 and L4 there are new hydrogen bond networks formed, which can stabilize their binding. This work provides a theoretical contribution to design more promising SMO inhibitors.

1. Introduction Globally, malignant tumors are still one of the most serious diseases that endanger human health. There were about 10.86 million new cases of malignant tumors every year in the world, and about 6.72 million patients died from carcinoma annually [1]. The occurrence of the cancer is a multi-factor process involving the abnormal activation of a series of proto-oncogenes and the inactivation of tumor-suppressor genes due to mutation, which in turn leads to abnormalities in some signaling pathways. In particular, recent studies on hedgehog (Hh) signaling pathways have shown that the aberrant activation of Hh signaling pathways plays an important role in the development of numerous tumors including basal-cell carcinoma (BCC) [2–4], breast cancer [5], medulloblastoma (MB) [6], and prostate cancer [7,8] etc. Recent studies have also shown that Hh signaling is critical for stem cell renewal, tissue repair and regeneration [9]. Thus, it is not surprising that Hh signaling dysfunction can dedicate to multiple human diseases as described above. The Hh signaling pathway is largely dormant in the

healthy adults, and it is appropriately triggered through binding Hh ligands (Sonic hedgehog, Indian hedgehog or Desert hedgehog) [10,11] to the receptor Patched [12–14], a 12-transmembrane transporter that suppresses the activity of the smoothened receptor (SMO) in the absence of Hh ligands [15,16], which in turn brings the transposal of the SMO into cilia [17]. Then the activated SMO initiates a downstream signaling inducing the translocation of Gli transcription factors into the nucleus, thereby promoting the expression of Gli-targeted genes [18]. The SMO is a member of the family F of G protein-coupled receptor (GPCR) [19–21] which is a large family of seven transmembrane proteins involving in cellular signaling pathways. Besides, the sequence similarity between the SMO and frizzled (FZD) receptors [22–24] regulating the WNT signaling pathway [25] is high. It is reported that the function of the SMO acting as a GPCR can be regulated by the natural or synthetic antagonists and agonists [26,27]. The SMO is a key transducer in the Hh signaling pathway, and the inhibition of SMO activity can prevent Gli transcription factors from activating downstream signaling pathways, which mitigate the occurrence of tumors. Therefore,

⁎ Corresponding authors at: Shanghai Engineering Research Center of Molecular Therapeutics & New Drug Development, School of Chemistry and Molecular Engineering, East China Normal University, Shanghai, China. E-mail addresses: [email protected] (Y. Gao), [email protected] (T. Zhu).

https://doi.org/10.1016/j.cplett.2020.137126 Received 23 December 2019; Received in revised form 15 January 2020; Accepted 16 January 2020 Available online 17 January 2020 0009-2614/ © 2020 Elsevier B.V. All rights reserved.

Chemical Physics Letters 741 (2020) 137126

Z. Gao, et al.

researches on the inhibitors targeting the SMO have received extensive attention. Cyclopamine as the first discovered SMO inhibitor, a natural steroidal alkaloid derived from the veratrine with teratogenic properties, hinders the Hh signaling pathway [28–31]. However, poor oral solubility, acid sensitivity and severe by-effect in mice have led to the elimination of its clinical studies [32]. In recent years, researchers have made great efforts to develop inhibitors targeting the SMO primarily. Notably, vismodegib (GDC-0449, which is named GDC in the following) is the first Hh signaling pathway inhibitor to be marketed and approved by the US Food and Drug Administration (FDA) in 2012, which was used for the treatment of metastatic and advanced BCC [33]. Hereafter, sonidegib [34] (LDE-225), the second Hh signaling pathway inhibitor, received the FDA approval in 2015 for treating the BCC. Unfortunately, GDC and LDE-225 were found to show unexpected drug resistance in clinical use [35,36]. It was reported that the D473H-SMO mutant from a recurrent patient tumor tissue treated with GDC abolished drug binding in vitro [37]. LY2940680 (LY) is a representative SMO inhibitor for the treatment of advanced solid tumors in phase I and II trials [38]. Excitedly, LY also has a strong inhibitory effect on the GDC-resistant D473H-SMO mutant [39], highlighting its potential in the clinical treatment of GDC-resistant tumors. Recently, through modifying the 4methylamino-piperidine linker and the 1-methyl-1H-pyrazole group of LY [40], we designed and synthesized a novel candidate inhibitor L4, which strongly block the Hh signaling pathway by inhibiting both the wild-type (WT) and D473H mutated SMO. In addition, L4 showed good tolerance in the ICR mice toxicity test. In order to design more promising SMO inhibitors rationally, the resistance mechanisms of inhibitors induced by the D473H mutant should be explored in detail. To this end, comprehensive computational studies including molecular docking, MD simulations and free energy calculations were conducted in this work. This paper is organized as follows. Section 2 provided the detailed computational methods. Section 3 presented the results and analysis. Finally, a brief summary is given in Section 4.

Fig. 1. The model of the SMO receptor. (A) Superposition of the structure of human Smoothened receptors from two experiments (4JKV, gray; 5L7I, green). Four disulfide bonds are also shown. (B) One of the initial model for the MD simulation.

dynamics (MD) simulations.

2.3. Molecular dynamics simulations To perform simulations in a more realistic environment, the CHARMM-GUI [43,44] was used to build the lipid bilayer environment. The protein orientation in the membrane was defined according to the atomic coordinates of the crystal structure in the OPM (Orientations of Proteins in Membranes) database [45]. The lipid membrane structure was constructed using 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphatidylcholine (POPC) molecules which was placed in the rectangular TIP3P water tank with 0.15 M NaCl, the size of the water boxes are around 90 Å × 90 Å × 110 Å. The entire system containing the SMO receptor, inhibitor, phospholipids, water and ions was shown in Fig. 1B. MD simulations were performed by using the Amber16 software [46]. The force field parameters for three ligands were generated with the general AMBER force field (GAFF) [47], and the AM1-BCC method [48] was used to calculated the atomic point charges. The Amber ff14SB force field [49] and Lipid14 force field [50] were used for the protein and lipids, respectively. The disulfide bonds were constructed by the LEAP module. To remove the bad contacts in the initial structures, structural minimizations were performed by the steepest descent method, followed by the conjugate gradient method through three steps: (i) only the phospholipid tail moved freely; (ii) optimized the solvent while the complex and lipids were restricted with a force constant of 10 2 kcal/(mol·Å ) ; (iii) all atoms were relaxed. Subsequently, the system was gradually heated to 303 K by two steps in a Langevin thermostat with a collision frequency of 1.0 ps−1. Firstly, heating from 0 K to 100 K under the NVT ensemble for 100 ps, while restraining the complex and lipids 2 2 with a force constant of 5 kcal/(mol·Å ) and 2.5 kcal/(mol·Å ) , respectively; secondly, heating from 100 K to 303 K under the NPT ensemble 2 for 500 ps, during which restraints of 2.5 kcal/(mol·Å ) and 1.0 2 kcal/(mol·Å ) were added to the complex and lipids, respectively. Then, 500 ps simulations were performed to make the temperature stable. Thereafter, a 5 ns equilibrium simulation was performed without applying any constraints to each system. Finally, 50 ns MD simulations were carried out at 303 K and 1.0 bar pressure in the NPT ensemble with a time step of 2 fs. During the MD simulations, the particle-mesh Ewald (PME) approach [51,52] was applied to treat the long-range electrostatic interactions using periodic boundary conditions (PBC) with a 10.0 Å cut-off value. Bonds involving hydrogen atoms were confined by SHAKE algorithm [53]. Three independent replicas with different random seed numbers were generated for each system to improve sampling and to check the convergence of simulation.

2. Methods 2.1. System construction The crystal structures of human SMO complex bound with LY2940680 (PDB code: 4JKV) [41] and Vismodegib (PDB code: 5L7I) [42] were obtained from the Protein Data Bank. The subunit B of the crystal structures was used to build the receptor model. Besides, to further save the computational resources, the extracellular cysteine-rich domain (CRD) was removed, leaving the extracellular linker domain (ECD), the long extracellular loops (ECLs), and the heptahelical transmembrane domain (TMD) retained, which were stabilized by four disulfide bonds (Fig. 1A). On the basis of the WT model, the D473H mutant was constructed by replacing the specific residue ASP473 with histidine. 2.2. Molecular docking Since there is no experimental structure the molecular docking was adoped to construct the structure of the WT/Mutant SMO-L4 complex. The Glide module in Schrödinger 2013 software package was used. WT and D473H mutant LY receptors that have been built above were prepared by using the Protein Preparation Wizard. All water molecules (except those coordinated to the ligand or residues) were abandoned and missing hydrogen atoms were added. Minimization with restraints for the protein was conducted with the OPLS2005 force field to redirect side-chain hydroxyl groups and mitigate potential space collisions. The L4 molecule was prepared by using the LigPrep module. Ligand docking was conducted in the extra precision (XP) mode with the default parameters in the RRD (Receptor Rigid Docking) treaty. The binding pose which has the highest docking score was used for the molecular 2

Chemical Physics Letters 741 (2020) 137126

Z. Gao, et al.

2.4. Free energy calculation by MM-PBSA In this work, MM-PBSA (Molecular Mechanics Poisson-Boltzmann Surface Area) method was applied to calculate the binding free energy (ΔGbind ) of each system. The ΔGbind is shown:

ΔGbind

= Gcomplex − Greceptor − Gligand = ΔEvdw + ΔEele + ΔGpol + ΔGnopol − TΔS

(1)

in which the terms Gcomplex , Greceptor , and Gligand are the free energies of the complex, protein and ligand, respectively. The ΔEvdw and ΔEele denote the van der Waals and electrostatic energy change unpon binding. In this work, the implicit membrane model [54] was employed to calculate the solvation free energy including the polar (ΔGpol ) and the nopolar terms (ΔGnopol ). The ΔGpol is obtained by solving Poisson-Boltzmann equation, and the ΔGnopol is considered to be linearly proportional to the change of the solvent accessible surface area.

ΔGnopol = γ ΔSASA + β

(2)

The surface-tension (γ ) and the offset (β ) are empirical constants fitted according to experiment [55]. In this work, the two values are set 2 2 to 0.005 kcal/(mol·Å ) and 0 kcal/(mol·Å ) , repectively. This research primarily focused on the qualitative analysis of the binding modes, without considering the calculation of absolute free energy. Given that the entropy changes are almost the same due to the similarity of the three inhibitors and the calculation of the entropy using such as Nmode method is also inaccurate, the contribution of entropy is ignored in this work. The calculation of the enthalpy change and the energy decomposition are based on the 100 snapshots of the last 10 ns MD trajectories. 3. Results and discussion 3.1. Molecular docking Fig. 2. Three-dimensional interaction diagram of inhibitor L4 bound to WT (A) and D473H (B) SMO. All disrtances are in Å.

The initial binding poses of L4 in the active sites of WT/D473H SMO have been obtained by molecular docking with the docking scores of −8.936 and −9.643, respectively. There is no obvious difference for the docking score between the WT and mutant system, which is consistent with the corresponding experimental PIC50 values (Table 1). As shown in Fig. 2, the L4 molecule is wholly coordinated in the active sites of the WT/D473H SMO. Specifically, there is one hydrogen bond formed between the amide oxygen atom of L4 and the ASN219 residue. The 4,5-dimethylpyridazine group also forms two hydrogen bonds with the side chain of ARG400. These hydrogen bonds are similar to the crystal structure of SMO in complex with LY reported by Wang et al. [41]. The 1-fluoro-3(trifluoromethyl)-benzene group also makes a hydrogen bond with LYS395. In addition, the electron-rich benzene rings of PHE484, TYR394, and TRP281 form hydrophobic interactions with the benzene ring, the pyridazine ring, and the pyrazole ring of L4 through π-π stacking. As there is no derectly interaction between ASP473 with the ligand, these interactions do not change obviously after the D473H mutation.

3.2. MD simulations 3.2.1. Analysis of the stability of simulated systems The root mean square deviations (RMSDs) were tracked to characterize the stability of each system during the MD simulations. As shown in Figs. S1–S3, the RMSD values of protein backbone atoms for all simulated systems vary from 1 Å to 2 Å, implying the equilibration of all simulations. In addition, the time evolutions of the pocket RMSDs (The RMSD of Cα atoms of residues within 5 Å of ligands.) are plotted. Apparently, for all systems the values are very small and stable in the whole simulations. The RMSDs of the heavy atoms in ligands are also shown. Although a sharp rise is observed for the WT-GDC system in the beginning of the simulation, the RMSDs reach convergence quickly. The structural fluctuation is because the highly electronegative Cl atom on the GDC tends to approach the amide group in the side chain of the GLN477. The RMSD results indicate that the modeled systems have

Table 1 The structures of ligands and the corresponding experimental PIC50 values. Compound

GDC

Complex PIC50 (exp.)

WT-GDC 7.66 ± 0.11

LY

D473H-GDC 6.48 ± 0.16

L4

WT-LY 6.53 ± 0.10

The experimental data [40,56] was calculated by the equation: PIC50 = − log10 IC50 . 3

D473H-LY 6.22 ± 0.24

WT-L4 7.63 ± 0.0

D473H-L4 7.61 ± 0.0

Chemical Physics Letters 741 (2020) 137126

Z. Gao, et al.

Table 2 Average binding free energies of inhibitors bound to WT/D473H SMO. All values are in kcal/mol. The signs ± represented standard errors of the mean. Complex

Contribution

GDC

LY

L4

WT

ΔEvdw ΔEele ΔGpol

−54.72 ± 0.26 −38.02 ± 0.51 66.41 ± 0.37

−66.16 ± 0.29 −31.95 ± 0.36 56.91 ± 0.38

−63.50 ± 0.29 −35.86 ± 0.46 56.16 ± 0.40

D473H

ΔGnopol

−4.67 ± 0.01

−5.66 ± 0.01

−5.47 ± 0.01

ΔGbind

−31.00 ± 0.45

−46.86 ± 0.44

−48.67 ± 0.42

ΔEvdw ΔEele ΔGpol

−53.41 ± 0.32 −34.09 ± 0.47 66.38 ± 0.55

−65.43 ± 0.30 −38.02 ± 0.43 60.31 ± 0.48

−63.79 ± 0.28 −36.14 ± 0.45 54.52 ± 0.35

ΔGnopol

−4.83 ± 0.02

−5.64 ± 0.01

−5.41 ± 0.01

ΔGbind

−25.95 ± 0.43

−48.78 ± 0.39

−50.82 ± 0.38

contributed by the electrostatic interaction (Fig. 3E&F). It worth to mention that, ASP473 itself makes negative contributions to the binding for all WT systems. As shown in Fig. 3, this negative contribution is determined by the electrostatic interaction. By mutating the negative charged ASP to the neutral HIS, its electrostatic interaction with the ligands has almost vanished, leaving only the slightly beneficial van der Waals interaction. Therefore, the binding of mutated protein with LY and L4 is even stronger. However, the contribution of ASP473 to the binding in the D473H-GDC system is unable to compensate for the loss of binding affinities from ARG400 and GLN477.

reached equilibrium during the simulation, and it is reasonable to calculate the binding free energy based on the last 10 ns MD trajectories. 3.2.2. The binding free energies In order to obtain the specific energy information of the drug resistance caused by mutation, the binding free energies were calculated by the MM-PBSA method. Energy calculations were performed for three replicate simulations of each system (Tables S1–S3). As shown in Table 2, the rank of the predicted binding affinities agree very well with the experimental data (Table 1). It can be seen that there is a significant decline in the binding affinity for the GDC molecule after mutation (from −31.00 kcal/mol to −25.95 kcal/mol). However, the binding affinity of LY and L4 combined with WT/D473H SMO exhibites absent decrease, but increasing slightly by 1.92 and 2.15 kcal/mol, respectively. Comparing the individual energy terms between the WT and D473H SMO can bring deeper insights into the drug resistance mechanism. According to Table 2, the Van der Waals interactions (ΔEvdw ) between GDC with D473H-SMO are weakened by 1.31 kcal/mol. The electrostatic interactions (ΔEele ) offers the greatest contribution to the drug resistance with a decrease of 3.93 kcal/mol in the binding affinity relative to the WT-SMO. The polar solvation energy (ΔGpol ) hardly change, which is reduced by 0.03 kcal/mol compared to the WT-SMO. The non-polar solvation energies (ΔGnopol ) also produce barely contribution to the drug resistance. As for the LY and L4, the non-polar interactions, including the Van der Waals interactions (ΔEvdw ) and non-polar solvation energies (ΔGnopol ) almost don’t change after the mutation. Furthermore, the polar interactions of LY and L4 with the mutated SMO, consisting of the electrostatic interactions (ΔEele ) and the polar solvation free energies (ΔGpol ), are decreased moderately by 2.67 and 1.92 kcal/mol, respectively, which means that the binding of these two molecules with the mutated SMO are even stronger.

3.2.4. Structural comparison of the WT/D473H SMO in complex with inhibitors In order to further explain the energy changes discussed above, a comparative structural analysis were conducted. As shown in Table 3 and Fig. 4A, it is obvious that there is a hydrogen-bonding network in the WT-SMO. Residue ARG400 forms three hydrogen bonds with GLU518 with the occupancies 97.88%, 96.74% and 64.9% respectively. Besides, ASP473 forms two hydrogen bonds with residues GLN477 and ARG400, in which the occupancies are both 99.96%. It is important to note that the pyridine of GDC forms a hydrogen bond with ARG400 with 86.16% occupancy. However, in the D473H-SMO, the HIS473 cannot form these two hydrogen bonds. At the same time, the three hydrogen bonds between ARG400 and GLU518 are almost disrupted, their occupancies are reduced to 32.94%, 44.10%, and 41.92%, respectively. The disruption of the hydrogen bond network can cause the instability of the binding pocket. As shown in Fig. 4B, it seems that the side chain of ARG400 is far away from the GDC, so the hydrogen bond between ARG400 and the pyridine (The N1 atom) of GDC is disrupted. The distribution of the distance between NH2@ARG400 and N1@GDC is shown in Fig. 5. In the wild type, this distance mainly varies from 2.8 Å to 4.0 Å with the peak at 3.2 Å, suggesting that it is very stable over time. However, in the mutant system, the distance fluctuates greatly against time with an average distance of 5.0 Å. These results can explain that the electrostatic contribution of ARG400 is reduced sharply as shown in the Fig. 3B. Surprisingly, in the D473H-SMO, GLN477 forms a hydrogen bond with GLY393 with a 94.92% occupancy (Table 3). Since the residue GLY393 locates on the loop domain which is away from the GDC, the side chain of GLN477 cannot interact effectively with the electronegative Cl atom on the GDC (Figs. 4B & 5B), which agrees well with the decrease of the electrostatic contribution from GLN477 to the binding of GDC discussed above (Fig. 3B). For the LY, the hydrogen bond networks are both present in the binding pockets of the WT-SMO and D473H-SMO so that the residues ARG400 and GLN477 are well stabilized around LY. The average distances between the NH2 atom of ARG400 and the N4 atom of LY are 3.1 Å and 3.0 Å in the WT/D473H-SMO, respectively (Fig. 5C). The distributions of the distance between the NE2 atom of GLN477 and the C atom on trifluoromethyl of LY are highly overlapping in the wild type and mutated system (Fig. 5D). As shown in Table S4 as well as Fig. 4C and D, in the WT-SMO, the residue ARG400 forms two hydrogen bonds

3.2.3. Changes of ligand-residue interactions caused by D473H mutant To further illustrate the cause of the drug resistance induced by mutation, the decomposition of the free energy with each residue in the binding pocket was calculated. For GDC, the results reveal that the largest change in the binding affinity after mutation is contributed by the residue ARG400 and GLN477, which is reduced by 5.02 and 4.07 kcal/mol, respectively (Fig. 3A). The electrostatic interactions of ARG400 and GLN477 with GDC are weakened by 3.89 and 3.63 kcal/ mol by the mutation, while their van der Waals interactions with GDC are slightly weakened by 1.13 and 0.44 kcal/mol (Fig. 3B). These results indicate that the electrostatic interaction is the main driving force of the drug resistance. As shown in Fig. 3C, the contribution of ARG400 to the binding of LY almost has no alteration in the wild type and mutant systems. While the contribution of the GLN477 is decreased slightly, which is mainly due to the change of the electrostatic interaction (Fig. 3D). As for L4, the contribution of ARG400 to the binding is enhanced about 1 kcal/mol in the mutated system, which is mainly 4

Chemical Physics Letters 741 (2020) 137126

Z. Gao, et al.

Fig. 3. The ligand-residue interaction spectra for the corresponding complexes WT/D473H-GDC (A, B), WT/D473H-LY (C, D) and WT/D473H-L4 (E, F). Only the three residues with the most energy change by the mutation are given. The detailed contributions of electrostatic and van der Waals interactions are also shown.

Table 3 The hydrogen bonds for the GDC bound to the WT/D473H-SMO. The atom names of inhibors can be found in Fig. S4. Complex

Donor

Acceptor

Distance(Å)

Angle(°)

Occupancy(%)

WT-GDC

ARG400(NH1-HH11) ARG400(NH2-HH22) ARG400(NH2-HH22) ARG400(NH1-HH12) GLN477(NE2-HE22) ARG400(NH2-HH21)

ASP473(OD1) GLU518(OE2) GLU518(OE1) GLU518(OE1) ASP473(OD1) GDC(N1)

2.77 2.94 2.98 2.77 2.78 3.16

160.00 148.02 145.08 156.53 160.14 150.02

99.96 97.88 96.74 64.90 99.96 86.16

D473H-GDC

ARG400(NH2-HH21) ARG400(NH2-HH22) ARG400(NH2-HH22) ARG400(NH1-HH12) GLN477(NE2-HE21) ARG400(NH2-HH22)

HIS473(ND1) GLU518(OE2) GLU518(OE1) GLU518(OE1) GLY393(O) GDC(N1)

2.93 3.18 2.85 2.89 2.92 3.43

144.82 146.38 151.08 146.43 147.66 111.86

30.32 32.94 44.10 41.92 94.92 0.02

5

Chemical Physics Letters 741 (2020) 137126

Z. Gao, et al.

Fig. 4. The hydrogen-bonding network surrounding ARG400 for the corresponding inhibitors GDC (top), LY (middle) and L4 (bottom) in the WT-SMO (left) and D473H-SMO (right), respectively. All disrtances are in Å.

described by electrostatic interactions, thus these results are good explanations of the calculated energies discussed above (Fig. 3C and D). For the L4, the average distances between the NH2 atom of ARG400 and the N4 atom of L4 are 3.3 Å and 3.1 Å in the WT-SMO and D473HSMO, respectively (Fig. 5E); and the distribution of distances between the NE2 atom of GLN477 and the C atom on trifluoromethyl of L4 are also similar in the two systems (Fig. 5F). This phenomenon is attributed to the presence of hydrogen bond networks in the binding pockets of WT/D473H-SMO. According to the Table S5, Fig. 4E and F, in the WTSMO, there exist three hydrogen bonds between the residue ARG400 and L4 (96.0%, 86.22%, and 71.02% occupancy, respectively). In the mutated system, these three hydrogen bonds still exist, for which the occupancies are 90.9%, 94.08% and 72.06%, respectively. The residue ASP473 forms two hydrogen bonds in the WT-SMO: one is with GLN477

with LY, which the occupancies are 99.58% and 98.18%, respectively. There also exists a hydrogen bond between ASP473 and GLN477 (99.92% occupancy). These three hydrogen bonds are still stable in the D473H-SMO (The occupancies are 96.62%, 86.22%, and 98.78%, respectively). In addition, the residue ARG400 forms three hydrogen bonds with other residues in the WT-SMO: two with ASP473 (with 99.08% and 64.14% occupancy, respectively) and one with GLU518 (with 94.60% occupancy). Although these three hydrogen bonds are disrupted in the mutant system, two new hydrogen bonds are generated: one H-bond is between HIS473 and GLU518 (98.24% occupancy), and the other one is between ARG400 and HIS470 (96.24% occupancy). The existence of these hydrogen bond networks can maintain the stability of the binding pocket. It should be mentioned that in the MD simulation with Amber ff14SB force field, the hydrogen bond are 6

Chemical Physics Letters 741 (2020) 137126

Z. Gao, et al.

Fig. 5. Distance distribution of (A), NH2@ARG400–N1@GDC, (B), NE2@GLN477–Cl@GDC, (C), NH2@ARG400–N4@LY, (D), NE2@GLN477–C21@LY, (E), NH2@ARG400–N4@L4, and (F), NE2@GLN477–C23@L4. The atom names of inhibors can be found in Fig. S4.

because that these is a new hydrogen bond formed beteen GLN477 and GLY393, so the side chain of GLN477 is farther apart from the Cl atom on the GDC. But for the LY and L4, there both exist hydrogen bond networks around ARG400 in the WT-SMO and D473H-SMO, which can maintain the stability of the binding pocket and the interaction between ligand and protein. These results and analysis can provide a comprehensive explanation about the molecular mechanisms of drug resistance toward SMO, with the aim of making a theoretical contribution to design more promising SMO inhibitors.

(99.94% occupancy), and the other one is with ARG400 (96% occupancy). The residue ARG400 also forms a H-bond with GLU518 (67.68% occupancy). In the D473H-SMO, although only one H-bond between HIS473 and GLN477 (99.96% occupancy) is observed, three new H-bonds are formed: HIS473 forms two H-bonds with the residues HIS470 and GLU518 (72.34% and 69.08% occupancy, respectively) and ARG400 forms a hydrogen bond with HIS470 (96.98% occupancy). To some extent, the hydrogen bond networks in both systems stabilize the binding pocket. This result agrees with the minor difference in the contribution of ARG400 and GLN477 to the binding of L4 with the WTSMO and D473H-SMO (Fig. 3E and F). In summary, a hydrogen bond network is formed around ARG400 for all three inhibitors in the WT-SMO. While although this hydrogen bond network is destroyed by the mutation, new hydrogen bonds can be formed in the D473H-LY and D473H-L4 systems, which can also stabilize the binding pocket. However, for GDC, the hydrogen bond network destroyed by the mutation cannot be compensated by new hydrogen bonds, which may induce the drug resistance.

CRediT authorship contribution statement Zhenmei Gao: Writing - original draft. Jingxiao Bao: Methodology, Software. Shuhua Shi: Validation. Xiongwen Zhang: Resources. Ya Gao: Methodology, Writing - original draft. Tong Zhu: Conceptualization, Methodology, Writing - review & editing. Declaration of Competing Interest

4. Conclusions

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

In this paper, the drug resistance mechanisms of SMO to inhibitors induced by the D473H mutant are explored based on molecular docking, MD simulations and free energy calculations. The calculated results indicate that there is a obvious decline in the binding affinity of D473H-GDC with a decrease of 5.05 kcal/mol relative to the WT-GDC. While the binding affinities of LY and L4 bound to WT/D473H SMO show no obvious differences. In detail, the difference of interactions causing drug resistance is primarily reflected on the electrostatic and Van der Waals terms. Comparisons of the energy contribution of key residues between WT-GDC and D473H-GDC reveals that the largest changes in binding free energies are contributed by the residue ARG400 and GLN477. The further analysis of hydrogen bonds shows that the decrease of contributions of ARG400 to the binding of GDC in the mutated system is because that the D473H mutation disrupts the hydrogen-bonding network surrounding ARG400. While the decrease of contributions of GLN477 to the binding of GDC in the mutated system is

Acknowledgments This work was supported by the Ministry of Science and Technology of China (Grant No. 2016YFA0501700), the National Natural Science Foundation of China (Grants No.91641116), and the Key Research and Development Project of Shandong Province (No. 2018GSF121014). We thank the ECNU Multifunctional Platform for Innovation (001) for providing us with computer time. Appendix A. Supplementary material Supplementary data to this article can be found online at https:// doi.org/10.1016/j.cplett.2020.137126. 7

Chemical Physics Letters 741 (2020) 137126

Z. Gao, et al.

References

[33] A. Sekulic, M.R. Migden, A.E. Oro, L. Dirix, K.D. Lewis, J.D. Hainsworth, J.A. Solomon, S. Yoo, S.T. Arron, P.A. Friedlander, E. Marmur, C.M. Rudin, A.L.S. Chang, J.A. Low, H.M. Mackey, R.L. Yauch, R.A. Graham, J.C. Reddy, A. Hauschild, N. Engl, J. Med. 366 (2012) 2171. [34] S. Pan, X. Wu, J. Jiang, W. Gao, Y. Wan, D. Cheng, D. Han, J. Liu, N.P. Englund, Y. Wang, A.C.S. Med, Chem. Lett. 1 (2010) 130. [35] C.M. Rudin, C.L. Hann, J. Laterra, R.L. Yauch, C.A. Callahan, L. Fu, T. Holcomb, J. Stinson, S.E. Gould, B. Coleman, N. Engl, J. Med. 361 (2009) 1173. [36] S. Buonamici, J. Williams, M. Morrissey, A. Wang, R. Guo, A. Vattay, K. Hsiao, J. Yuan, J. Green, B. Ospina, Sci. Transl. Med. 2 (2010) 51ra70. [37] R.L. Yauch, G.J. Dijkgraaf, B. Alicke, T. Januario, C.P. Ahn, T. Holcomb, K. Pujara, J. Stinson, C.A. Callahan, T. Tang, Science 326 (2009) 572. [38] J. Bendell, V. Andre, A. Ho, R. Kudchadkar, M. Migden, J. Infante, R.V. Tiu, C. Pitou, T. Tucker, L. Brail, D. Von Hoff, Clin. Cancer Res. 24 (2018) 2082. [39] M.H. Bender, P.A. Hipskind, A.R. Capen, M. Cockman, K.M. Credille, H. Gao, J.A. Bastian, J.M. Clay, K.L. Lobb, D.J. Sall, M.L. Thompson, T. Wilson, G.N. Wishart, B.K.R. Patel, Cancer Res. 71 (2011) 1. [40] M.F. Zhu, H. Wang, C.L. Wang, Y.F. Fang, T. Zhu, W.L. Zhao, X.C. Dong, X.W. Zhang, Front. Pharmacol. 10 (2019) 13. [41] C. Wang, H. Wu, V. Katritch, G.W. Han, X.-P. Huang, W. Liu, F.Y. Siu, B.L. Roth, V. Cherezov, R.C. Stevens, Nature 497 (2013) 338. [42] E.F.X. Byrne, R. Sircar, P.S. Miller, G. Hedger, G. Luchetti, S. Nachtergaele, M.D. Tully, L. Mydock-McGrane, D.F. Covey, R.P. Rambo, M.S.P. Sansom, S. Newstead, R. Rohatgi, C. Siebold, Nature 535 (2016) 517. [43] S. Jo, T. Kim, W. Im, PLoS One 2 (2007) 9. [44] S. Jo, T. Kim, V.G. Iyer, W. Im, J. Comput. Chem. 29 (2008) 1859. [45] M.A. Lomize, I.D. Pogozheva, H. Joo, H.I. Mosberg, A.L. Lomize, Nucleic Acids Res. 40 (2012) D370. [46] D.A. Pearlman, D.A. Case, J.W. Caldwell, W.S. Ross, T.E. Cheatham III, S. DeBolt, D. Ferguson, G. Seibel, P. Kollman, Comput. Phys. Commun. (Netherlands) 91 (1995) 1. [47] J. Wang, R.M. Wolf, J.W. Caldwell, P.A. Kollman, D.A. Case, J. Comput. Chem. 25 (2004) 1157. [48] A. Jakalian, B.L. Bush, D.B. Jack, C.I. Bayly, J. Comput. Chem. 21 (2000) 132. [49] J.A. Maier, C. Martinez, K. Kasavajhala, L. Wickstrom, K.E. Hauser, C. Simmerling, J. Chem. Theory Comput. 11 (2015) 3696. [50] C.J. Dickson, B.D. Madej, A.A. Skjevik, R.M. Betz, K. Teigen, I.R. Gould, R.C. Walker, J. Chem. Theory Comput. 10 (2014) 865. [51] T. Darden, D. York, L. Pedersen, J. Chem. Phys. (USA) 98 (1993) 10089. [52] U. Essmann, L. Perera, M.L. Berkowitz, T. Darden, L. Hsing, L.G. Pedersen, J. Chem. Phys. (USA) 103 (1995) 8577. [53] J.P. Ryckaert, G. Ciccotti, H.J.C. Berendsen, J. Comput. Phys. (USA) 23 (1977) 327. [54] W.M. Botello-Smith, X.P. Liu, Q. Cai, Z.L. Li, H.K. Zhao, R. Luo, Chem. Phys. Lett. 555 (2013) 274. [55] D. Sitkoff, K.A. Sharp, B. Honig, J. Chem. Phys. 98 (1994) 1978. [56] E. Lauressergues, P. Heusler, F. Lestienne, D. Troulier, I. Rauly-Lestienne, A. Tourette, M.C. Ailhaud, C. Cathala, S. Tardif, D. Denais-Lalieve, M.T. Calmettes, A.D. Degryse, A. Dumoulin, L. De Vries, D. Cussac, Pharmacol. Res. Perspect. 4 (2016) 14.

[1] D.M. Parkin, F. Bray, J. Ferlay, P. Pisani, CA: Cancer J. Clinicians 55 (2005) 74. [2] A.E. Oro, K.M. Higgins, Z.L. Hu, J.M. Bonifas, E.H. Epstein, M.P. Scott, Science 276 (1997) 817. [3] J.W. Xie, M. Murone, S.M. Luoh, A. Ryan, Q.M. Gu, C.H. Zhang, J.M. Bonifas, C.W. Lam, M. Hynes, A. Goddard, A. Rosenthal, E.H. Epstein, F.J. de Sauvage, Nature 391 (1998) 90. [4] E.H. Epstein, Nat. Rev. Cancer 8 (2008) 743. [5] A. Flemban, D. Qualtrough, Cancers 7 (2015) 1863. [6] J.T. Romer, H. Kimura, S. Magdaleno, K. Sasai, C. Fuller, H. Baines, M. Connelly, C.F. Stewart, S. Gould, L.L. Rubin, T. Curran, Cancer Cell 6 (2004) 229. [7] P.C. Walsh, J. Urol. 173 (2005) 1169. [8] J. Xie, Future Oncology (London, England) 1 (2005) 331. [9] P.A. Beachy, S.S. Karhadkar, D.M. Berman, Nature 432 (2004) 324. [10] Y. Echelard, D.J. Epstein, B. St-Jacques, L. Shen, J. Mohler, J.A. McMahon, A.P. McMahon, Cell 75 (1993) 1417. [11] P.W. Ingham, A.P. McMahon, Genes Dev. 15 (2001) 3059. [12] N. Fuse, T. Maiti, B. Wang, J.A. Porter, T.M.T. Hall, D.J. Leahy, P.A. Beachy, Proc. Natl. Acad. Sci. USA 96 (1999) 10992. [13] V. Marigo, R.A. Davey, Y. Zuo, J.M. Cunningham, C.J. Tabin, Nature 384 (1996) 176. [14] D.M. Stone, M. Hynes, M. Armanini, T.A. Swanson, Q. Gu, R.L. Johnson, M.P. Scott, D. Pennica, A. Goddard, H. Phillips, Nature 384 (1996) 129. [15] J. Taipale, J.K. Chen, M.K. Cooper, B.L. Wang, R.K. Mann, L. Milenkovic, M.P. Scott, P.A. Beachy, Nature 406 (2000) 1005. [16] J. Taipale, M. Cooper, T. Maiti, P. Beachy, Nature 418 (2002) 892. [17] K.C. Corbit, P. Aanstad, V. Singla, A.R. Norman, D.Y.R. Stainier, J.F. Reiter, Nature 437 (2005) 1018. [18] S.J. Scales, F.J. de Sauvage, Trends Pharmacol. Sci. 30 (2009) 303. [19] S.S.G. Ferguson, L.S. Barak, J. Zhang, M.G. Caron, Can. J. Physiol. Pharmacol. 74 (1996) 1095. [20] R.T. Premont, R.R. Gainetdinov, Annu. Rev. Physiol. 69 (2007) 511. [21] M. Zalewska, M. Siara, W. Sajewicz, Acta Pol. Pharm. 71 (2014) 229. [22] P. Bhanot, M. Brink, C.H. Samos, J.C. Hsieh, Y.S. Wang, J.P. Macke, D. Andrew, J. Nathans, R. Nusse, Nature 382 (1996) 225. [23] Y.S. Wang, J.P. Macke, B.S. Abella, K. Andreasson, P. Worley, D.J. Gilbert, N.G. Copeland, N.A. Jenkins, J. Nathans, J. Biol. Chem. 271 (1996) 4468. [24] X. He, J.P. SaintJeannet, Y.S. Wang, J. Nathans, I. Dawid, H. Varmus, Science 275 (1997) 1652. [25] A. Klaus, W. Birchmeier, Nat. Rev. Cancer 8 (2008) 387. [26] K.L. Ayers, P.P. Therond, Trends Cell Biol. 20 (2010) 287. [27] P. Heretsch, L. Tzagkaroulaki, A. Giannis, Bioorg. Med. Chem. 18 (2010) 6613. [28] M.K. Cooper, J.A. Porter, K.E. Young, P.A. Beachy, Science 280 (1998) 1603. [29] J.P. Incardona, W. Gaffield, R.P. Kapur, H. Roelink, Development 125 (1998) 3553. [30] R.F. Keeler, Proc. Soc. Exp. Biol. Med. 149 (1975) 302. [31] R.F. Keller, Phytochemistry 7 (1968) 303. [32] M.R. Tremblay, A. Lescarbeau, M.J. Grogan, E. Tan, G. Lin, B.C. Austad, L.-C. Yu, M.L. Behnke, S.J. Nair, M. Hagel, J. Med. Chem. 52 (2009) 4400.

8