Targeting programmed death ligand 1 for anticancer therapy using computational drug repurposing and molecular simulations

Abstract

Discovering new drug candidates for complex diseases like cancer is a significant challenge in modern drug discovery. Drug repurposing provides a cost-effective and time-efficient strategy to identify existing drugs for novel therapeutic targets. Here, we exploited an integrated in-silico approach to identify repurposed drugs that could inhibit programmed death-ligand 1 (PD-L1). PD-L1 is a crucial protein that plays a pivotal role in immune checkpoint regulation, making it a potential target for cancer treatment. Using a drug repurposing approach, we combined molecular docking and molecular dynamics (MD) simulations to study the binding efficiency of FDA-approved drug molecules targeting PD-L1. From the binding affinities and interaction analysis of the first screening, several molecules emerged as PD-L1 binders. Two of them, Lumacaftor and Vedaprofen, showed appropriate drug profiles and biological activities and stood out as highly potent binding partners of the PD-L1. MD simulation was performed for 500 ns to assess the conformational and stability changes of PD-L1-Lumacaftor and PD-L1-Vedaprofen complexes. The simulations revealed sustained structural integrity and stable binding of both complexes throughout the 500 ns trajectories, supporting their potential as PD-L1 inhibitors. While the findings are promising, they remain computational and require experimental validation to confirm biological efficacy and specificity. This study also emphasizes the role of bioinformatics approaches in drug repurposing that can help in the identification of novel anticancer agents.

Introduction

The conventional model of drug discovery, where the emphasis is made on the identification of small organic molecules with distinct structures and modes of action, has been associated with various problems, including high levels of failure, long development cycles, and increasing costs¹. To avoid these challenges, drug repurposing has been considered as one of the effective solutions². This strategy builds on the fact that the safety, pharmacokinetic, and pharmacodynamic properties of approved drugs are already known and hence transforms them into new therapeutic uses³. It also enhances the conversion of potential molecules into useful therapies and minimizes resource wastage in addressing unmet medical needs more effectively. Computational biology and bioinformatics have become an integral part of the drug discovery process⁴. The application of bioinformatics approaches, including molecular docking and molecular dynamics (MD) simulation, has provided a solid groundwork for the structure-guided drug repurposing strategies⁵.

Programmed death-ligand 1 (PD-L1) is a critical immune checkpoint protein that helps cancer cells evade immune surveillance⁶. PD-L1 interacts with its receptor PD-1 on T-cells and prevents T-cell proliferation and cytokine release, leading to immune tolerance of tumor cells⁷. This has been achieved by the blockade of the interaction between PD-1/PD-L1, an aspect that has been proven to boost the capacity of the immune system to detect cancer cells, enhancing PD-L1 as a more viable anticancer therapy target⁸. Currently, available PD-L1 inhibitors like Atezolizumab and Avelumab have been proven to be effective in clinical settings, especially in cancers like melanoma, non-small cell lung cancer, and renal cell carcinoma⁹. However, these therapies are not so effective and can be accompanied by side effects and resistance mechanisms¹⁰. Thus, more effective and safer PD-L1 inhibitors are urgently needed.

In this context, the present work intends to fill the gap between bioinformatics approaches and therapeutic possibilities. To find and assess potential repurposed drugs as PD-L1 inhibitors, we used a rational and systematic bioinformatics strategy. Repurposing leverages the established safety and pharmacokinetic profiles of approved drugs, thus reducing clinical development time and cost. Herein, molecular docking and MD simulations are employed as computational tools to estimate the binding efficiency and interaction modes of repurposed drugs with the PD-L1 protein, as well as the dynamic behaviors of the protein-drug complexes. This strategy is intended to identify the molecules with the affinity to interact with the binding site of PD-L1 and affect its function in disease-associated pathways.

First, molecules were preselected according to the binding affinities and the results of interaction studies. Out of all the screened molecules, a few of them proved to be potential molecules due to high binding affinity and high binding efficiency to the PD-L1 binding pocket. We then proceeded to perform a detailed assessment of selected molecules for their ability to inhibit PD-L1, stressing their biological properties and the fact that they preferentially bind with the binding site of PD-L1. To supplement the structural information, we used the MD simulation approach to dissect the dynamic aspects of PD-L1 and its docked complexes with the screened drugs. These simulations, performed up to a 500 ns trajectory, provided information on the stability and the interactive behavior of the complexes. Besides, the present work has also emphasized the role of bioinformatics techniques in drug discovery. Overall, the study highlights repurposed drugs as potential PD-L1 inhibitors that have significance for therapeutic applications in multifaceted diseases like cancer.

Materials and methods

Data preparation and bioinformatics resources

