Deciphering the Allosteric Effect of Antagonist Vismodegib on Smoothened Receptor Deactivation Using Metadynamics Simulation

The smoothened receptor (Smo) plays a key role in Hedgehog (Hh) signaling pathway and it has been regarded as an efficacious therapeutic target for basal cell carcinoma (BCC) and medulloblastoma (MB). Nevertheless, the resistance mutation and active mutants of Smo have put forward the requirement of finding more effective inhibitors. Herein, we performed metadynamics simulations on Smo bound with vismodegib (Smo-Vismod) and with cholesterol (Smo-CLR), respectively, to explore the inhibition mechanism of vismodegib. The simulation results indicated that vismodegib-induced shifts of TM5, TM6, and TM7, which permitted the extracellular extension of TM6 and extracellular loop3 (ECL3) to enter the extracellular cysteine-rich domain (CRD) groove. Therefore, an open CRD groove that has not been noticed previously was observed in Smo-Vismod complex. As a consequence, the occupied CRD groove prevents the binding of cholesterol. In addition, the HD and ECLs play crucial roles in the interaction of CRD and TMD. These results reveal that TM5, TM6, and TM7 play important roles in allosteric inhibition the activation of Smo and disrupting cholesterol binding by vismodegib binding. Our results are expected to contribute to understanding the allosteric inhibition mechanism of Smo by vismodegib. Moreover, the detailed conformational changes contribute to the development of novel Smo inhibitors against resistance mutation and active mutants of Smo.


INTRODUCTION
The Hedgehog (Hh) signaling pathway plays central roles in the animal development and stem-cell function, linking to cell growth, and differentiation, with normal roles in embryonic pattern formation and adult tissue homeostasis and pathological roles in tumor initiation and growth (Pasca di Magliano and Hebrok, 2003;Lum and Beachy, 2004;Rohatgi and Scott, 2007). The smoothened receptor (Smo) serves as one of the key proteins in this signaling pathway, the phosphorylation, stabilization, and accumulate of Smo cause signal to intracellular effectors (Denef et al., 2000;Kalderon, 2000).
Smo has been regarded as an efficacious therapeutic target for basal cell carcinoma (BCC) and medulloblastoma (MB) (Kim et al., 2013). Nevertheless, therapeutic challenges remain in the tumors that acquire resistance to Smo antagonists (Rudin et al., 2009;Yauch et al., 2009;Dijkgraaf et al., 2011;Sharpe et al., 2015a), also in cases of active mutants of Smo that were reported as oncogenic drivers (Reifenberger et al., 1998;Sweeney et al., 2014). So the emergence of resistant mutations and active mutations has raised the need for the discovery of novel Smo inhibitors or novel target sites of Smo.
Previous biochemical studies have shown that sterols binding to CRD groove site activate Hh signaling pathway. Binding of an agonist or antagonist to the TMD-site can activate or inactivate Hh signaling pathway. In addition, the vismodegib binding to the TMD-site results in loss of cholesterol from the CRD-linker domain-TMD interface (Byrne et al., 2016). Zhang et al. have demonstrated that TM6, extracellular loop 3 (ECL3) and the HD play a central role in signal transmission, and their structures reveal a precise arrangement of TMD, HD, and CRD. This structure enables allosteric interactions between the three domains that are important to ligand recognition and receptor activation (Zhang et al., 2017). Yet, the detailed mechanisms of how vismodegib allosterically inhibits the activation of Smo and binding of cholesterol remain unknown.
The crystal structures of Smo binding with vismodegib and cholesterol have similar conformation, both in inactive states. We have little known from the static crystal structures about the conformational variation in Smo upon binding different ligands. The computational methods such as molecular dynamics (MD) simulations can provide the information about the dynamic process of conformational changes at the atomic level upon agonist and antagonist binding to GPCRs (McRobb et al., 2016;Miao and McCammon, 2016;Latorraca et al., 2017). The enhanced sampling method, metadynamics (Laio and Parrinello, 2002;Micheletti et al., 2004), supplies an effective and reliable way to explore the binding and unbinding of ligand from GPCRs, and conformational dynamics of GPCRs binding different ligands (Li et al., 2013;Schneider et al., 2015;Saleh et al., 2017). Bai et al performed metadynamics simulation on Smovismodegib complex to explore the binding mechanism between the vismodegib and Smo (Bai et al., 2014).
Herein we performed metadynamics simulation to shed light on the mechanism that vismodegib allosterically inhibits the activation of human

