Discovery of Novel Inhibitors From Medicinal Plants for V-Domain Ig Suppressor of T-Cell Activation

V-domain Ig suppressor of T cell activation (VISTA) is an immune checkpoint and is a type I transmembrane protein. VISTA is linked to immunotherapy resistance, and it is a potential immune therapeutic target, especially for triple-negative breast cancer. It expresses at a high concentration in regulatory T cells and myeloid-derived suppressor cells, and its functional blockade is found to delay tumor growth. A useful medicinal plant database for drug designing (MPD3), which is a collection of phytochemicals from diverse plant families, was employed in virtual screening against VISTA to prioritize natural inhibitors against VISTA. Three compounds, Paratocarpin K (PubChem ID: 14187087), 3-(1H-Indol-3-yl)-2-(trimethylazaniumyl)propanoate (PubChem ID: 3861164), and 2-[(5-Benzyl-4-ethyl-1,2,4-triazol-3-yl)sulfanylmethyl]-5-methyl-1,3,4-oxadiazole (PubChem ID: 6494266), having binding energies stronger than −6 kcal/mol were found to have two common hydrogen bond interactions with VISTA active site residues: Arg54 and Arg127. The dynamics of the compound–VISTA complexes were further explored to infer binding stability of the systems. Results revealed that the compound 14187087 and 6494266 systems are highly stable with an average RMSD of 1.31 Å. Further affirmation on the results was achieved by running MM-GBSA on the MD simulation trajectories, which re-ranked 14187087 as the top-binder with a net binding energy value of −33.33 kcal/mol. In conclusion, the present study successfully predicted natural compounds that have the potential to block the function of VISTA and therefore can be utilized further in experimental studies to validate their real anti-VISTA activity.


INTRODUCTION
Immunotherapy has turned into an important pillar of cancer treatment due to the successful blocking of the programmed cell death protein 1 (PD-1) and its ligand-programmed death-ligand 1 (PD-L1) immune checkpoints. Immune checkpoint receptors control the duration and intensity of immune response by inhibiting T cell activation (Tang et al., 2018). Several immune checkpoint proteins have been discovered, such as PD-1/PD-L1, TIGIT, VISTA, cytotoxic T lymphocyte antigen-4 (CTLA-4), TIM3, BTLA, and LAG3. PD-1 inhibitors such as nivolumab, pembrolizumab, and cemiplimab and the human IgG1 k anti-CTLA-4 monoclonal antibody ipilimumab have been approved by the Food and Drug Administration (FDA). These approved drugs have become successful cancer therapies. However, the relatively low response rate of current immunotherapy drugs (less than 30%) is still a serious challenge, and efforts are needed to identify and overcome other immunosuppressive pathways (Ventola, 2017).
In the ever-expanding list of immune checkpoints, VISTA (V-domain immunoglobulin suppressor of T cell activation) is considered to be an important regulator of the immune system. VISTA immune checkpoint protein is a type 1 transmembrane protein that is encoded by the C10orf54 gene . VISTA is part of the B7 family consisting of a single extracellular N-terminal Ig-V domain, a stalk with approximately 30 amino acids, a transmembrane domain, and a cytoplasmic domain (Flies et al., 2011). The closest homolog of VISTA in the B7 family is PD-L1, sharing 23% sequence identity. VISTA is highly expressed in tumor-infiltrating lymphocytes. VISTA is also expressed in CD4 + and CD8 + cells, where it negatively regulates T cell responses (Borggrewe et al., 2018;ElTanbouly et al., 2019). It has also been observed that VISTA is highly expressed in breast cancer as compared to other cancer types, indicating that targeting VISTA may benefit breast cancer immunotherapy (Xie et al., 2020). Interestingly, expression of VISTA has also been observed in different cancer types such as breast invasive carcinoma (BRCA), invasive ductal carcinoma (IDC), bladder urothelial carcinoma (BLCA), colon adenocarcinoma (COAD), kidney chromophobe (KICH), lung squamous cell carcinoma (LUSC), uterine carcinosarcoma (UCS), and skin cutaneous melanoma (SKCM). Recently, it has been reported that VISTA is the acidic pH selective ligand of PSGL-1, which means that it may engender resistance to antitumor immune response (Johnston et al., 2019). Research on a variety of clinical samples, autoimmune disease models, and tumor models has shown that VISTA has a key regulatory effect on the immune system and has the potential to be used as a therapeutic or combined drug target. These findings indicated that the high expression of VISTA on tumor cells in about 20% of NSCLC specimens can prove the feasibility of targeting VISTA for cancer therapy (Cuzick et al., 2015). Clinical studies have shown that the expression of VISTA is upregulated in oral squamous cell carcinoma and gastric cancer (Böger et al., 2017;Wu et al., 2017). After ipilimumab therapy, the VISTA immune checkpoint has also increased in patients with prostate cancer (Gao et al., 2017;Kakavand et al., 2017). In addition, previous studies have shown that VISTA is highly expressed in the immune cell subsets of human pancreatic cancer patients (Xie et al., 2018b).
Currently, compound CA-170 is undergoing phase I clinical trial for advanced tumors and lymphoma. CA-170 exhibits powerful activity to stop the lymphocyte proliferation and effector functions inhibited by VISTA proteins. CA-170 also exhibits selectivity for other immune checkpoints such as CTLA4, BTLA, and LAG3. These nonclinical data provide a strong basis for the clinical development of CA-170 Wang et al., 2018;Blevins et al., 2019). In this study, we performed molecular docking to select natural drugable molecules from medicinal plants which may act as antagonists against VISTA. Molecular dynamics (MD) studies were carried out to further verify the binding of natural leads with VISTA protein.