To identify candidate inhibitors for PD-L1, we used an integrated in-silico pipeline that entails several computational resources. These were MGL AutoDock Tools version 1.5.7 (https://ccsb.scripps.edu/mgltools)¹¹ InstaDock 1.2 (https://hassanlab.org/instadock/)¹² Discovery Studio Visualizer version 2023 (https://discover.3ds.com/discovery-studio-visualizer)¹³ and GROMACS 2020β (https://www.gromacs.org/)¹⁴ for virtual screening and MD simulation applications. The PD-L1 structure was downloaded from the RCSB Protein Data Bank (PDB ID: 8ALX)¹⁵. The structure was inspected for missing atoms, water molecules, and heteroatoms. The next step of preparation included the removal of water molecules and other unnecessary non-protein coordinates. Further, energy minimization was performed using the Swiss-PDB viewer tool (https://spdbv.unil.ch/)¹⁶. We appended hydrogens and set up the protonation states to the physiological pH using MGL AutoDock Tools. The drug data set used in our study included all 3,648 FDA-approved drugs that were retrieved from DrugBank¹⁷. To this dataset, we applied additional filtering to remove biologics and peptides to arrive at a list of FDA-approved drug molecules. All the molecules were obtained in the SDF 3D structure format for the next phase of molecular docking. The conformations of the chosen FDA-approved molecules were submitted to energy minimization and all-atom simulations with PPD-L1 through the GROMOS 54A7 force field¹⁸ in GROMACS¹⁴.

Molecular docking screening

Molecular docking screening was used to filter molecules based on their binding affinities to PD-L1. In the preparation of the protein structure, steps such as the addition of polar hydrogens and the charging of the structure were done. The InstaDock v1.2 graphical user interface program was used for the docking screening¹². We employed blind docking dimensions to generate a grid that extended over the whole structure of PD-L1. The dimensions of the grid were set at X = 55 Å, Y = 63 Å, and Z = 48 Å. The center coordinates were set as X = − 14.656 Å, Y = − 6.229 Å, and Z = 10.739 Å. We used a grid spacing of 1 Å across the entire surface of the structure. As for the other parameters for the molecular docking procedure, all of them were left unchanged from the default settings. To validate the molecular docking protocol, a retrospective redocking analysis was conducted. The co-crystallized ligand 8HW (4-[[4-[[3-(2,3-dihydro-1,4-benzodioxin-6-yl)−2-methyl-phenyl]methoxy]−2,5-bis(fluoranyl)phenyl]methylamino]−3-oxidanylidene-butanoic acid) was removed from the PD-L1 binding site (PDB ID: 5N2F) and subsequently redocked into the same site. The resulting pose closely aligned with the original crystallographic conformation, yielding a root-mean-square deviation (RMSD) of only 0.172 Å (Supplementary Figure S1). This strong agreement validates the accuracy and reliability of the docking protocol in reproducing experimentally determined ligand orientations within the PD-L1 binding pocket.

Scoring and interaction analysis

The results of the docking process were analyzed based on the docking scores and ligand efficiency, which is the combined measure of multiple energy terms used to establish the binding affinity¹⁹. The molecules from the top with higher docking scores were further subjected to interaction analysis. The most promising candidates were compared for the interactions with the key residues in the PD-L1 binding site using PyMOL (version 3.1; https://www.pymol.org/)²⁰ and Discovery Studio Visualizer¹³. After these analyses, we selected only the FDA-approved drugs that had the most favorable interactions with PD-L1.

PASS analysis

The molecules that were selected in the previous step were further evaluated for their drug profiles and pharmacological activities. This step sought to establish whether they could prevent or treat diseases related to PD-L1. The biological properties of the selected molecules were analyzed with the help of the PASS web server²¹. This tool helps forecast possible biological activities resulting from the chemical structure of compounds. It uses 2D molecular fragment descriptors linking the biological effects of a compound and its molecular structure. The PASS web server calculates a predictive score for several biological properties, which is defined as the ratio of the probability to be active (Pa) to the probability to be inactive (Pi). It is thus easier to predict that a compound with a higher Pa value will have a certain biological activity.

Molecular dynamics simulations

All-atom MD simulations were performed on the PD-L1 protein and its complexes with the selected molecules. This simulation was performed at 300 K employing the GROMOS 54A7 force field in GROMACS 2020β¹⁴. The ligand’s molecular topology parameters were prepared using the PRODRG server and then incorporated into protein topology files to build the protein-ligand complexes. Every system was placed in a periodic water box before equilibration, and counterions were added to balance the system’s charges. All three systems were placed in cubic simulation boxes in the presence of SPC216 solvent²². Energy minimization was done with 1500 steps of the steepest descent method²³. Bonds involving hydrogen were constrained using the LINCS algorithm, allowing a 2 fs integration time step. Nonbonded interactions were computed with a 1.0 nm cutoff, and long-range electrostatics were treated with the Particle-Mesh Ewald (PME) method. The equilibration phase was carried out at constant volume, with periodic boundary conditions at a pressure of 1 bar, at which the temperature of the system was gradually raised from 0 to 300 K in 1000 ps. After equilibration, all three systems were followed by 500 ns production MD simulations. The trajectories obtained were further analyzed with the help of the tools available in GROMACS to assess the stability of the PD-L1-ligand complexes during the simulation. In this analysis, attention was paid to the change and dynamics of load-bearing structures of complexes and their interaction with the environment throughout the simulation time. For reproducibility, an additional independent 500 ns MD simulation (replicate) was conducted for each system using a different random seed.

Principal component analysis and free energy landscapes

To analyze the conformational flexibility, atomic movements, and structural stability of PD-L1 and its complexes, we have used principal component analysis (PCA) by employing essential dynamics²⁴. We also generated free energy landscapes (FELs) of PD-L1 and its complexes with the screened molecules using a conformational sampling approach²⁵. It enabled us to sample the protein conformations close to the native state conformations. The FELs were generated to evaluate the stability and native configurations of PD-L1 before and after binding with the elucidated compounds.

Results and discussion

Molecular docking screening

The study was initiated with a structure-based molecular docking screening of 3648 FDA-approved drug molecules against PD-L1. While exploiting this molecular docking procedure, we were able to rapidly assess the entire dataset and select the subset of compounds with reasonably favorable binding affinities to PD-L1. Then, we evaluated all the top-ranked molecules for ligand efficiency to select the most efficient PD-L1 binders. From this analysis, five molecules were chosen from the initial pool with high binding scores varying between − 10.2 kcal/mol and − 10.7 kcal/mol, a ligand efficiency of > 0.3 kcal/mol/non-H atom to PD-L1 (Table 1). Each of these candidates demonstrated stronger binding affinities toward PD-L1 compared to the co-crystallized reference inhibitor 8HW²⁶which had a binding score of − 6.0 kcal/mol. The results indicate that these molecules might have great potential to develop high-affinity PD-L1 inhibitors.

Table 1. Top five docked FDA-approved molecules targeting PD-L1. Docking results include binding affinity (in kcal/mol), predicted pKi, ligand efficiency (normalized per non-hydrogen atom), and torsional energy. All scores were computed using InstaDock v1.2 against the PD-L1 crystal structure (PDB ID: 8ALX). The reference inhibitor 8HW (PDB ID: 5N2F) is listed for comparison.

S. no.	Drug molecule	Binding affinity (kcal/mol)	pKi	Ligand efficiency (kcal/mol/non-H atom)	Torsional energy
1.	Lumacaftor	−10.7	7.85	0.3242	1.8678
2.	Vedaprofen	−10.6	7.77	0.5048	1.2452
3.	Androstenedione	−10.4	7.63	0.4952	0
4.	Meclizine	−10.3	7.55	0.3679	1.5565
5.	Etoricoxib	−10.2	7.48	0.425	0.9339
6.	8HW	−6.0	4.4	0.1667	3.113

Drug profiling and biological activity evaluation

The hit molecules from the molecular docking screening were investigated for their drug profiles and possible biological activities concerning diseases associated with PD-L1. The study pinpointed two molecules, Lumacaftor and Vedaprofen, based on their drug profiles and PASS analysis, where they showed anticancer activities with similar attributes. Lumacaftor and Vedaprofen were subjected to a detailed pharmacological analysis to determine their suitability for PD-L1 inhibition; they possess drug-like properties appropriate for therapeutic use. Both molecules exhibited the expected antineoplastic effect, Pa varying from 0.971 to 0.199, with higher values of Pa meaning a higher chance of activity (Table 2). The reference compound, 8HW, further corroborated the inhibitory activity against PD-L1. This analysis highlights the high potential of the selected drug molecules in drug repurposing against cancer.

Table 2. Predicted biological activities of selected drug candidates via PASS analysis. PASS output includes the probability of activity (Pa) and inactivity (Pi) for predicted Pharmacological functions.

S. no.	Drug	Pa	Pi	Biological activity
1.	Lumacaftor	0,902	0,001	Cystic fibrosis treatment
		0,681	0,003	Protein kinase (Mos, Mil/Raf, MEKK, RIPK, TESK, LIMK, IRAK, ILK, Activin/TGF-beta) inhibitor
		0,398	0,009	Antineoplastic (brain cancer)
		0,379	0,042	Angiogenesis inhibitor
		0,362	0,116	Antiinflammatory
2.	Vedaprofen	0,853	0,005	Antiinflammatory
		0,533	0,024	MMP9 expression inhibitor
		0,479	0,004	Antiinflammatory, ophthalmic
		0,494	0,028	AR expression inhibitor
		0,530	0,076	TP53 expression enhancer
3.	Androstenedione	0,985	0,001	Testosterone 17beta-dehydrogenase (NADP+) inhibitor
		0,937	0,002	Lysase inhibitor
		0,830	0,007	Respiratory analeptic
		0,691	0,011	Antisecretoric
		0,682	0,004	Contraceptive
4.	Meclizine	0,917	0,004	Antieczematic
		0,721	0,008	Analgesic
		0,610	0,011	Antipsychotic
		0,478	0,024	Antiobesity
		0,481	0,032	Antiallergic
5.	Etoricoxib	0,825	0,006	Antiarthritic
		0,670	0,020	Antiinflammatory
		0,568	0,004	Diabetic neuropathy treatment
		0,573	0,016	Analgesic, non-opioid
		0,409	0,027	Alzheimer’s disease treatment
6.	8HW	0,314	0,004	Hydroxymethylglutaryl-CoA reductase inhibitor
		0,410	0,182	Phosphatase inhibitor
		0,338	0,130	Antiinflammatory
		0,141	0,051	Antineoplastic (carcinoma)
		0,130	0,073	Antineoplastic (cervical cancer)

Interaction analysis

The molecules selected from the PASS analysis based on the drug profiles and biological activities were studied for their interaction patterns with PD-L1. While exploiting PyMOL and Discovery Studio Visualizer, we examined detailed interactions of Lumacaftor and Vedaprofen with PD-L1. The docked conformations from these hits and the reference PD-L1 inhibitor were selected based on their specific interaction with the binding site residues of PD-L1. Lumacaftor and Vedaprofen were found to interact with key residues within the binding pocket of PD-L1 through specific hydrogen bonding and hydrophobic interactions (Fig. 1). The specific interactions of Lumacaftor and Vedaprofen with PD-L1 are depicted in Fig. 1A. Both molecules showed high affinity to bind to the binding site of PD-L1, where most of the co-crystalized inhibitors are bound²⁶. Both molecules showed various common residues and interacted directly with the binding site residues as the reference inhibitor (Fig. 1B). They complemented PD-L1’s binding pocket well and had favorable structural overlap (Fig. 1C). Considering the orientation of Lumacaftor and Vedaprofen with 8HW, it can be concluded that these molecules fit into the PD-L1 binding site and can be developed as repurposed high-affinity binding inhibitors against PD-L1.

Fig. 1 [Images not available. See PDF.]

3D binding orientation of docked ligands within the PD-L1 binding site. (A) PD-L1 in secondary structure color code bound with Lumacaftor (orange) and Vedaprofen (magenta) along with the reference inhibitor 8HW (cyan). (B) Binding pocket residues of PD-L1 occupied by the selected molecules and 8HW. (C) Surface view of PD-L1’s ligand-binding site indicating spatial occupancy and complementarity of the ligands. Structural visualization was performed in PyMOL.

Further, we examined the detailed molecular interactions of these ligands bound to the binding pocket of PD-L1, especially with the functionally important amino acids of the protein. Lumacaftor binds to critical residues in the binding pocket of PD-L1, including hydrogen bonds with Asn63, Gln66, and Asp73, together with other binding pocket residues (Fig. 2A). Likewise, Vedaprofen interacts with Tyr56 as Pi-Pi stacked, Ile54 and Met115 as Alkyl bonds, along with the crucial residues of the PD-L1 binding pocket (Fig. 2B). Both molecules were found to have similar binding patterns with the reference inhibitor 8HW, as shown in Fig. 2C. Considering these findings, it can be concluded that the binding pocket Lumacaftor and Vedaprofen to the crucial residues of the PD-L1 binding site may result in inhibition of its function and can be exploited for repurposed drug development.

Fig. 2 [Images not available. See PDF.]

2D interaction diagrams of ligand-PD-L1 complexes. (A) Lumacaftor–PD-L1 interaction map showing hydrogen bonding with residues Asn63, Gln66, Asp73 and hydrophobic contacts. (B) Vedaprofen–PD-L1 interaction showing π-π stacking with Tyr56 and alkyl interactions with Ile54 and Met115. (C) Reference compound 8HW’s interaction profile is also shown for comparison. Diagrams were generated using Discovery Studio Visualizer to highlight key residues and interaction types.

MD simulations

MD simulations serve as a powerful tool in drug discovery, which provides detailed insights into the dynamic behavior of proteins and protein-ligand complexes^27,28. In this study, we aimed to analyze the stability and dynamics of PD-L1 and the PD-L1-Lumacaftor and PD-L1-Vedaprofen complexes and thus performed all-atom MD simulations for 500 ns (Table 3). To ensure the robustness of the results, we also ran an independent replicate simulation (500 ns) for each system with a different random seed for initial velocities. The replicate runs produced comparable structural stability and convergence (Supplementary Figures S2-S3). A detailed description of the time-evolution analysis of different structural dynamics and stability parameters of PD-L1 and its complexes with Lumacaftor and Vedaprofen is given below.

Structural deviations and compactness

To evaluate structural changes and dynamics at the protein level, we employed the RMSD criterion²⁹. We analyzed time-dependent RMSD fluctuations for PD-L1, PD-L1-Lumacaftor, and PD-L1-Vedaprofen complexes and calculated their average values: 0.14 nm, 0.16 nm, and 0.18 nm, respectively. The average RMSD values of both PD-L1-Lumacaftor and PD-L1-Vedaprofen complexes were slightly increased and were stable in the trajectory as compared to free PD-L1. RMSD analysis (Fig. 3A) revealed minor fluctuations (~ 0.25 nm), indicating structural convergence and overall stability of the complexes. A similar trend was also observed in the replica simulations with some minor variations (Supplementary Figure S2A). We also analyzed the residual flexibility of both PD-L1 and its complexes with Lumacaftor and Vedaprofen by calculating the root-mean-square fluctuation (RMSF) of all residues (Fig. 3B). The levels of PD-L1-ligand complexes fluctuated in the same way as free PD-L1 from the N- to C-termini. The RMSF analysis indicated stable residual fluctuations after ligand binding, with Lumacaftor showing slightly increased vibrations in some loop regions. The RMSF patterns of the docked complexes were found to have a similar trend as for the free protein, but with some increased fluctuations. Overall, PD-L1-Vedaprofen exhibited higher levels compared to PD-L1-Lumacaftor and free PD-L1. A similar trend was also observed in the replica simulations with some increased fluctuation in PD-L1 around 100–110 amino acid residues (Supplementary Figure S2B).

Fig. 3 [Images not available. See PDF.]

Backbone stability and flexibility of PD-L1 and complexes. (A) RMSD (Root Mean Square Deviation) profiles over 500 ns simulations for PD-L1 (black), PD-L1–Lumacaftor (orange), and PD-L1–Vedaprofen (aqua). (B) RMSF (Root Mean Square Fluctuation) per residue indicating flexibility variations across protein regions. Lower panels show the corresponding probability distribution functions (PDFs) indicating stability convergence. Data analyzed from 500 ns MD trajectories using Grace (GRaphing, Advanced Computation and Exploration), version 5.1.25 http://plasma-gate.weizmann.ac.il/Grace/.

Table 3. Summary of average molecular dynamics parameters over 500 ns. Reported metrics include RMSD, RMSF, radius of gyration (Rg), solvent accessible surface area (SASA), and average number of intramolecular hydrogen bonds. Simulations were conducted at 300 K using GROMOS 54A7 force field with the GROMACS 2020 suite.

System	RMSD (nm)	RMSF (nm)	Rg (nm)	SASA (nm2)	#H-bonds
PD-L1	0.14	0.08	1.42	70.27	77
PD-L1-Lumacaftor	0.16	0.09	1.41	69.40	76
PD-L1-Vedaprofen	0.18	0.10	1.40	70.00	74

Further, we investigated the conformational stability of PD-L1 with Lumacaftor and Vedaprofen by determining the radius of gyration (Rg) in both cases. The obtained average Rg values for PD-L1 in its free form and complexed with Lumacaftor and Vedaprofen were in the range of 1.35 nm to 2.45 nm. The time-dependent average Rg values for PD-L1, PD-L1-Lumacaftor, and PD-L1-Vedaprofen complexes were found to be 1.42 nm, 1.41 nm, and 1.40 nm, respectively. The Rg plot analysis suggested that there was a small reduction in the Rg of PD-L1 when bound to Lumacaftor and Vedaprofen. Small deviations of Rg that could be detected during the MD trajectories probably indicate the compact packing rearrangements of PD-L1. The Rg values became constant later, pointing towards the stability during the whole simulation period and the stability of the complexes (Fig. 4A). The replica simulation also showed a similar trend in the Rg values with some minor change (Supplementary Figure S2C). Solvent accessible surface area (SASA) of PD-L1, PD-L1-Lumacaftor, and PD-L1-Vedaprofen complexes was also evaluated for conformational analysis. The average SASA was calculated to be 70.27 nm²69.40 nm²and 70.00 nm²respectively. The equilibration pattern of all the systems in the SASA plot was similar, with higher compactness of the protein-ligand complexes than the free protein (Fig. 4B). A small reduction in average SASA values for PD-L1-Lumacaftor and PD-L1-Vedaprofen indicates that PD-L1 becomes more compact upon Lumacaftor and Vedaprofen binding and might have stronger contacts. SASA was distributed with a similar trend in the replica simulations with some minor variations (Supplementary Figure S2D).

Fig. 4 [Images not available. See PDF.]

Compactness and solvent exposure of PD-L1 complexes. (A) Radius of gyration (Rg) analysis showing conformational compactness of PD-L1 and its ligand complexes. (B) Solvent Accessible Surface Area (SASA) plots showing changes in hydration shell and surface exposure due to ligand binding. Lower panels display corresponding PDFs. Simulations span 500 ns.

Dynamics of hydrogen bonds

Hydrogen bonds play a significant role in maintaining the structural integrity of proteins and in mediating protein-ligand binding³⁰. The dynamics of hydrogen bonds in MD simulations are crucial for understanding protein stability, protein-ligand interactions, and the overall behavior of biomolecular systems³¹. To have a clear view of the stability of PD-L1 with Lumacaftor and Vedaprofen, we studied the temporal dynamics of hydrogen bonds formed within 0.35 nm during simulations (Fig. 5). At first, there were 77 intramolecular hydrogen bonds in PD-L1, which changed slightly to 76 and 74 after corresponding interaction with Lumacaftor and Vedaprofen (Fig. 5A). This small decrease in the hydrogen bonding is probably the result of the intramolecular space occupancy in PD-L1 by the ligands. To examine the hydrogen bond dynamics in more detail, we plotted probability distribution functions (PDFs) of the probability of hydrogen bonds for both cases (Fig. 5B). The results showed a similar distribution of hydrogen bonding before and after ligand binding. The intramolecular hydrogen bond distribution showed a similar trend in the replica simulations with some minor variations (Supplementary Figure S3A-B).

Fig. 5 [Images not available. See PDF.]

Intramolecular hydrogen bond dynamics within PD-L1. (A) Number of internal H-bonds over time for PD-L1 (black), PD-L1–Lumacaftor (orange), and PD-L1–Vedaprofen (aqua). (B) Corresponding PDF plots showing frequency distribution of hydrogen bond numbers across the simulation. Analysis reveals minimal loss in internal H-bonds post-ligand binding, indicating structural preservation.

We also analyzed the occurrence of the intermolecular hydrogen bonds formed between the protein-ligand complexes (Fig. 6). In the current study, we employed the gmx hbond utility in the GROMACS package to study the time-dependent formation of these bonds. The study showed that both PD-L1-Lumacaftor and PD-L1-Vedaprofen complexes were stable, having an average of 1–2 bonds (Fig. 6A-B, upper panels). The study showed that Lumacaftor forms 1–9 bonds to PD-L1 with higher fluctuations and 1–2 bonds with quite stability throughout the simulation (Fig. 6A). At the same time, Vedaprofen forms 1–8 bonds to PD-L1 with higher fluctuations and 1–2 bonds with quite stability throughout the simulation (Fig. 6A). This was also evident in the PDF values, which revealed that one hydrogen bond was present with the highest frequency in these complexes, as explained in Figs. 6A-B, lower panels. The intramolecular hydrogen bond distribution showed a similar trend in the replica simulations with some minor variations (Supplementary FigureS3C-D). These results point out the importance of hydrogen bonds in stabilizing the complexes of PD-L1-Lumacaftor and PD-L1-Vedaprofen and their structural organization. The outcomes indicate that the docking positions of Lumacaftor and Vedaprofen did not change during the simulation, hence providing more supporting data on complex stability. In summary, these findings provide useful information on the molecular recognition of protein-ligand binding and improve the knowledge of stable complex formation.

Fig. 6 [Images not available. See PDF.]

Intermolecular hydrogen bonding between PD-L1 and ligands. (A) Time-dependent intermolecular H-bond count for PD-L1–Lumacaftor, with fluctuations between 1–9 bonds. (B) PD-L1–Vedaprofen complex forming 1–8 H-bonds during the simulation. PDF panels (below) highlight that 1–2 persistent hydrogen bonds dominate ligand stabilization during the 500 ns trajectories.

Principal component analysis

PCA is an important technique for the identification of primary motion patterns of proteins and the investigation of dynamics using configurational space and protein motion based on essential degrees of freedom³². In this study, PCA was used to analyze the conformational flexibility of free PD-L1, the PD-L1-Lumacaftor, and the PD-L1-Vedaprofen. This made it possible for us to study their collective movements in terms of what might be referred to as the basic dynamics. The conformational complexity of these complexes was well captured within the essential subspace, as illustrated in Fig. 7; the conformational space of PD-L1 correlated with eigenvectors (EV-1 and EV-2) Cα. The data showed that PD-L1-Lumacaftor and PD-L1-Vedaprofen explore the same conformational subspace as free PD-L1 (Fig. 7A). Both EVs showed a similar pattern of conformational projections without any significant drift in the protein motional modes. This again indicates that PD-L1 complexes with Lumacaftor and Vedaprofen were quite stable during the simulations.

Fig. 7 [Images not available. See PDF.]

Principal component analysis (PCA) of essential dynamics. 2D projection of first two eigenvectors (EV1 vs. EV2) depicting conformational space sampled by PD-L1 (black), PD-L1–Lumacaftor (orange), and PD-L1–Vedaprofen (aqua). Similar clustering patterns support convergence and structural stability of the ligand-bound complexes relative to the apo form.

Free energy landscapes

The FEL analysis delivers detailed information about the stability of protein and the protein-ligand binding systems, the prediction of the binding transition states, and metastable conformations that are important for drug design²⁵. To assess the folding dynamics of free PD-L1, PD-L1-Lumacaftor, and PD-L1-Vedaprofen, FELs were built based on the first two principal components (PCs) (Fig. 8). The darker the blue, the lower the energy level, which indicates the native structure is closer to global minima. The FEL analysis yielded distinct findings: For PD-L1, there is only one global minimum, which is surrounded by two local minima (Fig. 8A). On the other hand, when Lumacaftor binds with PD-L1, the new metastates appear to be much larger but with well-defined global minima, thus suggesting that the conformational changes may be rather different with 4–5 discrete local basins (Fig. 8B). At the time, PD-L1 with Vedaprofen showed a well-defined global minimum, indicating 1–2 basins with stable conformational states (Fig. 8C). Taken together, MD simulations along with the PCA and FELs analysis of PD-L1 with Lumacaftor and Vedaprofen showed their stable interaction and offer indications of the Lumacaftor and Vedaprofen applicability as repurposed molecules as PD-L1 inhibitors.

Fig. 8 [Images not available. See PDF.]

Free energy landscapes (FELs) derived from MD simulations. (A) FEL of unbound PD-L1 showing a single global minimum and two local minima. (B) PD-L1–Lumacaftor complex exhibits broader conformational sampling with multiple local energy basins. (C) PD-L1–Vedaprofen complex shows 1–2 discrete metastable states. Energy minima are represented in darker blue. FELs were calculated based on projections onto the first two principal components.

Taken together, our integrated computational study identifies Lumacaftor and Vedaprofen as two FDA-approved drugs that bind stably to the PD-L1 pocket. These results suggest they could potentially inhibit PD-L1 and restore immune detection of cancer cells. However, this study has important limitations. It is entirely based on in silico modeling; experimental validation is needed to confirm that these compounds inhibit PD-L1 function. Moreover, Lumacaftor and Vedaprofen have known primary targets outside of PD-L1. Lumacaftor is designed to stabilize the F508del mutant of the cystic fibrosis transmembrane conductance regulator (CFTR) channel³³and Vedaprofen is an NSAID indicated for pain and inflammation³⁴. These differing targets mean off-target effects are possible if these drugs are used for PD-L1. Finally, MD simulations rely on force-field approximations and finite sampling, which cannot capture all aspects of biological systems. Future work should address these limitations with experimental tests and broader pharmacological profiling of Lumacaftor and Vedaprofen for human use. Such studies will be needed to confirm their specificity and effectiveness as repurposed PD-L1 inhibitors.

Conclusions

Among the various immune checkpoint molecules, the protein PD-L1 stands out as a potential drug target. This study explored the binding and inhibitory potential of repurposed drug molecules to the PD-L1 binding pocket. An integrated bioinformatics study involving molecular docking, drug profiling, MD simulations, and essential dynamics was performed to reveal the binding kinetics, thermodynamics, and dynamics of the molecules within the PD-L1 binding pocket. The study pinpointed Lumacaftor and Vedaprofen as two FDA-approved drugs with strong binding affinities and specific interactions with residues in the PD-L1 protein. Extended MD simulations confirmed the structural stability and favorable interaction dynamics of PD-L1 with Lumacaftor and Vedaprofen. Based on structural analyses such as RMSD, RMSF, Rg, and SASA, it was found that both complexes were stable at the end of the simulation period, and the interactions were stable as well. Moreover, the hydrogen bonding studies suggested that there was regularity in the interactions that are often required for the stabilization of complexes formed between PD-L1 and its ligands. Furthermore, PCA and FEL analysis offered an additional understanding of the conformational behavior and energy profiles of PD-L1 in complex with Lumacaftor and Vedaprofen, which backed up the stability and preferred binding conformations. Nevertheless, this in silico study does not account for off-target effects, in vivo pharmacokinetics, or therapeutic efficacy, which should be addressed through experimental and clinical validation. Altogether, the present results indicate that Lumacaftor and Vedaprofen can be considered potential anti-PD-L1 repurposed molecules for cancer treatment after further experimental and clinical research.

Acknowledgements

Mohd Shahnawaz Khan extend his appreciation to Ongoing Research Funding Program (ORF-2025-352), King Saud University, Riyadh, Saudi Arabia for funding this research. The authors are also grateful to Ajman University, UAE for supporting the publication.

Author contributions

M.S.K., A.S. and D.K.Y. wrote the main manuscript text and A.S., M.S.K, K.D. and M.S. prepared all the figures. All authors reviewed the manuscript and approved the revised version of the manuscript.

Data availability

All data supporting the findings of this study are available within the paper.

Declarations

Conflict of interest

The authors declare no competing interests.

Publisher’s note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

References

1. Khanna, I. Drug discovery in pharmaceutical industry: productivity challenges and trends. Drug Discovery Today; 2012; 17, pp. 1088-1102. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/22627006]

2. Parvathaneni, V; Kulkarni, NS; Muth, A; Gupta, V. Drug repurposing: a promising tool to accelerate the drug discovery process. Drug Discovery Today; 2019; 24, pp. 2076-2085.1:CAS:528:DC%2BC1MXht1ers7nL [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31238113][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC11920972]

3. Park, K. A review of computational drug repurposing. Translational Clin. Pharmacol.; 2019; 27, pp. 59-63.

4. Shamsi, A. et al. Targeting PDE4A for therapeutic potential: exploiting drug repurposing approach through virtual screening and molecular dynamics. Journal Biomol. Struct. Dynamics, 43(11), 5423–5435 (2024).

5. Shamsi, A; Khan, MS; Yadav, DK; Shahwan, M. Structure-based screening of FDA-approved drugs identifies potential histone deacetylase 3 repurposed inhibitor: molecular Docking and molecular dynamic simulation approaches. Front. Pharmacol.; 2024; 15, 1424175.1:CAS:528:DC%2BB2cXisFOhsLjL [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/39005934][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC11239971]

6. Duan, J et al. Use of immunotherapy with programmed cell death 1 vs programmed cell death ligand 1 inhibitors in patients with cancer: a systematic review and meta-analysis. JAMA Oncol.; 2020; 6, pp. 375-384. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31876895]