Preparation of Simulation Systems
As a starting point in the simulations, and as a reference conformation to analyze the results, we used previously determined structures of the human Smo bound to vismodegib and cholesterol (Byrne et al., 2016). The crystal structures were obtained from the PDB database (PDB ID: 5L7I, 5L7D  Sastry et al., 2013). Side chain ionization states were modeled with the PROPKA tool (Søndergaard et al., 2011). The membrane around the transmembrane domain of Smo was built by 85 × 85 Å POPC: cholesterols with 9:1 using CHARMM-GUI webserver (Lee et al., 2016), the receptor crystal structure pre-aligned in the OPM (Orientations of Proteins in Membranes) database (Lomize et al., 2006). Each system was solvated by 12 Å with a truncated rectangular box of TIP3P waters (Jorgensen et al., 1983) and neutralized to a concentration of 0.15 M NaCl.
The proteins were modeled using the AMBER FF99SB force field (Hornak et al., 2006), the ligands were modeled using the generalized AMBER force field (GAFF) (Wang et al., 2004). Geometry optimization and the electrostatic potential calculations on the ligands were performed at the HF/6-31G * level in the Gaussian09 software (Frisch et al., 2009), and the partial charges were calculated with the RESP (Fox and Kollman, 1998). The force field parameters for the ligands were created by the Antechamber package.
Before metadynamics simulation, the energy minimization and equilibration were conducted by NAMD 2.9 simulation package (Phillips et al., 2005) in order to equilibrate the systems. Firstly, to remove bad contacts in the initial structures, steepest descent was carried out. After energy minimization, each system was gradually heated in NVT ensemble from 0 to 300 K in 300 ps. Subsequently, constant temperature equilibration at 300 K for a total of 5 ns was performed to adjusting the solvent density. Finally, 20 ns conventional molecular dynamic simulations were carried out for each system in NPT ensemble with periodic boundary conditions; an integration step of 2 fs was used. The particle mesh Ewald (PME) algorithm (Darden et al., 1993) was employed to treat long-range electrostatic interactions, while the non-bonded interactions were calculated based on a cutoff of 12 Å. The SHAKE algorithm (Ryckaert et al., 1977) was applied to constrain all covalent bonds involving hydrogen atoms.

Metadynamics Simulations
Metadynamics is an efficient enhanced sampling method, allows the system to escape from local minima in the free energy surface (FES) to explore the conformational space by filling the minima with an external history-dependent bias potential, and permits an accurate determination of the FES (Laio and Parrinello, 2002;Laio and Gervasio, 2008). This bias potential is built as a sum of Gaussians deposited along the trajectory in the pre-defined collective coordinates (CVs) space.
where V − → s , t is the bias potential added to the system, τ is the Gaussian deposition stride, σ i is the width of the Gaussian for the ith CV, and W kτ is the height of the Gaussian. The effect of the metadynamics biased potential is to push the system away from local minima into visiting new regions of the phase space. Furthermore, in the long time limit, the bias potential converges to minus the free energy as the function of the CVs: In standard metadynamics, Gaussians of constant height is added for the entire course of a simulation. As a result, the system is eventually pushed to explore high free-energy regions and the estimate of the free energy calculated from the bias potential oscillates around the real value.
In this work, we carried out 100 ns metadynamics simulations for each system, biasing the potential along the following two CVs: the distance between the center mass of W109 and R161, as well as the distance between the center mass of P263 and W535, ignoring hydrogen atoms and labeled as d1, d2, respectively. A Gaussian width of 0.15 Å was used for both CVs, and a Gaussian deposition rate of 0.1 kcalmolps was used.

