Graph Neural Networks Bootstrapped for Synthetic Selection and Validation of Small Molecule Immunomodulators

Prageeth R. Wijewardhane,#,1 Krupal P. Jethava,#,1 Jonathan A. Fine,1 Gaurav Chopra*,1,2 1Department of Chemistry, Purdue University, 560 Oval Drive, West Lafayette, IN, USA 2Purdue Institute for Drug Discovery, Integrative Data Science Institute, Purdue Center for Cancer Research, Purdue Institute for Inflammation, Immunology and Infectious Disease, Purdue Institute for Integrative Neuroscience, West Lafayette, IN, USA


INTRODUCTION
Programmed cell death protein 1 (PD-1) is an immune checkpoint receptor implicated for the creation of new cancer therapeutics. 1 The prolonged interaction between the T-cell receptor and the major histocompatibility complex (MHC) leads to upregulation of PD-1 on the activated T-cell surface. 2 Activated T cells produce cytokines, such as Interferon-γ, which in turn cause tumor cells to express programmed death ligand 1 (PD-L1) on the their cell surface. 2 Tumors escape the action of immune system by utilizing the interaction between PD-1 and ligand PD-L1 resulting in lower effector T-cell function and survival, as such resulting in a suppressive immune response in the tumor microenvironment. 2 The inhibition of the PD-1/PD-L1 interaction can enhance anti-tumor immunity and a large amount of work has been done to develop monoclonal antibodies as inhibitors of PD-1/PD-L1 interaction inhibitors. 3,4 For example, pembrolizumab, cemiplimab, and nivolumab are three FDA approved anti-PD-1 monoclonal antibodies. 4 The discovery of small-molecule inhibitors would be an advantageous over monoclonal antibodies, such as being fast-acting, simple for in vivo administration, ability to penetrate through cell membranes and interact with the cytoplasmic domains of cell surface receptors. 5 Since a few years, there has been significant development in designing PD-1/PD-L1 inhibitors. 6,7 Specifically, Bristol-Myers Squibb (BMS) discovered a set of potent PD-1/PD-L1 small molecule inhibitors based on the peptidomimetic molecules and non-peptidic small molecules. 6,7 In particular, BMS revealed a 2-methyl-3-biphenyl-methanol scaffold containing chemical libraries. Later, Holak et al. studied the interaction of BMS molecules with PD-L1 suggesting that BMS molecules induce PD-L1 dimerization and also reported crystal structures of compounds with dimeric PD-L1. 8,9 Based on these findings, we envisioned to develop a machine learning (ML) framework for selecting and testing new PD-1/PD-L1 inhibitors.
Traditionally, the development of small-molecule inhibitors requires high throughput screening of a large library of diverse drug-like compounds 10 or a medicinal chemist iterating over a scaffold with weak receptor activity to enhance potency. 11 This entire process is -(i) time consuming; (ii) needs expensive instrumentation and robotics; (iii) based on trial-and-error; and (iv) highly inefficient to identify several new scaffolds rapidly. 12 In addition, virtual screening using docking methods have been developed to improve this process but with limited success. 13 Further, ML architectures such as Support Vector Machine (SVM) [14][15][16] , Random Forest (RF) [17][18][19] , Graph Convolution Network 20 , and Graph Neural Networks (GNN) 21,22 have been used for drug design and predicting drug-target interactions 23,24 . Recently, new architectures utilizing a combination of graph features in the binding site of a protein have shown great promise for calculating binding affinities and determining whether a compound will bind to a target. 20,22 Several new neural network-based architectures have also been proposed which promise to identify potent scaffolds, but many have not been tested experimentally, 15,16,[25][26][27][28] and developments in the ability to mine and characterize protein crystallography data hopes to drive the creation of these models. 29 Recently, it has been shown that molecular sub-graph features incorporated through a GNN and protein features encoded by their sequence can be combined to predict if a compound can target a given protein. 24 Inspired by this work and based on our interest in developing methods for drug design and immunology [29][30][31][32][33][34][35][36] , we have developed a new machine learning model to predict if a compound can inhibit the PD-1/PD-L1 interaction. Our method replaces the protein sequence features with docking scores representing the free energy of binding and due to this global energetic interaction of the small molecule in the binding pocket, we have termed this model as an "Energy Graph Neural Network" (EGNN). The three-dimensional atomic interaction energetic scores are calculated using CANDOCK 31 ( Figure 1B) and are combined with local molecular graph features ( Figure 1A) using an end-to-end training methodology (Figure 1C-D). In this work, we use this EGNN model to select designs for synthesis, and experimentally test a curated list of compounds from these predictions to prospectively identify potent PD-L1 small molecule inhibitors using the Homogenous Time-Resolved Fluorescence (HTRF) assay. We also tested negative predictions suggesting the utility of the model to be used for selecting potent leads as PD-1/PD-L1 inhibitors.

