Molecular mechanics and dynamics simulation of CD-47/SIRPα blockade study: A computational study on overcoming immunotherapeutic resistance in pancreatic ductal adenocarcinoma

2.1 Protein preparation

The three-dimensional X-ray crystallographic structure of a complex of human signal regulatory protein SIRPα complex with CD47 (PDB ID: 2JJS) was chosen. The protein was prepared using the Schrodinger suite's protein preparation package. The preparation included the assignment of hydrogen bonds, bond orders, hydrogen addition, optimization, protein minimization, removal of all chains excluding chain C (CD47 chain), and deletion of waters (Newman et al., 2003).

2.2 Grid generation

The receptor grid was created using the amino acids found on the surface of CD47 that are important for binding to SIRPα (Sandeep and Sobhia, 2018). GLU 35, THR 99, ARG 103, GLU 97, THR 102, GLU 100, GLU 104, LEU 101 and GLU 106 are among the amino residues. The coordinates for the X, Y, and Z axis are 28.98, −14.73, and 35.32 respectively.

2.3 Ligand preparation

Flavonoids with cancer-inhibitory properties from medicinal plants were retrived from literature. The flavonoids' three-dimensional structure was downloaded in Sdf format from Pubchem database. The compounds were prepared with the Maestro Schrodinger software's LigPrep module . Fig. 2 shows the initially screened top ranking ligand structures as well as their chemical names.

Close

Fig. 2

2D Structure and name of the selected flavonoids.

Export to PPT

2.4 Structure-based screening

GLIDE structure-based screening includes three precision methods: HTVS, standard precision (SP) docking, and XP docking. This study used two docking precisions to obtain a potential lead molecule quickly. HTVS can swiftly screen a vast number of molecules; however, the sampling techniques were constrained, making it difficult to interpret the results. As a result, SP was used to dock all the ligands (16) with high glide score from the HTVS precision analysis, which chooses an appropriate binding pose from a broad pool of ligands.

2.5 Generation of pharmacophore hypothesis

The standard drug (NCGC00138783) was used to generate an energy optimized pharmacophore model for the crystal structure of CD47. The model was generated from protein-ligand option of the phase develop pharmacophore tool of schrodinger suite (2017–1) (Omoboyowa et al., 2023). The pharmacophore model is shown in Fig. 3.

Close

Fig. 3

The pharmacophore hypothesis generated with the standard drug.

Export to PPT

The virtual screening base on E-pharmacophore was performed with the five (5) flavonoids with top docking scores after preparation using macro model minimization. The pharmacophore-based analysis was carried out with phase module to generate a subset of molecules having chemical features for binding to CD47 according to the generated model. The fitness scores were used to justify the best hits.

2.6 Development of automated QSAR model