7. Jiang, X et al. Role of the tumor microenvironment in PD-L1/PD-1-mediated tumor immune escape. Mol. Cancer; 2019; 18, pp. 1-17.1:CAS:528:DC%2BC1MXhtl2qt7nE

8. Chen, X et al. Epigenetic strategies synergize with PD-L1/PD-1 targeted cancer immunotherapies to enhance antitumor responses. Acta Pharm. Sinica B; 2020; 10, pp. 723-733.1:CAS:528:DC%2BB3MXnslaisg%3D%3D

9. Akinleye, A; Rasool, Z. Immune checkpoint inhibitors of PD-L1 as cancer therapeutics. J. Hemato Onco; 2019; 12, 92.

10. Sun, JY et al. Resistance to PD-1/PD-L1 Blockade cancer immunotherapy: mechanisms, predictive factors, and future perspectives. Biomark. Res.; 2020; 8, pp. 1-10.1:CAS:528:DC%2BB3cXitF2ns7vK

11. Huey, R; Morris, GM; Forli, S. Using AutoDock 4 and AutoDock Vina with autodocktools: a tutorial. Scripps Res. Inst. Mol. Graphics Lab.; 2012; 10550, 1000.

12. Mohammad, T; Mathur, Y; Hassan, MI; InstaDock,. A single-click graphical user interface for molecular docking-based virtual high-throughput screening. Brief. Bioinform.; 2021; 22, bbaa279. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33105480]