Patent Data for Training the EGNN Model
We used PD-1/PD-L1 small molecule inhibition data for 762 compounds from four patents to train our models: WO 2015/034820 A1 7 and WO 2015/160641 A2 37 by BMS (674 compounds),  A homogeneous time-resolved fluorescence (HTRF) binding assay was used to show activity against PD-1/PD-L1 interaction in the patents. However, the patents did not list individual IC50 values for all compounds but provided a range of inhibition with different molecules. Therefore, we trained a binary classifier with cutoffs for both datasets to treat a molecule as "High potency" or "Low potency" (Figure 1). If the reported IC50 of a molecule is less than or equal to 100 nM in the patent it was considered as a "High potency" molecule, otherwise it was considered as a "Low potency" molecule. This threshold was selected as it is the only common threshold value among four patents (Table S4). It should be noted that the actual value of IC50 should not be considered here as our experiments with multiple replicates were not able to obtain exactly reported results for some molecules in the patents (see IC50 value of compound 4a in Table 2, BMS-1 annotated with 6-100 nM in the WO 2015/034820 A1 patent 7 ) since our experimental conditions were different from the conditions reported in patents (See HTRF assay part in the experimental section). Therefore, we consider positive prediction (high potency) based on our experimental IC50 value as compared to the upper limit of a BMS control molecule (compound 4a/BMS-1) in WO 2015/034820 A1 patent 7 . The training dataset of 762 small molecules with the BMS or Incyte annotation is shown in Supporting Information File (TrainingData.xlsx).
We selected BMS and Incyte patents to include chemical diversity of the molecules in the training data set. Figure 2A shows the distribution of low and high potency molecules and general scaffolds in the BMS and Incyte patents. The BMS patents have 372 high potency compounds and 302 low potency compounds while the Incyte patents have 47 high potency compounds and 41 low potency compounds respectively. The BMS patent scaffolds contains 417 derivatives of (2-methyl-3-biphenylyl)methanol and 257 derivatives of [3-(2,3-dihydro-1,4-benzodioxin-6-yl)-2methylphenyl]methanol shown in Figure 2A (bottom -left) with R groups as CN, Cl, Br, and CH3.
On the other hand, Incyte patent scaffolds have distinct sub-scaffolds, denoted as A and B in Figure 2A (bottom -right). For Incyte scaffolds, X denotes for either N or C-R groups (R: Alkyl groups). These scaffolds suggest that the chemical diversity of Incyte compounds is higher than that of the BMS compounds because the general structures of the compounds in Incyte patents have more structural diversity for sub-scaffolds and atoms. We validated this observation using pairwise Tanimoto similarity scores of BMS and Incyte compounds as shown as heatmaps in Figure 2B and 2C, respectively. Morgan fingerprints with radius of 2 and bit length of 1024 were used to calculate pairwise Tanimoto similarities. High red color areas in the BMS heatmap indicates that the molecular pairs are structurally similar to each other. Low red areas in the Incyte heatmap suggests it has more chemical diversity in molecular structures. Furthermore, the average pairwise Tanimoto 40 similarity score of all BMS compounds was found to be 0.4434 and 0.3920 for all Incyte compounds, confirming higher chemical diversity in Incyte compounds when compared to BMS compounds.