Docking Studies
The human VISTA extracellular domain crystal structure present in the RCSB PDB database with the PDB ID: 6OIL was retrieved and processed in UCSF Chimera (Pettersen et al., 2004) for the molecular docking process. The structure was prepared first by removing co-crystalized water molecules and the NAG ligand, and then missing hydrogen atoms were added and minimized for energy via two algorithms, conjugate gradient and steepest descent, keeping the step size of 0.02 Å. Autodock Vina (Trott and Olson, 2010)

Analysis of Complex Dynamics Using MD Simulations
The FF14SB force field of the AMBER 18 (Case et al., 2012) molecular dynamics (MD) simulation package was used for preparation of protein parameters, while its GAFF force field was used for generating ligand parameters (Wang et al., 2004).
The whole system was solvated in the water box (TIP3P), considering the padding distance of 12 Å (Jorgensen et al., 1983). Particle mesh Ewald (PME) was employed for processing long-range electrostatic interactions (Darden et al., 1993), and for the nonbonded interactions, the distance cutoff was allowed to be 10 Å. The SHAKE algorithm was used to constrain the bonds involving hydrogen (Ryckaert et al., 1977). All the systems were subjected to energy minimization by running 1,000 steps of the steepest descent and conjugate gradient algorithms. Temperature of each system was equilibrated to 300 K using NVT for a time period of 20 picoseconds (ps), gradually. Afterward, the system equilibration was achieved using NPT ensemble. Finally, a production run of 50 ns was performed, and each trajectory was saved after every 2 fs. Root mean square deviation (RMSD) and root mean square fluctuation (RMSF) analysis of all trajectories was performed to check the system stability by using module CPPTRAJ (Roe and Cheatham, 2013).

Free Energies Estimation by AMBER MMPBSA.py
The MM-GBSA method in AMBER 18 was used to estimate free energies binding for complexes (Miller et al., 2012). 100 snapshots separated at equal intervals were collected from MD trajectories to carry out the binding free energy calculations. In MM-GBSA, estimation of the net binding free energy (ΔG bind ) is done as follows: In Equation 1, ΔG complex is complex free energy, ΔG receptor is receptor free energy, and ΔG ligand is ligand free energy. The free energy of the above terms can be gained by using the equations given below: In Equation 2, ΔG is the free energy. TΔS corresponds to entropy energy. In Eq. 3, the electrostatic interaction energy (ΔE elec ) and van der Waals interaction energy (ΔE vdw ) collectively correspond to the molecular mechanics energy in the gas phase (ΔG gas ). The polar contribution (ΔG GB ) and the nonpolar contribution (ΔG SA ) result in solvation free energy (ΔG sol ). The MM-PBSA method takes more time than MM-GBSA. Hou T et al. described that to calculate the relative ΔG bind , MM-GBSA is better in terms of result accuracy than MM-PBSA (Gohlke et al., 2004;Hou et al., 2011;Jyrkka€Rinne et al., 2012). This approach has been extensively employed in protein-protein interaction and protein-ligand binding studies (Alamri et al., 2021;Tahir ul Qamar et al., 2021).

Computational Prediction of Compound Pharmacokinetics
The selected compounds were also subjected to different predictions such as drug-likeness, lead-likeness, pharmacokinetics, medicinal chemistry, and toxicity to guide synthetic chemists in optimizing the structure to be successful in clinical studies. Computational predictions of the compound parameters as discussed above were done using the SWISSADME server (Daina et al., 2017).