13. Visualizer, D. Discovery Studio Visualizer. 2. Accelrys software inc (2005).

14. Van Der Spoel, D et al. GROMACS: fast, flexible, and free. J. Comput. Chem.; 2005; 26, pp. 1701-1718. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/16211538]

15. Rodriguez, I et al. Structural and biological characterization of pAC65, a macrocyclic peptide that blocks PD-L1 with equivalent potency to the FDA-approved antibodies. Mol. Cancer; 2023; 22, 150.1:CAS:528:DC%2BB3sXhvVOnsr7P [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/37679783][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC10483858]

16. Kaplan, W; Littlejohn, TG. Swiss-PDB viewer (deep view). Brief. Bioinform.; 2001; 2, pp. 195-197.1:STN:280:DC%2BD3Mvhs1Ggtw%3D%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/11465736]

17. Knox, C et al. DrugBank 6.0: the drugbank knowledgebase for 2024. Nucleic Acids Res.; 2024; 52, pp. D1265-D1275.1:CAS:528:DC%2BB2cXivVamt7vP [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/37953279]

18. Schuler, LD; Daura, X; Van Gunsteren, WF. An improved GROMOS96 force field for aliphatic hydrocarbons in the condensed phase. J. Comput. Chem.; 2001; 22, pp. 1205-1218.1:CAS:528:DC%2BD3MXltVelurs%3D