PD-L1 homodimer and PD-1/PD-L1 Crystal Structures Reveals a Binding Site for Docking
It has been shown previously that BMS compounds inhibit the PD-1/PD-L1 interaction by inducing dimerization of PD-L1. 8,9 Therefore, a PD-L1 homodimer crystal structure (PDB ID: 5N2F) was selected for docking all the compounds in this manuscript. A PD-1/PD-L1 crystal structure (PDB ID: 4ZQK) was also used to check whether the binding site location of PD-L1 in the homodimer crystal structure (5N2F) overlapped and aligned with each other using the PyMol software package 41 ( Figure 3A). In Figure 3B, the selected binding site of the PD-L1 homodimer on the overlapped and aligned crystal structures is shown to indicate that the formation of the homodimer of PD-L1 with small molecules blocks the PD-1/PD-L1 interaction. A known inhibitor of the PD-1/PD-L1 interaction (ligand ID: 8HW) 8 in the selected binding site ( Figure 3C) suggests that the selected binding site corresponding to PD-L1 homodimers is relevant to develop PD-1/PD-L1 inhibitors. Therefore, the docking interactions of the PD-L1 homodimer will be relevant towards identifying PD-1/PD-L1 inhibitors. Also, direct docking with the PD-1/PD-L1 was not carried out since the binding site in between the PD-1 and the PD-L1 is filled with interacting amino acid residues from both proteins. Therefore, there is no space to dock a small compound with the PD-1/PD-L1 complex.   Table S1). The scoring function, radial cumulative complete 15 (RCC15) acquired the highest Cohen's kappa score of 0.41447. However, RCC15 scores were not able to clearly separate all the high and low potent classified molecules in the training data (see Violin plots in Figure S1). Using only one scoring function is not sufficient to capture the different states of PD-1/PD-L1 inhibition with small molecules. Therefore, we developed an EGNN  Training the EGNN Model). We calculated variation in the average F1 score (over five crossvalidated folds) versus the number of epochs for different hyperparameters ( Figure S2). Also, model's performance with different datasets with changing hyperparameters are shown in Table   S6. Optimal hyperparameters were selected to avoid overfitting and underfitting for EGNN.

EGNN Model with Hyperparameter Optimization Outperforms GNN and Other Baseline Models
Dimension of the hidden molecular vector (dim) = 10, sub-graph radius = 2, and number of hidden layers = 1 were selected as optimum hyperparameters for the EGNN model (see Experimental

Section on EGNN Training and Hyperparameter Optimization).
The EGNN and GNN models were trained with different training sets to examine the effect of chemical diversity on model performance for classification of high and low potency molecules.
Two datasets (BMS and Incyte) were used separately and in combination to train the EGNN model and determine the best dataset to predict PD-1/PD-L1 inhibitors. Splitting of the dataset into training-validation set and test set (4:1) were carried out using two different methods: (1) using a random splitter on shuffled data and (2) using a scaffold splitting method by DeepChem library. 43 Then training was carried out with fivefold cross validation and test sets were used to evaluate the models' prediction ability. Here, Cohen's kappa, F1 score and Area Under the Receiver Operator Characteristic Curves (AUROC) were measured to compare three models trained with BMS data only, Incyte data only and BMS-Incyte combined data. Further, as a separate experiment, all the measures were obtained for EGNN and GNN models trained only on BMS data while predicting for Incyte data, and vice versa as well. Figure 4A shows how data sets were split within the train-validation-test set scheme.
Initial dataset was split into two sets with the 4:1 ratio based on the scaffold splitting or random shuffled splitting. Then the 80% dataset was used as the training and validation dataset while the 20% dataset was used as the hold out test set to evaluate model performances.
GNN models with scaffold splitting appeared to generate comparable results with the EGNN ( Figure S4). However, this was expected since the graph neural network uses the twodimensional molecular framework/topology in training. When the framework distributions of the compounds are similar in the train-validation and test sets, GNN performs well. However, our intension is to develop a model which could be used to screen a large compound library which would not be necessary to share the same distribution of scaffolds with training sets (i.e. Incyte or BMS). Hence, we selected the random splitter with shuffling to create the test set for performance evaluations to develop a more generalized model.
Cohen's kappa scores of different test sets (hold out test set is based on random splitting) for both models trained with BMS compounds, Incyte compounds, and the union of these sets are shown in Figure 4B. The kappa scores of the EGNN and GNN models trained with Incyte data and tested on the hold out test set were 0.8861 and 0.4304, respectively ( Figure 4B). This result suggests that the EGNN trained model with Incyte data that contains diverse chemical scaffolds ( Figure 2C) performs much better than the GNN trained with the same data set. However, when the same test was done with only BMS compounds with lower chemical diversity than Incyte, the Cohen's kappa score is comparable for both models with 0.6416 for the EGNN model and 0.7164 for the GNN model. This suggests that the GNN model performs well with smaller chemical diversity in the training and test data as compared to larger chemical diversity. Both models show comparable performances with the combined datasets as well. When both BMS and Incyte datasets were combined, the kappa score for the hold out test set of the EGNN model was 0.6072 and 0.6729 for the GNN model. A similar trend is observed for the F1 scores for the three different training set comparisons ( Figure 4C). These results suggest that the EGNN model outperforms the GNN model for chemically diverse data sets such as Incyte data. We believe this to be due to the addition of 'global' energy features captured by the docking scores of PD-L1 homodimers as training data in EGNN compared to only the 'local' structural features of small molecules in the training data for the GNN model.
We also investigated the ability of the EGNN and GNN models trained on one compound set to predict high and low potency inhibitors of PD-1/PD-L1 in the other compound set. These results are represented in Figure 4B and 4C (Kappa scores and F1 scores respectively) with different bar patterns to represent different test sets. Tanimoto similarities between Incyte and BMS compounds are also shown in the Figure 4D heat map. The average pairwise Tanimoto similarity score of 0.3044 shows that compounds in these two datasets are very dissimilar to each other. When the EGNN and GNN models are trained on BMS compounds and used to predict Incyte compounds, a Cohen's kappa score of 0.1505 for the EGNN and 0.1200 for the GNN was observed and F1 score of 0.3810 and 0.2264 were observed, respectively. On the other hand, both F1 and kappa scores for both models improved when they were trained with Incyte data and used to predict the BMS compounds (kappa score of EGNN = 0.3852, GNN = 0.3196 and F1 score of EGNN = 0.7400, GNN = 0.6958). These results show that there is a marked improvement in F1 scores and Cohen's kappa scores for both EGNN and GNN models when trained on Incyte data and tested on BMS. However, AUROC score cannot distinguish these models correctly  Table 1 for all four models. The SVM model was trained using the 'svm' package in scikit-learn library 46 with the "linear" kernel and the RF was trained using the 'RandomForestClassifier' in scikit-learn library 46 with 500 trees. The 'metrics' module in scikitlearn package 46 was used for statistics AUROC, precision, recall, F1 score and Cohen's kappa.
Precision-recall curves for models and AUPRC values were obtained using the 'precrec' library 47 in R programming language. The EGNN model outperforms all the other models with values of 0.9250, 0.9212, 0.9091, 1.0000, 0.9524, and 0.8861 for AUROC, AUPRC, precision, recall, F1 score, and Cohen's kappa respectively ( Table 1). Comparing precision-recall curves of these four models ( Figure 4E) also confirmed that the EGNN model outperforms all three other models. Taken together, the combined local and global features in EGNN gives the best performance with the Incyte dataset. inhibitor predictions compared to other baseline models, such as, Random Forest, SVM, and GNN models.
All models were trained on the Incyte dataset and evaluated based on the same hold out test set.