Determination of Conformational States From Free Energy Surface
W109 and R161 of CRD groove located in the opposite place are used to characterize the conformational dynamics of the groove site, and they are involved in binding of sterols (Rana et al., 2013). P263 located in the C-terminal of ICL1 and W535 located in the intercellular tip of TM7. The communication between ICL1 and W535 may trigger Smo activation (Arensdorf et al., 2016), since the activating Smo mutation W535L has been already known in BCCs (Xie et al., 1998;Lam et al., 1999). The free energy surfaces (FES) are shown in Figure 1. The FES along  Frontiers in Chemistry | www.frontiersin.org d1 and d2 achieves convergence in last 10 ns for each system (Supplementary Figure 1). As seen in Figure 1A, the FES of Smo-Vismod in the spatial distribution is separated obviously by four energy basins, marked as 1, 2, 3, and 4. While the FES of Smo-CLR owns two energy basins and a minor energy basin, marked as 1, 2, and 3 ( Figure 1B).  Frontiers in Chemistry | www.frontiersin.org Rana et al., 2013;Wang et al., 2014;Weierstall et al., 2014;Byrne et al., 2016;Huang et al., 2016;Zhang et al., 2017) (Supplementary Table 1), we speculated that the CRD groove of the inhibited state is open. ICL1 and TM7 are away from each other, with the corresponding d1 is 23 ± 3 Å and d2 is 10 ± 2 Å. The d1 and d2 shifting to 12 ± 2 Å and 5 ± 1 Å are considered as activated state of Smo, the CRD groove is in a closed conformation where can accommodate cholesterol, and ICL1 and TM7 remain communicating. The d1 varying between 17 ± 3 or 7 ± 2 Å indicates the intermediate states of Smo, meanwhile, the d2 varies between 7 ± 1 Å ( Table 1) (Figures 2b-e). Firstly, we studied the cholesterol binding site, a hydrophobic pocket formed by the hydrophobic residues in CRD (residues V107, L108, L112, I156, V157), HD (residue V210), TM6 (residues V488, L489) and ECL3 (residues V494, I496) (Supplementary Figure S2B). And the hydrophobic pocket is still maintained by cholesterol binding in Smo-CLR (Figure 2b). However, in the stable state of Smo-Vismod, the extracellular extension of TM6 and ECL3 occupy the CRD groove (Figure 2c) and the CRD groove site collapses. The CRD groove, HD, TM6 and ECL3 form strong hydrophobic interactions. Therefore, TM6 and ECL3 occupy the CRD groove in Smo-Vismod, which hinders the cholesterol binding to groove site.
The  Table 1) (Huang et al., 2016), exists neither in other CRD crystal structures nor in the multiple-domain crystal structures that resolved previously. Obviously, the CRD groove of inactivated Smo takes the open conformation, which accommodates TM6 and ECL3 by forming strong hydrophobic interactions. In the TMD-site (Figure 2e), the notable variations are the movements of TM5, TM6, and TM7. Compared with Smo-CLR, the TM6 of Smo-Vismod moves outward from the helix bundle, while TM5 and TM7 move inward (Figure 2e).