19. Tao, P; Lai, L. Protein ligand Docking based on empirical method for binding affinity Estimation. J. Comput. Aided Mol. Des.; 2001; 15, pp. 429-446.2001JCAMD.15.429T1:CAS:528:DC%2BD3MXktVSltr4%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/11394737]

20. DeLano, WL; Pymol,. An open-source molecular graphics tool. CCP4 Newsl. Protein Crystallogr.; 2002; 40, pp. 82-92.

21. Filimonov, D et al. Prediction of the biological activity spectra of organic compounds using the PASS online web resource. Chem. Heterocycl. Compd.; 2014; 50, pp. 444-457.1:CAS:528:DC%2BC2cXovVaitbs%3D

22. Mark, P; Nilsson, L. Structure and dynamics of the TIP3P, SPC, and SPC/E water models at 298 K. J. Phys. Chem. A; 2001; 105, pp. 9954-9960.1:CAS:528:DC%2BD3MXntlWrurs%3D

23. Jaidhan, B; Rao, PS; Apparao, A. Energy minimization and conformation analysis of molecules using steepest descent method. Int. J. Comput. Sci. Inf. Technol.; 2014; 5, pp. 3525-3528.

24. Sittel, F., Jain, A. & Stock, G. Principal component analysis of molecular dynamics: on the use of cartesian vs. internal coordinates. The J. Chem. Physics141, 014111 (2014).

