key: cord-0032038-pvofoi0h
authors: Weng, Ye; Pan, Chenghao; Shen, Zheyuan; Chen, Sikang; Xu, Lei; Dong, Xiaowu; Chen, Jing
title: Identification of Potential WSB1 Inhibitors by AlphaFold Modeling, Virtual Screening, and Molecular Dynamics Simulation Studies
date: 2022-05-13
journal: Evid Based Complement Alternat Med
DOI: 10.1155/2022/4629392
sha: fbd915837b652ab8f0b64a43d3b949b479f3b77c
doc_id: 32038
cord_uid: pvofoi0h

WD40 repeat and SOCS box containing 1 (WSB1) consists of seven WD40 repeat structural domains at the N-terminal end and one SOCS box structural domain at the C-terminal end. WSB1 promotes cancer progression by affecting the Von Hippel–Lindau tumor suppressor protein (pVHL) and upregulating hypoxia inducible factor-1α (HIF-1α) target gene expression. However, the crystal structure of WSB1 has not been reported, which is not beneficial to the research on WSB1 inhibitors. Therefore, we focused on specific small molecule inhibitors of WSB1. This study applied virtual screening and molecular dynamics simulations; finally, 20 compounds were obtained. Among them, compound G490-0341 showed the best stable structure and was a promising composite for further development of WSB1 inhibitors.

Metastasis is one of the causes of cancer death [1] [2] [3] . However, due to the lack of effective interventions, the study of new antimetastatic targets has become a popular research topic in oncology-related research areas. Tumor-derived blood vessels are unevenly distributed, and their function is abnormally compared to normal blood vessels, leading to the persistence of a hypoxic microenvironment in the tumor, which shows that tumor metastasis is closely related to hypoxia [4, 5] . HIF-l has been reported to regulate various related signaling pathways during the adaptation of tumor cells to the hypoxic environment and plays a crucial role in tumor cell proliferation, angiogenesis, malignant invasion, and metastasis [6] [7] [8] [9] . For example, high expression of HIF-l has been associated with local infiltration and distal metastasis of tumors, including colon, prostate, and lung cancers [10, 11] .

HIF-1 regulates WSB1 expression. WSB1 increases under hypoxic conditions, and it has been identified to be dependent on HIF-1 [12, 13] . WSB1 has been shown to regulate pVHL protein stability not only under hypoxia but also under normoxia. WSB1 can upregulate HIF-1 through ubiquitination of pVHL [14] [15] [16] . us, HIF-1 and WSB1 form a positive feedback loop, which provides a strong activation of HIF-1. Under hypoxic conditions, the HIF-l protein transcriptionally activates the E3 ubiquitin ligase family protein WSB1 and induces an increase in its protein expression.

ere are many research studies about WSB1 regulating tumor progression [17] [18] [19] [20] [21] . WSB1 was found to drive the metastatic potential in osteosarcoma cells which correlated with the pulmonary metastatic potential. WSB1 plays a role in neuroblastoma cell growth and involved in pancreatic cancer progression. Besides, it has also been reported to participate in the carcinogenesis of lung cancer.

Currently, studies on the crystal structure of WSB1 and the inhibition of WSB1 activity by gene silencing and gene mutation techniques are still in progress. ere are still many challenges for tumor control, and it is even more important for us to investigate small molecule inhibitors to support the development of subsequent new drugs.

However, in cancer drug discovery, the experimental tests used to examine small molecules are often expensive and laborious [22] . In recent years, artificial intelligence (AI) has provided new opportunities for drug discovery [23] [24] [25] .

In this study, we applied AlphaFold2 to predict the protein structure of WSB1. After that, molecular dynamics simulations were applied to optimize the structure as well as software to validate the accuracy of the modeled structure, followed by peptide-protein docking and structure-based virtual screening including AutoDock-GPU and Glide. Virtual screening (VS) is a powerful drug discovery tool that takes advantage of high-performance computers to filter compounds by using ligand or structure-based methods [26] [27] [28] . Finally, four compounds with different compound scaffolds were selected, namely, G490-0341, G610-0188, Y043-6168, and Y044-5019 compounds. Moreover, the binding mode of compound G490-0341 was investigated, which provided important information for further structural modifications.