The Conformational Dynamics of the CRD Groove Site
The Smo-Vismod undergoes significant conformational changes during the simulation (Figures 1, 2). We analyzed the conformational dynamics process of Smo-Vismod from the FES (Figure 1). Based on the hydrophobicity of the CRD groove site, and the interaction among the hydrophobic residues of CRD groove, HD, TM6, and ECL3, the solvent-accessible surface (SASA) of the hydrophobic pocket composed of them is observed over time (Figures 3C,D), along with the evolution of d1 (Figures 3A,B) to study the conformational dynamics of the ECD of Smo-Vismod.
At the ∼30 ns of the start of the simulation, the d1 is <15 Å (Figure 3A), and the SASA of the hydrophobic pocket during the simulation is not significant change (Figure 3C), indicating that the CRD groove is still in the closed state (Supplementary Table 1). The corresponding energy basin is 4 ( Figure 1A). The representative conformation of the energy basin 4 was superimposed with the crystal structure of the Smo-Vismod ( Figure 4A). As seen that the CRD groove is in a closed state as in the initial state. The extracellular extension of TM6 and ECL3 do not occupy the CRD groove, nor form hydrophobic interaction. Compared to the initial structure, TM6 and ECL3 have begun to approach the groove.
The second obvious fluctuation is from ∼30 to ∼50 ns, the d1 increases rapidly and fluctuates in the range of ∼15 to ∼20 Å. The corresponding SASA of the hydrophobic pocket decreases rapidly, indicating that the CRD groove is beginning to open. The appropriate energy basin is 3 ( Figure 1A). As shown in Figure 4A, we found that TM6 and ECL3 enter into the CRD groove. The hydrophobic pockets formed among the hydrophobic residues of CRD, HD, TM6, and ECL3 begin to be closed. Several hydrophilic residues (E158-G162) of the side of CRD groove begin to leave the groove site and hence the CRD groove gets gradually open.
At the end of 50 ns, d1 increases and fluctuates between 20 and 25 Å during the subsequent simulation ( Figure 3A). The SASA of the hydrophobic pocket also fluctuates at a steady level in the last 50 ns (Figure 3C), indicating that the hydrophobic pocket maintains in compact contact. The representation of energy basins are 1 and 2 ( Figure 1A). Compared the representative structures of the two energy basins (Figure 4B), the TM6 and ECL3 have completely occupied the CRD groove and form strong hydrophobic interactions with the CRD groove. CRD groove presents a significant open conformation. With the TM6 and ECL3 of Smo-Vismod gradually approach to CRD groove, CRD is away from the membrane plane (Figure 4C), the HD moves toward the direction of TM6 (Figure 4D). The synergistic movement of the extracellular domain of Smo-Vismod contributes to stabilize CRD in inactivated conformation.
In contrast, the d1 smoothly fluctuates at the beginning of ∼25 ns; after ∼25ns, mainly fluctuates between 10 and 15 Å ( Figure 3B); the SASA shows increased at beginning ∼ 20 ns ( Figure 3D) in Smo-CLR. We compared the representative structures of Smo-CLR in the energy basins (Supplementary Figure S3). It is clear that CRD tilts toward the membrane plane during the simulations, the TM6 shifts toward the flank of CRD, simultaneously, the ECL3 deviates from CRD. And the HD moves away from TM6 (Supplementary Figure S3A). However, the CRD groove does not change significantly ( Figure S3B). Noticeably, the activated conformational dynamics of ECD is in contrast to the inhibited conformation of Smo-Vismod.
At this point, we have proved that the movement of TM6 and the hydrophobic residues of CRD groove, HD, TM6, and ECL3 play crucial roles in the deactivation of Smo by vismodegib binding. Indeed, the TM6 and ECL3 occupy the CRD groove in Smo-Vismod, which hinders the cholesterol binding to groove site.