Synthetic Selection and Validation of EGNN Predictions for PD-1/PD-L1 Inhibition
The EGNN model trained with optimum hyperparameters and the Incyte dataset was used to get predictions for an in-house database of small molecular designs. We developed a bootstrapped  (Table S5). Bootstrapping is an essential statistical technique that can be used to select confident molecules for synthesis and experimental validation based on agreements among multiple models. The bootstrapped EGNN model identified high and low potency small molecules as PD-1/PD-L1 inhibitors that were synthesized and then experimentally verified with HTRF binding assay (see Table 2 for summary). Specifically, we selected 4 molecules predicted to be high or low potent for PD-1/PD-L1 inhibition for testing based on bootstrapped EGNN SoftMax average scores and standard deviation (see EGNN SoftMax scores in Table 2).
Out of EGNN bootstrapped predictions, we have selected 1 molecule as highly potent   The D-serine end of the 4b compound forms hydrogen bonds with AThr20 and AAla121 along with a plausible hydrogen bond formation between backbone NH of ATyr123 and the oxygen in one of the two methoxy groups of the 4b molecule ( Figure 5C). These results suggest favorable interactions of compound 4b that could dimerize PD-L1 will result in PD-1/PD-L1 inhibition.
The HTRF assay confirmed that compound 4b has an IC50 of 339.9 nM (see Experimental Section for details) to inhibit PD-1/PD-L1 interaction ( Figure 5D). This is comparatively better  Table 2). The IC50 plots for each compound tested ( Figure S6) as well as the 13 C and 1 H NMR spectra are provided as Supporting Information. Pairwise Tanimoto similarity scores between these 4a-4e (   Our EGNN methodology can be further developed with the addition of more chemically diverse data, and incorporating reinforcement iterative learning with experiments performed in each step for developing a library of structurally diverse small molecule inhibiting PD-1/PD-L1 interaction to guide structure-activity relationships. Given the general nature of the machine learning model and docking methodology that is readily available for use, this approach can be adapted to identify small molecule immunomodulators by targeting other immune checkpoints, as well as, generally used to include local and global features for target-based drug design.

Compounds
Inhibition of PD-1/PD-L1 interaction was tested for 4 high and low potent predicted compounds using the PD1/PD-L1 HTRF assay kit from Cisbio US, Inc. The assay protocol was used as mentioned in the kit for each predicted compound (4b, 4c, 4d and 4e) and the BMS control compound (4a). Briefly, 2 µL of the compound, 4 µL from a 25 nM Tag1-PD-L1 protein solution and 4 µL from a 250 nM Tag2-PD1 protein were added into a Cisbio's HTRF 96-well low volume white plate. Then, the plate was incubated for 15 minutes at room temperature. Next, 10 µL from pre-mixed anti-tag detection reagents (5 µL from 1X anti-Tag1-Eu 3+ and 5 µL from 1X anti-Tag2-XL665) were added and the sealed plate was incubated for 2 hours at room temperature. Finally, the plate sealer was removed, and measurements were taken using a HTRF ® compatible reader. To calculate ∆ /∆ , first the HTRF ratio is calculated as follows; = 665 620 × 10000 A multiplication factor of 10000 factor was used to not deal with decimal values that improves data accuracy during calculation. The ΔR ratio indicating "specific signal" of the compound disrupting the PD-1/PD-L1 interaction was calculated by subtracting background HTRF ratio (negative DMSO control in our work) from each compound (sample) HTRF ratio as follows; Next, data normalization was done to minimize variation in values on different days, different plate reader instruments, or if the assay was done by different individuals. The normalization was done with respect to the background HTRF ratio and was calculated as follows; Finally, the ∆ /∆ ratio was calculated to enable comparison of values between multiple experiments.
where ∆ is taken as the ∆ of the positive DMSO control in the assay.

Calculation of IC50 values
The IC50 value for PD-1/PD-L1 inhibition was determined by analyzing the log of the concentration−response curves to fit a sigmoid curve with four-parameter logistic (4PL) regression using the GraphPad Prism Software version 8.3.0 for Windows, GraphPad Software, La Jolla California USA, www.graphpad.com. The IC50 values are provided in Table 2.

Machine Learning Architecture of the EGNN model
The EGNN model was developed using PyTorch 52 . All scripts for implementing the machine learning model and results are provided on GitHub at https://github.com/chopralab/egnn. The Figure 1 shows the overview of the EGNN machine learning architecture. We implemented the Graph Neural Networks for the molecular graph by Tsubaki and coworkers. 24 . Briefly, the molecular structures were converted into SMILES strings using ChemAxon MolConverter 53 software. Then RDKit 54 software package and the Weisfeiler-Lehman algorithm was used to extract r-radius subgraphs graphs for molecules ( Figure 1A). The following sections include details of the EGNN architecture.

Graph Neural Network for Molecular Graphs in EGNN
The following equations and notations with details for molecular GNN  Randomly initialized embeddings (Figure 1) are assigned to each r-radius edge /: (9) and vertex ( / (9) ) based on the type. Backpropagation has been used to train these random embeddings.

Vertex Transition Function
Say / (<) ∈ ℝ 5 is the embedded vector for the th vertex of a given molecular graph at time step . Then the updated / (<=;) ∈ ℝ 5 vector can be written as follows; Here, is the Rectified Linear Units (ReLU), a non-linear activation function such ( ) = (0, ). 83/%AB$9 ∈ ℝ 5×C5 and 83/%AB$9 ∈ ℝ 5×C5 are the weight matrix and the bias vector respectively. The vector between the th and th atoms (vertices) of the molecular graph after the time step is defined as /: (<) .
Moreover, / (<) and : (<) are added, because there is no direction for edges in molecular graphs.

Molecular Vector Output of Molecular GNN
The transition function generates an updated set of atom (vertex) vectors p = s ; (<) , C (<) , . . . , |E| (<) ur. Then the output function uses this set of atom vectors to obtain an unique molecular vector 6$03FG03 ∈ ℝ 5 (Figure 1A), which is defined as follows; Here, the total number of vertices in the full molecular graph is denoted by the | |.

Generation of Energy Features with Docking and Energy Vector (E) in EGNN
First, all the reported molecules were carefully drawn using MarvinSketch 56 software. Then, all the drawn molecules were cleaned in 3D and converted into a sybyl.mol2 file, which was used for docking with our in-house CANDOCK 31  Cohen's kappa score were selected and top in each class was selected for the EGNN model. Thus, the normalized docking score vector for each molecule in the EGNN model is represented using RCR15 and RCC15 normalized potential energy scoring functions as 3839%J ∈ ℝ C (Figure 1B).

Output of EGNN
As represented in Figure 1C, the normalized docking energy score vector ( 3839%J ) is concatenated with the molecular vector output of the GNN ( 6$03FG03 ). Then, the concatenated long vector ( 6$03FG03 ⊕ 3839%J ) ∈ ℝ (5=C) was used for the training as follows to obtain an output vector $G<2G< ∈ ℝ C ; Here ⊕ denotes concatenation, $G<2G< ∈ ℝ C×(5=C) denotes the weight matrix and the $G<2G< ∈ ℝ C denotes the bias vector. Then, a SoftMax classifier ( Figure 1D) is added on to the top of the $G<2G< = [ ) , ; ] vector to get the high or low potency probabilities.
Here, ∈ {0,1}; 0 indicates low potency and 1 indicates high potency, and the < is the probability of the given < .