25. Papaleo, E; Mereghetti, P; Fantucci, P; Grandori, R; De Gioia, L. Free-energy landscape, principal component analysis, and structural clustering to identify representative conformations from molecular dynamics simulations: the myoglobin case. J. Mol. Graph. Model.; 2009; 27, pp. 889-899.1:CAS:528:DC%2BD1MXms1ygs7Y%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/19264523]

26. Guzik, K et al. Small-molecule inhibitors of the programmed cell death-1/programmed death-ligand 1 (PD-1/PD-L1) interaction via transiently induced protein States and dimerization of PD-L1. J. Med. Chem.; 2017; 60, pp. 5857-5867.1:CAS:528:DC%2BC2sXpvFeisLc%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/28613862]

27. Shukla, R. & Tripathi, T. Molecular dynamics simulation of protein and protein–ligand complexes. Computer-aided Drug Design, 133–161 (2020) Springer, Singapore. https://doi.org/10.1007/978-981-15-6815-2_7.

28. Shukla, R; Munjal, NS; Singh, TR. Identification of novel small molecules against GSK3β for alzheimer’s disease using chemoinformatics approach. J. Mol. Graph. Model.; 2019; 91, pp. 91-104.1:CAS:528:DC%2BC1MXhtFOktL%2FF [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31202091]

29. Maiorov, VN; Crippen, GM. Significance of root-mean-square deviation in comparing three-dimensional structures of globular proteins. J. Mol. Biol.; 1994; 235, pp. 625-634.1:CAS:528:DyaK2cXhsFyisr0%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/8289285]

30. Hubbard, R. E. & Haider, M. K. Hydrogen bonds in proteins: role and strength. eLS (2010).

31. Bitencourt-Ferreira, G. & Veit-Acosta, M. & de Azevedo, W. F. Hydrogen bonds in protein-ligand complexes. Docking Screens Drug Discovery, 2053, 93–107 (2019).

32. Tharwat, A. Principal component analysis-a tutorial. Int. J. Appl. Pattern Recognit.; 2016; 3, pp. 197-240.

33. Connett, G. Lumacaftor-ivacaftor in the treatment of cystic fibrosis: design, development and place in therapy. Drug Des. Devel Ther, 13, 2405–2412 (2019).

34. Lees, M; Hoeijmakers, C; Rens,. A pharmacodynamic and Pharmacokinetic study with Vedaprofen in an equine model of acute nonimmune inflammation. J. Vet. Pharmacol. Ther.; 1999; 22, pp. 96-106.1:CAS:528:DyaK1MXjvFOmtL4%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/10372594]