The Conformational Dynamics of the TMD-Site
As observed above, the movement of TM6 is crucial in the conformational dynamics of the CRD groove site. Considering the binding mode of vismodegib that directly interacts with TM5, TM6, and TM7, we firstly studied the interaction between vismodegib and TMD-site. As shown in Figure 5A, we observed that the chlorophenyl-methylsulfone moiety of vismodegib overturns nearly 90 • in the stable inactivated state and closes to HD-ECL2 compared with the crystal structure. The atom O3 of the methylsulfone moiety forms H-bond with atom ND2 of side-chain of N219. The RMSD of vismodegib was calculated throughout the simulation (Figure 5B). We caught sight of the obvious remolding, that is, the RMSD suddenly increases from ∼0.4 to ∼1.4 Å at the beginning ∼5 ns, and fluctuates stably between 1.2 and 1.6 Å, indicating that the binding mode of vismodegib remodels immediately during the simulation. We also monitored the distance between atom O3 of the methylsulfone moiety and atom ND2 of N219 ( Figure 5C). The distance rapidly drops to ∼3 Å from around 5.2 Å at the beginning simulation of ∼5 ns and eventually runs aground. This means that the vismodegib rapidly leaves the initial binding mode shortly after the start of the simulation, afterwards, stabilizes in the new binding mode, forming steady H-bond with side-chain of N219, which also indicates the importance of HD.
At the other side of TMD-site, the charged residue R400 of TM5, D473 of TM6, and E518 of TM7 contribute to the binding of vismodegib Byrne et al., 2016). The movement of TM5, TM6, and TM7 is obvious in the stable states of Smo-Vismod and Smo-CLR (Figure 2e). Therefore, we analyzed the conformational dynamics of the three polar residues (Figure 6). As seen in Figure 6A, the distance between R400 and   E518 of Smo-Vismod is much smaller than Smo-CLR throughout the simulation. This indicates that the two residues going close to each other at the beginning of the simulation and form stable electrostatic interaction Smo-Vismod, which results in TM5 and TM7 moving inward (Figures 6B,F). TM6 is extruded from the helix bundle (Figure 6F). At the same time, the two residues going close to each other results in a remolding binding mode of vismodegib. Since in the crystal structure, the negatively charged E518 is near the amide oxygen atom of vismodegib with the distance of 3.1 Å from OE2 of E518 (Supplementary Figure S4), the electrostatic repulsion causes overturning of vismodegib.
In addition, due to R400 and E518 going close to each other, the electrostatic repulsion of the E518 and the electrostatic attraction of the R400 result in the D473 deflected (Figures 6C,E). The dihedral angles of the atoms C, CA, CB, and CG of the D473 of the Smo-Vismod is far less than the Smo-CLR (Figure 6C). D473 points toward R400 in the representative state of Smo-Vismod ( Figure 6E) instead of pointing toward E518 in the crystal structure (Supplementary Figure S4).
From these results, the remodeled binding model of vismodegib is stabilized by forming H-bonds with N219 of HD, D384, and Y394 of ECL2, respectively. Additionally, R400, E518, and D473 constitute the electrostatic interface to form polar interaction with amide linker of vismodegib (Figures 5A, 6E). R400, D473, and E518 play a vital role in the conformational dynamics of helices. Moreover, the outward movement of TM6 makes the extracellular extension of TM6 close to the CRD groove, and thus along with ECL3 interact with hydrophobic residues of CRD groove and HD, which hinders the binding of the cholesterol. On the contrary, cholesterol occupied the groove to push the ECL3 away from the groove and keep the TM6 out of the groove, allowing the CRD to tilt toward membrane plane.

The Conformational Dynamics of the ICD
Compared with d1, the fluctuation of d2 is not significant throughout the simulation in both of systems (Figures 3A,B). Nevertheless, three major energy basins of Smo-Vismod along d2, termed to 4, 2-3 and 1, were obtained ( Figure 1A). We The vismodegib binding to TMD site leads to the movement of TM6, which promotes the extracellular extension of TM6 and ECL3 to occupy the groove by forming strong hydrophobic interactions (Figure 2c). Hence the vismodegib binding blocks the binding of cholesterol. At the intracellular side, the ICL2 shows closed owing to the shifts of TM3. ICL3 is open caused by the shifts of TM5 and TM6. The ICL1 and W535 are away from each other due to the TM7 shifts outward. The interaction of the multidomain of Smo induced by vismodegib hinders the binding of the cholesterol, and destroys the intracellular interactional interface of the effector.

The Communication of Multi-Domain of Smo
We have made it clear so far that the antagonist vismodegib binding leads to the movements of TM5, TM6, and TM7, and then stabilizes Smo in an inactive state through coordinated movement between multiple domains. We further carried out the dynamical network analysis (Sethi et al., 2009)

CONCLUSIONS
In this work, the mechanism of the vismodegib allosterically inhibits Smo activation and hinders cholesterol binding is revealed by metadynamics simulations. We revealed that the vismodegib binding leads movements of TM5, TM6, and TM7, and the shift of TM6 triggers the entrance extracellular extension of TM6 and ECL3. Therefore, the TM5, TM6, and TM7 are key factors to the deactivation of Smo upon vismodegib binding. Moreover, we also found an inhibited open conformation of CRD groove, which is not shown in crystal structures. The open CRD groove accommodates the TM6 and ECL3 so that hinders cholesterol to bind and holds the CRD stacked atop the TMD. The strong hydrophobic interaction of CRD groove, HD, TM6 and ECL3 stabilizes the interaction of them. Furthermore, the HD and ECLs play a key role in the coordinated interaction of CRD, TMD and themselves. Therefore, blocking the coordinated movement of CRD, HD, and ECLs may potentially inhibit the activation of Smo. And we observed a remolding of vismodegib binding, stabilized by hydrogen bond formed with N219, D384 and Y394. These results can be taken into account for the design and discovery of novel Smo inhibitors, or providing structural information for discovering potential active sites.