Bootstrapping the EGNN model
Bootstrapping which carried out random sampling with replacement was used with the final model to get predictions. One hundred different models with different sampled training data (sample size of 60) were trained and predictions were obtained for an in-house molecular designs test set. Averaged SoftMax scores were used as the final prediction results of the bootstrapped model. Finally, molecules in the synthetic test set were classified as high Potency or low potency based on the averaged SoftMax score. If it is greater or equal to 0.5, it was considered as high Potency, else low potency (Figure 1). Thus, the EGNN model was trained with back propagation with given SMILES strings, the vectors of RCR15 and RCC15 scores generated by CANDOCK 31 and their high potency or low potency status with the PD-L1 protein. The trained model can be used to predict the probability of a given molecule to be a high or low potent molecule towards the PD-L1 protein.

EGNN Training and Hyperparameter Optimization
The model takes a SMILES string and a docking energy score string for a given molecule as inputs.
Hyperparameters of the model were optimized before using it for predictions. Dimension of the GNN hidden vector (dim), number of hidden layers of the GNN, and sub-graph radius were optimized by considering the five-fold cross validated F1 score. Three values were used for the dimension of the GNN hidden molecular vector output (i.e. dim = 5, 10 and 15). Numbers 1, 2, and 3 were used to check for the optimum number of hidden layers in the GNN. Finally, the optimum sub-graph radius for the model was selected out of radius = 1, 2 and 3.