Retrieval of Lead Compounds From MPD3 Database
The proposed research involves virtual screening of the MPD3 database against VISTA protein, followed by MD simulations and MM-GBSA methods. MPD3 is a collection of uniquely retrieved phytochemical compounds with reported therapeutic potential. The natural compounds were preferred because they are safer, possess better pharmacokinetics, and are easy to test in further experimental studies (Riaz et al., 2017). The lead-like compounds from MPD3 were considered to be therapeutically useful in the drug discovery process, as such compounds have improved selectivity, potency, and medicinal chemistry parameters (Hughes et al., 2011). Additionally, such compounds' structures can be easily optimized to get the desired biological activity. Previously, only Li et al. (2020) and Gabr and Gambhir (2020) reported small-molecule inhibitors against VISTA. Thus, lead-like natural compounds from MPD3 were retrieved ( Figure 1). In total, 1,634 molecules were able to fulfill the criteria of lead-like compounds. Theses 1,634 compounds were used for subsequent docking studies.

Molecular Docking of CA-170 Into VISTA Immune Checkpoint
The CA-170 small molecule has been reported as a dual inhibitor of PDL1/L2 and VISTA in order to treat advanced solid tumors and lymphomas. CA-170 is under phase II clinical trials for head and neck/oral cavity cancer, MSI-H positive cancers, lung cancer, and Hodgkin lymphoma in India. Its exact chemical structure has not been disclosed; however, some studies suggested that CA-170 is a peptidomimetic compound, composed of D-asparagine, L-serine, and L-threonine Musielak et al., 2019). Recently, the X-ray structure of the human VISTA extracellular domain has been solved at a resolution of 1.85 Å (PDB ID: 6OIL). VISTA is implicated in different cancers, including breast (Zong et al., 2020), skin (melanoma) (Kakavand et al., 2017), prostate (Gao et al., 2017), colon (Xie et al., 2018a), pancreatic (Xie et al., 2018b), ovarian (Mulati et al., 2019), and lung cancer (Villarroel-Espindola et al., 2018). A single-point mutation study to find the essential residues involved in the interaction of anti-VISTA antibody VSTB showed that three residues, Arg54, Phe62, and Gln63, are essential for the binding of VISTA to VSTB. The latter suggested that targeting these residues would be a valuable approach to inhibiting the VISTA immune checkpoint. In order to predict the binding pocket of CA-170 within the VISTA immune checkpoint, a flexible structure-based docking of CA-170 (PubChem ID: 123843830) using Autodock vina software was performed, following the same protocol as mentioned in the methodology. The grid box which represents the docking search area was centered to cover three key residues (Arg54, Phe62, and Gln63). Interestingly, the top pose of CA-170 with the lowest binding energy was forming hydrogen bonds with the Tyr41, Tyr37, Cys51, Ser52, and Arg127, including two crucial residues, Arg54 and Gln63 (Figure 2). Previous study indicated that a single-point mutation of Arg54 into Ala led to the abolition of the binding of anti-VISTA antibody VSTB to VISTA (Mehta et al., 2019). These results validated the docking protocol being applied in this study.

Virtual Screening
The molecular docking approach, one of the reliable approaches in the drug discovery process, was used to determine the natural inhibitors of VISTA protein. Docked ligands were graded based on their binding energy scores. The pose with the lowest score compared to CA-170 was regarded as the stable binding mode of the ligand. The top 20 compounds were visually analyzed using PyMol; out of those 20, three compounds were selected based on the binding conformation and interactions with the active site key  residues. These top selected natural ligands were successfully docked in the target active site. The ligand binding poses are depicted in Figure 3. All the compounds had at least one hydrogen bond with the critical active site residues. Among all three ligands' compounds, 14187087 has a greater number of hydrogen bonds with an energy value of −6.3 kcal/mol. It formed hydrogen bonds with Arg54, Gln63, His68, and Arg127 residues, out of which two residues (Arg54 and Gln63) are important active site residues. Compounds 3861164 and 6494266 formed two hydrogen bonds with Arg54 and Arg127 with the binding scores −6.8 kcal/mol and −6.7 kcal/mol, respectively ( Table 1). All the ligands have two common interactions with Arg54 and Arg127. As the compounds revealed favorable docking scores and good atomic-level chemical interactions, including hydrogen bonding with the VISTA, dynamics supported by binding free estimation were undertaken to further investigate the applicability of these compounds as effective VISTA inhibitors.