By blasting the FASTA sequence of the protein received from PDB, the protein inhibitors with their corresponding IC50 were extracted from the CHEMBL database (https://www.ebi.ac.uk/chembl/) with the chain C (CD47 chain) sequence (shown below) used for the blasting. The inhibitory chemicals were translated to sdf format using the Data-Warrior software (v.2) (Omoboyowa et al., 2022). The sdf format of the inhibitors was imported into the Schrodinger suite workspace and prepared using the macro model minimization tool. The QSAR model for the protein was developed. The best-projected rank was kpls molprint 24, with the model's prediction precision measured by the ranking score, RMSE, SD, Q2, and R2. The pIC50 of a lead compound was predicted using this approach.

2.7 CD47 FASTA sequence for CHEMBL blast

>2JJS_2|Chains C|LEUKOCYTE SURFACE ANTIGEN CD47|HOMO SAPIENS (9606).

QLLFNKTKSVEFTFGNDTVVIPCFVTNMEAQNTTEVYVKWKFKGRDIYTFDGALNKSTVPTDFSSAKIEVSQLLKGDASLKMDKSDAVSHTGNYTCEVTELTREGETIIELKYRVVSWSTRHHHHHH.

2.8 Free binding energy calculation (MM-GBSA)

The docking results showed that the selected ligands bind to the protein's active site via the receptor grid. However, it was unclear if this binding would be sufficient to cause a biological response as this depends primarily on the free binding energies of the protein–ligand complex. The binding free strengths of the top-ranking compounds and the reference drug were determined using the MM-GBSA module of the Schrodinger software.

2.9 Quantum chemical methods

Theoretical calculations have been widely studied through a powerful tool denoted as the Density Functional Theory method (DFT) due to its structural and spectral explanation of organic molecules. All calculations were performed and computed by spartan 14 software by wavefunction Inc on the top five hit compounds. The study was carried out with complete optimization of all geometrical variables via 6-31G* basis set, and this was accomplished with B3LYP density functionals (Huang et al., 2021). Frontier molecular orbitals (FMOs) energy was calculated, i.e., for the top five hit compounds identified, which by calculation deduce the energy band gap (Eq (1), which according to Koopman's theorem, predicts the interacting center.

2.10 Insilco ADMET prediction

We were interested in the drug-likeness properties of the top five hits identified through molecular docking. This was accomplished with the help of the QikProp module. Then, additional ADME properties were obtained from the online web server AdmetSar (https://lmmd.ecust.edu.cn/). The ligands' chemical structures were provided. The structures were translated to their canonical simplified molecular-input line-entry system (SMILE) to determine the ligand's physicochemical properties and pharmacokinetic models (Jensen, 2001). The drug-likeness of can be used to establish if a ligand is suitable for oral administration. The in silico prediction is based on Lipinski's rules (molecular weight, hydrogen bond donor, and hydrogen bond acceptor) (Cheng et al., 2012).

2.11 Molecular dynamic (MD) simulation

After the docking screening, the top ligand–protein complexes were further submitted to molecular dynamics simulation. Other approaches are needed for validation because the dynamics of the complexes were not considered during the molecular docking procedure. Using molecular dynamic simulations, the stability of docked Protein-Ligand complexes was evaluated. LiGRO (Yao et al., 2017) was used to set up the simulation system, and a GUI-based tool produced the system file needed to execute the simulation. The target was dissolved in a cubic box filled with TIP3P molecules of water with 150 mM NaCl ion concentration. The target was parameterized using the Amber99sb force field, and the ligand molecules were parameterized using LiGRO's ACPYPE module. The study by Omoboyowa et al. served as the guide for preparing all other simulation systems and parameters. GROMACS 5.1.5 was used to generate a run for each complex system that lasted 100 ns (Omoboyowa et al., 2022).

2.12 MD simulation trajectories analysis

Using conventional GROMACS tools, the trajectories of the MD simulation were evaluated for the root mean square fluctuation, radius of gyration (RG), H-bond mapping and root mean square deviation. The generated trajectories were shown in the PyMOL visualization graphics system (Version 2.0 Schrödinger, LLC.). Using the Molecular-dynamics-Interaction-plot tool, the interaction proportions of the target residues interacting with the compounds were computed (Jiang et al., 2021).

3.1 Molecular docking, MMGBSA and interaction profiling

The protein CD47 was docked by the selected flavonoid derivatives to predict binding energy and interaction with the protein. Before docking was carried out, the protein was analyzed for its residues, which were vital for interacting with its partner protein CD-47.

Pelargonidin and (+)-Gallocatechin exhibited the highest binding energies of −6.715 kcal/mol and −6.353 kcal/mol respectively compared to the standard drug, which has a binding energy of −3.445 kcal/mol (Fig. 4). Thus, the docking results were analyzed, and it was finally reported that among the top flavonoid compounds, Pelargonidin and (+)-Gallocatechin exhibit the best binding interaction, essential in identifying and developing new therapeutics targeting CD47. This will block its interaction with CD47 as a strategy to prevent cancer. The 2D interaction diagram shows the residues involved in the ligand's binding (Fig. 5) and the 3D binding complexes were presented in Fig. 6. Pelargonidin, the top-ranked ligand, formed H-bonds with THR 34, GLU 97, LEU 101, and GLU 104. While the other ligands including the standard drug, interacted via H-Bond with GLU 104, a critical amino residue in the CD47 binding site. Malvidin and Peonidin had another unique interaction with TYR 37, called Pi-Cation. VAL 36, TYR 37, ALA 53, and LEU 101 were the same interacting hydrophobic amino acids in the top-ranked ligand and the standard drug. Additionally, all reported ligands, including the standard drug, have similar interactions with the following key hydrophobic amino acids: TYR 37, ALA 53, and LEU 101 as presented in Table 1. The formation of H-bond and other hydrophobic interactions between small molecules and the amino acid residues at the binding domain of the target is vital for their inhibitory potential (Omoboyowa, 2022). The Schrodinger suite's Prime module's MM-GBSA approach was used to calculate the binding energies of the top compounds with the highest docking scores. The lower the score, the higher the binding energy, this method provide a reliable statistical post-docking examination of docked complexes. The relative free binding energies of pelargonidin, (+)-gallocatechin, malvidin, peonidin, and catechin are −29.25, −15.87, −29.59, −24.92, and −16.00, respectively. Pelargonidin, Malvidin, and Peonidin exhibit higher binding energies than the reference drug, per the results of the MM-GBSA (−22.31).

Close

Fig. 4

Graph depicting the docking and MM-GBSA scores for the top five hits including the standard drug.

Export to PPT

Close

Fig. 5

2D molecular interaction SIRPα binding pocket with hit compounds. (a) Pelargonidin (b) (+)-Gallocatechin (c) Malvidin (d) Peonidin (e) Catechin (f) NCGC00138783(Standard Drug).

Export to PPT

Close

Fig. 6

3D representation of hits with binding pocket of SIRPα.

Export to PPT

3.2 Virtual screening using E-pharmacophore model

Ligand-based pharmacophore approach is a vital computational model for drug design without macromolecular protein structure. This hypothesis is an ensemble of electronic and steric characters which are vital in ensuring molecular interactions with specific biological molecules and to stimulate or inhibit signaling pathways of such protein (Yang, 2010; Omoboyowa et al., 2023). In this study, Herein, the e-pharmacophore model was developed based on the standard drug-CD47 complex using four partitioning features. Fig. 3 showed the generated hypothesis with the standard drug. The model with the best fitness score has one hydrogen bond acceptor, two hydrogen bond donors, one aromatic ring and hydrophobic interaction.

The fitness scores of the five (5) hits and standard drug are shown in Table 2. All the hit compounds showed fitness score greater than 1.0 with (+)-gallocatechin and catechin having the highest score of 1.191 followed by Catechin (1.531). Although the standard drug showed higher fitness scores higher (1.530).

3.3 Auto QSAR analysis

As it uncovers the correlations between the structural characteristics of chemical compounds and their biological activities, the quantitative structure–activity relationship (QSAR) is a computational model significant in drug development (Omoboyowa et al., 2023). AutoQSAR is a machine-learning approach that generates streams of independent variable-building models with various topological and physiochemical descriptors (Dixon et al., 2016). The autoQSAR model divides the test compounds into a 25% test and a 75% train set, as shown in Table 5. The best predictive model from the experimental data was determined using the based partial least-squares regression (kpls) analysis: kpls molprint2D 29. The model parameters resulted in a standard deviation (S.D.) of 0.4354, R2 of 0.8774, RMSE of 0.4220, and Q2 of 0.8706 (Table 3). The experimental compounds' model showed more training sets (blue colour) than the test set (red colour) as observed in the scatter plot in Fig. 7. The distribution of the test and train set in Fig. 7 was consistent with the distribution observed in Table 5. The QSAR model was use to obtain the pIC50 of the hit compounds shown in Table 5, All the hit compounds and standard drug showed high pIC50 above 5.00 μM with pelargonidin and peonidin (6.001 μM) showing better pIC50 among the hit compounds comparable with the standard drug (6.841 μM).(See Table 4)

Close

Fig. 7

Scatter plot of the best model.

Export to PPT

3.4 Quantum calculations

The development of density functional theory (DFT) of small molecules is an important tool use to describe the molecular interactions between molecules and gives information concerning electron transfer within the molecule which is required to predict the chemical stability and reactivity of a molecule (Balogun et al., 2021). The energy of high occupied molecular orbital (E_HOMO) and low unoccupied molecular orbital (E_LUMO) and energy gap (Eg) estimated from the quantum mechanics calculations are vital in predicting molecular reactivity of bioactive compounds (Omoboyowa et al., 2023). From the results presented in table Fig. 8, the complete geometry optimization of the selected flavonoid derivatives at their low energy level, in which all analysis were computed on these least optimized flavonoid molecules.

Close

Fig. 8

HOMO and LUMO density for the five flavonoid derivatives.

Export to PPT

From the results presented in Table 6 and Fig. 8, the electron donation potential of the hit compounds was suggested by the E_HOMO values which range from −5.63 to −9.25 eV and the electron accepting potential was predicted by the E_LUMO values which range from −6.12 –n0.08 eV. Hence, the high E_HOMO value and lower E_LUMO value are necessary for the predicted reactivity of the compounds. Energy gap (Eg) is the difference between the HOMO and LUMO energies and has been reported by Uzzaman and Mahmud (2020), as a predictor of the molecule stability and reactivity, higher value of Eg denote greater stability with less bioavailability and low reactivity (Omoboyowa et al., 2023). The Eg of the hits were estimated between the range of 2.56 to 5.71 eV suggesting that the hit compounds are reactive and stable molecules.(See Table 7 and Table 8.

3.5 Prediction of physicochemical and ADMET-TOX properties

The in silico drug-likeness predictions are founded on Lipinski's rule of five; hydrogen bond donor, hydrogen bond acceptor and molecular weight (Cheng et al., 2012). The draggability of a molecule is generally based on this rule. To ascertain if substances could penetrate the central nervous system, the blood–brain barrier penetration was evaluated (Table 5) (Karelson et al., 1996). Out of the five predicted molecules, only (+)-Gallocatechin and Catechin are likely to penetrate the blood–brain barrier. In humans, many cytochrome P450s catalyze the metabolism of various substances, including xenobiotics and medicines.

Thus, inhibition of cytochrome P450 isoforms may result in drug-drug interactions in which co-administered drugs fail to be metabolized and accumulate to toxic levels. The predicted compounds are CYP2D6 inhibitors and CYP2C9 substrates. All are predicted to be CYP3A4 substrates except for pelargonidin. Adverse drug administration-related effects are referred to as acute oral toxicity (Nyandoro, et al., 2018). Fortunately, all substances proved negative for acute toxicity and mutagenicity in the AMES test. Based on their acute oral toxicity, compounds are categorized into four classes. LD50 values for substances in Category I are less than or equal to 50 mg/kg. LD50 values for substances in Category II are higher than 50 mg/kg but lower than 500 mg/kg. LD50 values for substances in Category III are higher than 500 mg/kg but lower than 5000 mg/kg. The LD50 values of the substances in Category IV were higher than 5000 mg/kg. It is predicted that all the compounds will fall into Category III or IV. All of the compounds are thought to be non-carcinogenic. Drug solubility has previously been defined as a critical property of evaluating pharmaceuticals in the drug development cycle. This is because it helps to determine the concentration of the drug in the systemic circulation, resulting in a maximal optimal response (Walum, 1998). As a result of the presence of hydroxyl groups in the compounds, they were water-soluble. None of the substances is P-glycoprotein substrates (P-GB). The potassium channel, the human Ether-a-go-go Related Gene (hERG), is required for cardiac excitability control and regular heartbeat (Savjani et al., 2012). Since none of the chemicals inhibits the hERG gene, they cannot induce proarrhythmia. Pelargonidin, the lead compound, has been shown to have antioxidant properties after reducing oxidative markers. A research found the antioxidant efficacy of Pelargonidin via an increase in the level of the natural antioxidant enzymes catalase and superoxide dismutase (Lamothe et al., 2016; Mirshekar et al., 2010). The therapeutic potential of the screened compounds was determined using Christopher Lipinski's proposed rules of five (ROF) with a molecular weight of 500, HB donors of 5, HB acceptors of 10, and an octanol/water partition coefficient (log p 5) (Table 3) (Cheng et al., 2012). Except for gallocatechin, which possesses more than five HB donors, all the compounds met all the criteria.

3.6 Molecular dynamics simulation

Following the results generated from the molecular docking campaign, we employed the use of the following molecular dynamics parameters; RMSD, RMSF, HBOND and ROG in a 100 ns simulation to evaluate the dynamical stability of the three selected hit druglike candidates (Pelargonidin, (+)-Gallocatechin and Malvidin) and the reference drug (NCGC00138783). Their behavioral kinetics was also evaluated in both bound and unbound states with the protein structure (2JJS). Each complex was confined to an environment closely related to a normal physiological condition of a natural cell in terms of temperature, solvent, pressure, and ions throughout the simulation. To this end, we evaluated and compared the 100 ns spectrums of each prospective drug candidate complex against the reference drug (NCGC00138783).

3.7 Root mean square deviation

Using RMSD analysis, it is possible to determine the conformational and structural alterations of the backbone atoms of the protein–ligand entities (2JJS-NCGC00138783, 2JJS-GALLOCATECHIN, 2JJS-MALVIDIN, 2JJS-PELARGONIDIN and 2JJS). By comparison of RMSD spectrum of the referenced inhibitor and the unbound protein with the proposed hit compounds, one could speculate if the binding of the chemical entities is stable and capable of binding stably at the binding pocket of the receptor (Roy et al., 2017). Hence, we calculated and plotted the RMSD of each simulated complex for the entire 100 ns MD production run (Fig. 9). Looking at the RMSD spectrum, Pelargonidin and Malvidin depicted almost the same graphical pattern with the standard compound throughout the entire simulation suggesting their potential to act as a putative binder of the protein target. This is evidenced with their close average RMSD of 0.1899 nm and 0.186 nm respectively, compared with the 0.175 nm of the referenced compound. In contrast, Gallocatechin showed a different RMSD pattern compared to the standard drug and the apoprotein (2JJS). However, with mean RMSD value of 0.30 nm, the ligand average RMSD value falls below 0.5 nm value which is acceptable for a considerable stable system (Adelusi et al., 2020). Overall, the RMSD result showed that the binding of Malvidin, Pelargonidin, and Gallocatechin with the protein target are stable and suggest their ability to disrupt the SIRPα-CD47 signaling cascade.

Close

Fig. 9

The RMSD of the backbone atom of the protein in top docking scored hit complexes throughout 100 ns simulation.

Export to PPT

3.8 Root mean square fluctuation of the top docking scored hit complexes

For examining the structural fluctuation of the amino acid residues of the bound and unbound complexes, as well as changes in the position of the ligand, the Root Mean Square Fluctuation (RMSF) is a commonly used metric. A higher average RMSF value denotes a system with more fluctuating residues while a lower average value corresponds to a system with less fluctuating residues (Roy et al., 2017). According to the result represented in Fig. 12, all the bounded system including the standard (2JJS-GALLOCATECHIN, 2JJS-MALVIDIN, 2JJS-PELARGONIDIN, and 2JJS-NCGC00138783) demonstrate similar RMSF pattern. In comparison with the mean RMSF of the unbound receptor (0.199 nm), 2JJS-GALLOCATECHIN, 2JJS-MALVIDIN, 2JJS-PELARGONIDIN, and 2JJS-NCGC00138783, have lower average RMSF values of 0.146 nm, 0.135 nm, 0.152 nm, and 0.138 nm respectively. Interestingly, the RMSF result is consistent with the RMSD analysis which revealed that the binding of the ligands do not disrupt the conformational dynamics of the protein by depicting a lower average RMSD (with the exception of Gallocatechin) and RMSF values compared with the unbound form of the receptor (2JJS) (Fig. 10).

Close

Fig. 10

The RMSF of the residues present in the protein (2JJS) of top docking scored hit complexes throughout the simulation period.

Export to PPT

3.9 Hydrogen bond mapping of top docking scored hit complexes

Intermolecular H-bond is an important type of interaction that exists between protein–ligand complexes. Higher H-bond could be responsible for the greater stability of a compound at the binding pocket of the receptor while lower H-bond may indicate a lower stable system (Lee et al., 2012). In the H-bond spectrum depicted in Fig. 11, Pelargonidin-one of the promising drug candidates showed thicker H-bond spectrum compared to all other binary complexes including the referenced inhibitor. Statistically, Pelargonidin, Malvidin, Gallocatechin, and the standard drug averaged 1.85, 0.44, 0.46, and 0.53 H-bonds respectively with their targeted receptor. Although Malvidin, and Gallocatechin have lower intermolecular H-bond, however, the difference between their values with the standard drug is not statistically significant and hence they could also be a good inhibitor of the protein using the H-bond metric.

Close

Fig. 11

H-bond distribution of the top docking scored hit complexes during 100 ns simulation run in the active site of the protein target.

Export to PPT

Close

Fig. 12

Radius of gyration spectrum of the top docking scored hit complexes during 100 ns simulation run.

Export to PPT

3.10 Radius of gyration

The compactness of the secondary structures of the bound and unbound form of a simulated receptor is measured by radius of gyration (ROG). From this, one could determine if a system is stably folded or not. Fig. 12 shows the ROG plots of the apoprotein (2JJS) and the complexes (2JJS-NCGC00138783, 2JJS-Gallocatechin, 2JJS-Malvidin, 2JJS-Pelargonidin) as a function of time. Like the ROG result, Pelargonidin and Malvidin depicted similar ROG pattern with the standard drug throughout the entire simulation time. Their average ROG values are 1.364 nm, 1.370 nm, and 1.369 nm respectively. In contrast, Gallocatechin demonstrated different ROG pattern with the standard drug, but it has the lowest ROG value (1.337 nm) when compared with all other bounded systems. Our overall analysis of ROG indicates that the binding of the ligands does not distort the structural compactness of the receptor as they all averaged lower ROG vales when compared with the apoprotein (1.376 nm) as shown in Fig. 12.