Calculation of the F1 Score and the Cohen's Kappa
Following terms were used to calculate the Cohen's Kappa and the F1 score. The number of compounds predicted to be high potency that experimentally reported to be high potency was   N-(2,6-dimethoxy-4-((2-methyl-[1,1'-biphenyl]-3-yl)methoxy) benzyl)-3,3,3-trifluoro-1-phenylpropan-1-amine: Compound 3b from scheme 1 (8 mg, 0.022 mmol), 3,3,3-trifluoro-1-phenylpropan-1-amine (16.7 mg, 0.088 mmol, 4 equiv), sodium cyanoborohydride (7.2 mg, 0.114 mmol, 5.2 equiv), were dissolved in DMF (0.5 mL) and then added acetic acid (1 drop). The reaction mixture was allowed to stir at 80 o C for 3 hours. The reaction was monitored by TLC. The crude was purified by 0-20% MeOH:DCM to afford desire product as oily product (68% yield Tibolone (156 mg, 0.5 mmol) was taken in a round bottom flask containing 10 mL of THF and 100 µL of water was added. Next, p-toluene sulfonic acid (85 mg, 0.5 mmol) was added to it and the mixture was refluxed at 80 o C for 48 hours and the progress of the reaction was monitored by TLC. The organic solvent was then evaporated to dryness to get the crude product, which was purified by flash column chromatography using 20% ethyl acetate in pet-ether solvent mixture as eluent to give off-white solid pure compound GCL2 (53% yield). 1   MeOH to afford desire product as oily product (43% yield). 1