Molecular Dynamics Simulation Analysis of the Docked Complexes
All atom MD simulation was conducted using the AMBER package to assess the validity of the docking data and results by analyzing the dynamics behavior of protein atoms and the stability of the compounds at the binding site. For a time scale of 50 ns, the systems were analyzed for structure stability using RMSD and RMSF. The CPPTRAJ module of AMBER 18 was used to calculate the RMSD values to determine the convergence of the trajectories. RMSF values were calculated to determine the structure flexibilities of protein. Even though docking studies have been used effectively for calculating the ligand binding pose for several proteins, they failed to assess the ligand binding affinity (Cheng et al., 2012). During docking, proteins are treated as rigid molecules which do not consider the conformational changes that occur due to the ligand binding (Heitz and Van Mau, 2002). These conformational changes can be studied using molecular dynamics simulations. MD simulations have been extensively used to study the conformational changes in the protein-ligand interactions .
Compound 6494266 fluctuated up to 3.5 Å during the first 5 ns, but later, after 15 ns, it reached equilibrium. Among all the complexes, the 14187087 compound was the most stable complex throughout the simulation with an average RMSD of 1.31 Å. However, the 3861164 complex kept oscillating throughout the simulation, indicating that this complex might be unstable among all the complexes. Thereafter, the 14187087 and 6494266 complexes were stabilized and showed steady state dynamic behavior, as shown in Figure 4.
The variability in the conformation of trajectories can be monitored by calculating RMSF for individual atoms. In order to investigate and explore the conformational variability of each trajectory, RMSF of residues was plotted with respect to the residue number to show the local conformational changes for all the simulated complexes ( Figure 5). Among all the docked  complexes, compound 3861164 and compound 6494266 showed high fluctuations as compared to other systems, which is also consistent with the RMSD results. It can be concluded that the apo-VISTA structure, despite one large peak (40-52 amino acids), is highly stable compared to the VISTA-compound complexes. Conformational rearrangement of the loops than the rest of the protein in the presence of compounds has been previously reported and is linked to greater flexibility (Streaker and Beckett, 1999;Danielson and Lill, 2012). As VISTA protein has higher loop percentage and has a small size, upon ligand binding it is highly likely that loops may behave more dynamically. However, these fluctuations did not disturb the ligand binding conformation and the chemical interaction network, which are key to the stable binding of the compounds throughout the simulation time.

Analysis of Intermolecular Binding Stability by MM-GBSA
MM-GBSA binding free energies of the complexes were estimated to validate the docking and simulation results.
Such MM-GBSA binding free energy is now regularly applied in drug-designing protocols as they are more reliable than conventional docking techniques and less computationally expensive (Alamri et al., 2020a;Alamri et al., 2020b). The binding energies of complexes are presented in Table 2. It was observed that van der Waals energy and electrostatic energy dominated chemical interactions between the compounds and VISTA protein and contributed majorly to the total energy. In the interaction of 14187087, the van der Waals and electrostatic energy values were −32.2723 and −49.3294 kcal/mol, respectively, suggesting that electrostatic interactions were the major forces in the binding of VISTA and compound-1.
In the case of 3861164, the contribution of van der Waals energy was −21.4642 kcal/mol and that of electrostatic energy was −16.8891 kcal/mol. In the case of complex 6494266, van der Waals and electrostatic energy was −27.7207 and −37.0189 kcal/mol, respectively. Among all three complexes, 14187087 had the minimum binding energy, indicating it as an effective inhibitor.

Computational Prediction of Compound Pharmacokinetics
SwissADME is an online server for calculating different physical and chemical indicators and predicting drug-like properties, ADME parameters, pharmacochemical friendliness, and pharmacokinetic properties to help drug discovery. Detailed results of all the compounds are listed in Table 3. The oral bioavailability radar of the compounds is shown in Figure 6.
The physicochemical properties of the compounds are within the scope of drug-likeness and do not violate any Lipinski rule parameter. In addition, the compounds have good lipophilicity, so they can be transported to the maximum extent and reach the target site (Arnott and Planey, 2012). The compounds were also demonstrated to fulfill all requirements of the prominent Lipinski (Lipinski, 2004), Egan (Egan et al., 2000), Muegge (Muegge et al., 2001), and Veber (Veber et al., 2002) drug-ability rules. The compounds were predicted to be soluble and thus can be good candidates for oral administration. All the compounds were also predicted to not contain pan-assay interference compounds (PAINS) alerts and will not interact nonspecifically with multiple biological targets. This analysis revealed that the screened hits are VISTA-specific and will not have off-target effects. The compounds also have good gastrointestinal absorption, thus indicating that the good concentration of the drugs can reach the target site for performing the required action. Also, the compounds have good synthetic accessibility scores; therefore, they are easy to synthesize for experimental studies.
The compounds' docked conformation dynamics study validated their stable binding nature and compounds remained intact at the active site by both hydrophobic and hydrophilic interactions with key active residues of the protein. Additionally, confirmation on the binding stability of the compounds was accomplished through the binding free energy approach, which also revealed consistent results with the docking and MD simulation outcomes. The study employed a comprehensive computational framework to identify anticancer molecules by targeting VISTA protein. Although each step is thoroughly validated and the results are investigated for accuracy via follow-up computational approaches, the study suffers from lack of experimental validation. Altogether, the study findings are promising and could be subjected to further experimental evaluation to disclose their anti-VISTA/cancer potency.