ese results provide relevant information for the study of WSB1 inhibitor drugs and may assist in future drug design.

We obtained the AlphaFold model from their web page (https://deepmind.com/research) as an initial machine learning-based model. Our latest improved protocol based on molecular dynamics simulations was applied to the protein model. is approach is an improved version of the protocol we used previously during CASP13. e repaired protein structure was further processed by the Protein Preparation Wizard for modules including hydrogenation, redistribution of bond levels, readdition of side-chain residues, recovery of selenomethionine, and specifying conditions such as the protonated state pH of the protein.

Finally, molecular dynamics simulations of 200 ns were performed to eliminate the irrational conformation in the structure. To fix the irrational factors in the confirmation, we performed longer (200 ns) simulations at 300 K. 

Molecular docking was performed in Maestro using GlideSP. First, the protein was prepared for docking using an MOE's quick prep tool, which included correcting structural issues, protonating the structure, removing unbound water molecules, and minimizing the structure to a specified gradient to make the pocket available for docking of the new molecule. e original ligand (D2) was used to define the binding site of the WSB1 active pocket.

e LigPrep module was used to prepare the WSB1 inhibitor. After preparation, the liganddocking module was used to perform the docking work. In the docking parameters, the maximum output confirmation number was set to 20 to improve the accuracy of docking. e molecular conformation with the highest docking score was selected for analysis. e results can be displayed in a ligand interaction diagram.

e ChemDiv database is a commercial small molecule database from ChemDiv Inc. (ChemDiv) containing over 1.6 million compounds as screening libraries.

In recent years, there has been an increasing interest in applying various methods to improve its accuracy in molecular docking and eventually to increase the discriminative ability of molecular docking to efficiency. Here, we report a virtual screening campaign for WSB1 inhibitors on the ChemDiv library through a novel selection strategy. e virtual screening workflow includes AutoDock-GPU, HTVS, SP, and XP models.

First, the first round of screening of ChemDiv was performed by AutoDock-GPU, and the top 10,000 molecules were selected for scoring. After that, the screening was performed by Glide-dock-SP with 5 docking conformation parameters, and 15,000 molecules were obtained for the next round of Glide SP accuracy screening. To improve the screening accuracy, the docking program was reprepared by LigPrep for better adaptation to Glide SP; docking screening was performed by Glide SP for compounds and WDR5 proteins with the highest output confirmation set to 20, and the remaining parameters were kept constant. e final 20 top scoring molecules were obtained, and the combinatorial pose metadynamics analysis was performed.

e combination of WSB1 and D2 was studied by using three 10 ns independent mild metadynamics simulations of Desmond 39, version 2.3 (Schrödinger, LLC). Metadynamics simulations are a widely used enhanced sampling method for sampling the free energy landscape. Schrödinger's binding pose metadynamics (BPMD) is a variant of metadynamics that samples the motion of ligands in and around their binding modes to use metadynamics methods to determine the relative stability of different binding modes. In this simulation, the biased collective variable is defined as the distance between the center of mass of the ligand molecule and the ligand-binding residues, i.e., Arg174, Arg315, and TYR218, used for WSB1 binding simulations. e initial Gaussian peak height and goodness-of-fit parameters were set to 0.1 and 2.4 kcal/mol, respectively.

In this experiment, we first induced fit docking of the 20 molecules obtained from the virtual screening and then performed BPMD with the 20 poses obtained from the induced fit docking. e stability of the resulting poses was assessed based on the PoseScore. e PoseScore is the root mean square value of the atomic coordinates of the ligand relative to the initial ligand weight. A Pose-Score ≤2Å was considered stable (this RMSD value is often used as a threshold to define the success of prospective docking simulations). e results were analyzed by PersScore, which indicates the strength of the hydrogen bond formed between the ligand and the protein residues. If 60% of the hydrogen bonds were retained during the simulation (e.g., PersScore ≥0.6), this was considered as a sign of good persistence. Finally, compound G490-0341 showed the best stable structure and was obtained as a promising compound for further development of WSB1 inhibitors. 

3.1. AlphaFold2 Protein Structure Prediction. Jinxin Che et al. found that 5,6-bis (4-methoxy-3-methyl phenyl) pyridine-2-amine acted as a degradation agent to inhibit cancer cell metastasis [18] . However, research on inhibitors targeting the WSB1 axis is still ongoing, so we put more emphasis on drug discovery specifically for WSB1 to identify potential inhibitors. Currently, the crystal structure of WSB1 has not been reported. In this study, we applied α-folding deep learning algorithms to predict the structure of WSB1. AlphaFold2 is artificial intelligence for protein structure prediction, which is a new neural network-based model that can predict protein structures with atomic-level accuracy [29, 30] . In recent years, artificial intelligence and machine learning techniques have played a crucial role in drug discovery and development [31] [32] [33] .

e neural network of AlphaFold2 can predict the structure of a typical protein in minutes, as well as larger proteins, such as a protein containing 2180 amino acids without a homologous structure. e model provides accurate predictions of the reliability of its predictions on a per amino acid basis. For this experiment, the prediction model was downloaded from the AlphaFold protein structure database [34] .

We performed molecular dynamics (Figure 1(b) ) to enable the repair of the irrational factors in the structure of WSB1 predicted by AlphaFold [35] . After 200 ns MD simulations, we obtained the repaired structure as well as RMSD and RMSF. Figure 1(b) shows the results of the RMSD analysis of the WSB1 interface. WSB1 reached equilibrium after 25 ns simulations and reached a stable conformation below 1Å, which was acceptable. Meanwhile, the residues of WSB1 were relatively unstable and prone to displacement (Figure 1(c) ). e secondary structure of WSB1 (Figure 1(d) ) has 6.67% α-helix and 34.64% β-strand. 

It was reported that thyroid-hormone-activating type 2-iodothyronine deiodinase (D2) expression was associated with activation of serum thyroid hormone, (de)ubiquitinase ubiquitin-specific peptidase 33, WD repeat sequence, and SOCS boxcontaining 1 (WSB1), correlated ions of cytokine expression, and inflammatory pathways [36] . To further explore inhibitors targeting WSB1, we hypothesized that D2 docking sites with WSB1 were very promising docking sites. We performed protein-peptide docking of WSB1 and D2 peptides using Maestro and simulated the binding conformation of WSB1 and D2, which is shown in Figure 1 (e). WSB1 is also the E3 ubiquitin ligase for D2 [37] . To position D2 relative to the ECSWSB-1 complex, a new loop of 18 amino acids is identified in D2 by comparing a three-dimensional model of D2 with nonubiquitinated D1 and D3 enzymes. It is shown that a better binding position do exist between the 18 amino acids of D2 and WSB1, and we also find further details of the binding between D2 and WSB1. e loops of 18 amino acids are critical for WSB1 to recognize D2 [38, 39] . D2 is integrated into the model by positioning the 18 amino acid loops proximal to the D-A/B-C side of the WSB1 propeller.

In other WD40 propeller E3 ubiquitin ligases, substrate recognition relies on a "supersite" for protein-protein interactions and most commonly the permeation of the second position of each blade, which is often occupied by an arginine residue that interacts with the phosphate group in the substrate. Arginine occupies three sites in WSB1, namely, Arg174, Arg315, and TYR218 [40, 41] . e study confirmed that the WSB-1 WD40 propeller supersite was essential for D2 recognition. e present study restored the details of WSB1 and D2 binding and provided us with more highdefinition views of the two bindings.

Screening. Based on the above analysis, we conducted virtual screening from the ChemDiv database to find the inhibitors targeting WSB1. e ChemDiv database is a commercially available small molecule database containing over 1.58 million compounds that serves as a screening library reference [42] . e virtual screening workflow based on the ChemDiv database is shown in Figure 1 (f ). Data are shown in Table 1 .

Initially, we screened the ChemDiv compound library of 1.58 million molecules by AutoDock-GPU [43] and selected TRP_38  THR_39  VAL_40  ALA_41  TRP_111  ALA_114  LYS_122  ARG_174  ASP_175  TYR_218  SER_219  VAL_260  ALA_261  CYS_262  TYR_276  ARG_315  SER_316  VAL_317  GLY_354  LEU_355  CYS_356 CYS_357 ALA_358

(c) (d) Evidence-Based Complementary and Alternative Medicine 127,105 molecules with a score of ≤−10.0 for the next round of Glide SP accuracy screening [44] . After screening, we retained 15,297 molecules with scores ≤−5.5 and then performed a conformational restriction screen test (the highest output confirmation was set to 10, and the remaining parameters were kept constant). Finally, we obtained 278 molecules and proceeded to the next screening stage. In the XP mode [45] , compounds with Glide SP scores ≤−7.0 were docked to WSB1.

is final docking procedure narrowed down our group of compounds to 20 (Glide XP score ≤−8.0), after which we subjected these molecules (Figure 2 ) to combinatorial pose metadynamics analysis.

A theoretical strategy able to probe the conformational profile of ligands in the enzyme active site is very important. It is well known that, from a theoretical standpoint, molecular dynamics simulations can be used to evaluate the molecular flexibility of ligands and receptors; however, it is worth mentioning that some conformational changes occur in the time scale of only dozens of nanoseconds, which could compromise the MD simulation viability for virtual screening studies, for instance. In this regard, a theoretical strategy to select promising configurations from the MD simulation is crucial to determine the theoretical accuracy. Hence, great computational effort is necessary to carry out this kind of simulation. Aiming, then, to reduce the number of frames of MD simulations to rationalize the theoretical findings without loss of the relevant information from the simulation, new methods based on the statistical inefficiency such as principal component analysis and wavelet analysis for selecting MD conformations had been developed. Based on the above analysis, we selected binding pose metadynamics [46, 47] . Metadynamics can enhance sampling by allowing efficient back and forth motions across large-free energy barriers, using a more realistic heap of computational resources to sample the relative stability of different binding conformations produced by IDFD while still maintaining full atomic resolution [48, 49] . In contrast, Schrödinger's binding pose metadynamics (BPMD), a variant of TRP_38  THR_39  VAL_40  ALA_41  PHE_42  LEU_113  ALA_114  GLU_121  LYS_122  ARG_174  ASP_175  LEU_176  THR_177  TYR_218  SER_219  CYS_220  ALA_221  VAL_260  ALA_261  CYS_262  ASP_263  TYR_276  ARG_315  SER_316  VAL_317  SER_318  PHE_318  LEU_355  CYS_356  CYS_357 ALA_358 PHE_359

(c) (d) Evidence-Based Complementary and Alternative Medicine metadynamics that samples the motion of ligands in and around their binding modes, aims to use metadynamics methods to determine the relative stability of different binding modes. Compared to molecular dynamics, BPMD can help select the dominant molecule more accurately and efficiently. In this experiment, three kinetic runs were performed for each of the 20 candidate poses. During the metadynamics calculations of the initial structure, the average RMSD of the 20 candidate poses increased from the beginning to the end of the simulation time and hydrogen bonds were present for finite time of the simulation run. Table 1 summarizes the scores of the 20 candidate poses obtained from the scheme, including PoseScore, PersScore, and ComScore. e stability of the poses is assessed based on the PoseScore, which is the root mean square value of the atomic coordinates of the ligand relative to the initial ligand weight, and a PoseScore ≤2Å is considered stable. PersScore indicates the strength of the hydrogen bond form between the ligand and the protein residues. If 60% of the total hydrogen bonds are retained during the simulation (e.g., PersScore ≥0.6), it is considered a sign of good persistence. e docking of the ligand to the WSB1 receptor in pose 6 provided a good example of metadynamics isolating a native-like pose (RMSD of 0.929Å). is pose was ranked 6th by PersScore, and none of the top 5 poses by PersScore were native. Figure 3 shows the average RMSD estimate versus simulation time for all 20 candidate poses.

To explore the binding modes between the compounds and WSB1 protein, we selected four compounds with more potential in the metadynamics results, including G490-0341, G610-0188, Y043-6168, and TRP_38  ARG_174  ASP_175  TYR_218  SER_219  CYS_220  ALA_221  VAL_260  ALA_261  CYS_262  ASP_263  TYR_276  ARG_315  TRP_313  VAL_314  SER_316  VAL_317  LEU_355  CYS_356 CYS_357 ALA_358

(c) (d) Evidence-Based Complementary and Alternative Medicine Y044-5019, and analyzed them using molecular dynamics (MD) simulations and interaction decomposition. e binding modes of the four complexes are shown in Figures 4-7 . Compound G490-0341 formed three hydrogen bonds with residues such as Arg315, Tyr276, and Asp175. An aromatic interaction with residue TRP38 was observed. e reported binding sites for WSB1 might include Arg174, Arg315, and Tyr218. Like G490-0341, G610-0188 ( Figure 5 ), Y043-6168 (Figure 6 ), and Y044-5019 (Figure 7 ) also had hydrogen-bonding interactions with residues Asp175 and Arg315. Compound G610-0188, identical to compound G490-0341, had hydrogen bonds with residues Arg315, Tyr276, and Asp175. In addition, both Y043-6168 and compound Y044-5019 formed hydrogen bonds with Ser316 and Y043-6168 formed π-π interactions with residue Tyr28. erefore, it was highly likely that the relatively advantageous compound G490-0341 in BPMD forminghydrogen bonds through hydrophobic residues of amino acids such as Arg315, Tyr276, and Asp175 might be more stable and contribute significantly to the formation of the complexes.

Considering that WSB1 plays an important role in tumor metastasis and tumor cell proliferation, clinically relevant drugs targeting the WSB1 axis as inhibitors are still being investigated. In this study, we predicted the protein structure of WSB1 by Alphafold2 and then used structure-based virtual screening, including AutoDock-GPU and Glide, to select compounds. rough structure-based virtual screening of WSB1 inhibitors, we finally found four compounds with different compound scaffolds such as G490-0341, G610-0188, Y043-6168, and Y044-5019. e binding pose meta-kinetics showed that compound G490-0341 binds tightly to residues Asp175 and Arg315 and is more stable than other compounds. is study will contribute to the further development of WSB1 inhibitors and provide some valuable information for understanding the structure of WSB1 inhibitors. To further study the medicinal properties of these compounds, the ADME/TOX properties of G490-0341 were calculated [50] . e detailed results for the pharmacokinetic parameters and toxicity analyses are shown in Table S1 . Computational pharmacokinetics and toxicology studies on G490-0341 suggest that it can be used as a good starting point for further developing and designing new derivatives.

e data used to support the findings of this study are available from the corresponding author upon request.

e authors declare that they have no conflicts of interest.

Metastatic colonization by circulating tumour cells

e challenge of targeting metastasis

Molecular basis of metastasis

Role of hypoxia in cancer therapy by regulating the tumor microenvironment

Clinical targets for anti-metastasis therapy

Hypoxia promotes dissemination of multiple myeloma through acquisition of epithelial to mesenchymal transition-like features

Hypoxia: a key regulator of angiogenesis in cancer

Oligodendrocyteencoded HIF function couples' postnatal myelination and white matter angiogenesis

Autophagy promotes paclitaxel resistance of cervical cancer cells: involvement of Warburg effect activated hypoxia-induced factor 1-α-mediated signaling

e hypoxic tumor microenvironment: driving the tumorigenesIs of non-smallcell lung cancer

Study on the inhibition of hyperthermic CO 2 pneumoperitoneum on the proliferation and migration of colon cancer cells and its mechanism

An integrative genomics approach identifies hypoxia inducible factor-1 (HIF-1)-target genes that form the core response to hypoxia

HIF1 regulates WSB-1 expression to promote hypoxia-induced chemoresistance in hepatocellular carcinoma cells

Deubiquitylase OTUD6B governs pVHL stability in an enzyme-independent manner and suppresses hepatocellular carcinoma metastasis

WSB1 promotes tumor metastasis by inducing pVHL degradation

VHL, the story of a tumour suppressor gene

Hypoxia-induced WSB1 promotes the metastatic potential of osteosarcoma cells

Potential role of WSB1 isoforms in growth and survival of neuroblastoma cells

e WSB1 gene is involved in pancreatic cancer progression

WD repeat and SOCS box containing protein 2 in the proliferation, cycle progression, and migration of melanoma cells

Discovery of 5,6-Bis(4-methoxy-3-methylphenyl) pyridin-2-amine as a WSB1 degrader to inhibit cancer cell metastasis

Personal muta-nomes meet modern oncology drug discovery and precision health

A structure-based drug discovery paradigm

Artificial intelligence in COVID-19 drug repurposing

Converging blockchain and next-generation artificial intelligence technologies to decentralize and accelerate biomedical research and healthcare

Recent applications of deep learning and machine intelligence on in silico drug discovery: methods, tools and databases

Machine learning-based virtual screening and its applications to alzheimer's drug discovery: a review

A comprehensive review on promising anti-viral therapeutic candidates identified against main protease from SARS-CoV-2 through various computational methods

It will change everything': DeepMind's AI makes gigantic leap in solving protein structures

High-accuracy protein structures by combining machine-learning with physics-based refinement

Combined SVMbased and docking-based virtual screening for retrieving novel inhibitors of c-Met

Artificial intelligence in drug design

Machine learning in protein structure prediction

Highly accurate protein structure prediction with AlphaFold

An integrated computational approach of molecular dynamics simulations, receptor binding studies and pharmacophore mapping analysis in search of potent inhibitors against tuberculosis

Skeletal muscle deiodinase type 2 regulation during illness in mice

Role of the type 2 iodothyronine deiodinase (D2) in the control of thyroid hormone signaling

e Hedgehog-inducible ubiquitin ligase subunit WSB-1 modulates thyroid hormone activation and PTHrP secretion in the developing growth plate

e iodothyronine selenodeiodinases are thioredoxin-fold family proteins containing a glycoside hydrolase-clan GH-A-like structure

Structural basis for phospho-dependent substrate selection and orientation by the SCFCdc4 ubiquitin ligase

A ubiquitin ligase complex essential for the NFkappa B, Wnt/Wingless, and Hedgehog signaling pathways

Identification of phelligridinbased compounds as novel human CD73 inhibitors

Accelerating autodock4 with GPUs and gradient-based local search

Multiple ligand docking by glide: implications for virtual second-site screening

Virtual fragment docking by glide: a validation study on 190 protein-fragment complexes

A minimization principle for transition paths of maximum flux for collective variables

e role of the oximes HI-6 and HS-6 inside human acetylcholinesterase inhibited with nerve agents: a computational study

Exploring ligand stability in protein crystal structures using binding pose metadynamics

Mechanism of small molecules inhibiting activator protein-1 DNA binding probed with induced fit docking and metadynamics simulations

ADMETlab 2.0:an integrated online platform for accurate and comprehensive predictions of ADMET properties

Acknowledgments e authors are thankful to Schrodinger, Bangalore, India, for providing software support for this research work. Table S1 .

e ADME/TOX prediction of G490-0341. (Supplementary Materials) 

Ye Weng was involved in data curation, investigation, review, and editing. Chenghao Pan took part in guidance, modification, review, and editing. Zheyuan Shen and Sikang Chen were responsible for guidance and modification. Lei Xu was involved in modification. Xiaowu Dong performed design, guidance, and revision. Jing Chen was responsible for guidance and revision.