- 1State Key Laboratory of Applied Organic Chemistry, Department of Chemistry, Lanzhou University, Lanzhou, China
- 2School of Basic Medical Science, Lanzhou University, Lanzhou, China
- 3Institute of Modern Physics of Chinese Academy of Sciences, Lanzhou, China
- 4School of Pharmacy, Lanzhou University, Lanzhou, China
- 5State Key Laboratory of Quality Research in Chinese Medicine, Macau Institute for Applied Research in Medicine and Health, Macau University of Science and Technology, Macau, China
Recently, small-molecule compounds have been reported to block the PD-1/PD-L1 interaction by inducing the dimerization of PD-L1. All these inhibitors had a common scaffold and interacted with the cavity formed by two PD-L1 monomers. This special interactive mode provided clues for the structure-based drug design, however, also showed limitations for the discovery of small-molecule inhibitors with new scaffolds. In this study, we revealed the structure-activity relationship of the current small-molecule inhibitors targeting dimerization of PD-L1 by predicting their binding and unbinding mechanism via conventional molecular dynamics and metadynamics simulation. During the binding process, the representative inhibitors (BMS-8 and BMS-1166) tended to have a more stable binding mode with one PD-L1 monomer than the other and the small-molecule inducing PD-L1 dimerization was further stabilized by the non-polar interaction of Ile54, Tyr56, Met115, Ala121, and Tyr123 on both monomers and the water bridges involved in ALys124. The unbinding process prediction showed that the PD-L1 dimerization kept stable upon the dissociation of ligands. It's indicated that the formation and stability of the small-molecule inducing PD-L1 dimerization was the key factor for the inhibitory activities of these ligands. The contact analysis, R-group based quantitative structure-activity relationship (QSAR) analysis and molecular docking further suggested that each attachment point on the core scaffold of ligands had a specific preference for pharmacophore elements when improving the inhibitory activities by structural modifications. Taken together, the results in this study could guide the structural optimization and the further discovery of novel small-molecule inhibitors targeting PD-L1.
Introduction
The blockage of the protein-protein interaction (PPI) between programmed cell death protein 1 (PD-1) and programmed cell death 1 ligand 1 (PD-L1) can reactivate the effector functions of T cell and eliminate tumor phenotypes with significant PD-L1 expression (Gatalica et al., 2014; Patel and Kurzrock, 2015; Sharma and Allison, 2015a,b). The crystal structures of PD-1/PD-L1 complex revealed the interface and hot-spot domains for both proteins (Zak et al., 2015; Pascolutti et al., 2016), which provided the structural basis for drug design. Ligands such as monoclonal antibodies (mAbs) (Lee et al., 2016, 2017; Liu K. et al., 2016; Tan et al., 2017; Zhang et al., 2017a,b), peptides (Chang et al., 2015; Magiera-Mularz et al., 2017), and small-molecule compounds (Abdel-Magid, 2015; Zak et al., 2016; Skalniak et al., 2017) had been discovered to interact with the PPI interface of PD-1 or PD-L1, showed obvious inhibitory activities against PD-1/PD-L1 signaling pathways. As the small-molecule inhibitors have better characteristic on aspects like production cost, drug-like property, immunogenic side effects, and half-life period (Liu K. et al., 2016) than peptides and monoclonal antibodies, the development of small-molecule inhibitor tended to be more promising. The crystal structures of small-molecule complex provided a good chance for the drug design of anti-cancer immunotherapy targeting on PD-1/PD-L1 immune checkpoint.
According the patents by Bristol-Myers Squibb (BMS) company, the compounds with (2-methyl-3-biphenylyl)methanol scaffold were privileged for inducing the dimerization of PD-L1 and interacted with the hydrophobic tunnel formed by two PD-L1 monomers (Zak et al., 2016; Guzik et al., 2017). Previously, George F. Gao's group resolved a dimeric interface of PD-L1 formed by B, C″, D, and E strands on each monomer, which was proved to be either a functional unit in immunological synapse formation or a revolution relics of B7 family (Tan et al., 2016). The crystal lattice analysis by Zak et al. also didn't suggest the spontaneous dimerization of PD-L1 (Zak et al., 2015), indicating that the interfacial interaction between two PD-L1 monomers was quite weak for dimerization process. As for the small molecule intervening PD-L1 dimerization, the interacting interface analysis showed that these ligands interacted with the G, F, C, C′ strands of PD-L1 in a competitive manner vs. PD-1 like mAbs or peptide inhibitors (Sharpe et al., 2011; Liu A. et al., 2016). Specially, the dimerized crystal structures tend to be a common pharmacodynamic characteristics for BMS small-molecule analogs despite of inhibitory activity difference from millimole to nanomole level (Abdel-Magid, 2015; Zak et al., 2016; Skalniak et al., 2017; Perry et al., 2019). Considering the potential relationship between the inhibitory activities of BMS small-molecule inhibitors and the stabilities of the dimerized complex systems, the dimerization process and the structure-activity relationship of small-molecule inhibitors need to be further elucidated. Besides, the broad, scattered and hydrophobic interface on PD-L1 makes it difficult for the discovery of novel small molecule ligands and also results in the strong hydrophobicity of BMS small-molecule inhibitors (Zarganes-Tzitzikas et al., 2016). Therefore, an understanding of the inhibitory mechanism of small-molecule ligands targeting PD-L1 such as key residues at the binding site, effect of the solvation and binding or unbinding process of small molecule inhibitors would help in the discovery of novel inhibitors and structural optimization of reported small-molecule inhibitors.
In this study, we aimed to reveal the detailed molecular mechanism of BMS small-molecule inhibitors from the formation and disassociation of PD-L1 dimers by multiple molecular modeling methods. Two representative compounds (BMS-8 and BMS-1166) with known inhibitory activities and complex crystal structures were selected to perform molecular dynamics simulations. During the formation process, both monomer and dimer systems of PD-L1 in complex with small-molecule ligands were applied to evaluate the stabilities of binding modes between ligands and PD-L1. The binding free energy calculation by MM-PBSA and MM-GBSA (Genheden and Ryde, 2015; Chen et al., 2018; Sun et al., 2018) were also used to analyze the energy contribution of the interfacial residues on PD-L1 dimers. During the disassociation process, metadynamics simulations (Bernardi et al., 2015) with specific collective variables (CVs) were performed to explore the key transition states along unbinding pathways. Based on the results of molecular modeling, an interplay mechanism of BMS small-molecule ligands with PD-L1 was proposed. Finally, R-group based QSAR analysis (Holliday et al., 2003; Hirons et al., 2005) and molecular docking were constructed on the reported BMS small-molecule inhibitors. The results of this study would provide a good guidance for the discovery of novel small-molecule inhibitors and structural modification of BMS small-molecule inhibitors targeting PD-L1.
Methods and Materials
The Conventional Molecular Dynamics Simulations
The complex crystal structures of BMS-8 and BMS-1166 were used to perform conventional molecular dynamics simulations. The Cartesian coordinates of the heavy atoms of PD-L1 (sequence 18–132) and small-molecule ligands were derived from the PDB database with accession number of 5J8O (Zak et al., 2016) and 5NIX (Skalniak et al., 2017). In order to eliminate the electrostatic effect of terminal residues, both monomers were capped with ACE and NME at two ends. The simulation details of the monomer systems and dimer systems were shown in Table 1. All the complex systems were firstly prepared through structural inspection and optimization in Schrödinger 2015 software suite (Schrödinger, LLC: New York, NY, 2015). Then, the complex proteins were solvated in a rectangular box of TIP3P waters and neutralized with Na+ ions. The periodic boundary conditions were setup with all the solvents at least 10 Å away from the complex. Then, the solvated systems were parameterized using the AMBER FF14SB force field (Case et al., 2014). The molecular dynamics simulations were performed in four steps. Firstly, energy minimization was performed to remove the local atomic collision in the systems. The energy minimization was conducted by both descent steepest method and conjugated gradient method with 5,000 steps. Then, the temperature of each system was gradually upgraded from 0 to 300 K in the NVT ensemble with all the solute atoms constrained with a force constant of 2.0 kcal mol−1·Å−2. After that, each system was equilibrated with the force constant decreasing from 2.0 to 0 kcal mol−1·Å−2 in a period of 1 ns. Finally, a production run of 150 ns was performed for each system in the NPT ensemble at 300 K and 1.0 atm condition. The snapshots for all the trajectories were saved every 2 ps.
The Binding Free Energy Calculation
For dimer systems of BMS-8 and BMS1166, two PD-L1 monomers were selected as the receptor, while small-molecule inhibitors were selected as the ligand. Both MM-PBSA and MM-GBSA methods were performed to calculate the binding free energy of BMS inhibitors according to the equation below:
Where < > represents the average value for all the snapshots used for MM-PBSA and MM-GBSA calculation. Different energy terms can be estimated as follows:
500 snapshots were extracted from the last 20 ns trajectories and used for MM-PBSA and MM-GBSA calculation. The parameter settings during MM-PBSA and MM-GBSA calculation were referred to the previous works published by our group (Xue et al., 2013). Then, the per-residue based decomposition was performed to identify the key residues in both dimer systems. Finally, the contribution of entropy change (–TΔS) was calculated by 100 snapshots from the last 20 ns trajectory.
The Calculation of Water Occupancies
The water molecules on the surface affected the conformational stability of proteins (Bellissent-Funel et al., 2016). By calculating the water occupancies on the surface of protein complex, water sites with a higher probability of finding a water molecule could be identified (Gauto et al., 2013). The water molecules at those sites were involved in the water bridges between protein and ligand and could enhance the stability of protein complex thermodynamically (Romero et al., 2016). To evaluate the effects of the water-mediated complex stability upon the binding of BMS inhibitors, the water occupancies and the water bridges were calculated over the last 20 ns trajectories for each dimer system using the “cpptraj” module of the AMBER14. All the trajectories were first imaged and fit to the first frame by the root mean square deviation (RMSD) of the heavy atoms of PD-L1 dimers. Then, the water occupancies were calculated using the “grid” command with a 0.5 Å *0.5 Å * 0.5 Å spacing over the whole box. And the water occupancies for both dimer systems were represented in the Chimera software (Pettersen et al., 2004).
Metadynamics Simulations
Metadynamics simulations have been widely used to predict the unbinding pathways and dissociation energy barrier of ligands for ligand-target systems (Cavalli et al., 2015). The sampling process of metadynamics simulations had an advantage of not requiring an initial estimate of the energy landscape to explore by periodically adding history-dependent biasing potential on selected collective variables (CVs) (Masetti et al., 2009; Barducci et al., 2011; Casasnovas et al., 2017; Sun et al., 2017). In this study, CV1 was selected as the distance between the mass center of the heavy atoms on ligand and the mass center of heavy atoms on key residues including Ile54, Tyr56, Met115, Ala121, Tyr123 in both chains; CV2 was selected as the angle between the Cα atom of Tyr56 and two carbon atoms that were the furthest away from each other on the core scaffold. The metadynamics simulations were performed for both dimer systems. The prepared topology files and coordinate files by AMBER ff14SB force field were further applied in the NAMD2.9 software (Kalé et al., 1999) implemented by PLUMED code (Bonomi et al., 2009). The initial structures were minimized for 5,000 steps with all the atoms on protein and ligand restrained with 5 kcal mol−1·Å−2 and all restraints released therewith. Then the temperature of systems were upgraded to 300 K in 30,000 steps. Afterward, all the systems were submitted to two short time NVT simulations (100,000 steps) to equilibrate the systems with restraining force constant of 5 kcal mol−1·Å−2 and all restraints released therewith. Finally, the equilibrated structures restarted from the NVT simulation were used for metadynamics simulations.
Metadynamics could reconstruct the free-energy surface as a function of specific collective variables (CVs). The general formalism of history-dependent Gaussian potential was shown as Equation (7). V represents the sum of the history-dependent Gaussian potential along the specific reactive coordinate (si) during time span (kτ). In this study, the deposition time (τ) was set as 1 ps to give enough dissociation time for ligands without adding biasing potential on the dissociation boundary. The Gaussian width (σ) of CV1 and CV2 were set to 0.8 Å and 0.02 rad, respectively. As for the well-tempered metadynamics, the height of the Gaussian potential (W) is affected by a parameter ΔT as Equation (8). The initial hill height (W0) of Gaussian potential was set to 0.6 kcal/mol·ps and the bias-factor (γ) was set to 10 with a temperature (T) of 300 K to control the decrease rate of the biasing potential as Equation (9).
R-Group QSAR Model and Molecular Docking of BMS Small-Molecule Inhibitors
The pharma R-group quantitative structure-activity relationship (RQSAR) models tended to be an effective approach for the SAR evaluation of the congeneric series of compounds (Adhikari et al., 2015). It was more suggestive than other approaches for the structural modification of small-molecule inhibitors by identifying the core scaffold and evaluating the effective pharma element at different attachment points (Kolarevic et al., 2018; Ts Mavrova et al., 2018). A total of 110 BMS small-molecule inhibitors with 2-methyl-3-(phenoxymethyl)-1,1'-biphenyl scaffold were collected from the patents of BMS company (Abdel-Magid, 2015; Table S1). All these small molecules had seven attachment points and diverse substitution groups, which were suitable to perform R-group QSAR analysis in the Canvas software of Schrödinger Suite (Duan et al., 2010; Sastry et al., 2010). The linear relationship between the substitutions and the activities (–log IC50) was analyzed and the importance of six key pharmacophore elements including hydrogen bond acceptor (A), hydrogen bond donor (D), hydrophobic group (H), negative ionic group (N), positive ionic group (P), and aromatic ring (A) were evaluated at each attachment point. During the process, the error and the importance were both set as 0.30. Eight representative small-molecule inhibitors (NO. of compound:4, 101, 102, 103, 104, 108, 109, 110) with substitutions on R1, R2, or R3 were selected to perform molecular docking to further study the binding modes. In order the compare the effect of R-groups, the core scaffold atoms with SMILES of “cOCc(c1C)cccc1c” were constrained with RMSD of 0.5 angstrom while other atoms were selected flexible. The standard precision (SP) docking score was used to evaluate the binding poses. The molecular docking was performed in Schrödinger 2015 software suite (Schrödinger, LLC: New York, NY, 2015).
Residue-Ligand Contact Analysis
In this study, we performed residue-ligand contact analysis to detect the surrounding residues around different substituent groups of BMS-8 and BMS-1166. It is assumed that the contacts exist between two groups as long as their distance was below a cutoff of 3.5 Å. The occupancy of each contact was calculated by the existence frequency in the 5,000 snapshots of the last 50 ns trajectories.
Results and Discussion
The Conformational Stabilities Between PD-L1 Monomer and BMS Small-Molecule Inhibitors
In order to explore the interactive process of BMS small-molecule inhibitors, we constructed two kinds of ligand-bound PD-L1 monomer systems as shown in Figure 1 and used molecular dynamics simulations to evaluate the stabilities of both binding modes by two replicas. As shown in Figure 2, the stabilities of both binding modes were evaluated by the RMSD of the heavy atoms of receptor, ligand and the core scaffold of ligands in two replicas. The core scaffold of BMS-8 tended to have a more stable contact with the conformation B than conformation A of PD-L1 according to the comparison of RMSD and the representative conformations of both binding modes. The detailed docking interactions diagram in Figures S1, S2 showed that the π-π stacking interaction between the biphenyl moiety and ATyr56 tended to be easily affected by the conformation of ATyr56 and unstable among three clusters, while the hydrophobic interaction between biphenyl moiety and residues on conformation B of PD-L1 tended to be stable among all three clusters. As for BMS-1166, both binding modes seemed to be quite stable, which probably accounting for the best inhibitory activities of BMS-1166 among the small-molecule inhibitors of BMS. The detailed docking interactions diagram in Figures S3, S4 showed that the biphenyl moiety had less conformational fluctuation and more stable hydrophobic interactions among the clusters of conformation A and B of PD-L1. The stabilities of the monomer complex of PD-L1 and ligand was affected by the hydrophobic interactions and turned out to be associated with the inhibitory activities of BMS small-molecule inhibitors.
Figure 1. The structural information of BMS-8 and BMS-1166. (A,B) The chemical formulas of BMS-8 and BMS-1166. The core scaffold is colored in red. (C,D) The conformational superposition of BMS-8 and BMS-1166 interacting with the monomer conformation A, B of PD-L1. (E,F) The surface of PD-L1 (A) and PD-L1 (B) interacting with BMS-8 and BMS-1166. The binding pockets formed by I54, V55, Y56, M115, I116, S117, A121, D122, and Y123 were colored in red.
Figure 2. The stability evaluation of the monomer systems. The RMSDs of the heavy atoms of PD-L1 monomer, ligands (BMS-8, BMS-1166) and the core scaffold of the ligand are shown in red, blue, and cyan lines, respectively. The representative conformations of three clusters for every monomer system was shown below the corresponding system. PD-L1 is shown in gray cartoon while the initial conformation of ligand is shown in orange sticks and the dynamics conformations of ligand is shown in green sticks.
The Interaction Stabilities Between PD-L1 Dimer and BMS Small-Molecule Inhibitors
The conformational stabilities of the dimer systems were evaluated by root mean square deviation (RMSD) and mean square root fluctuation (RMSF) as shown in Figure 3. The RMSDs of the complex, two PD-L1 monomers in complex systems showed that both dimer and monomers of PD-L1 had strong structural stabilities upon ligand binding. The conformational fluctuation of PD-L1 indicated that PD-L1 showed more flexibilities upon BMS-8 binding than BMS-1166. The comparison of the RMSDs of the core scaffold of two ligands showed that BMS-1166 had a more stable binding modes than BMS-8. The binding free energies were also calculated to evaluate the affinities of dimer systems. As shown in Table 2, the energy items of ΔΔGPB, ΔΔGGB, and ΔΔEexp by MM-PBSA and MM-GBSA methods could properly evaluate the difference of affinities of BMS-8 and BMS-1166, which showed the fact that the affinities between small-molecule inhibitors and PD-L1 dimer could reflect the inhibitory activities relatively. BMS-1166 had a stronger enthalpy contribution (ΔΔHPB, ΔΔHGB) and a worse entropy contribution (–TΔΔS) than BMS-8, which were consistent with the stability difference of BMS-8 and BMS-1166 dimer complex.
Figure 3. The RMSDs and RMSFs of the heavy atoms in dimer systems. (A) RMSD of complex, (B) RMSD of PD-L1 at conformation A, (C) RMSD of PD-L1 at conformation B, (D) RMSF of residues on PD-L1 at conformation A, B, (E) RMSD of the ligand, and (F) RMSD of the core scaffold.
Table 2. The binding free energies of BMS-8 and BMS-1166 evaluated by MM-PBSA, MM-GBSA, and metadynamics simulations.
The key residues on two PD-L1 monomers interacting with ligands were recognized by per-residue energy decomposition. The energy contribution for each residue were decomposed into the sidechain part and the backbone part, the non-polar part and the polar part as shown in Figures 4A–D. It could be seen that BMS-8 and BMS-1166 mainly formed non-polar interactions with the sidechain of the residues on PD-L1. With a cutoff value of −1.0 kcal/mol, the key residues in BMS-8 dimer system included ATyr56, AMet115, AAla121, ATyr123 and BIle54, BTyr56, BGln66, BMet115, BAla121 as shown in Figure 4E, while the key residues in BMS-1166 dimer system included AIle54, ATyr56, AMet115, AAla121, AAsp122, ATyr123, AArg125 and BIle54, BTyr56, BVal76, BMet115, BAla121, BAsp122 as shown in Figure 4F. Taken together, the interaction residues on conformation A and conformation B of PD-L1 were symmetrical both including Ile54, Tyr56, Met115, Ala121, and Tyr123. The hydrogen bond analysis in Table 3, Figures 5A,B showed that the protonated tertiary ammonium in BMS-8 formed a hydrogen bond with the side-chain oxygen of BGln66 with an occupancy of 57.21%, while the BMS-1166 dimer system also formed hydrogen bond between the ammonium group on BMS-1166 and the carboxyl group of AAsp122. The binding mode analysis of substitute groups on BMS-8 and BMS-1166 with the interfacial residues on PD-L1 indicated that the interaction with the peripheral residues including AIle54, AArg125, BVal76, and BAsp122 could significantly enhance the inhibitory activities of BMS small-molecule inhibitors.
Figure 4. The residue energy decomposition of the key residues for BMS-8 (A,C,E) and BMS-1166 (B,D,F) dimer systems. (A,B) The energy contribution of sidechain and backbone, (C,D) the energy contribution of polar and non-polar interaction, (E,F) the total energy contribution for each residue. The cutoff value of the energy contribution for key residues were set as −1 kcal/mol.
Figure 5. The binding modes and water occupancies in dimer systems. (A,B) The binding modes of BMS-8 and BMS-1166. The key residues were shown in green and cyan sticks and the ligand was shown in orange sticks. The hydrogen bonds were shown in red dash. (C,D) The water occupancies in BMS-8 and BMS-1166 dimer systems. The PD-L1 dimer is shown in cyan cartoon and the residues involved in water bridges are shown in green sticks. The water distributions are shown in red solid surface and the small molecule ligands are shown in orange sticks.
In order to analyze the effect of solvent on PD-L1 dimer complex, water occupancies and water bridges involved in receptor-ligand interaction were both calculated. As shown in Table 4, the residues or residue pairs involved in water bridges with an occupancy higher than 20% were extracted from both dimer systems. It can be seen in Figure 5 that three water bridges involved in AAsp122, ATyr123, ALys124, and BGln66 were stable in both dimer systems. Both ligands formed a strong water bridge with ALys124 with an occupancy higher than 90%, which indicated that ALys124 had a significant effect on the stabilities of the ligand conformations.
Table 4. The water bridge with occupancies higher than 20.00% in the BMS-8 and BMS-1166 dimer systems.
The Disassociation Process of BMS Small-Molecule Inhibitors
The free energy landscape of the unbinding processes of both BMS small-molecule inhibitors were constructed by CV1 and CV2. The distribution of minima in the landscapes showed that the most stable conformational state in the unbinding process was corresponding to the conformational states of the initial crystal structures as shown in Figures 6A,B. During the unbinding process, there were four different transition states for BMS-8 and three transition states for BMS-1166. In order to test the convergence of unbinding process, the free energies along both CVs were estimated. It can be seen that the free energy surface of CV1 (Figures 6C,D) and CV2 (Figures 6E,F) gradually came to a convergence along with the accumulation of time. As CV1 represented the distance between the ligand and the binding site of PD-L1 dimer and depicted the unbinding process better than CV2, the corresponding minimum points along CV1 were extracted from the unbinding trajectories. In BMS-8 complex systems, the minima along CV1 were 6.75 Å (−16.34 kcal/mol), 12.54 Å (−10.12 kcal/mol), and 14.97 Å (−7.31 kcal/mol). In BMS-1166 complex systems, the minima along CV1 were 5.67 Å (−28.79 kcal/mol), 10.95 Å (−8.50 kcal/mol), and 13.46 Å (−9.74 kcal/mol). The ultimate unbinding energy barriers of both small-molecule ligands estimated by CV1 and CV2 were shown in Table 1 and Figure 6, which were in good consistency with the inhibitory activities. Considering the difference between the binding free energies predicted by different methods, MM-PBSA and MM-GBSA calculated the binding free energies using implicit water models while metadynamics simulation considered the explicit water interaction between protein and ligand. Therefore, the results from metadynamics simulation tended to be more approximate to the experimental results.
Figure 6. The energy change during the unbinding process. (A,B) The free energy landscapes for the unbinding process of BMS-8 and BMS-1166. (C–F) The convergence of sampling process during the unbinding process of BMS-8 and BMS-1166. The history-dependent free energy surfaces along CV1 (C,D) and CV2 (E,F) are estimated by a segmented accumulation of simulation time. The unit for the free energy is kcal/mol.
From the free energy estimation of different conformational states, it can be seen that the conformational states of the crystal structures were much more stable than the other transition conformational states along the unbinding process. Therefore, the dissociation of small-molecule ligands of the initial conformational states tended to be the most important intermediate process for the unbinding of small-molecule ligands, which were corresponding to the minima of CV1 at 12.54 Å in BMS-8 dimer systems and the minima of CV1 at 10.95 Å in BMS-202 dimer systems. The corresponding transition states were extracted from the trajectories as shown in Figures 7A,B. The binding poses of BMS small-molecule ligands at the transition states were quite distinct from each other, which was probably owing to the difference of substituent groups. A common feature for both systems was that the ligands at transition states significantly lost the interaction with the chain A while the interaction with chain B were still compact and involved with a series of residues especially in BMS-1166 dimer systems. During the unbinding process, the core scaffold of ligands gradually divorced from the location of ATyr56 and got away from the pocket formed by PD-L1 dimer. In order to monitor the conformational change of the pocket formed by PD-L1 monomers, the distance between chain A and chain B were calculated by the distance between Ile54, Tyr56, Met115, Ala121, Tyr123 on each chain as shown in Figures 7C–F. The conformational fluctuation of PD-L1 monomers was reflected by the conformational change of the F-G loops on both PD-L1 monomers. It can be seen that the pockets in BMS-8 and BMS-1166 complex systems were quite stable with occasionally occurring conformational fluctuations. According to unbinding processes of BMS ligands, it can be seen that the dimer of PD-L1 had a large tendency to keep stable although accompanied with subtle conformational fluctuation of PD-L1 dimer.
Figure 7. The conformational change during the unbinding process. (A,B) The key transition conformational states of BMS-8 and BMS-1166 during the unbinding process with CV1 value of 12.54 Å and 10.95 Å, respectively. The chain A and chain B of PD-L1 dimers are shown in green and cyan surface, respectively. The key residues (green or cyan) and ligands (orange) are shown in sticks. (C,D) The distance between Ile54, Tyr56, Met115, Ala121, and Tyr123 on conformation A and conformation B of PD-L1 during the unbinding process of BMS-8 and BMS-1166. (E,F) The extracted conformational state (green cartoon) of each complex with the largest distance between two PD-L1 monomers was overlapped with the original crystal structures (magenta cartoon).
Taken together, the most possible deduction for the interaction mechanism of BMS small-molecule inhibitors with PD-L1 was depicted as shown in Figure 8. Firstly, all BMS small-molecule inhibitors with different activities tended to interacted with a monomer conformation B of PD-L1. As the PD-L1 dimer complex had strong conformational stability, the PD-L1 monomer complex further interacted with the other monomer of PD-L1 to form PD-L1 dimer complex. According to the results of metadynamics simulation, a complete dissociation for BMS inhibitors would probably be like that the small-molecule ligand was firstly unbound from the PD-L1 dimer and the rest receptor part was further depolymerized into monomer.
Figure 8. The complete binding and unbinding mechanism of BMS small-molecule inhibitors. The conformation A and conformation B of PD-L1 monomer, the small-molecule ligand are represented by green, light blue, orange objects, respectively.
The R-Group QSAR Model of BMS Small-Molecule Inhibitors
110 BMS small-molecule inhibitors with 2-methyl-3-(phenoxymethyl)-1,1′-biphenyl scaffolds were tested with diverse inhibitory activities with IC50 ranging from 9.492 μM to 1.4 nM. As shown in Figure 9A, there were 7 different of attachment points from R1 to R7 and the substituent groups of R6 and R7 had a relatively larger proportion than other attachment points. As shown in Figure 9B, the correlation coefficient between the predicted pIC50 and the experimental pIC50 was 0.7729. According to the evaluation of six key pharmacophore elements in Figure 9C, the substituent groups at R2, R4, R6, and R7 had obvious effect on the affinity of BMS small-molecule inhibitors. The substituent groups at R2, R4, R6, R7 of BMS-8 and BMS-116 as well as the interaction residues were recognized by the contact analysis as shown in Figures 9D,E. The contact analysis of BMS-1166 showed that the 1,4-benzodioxinyl group at R2 mainly was involved in the interaction with AIle54, ATyr56, BAsp122, BTyr123. The hydrogen bond acceptor and hydrophobic groups at R2 were favorable for BMS inhibitors such as the 2, 3-dihydro-1, 4-benzodioxinyl group on BMS-114, BMS-200, BMS-1001, and BMS-1166. The analysis of effect of solvent in dimer systems showed that the substituent groups at R4, R5, R6, and R7 were exposed to solvent environment. The hydrophobic groups at R4 were favorable for BMS inhibitors, which corresponded to the fact that the bromine atom on BMS-8 and the chlorine atom on BMS-1166 had a close contact with BIle54. The hydrophobic group at R6 was adverse while the negative ionic group was favorable. The substituent group at R6 of BMS-8 mainly interacted with BTyr56 and BGln66, however, that of BMS-1166 mainly interacted with AThr20 and AAsp122. Nevertheless, the substituent group at R6 of BMS-8 and BMS-1166 both formed hydrogen bonding with PD-L1. The positive ionic at R7 was adverse while the hydrogen bond acceptor and aromatic ring were favorable. The substituent group at R7 of BMS-1166 formed interaction with AAsp122, ATyr123, ALys124, AArg125, BTyr56, and BGln63. The comparison of the contact residues between BMS-8 and BMS-1166 indicated that the substituent group at R2 and R7 strongly strengthened the interactions with the conformation A of PD-L1, which was consistent with the stability of the monomer complex of conformation A and BMS-1166.
Figure 9. The R-group based QSAR for BMS small-molecule inhibitors. (A) The number of substituent groups at seven attachment points of BMS small molecule inhibitors. (B) The correlation validation between the predicted pIC50 and experimental pIC50. (C) R-group QSAR model for the BMS small molecule inhibitors. The case-insensitive alphabets A(a), D(d), H(h), N(n), P(p), and R(r), respectively represent the hydrogen bond acceptor, the hydrogen bond donor, the hydrophobic, the negative ionic, the positive ionic, aromatic ring. The significantly increase effect is colored in red while the significantly decrease effect is colored in blue. (D,E) The interacting residues of different substituent groups of BMS-8 and BMS-1166. The occupancies of the contact between each residues and substituent group were listed along with the black solid.
The further molecular docking study of eight representative small-molecule inhibitors showed that the docking scores had a good linear correlation with the experimental inhibitory activities (Figure S5). The further residue contribution comparison (Figure S6) and conformational analysis (Figure S7) of the residues within 5 angstroms showed that the residues interacting with R-group substituents had an obvious effect on the docking scores including BAsp122 (interacting with R1 to R3) and AAsp122, ALys124, BTyr56, BGln66 (interacting with R4 to R7). The binding mode analysis of novel series of [1,2,4]triazolo[4,3-a]pyridines designed by Qin et al. also revealed the retaining hydrophobic interaction with Tyr56, Met115, and Ala121 on both chain of PD-L1 and extra π-π stacking with the BTyr56 and π-anion interactions with AAsp122 (Qin et al., 2019). These interacting modes of [1,2,4]triazolo[4,3-a]pyridines inhibitors were consistent with the binding mode analysis of eight representative small-molecule inhibitors. It's suggested that the structure-activity relationship analysis of BMS small-molecule inhibitors was applicable for the further structure modifications.
Conclusions
In this study, we used multiple molecular modeling methods to study the detailed molecular mechanism of the interaction between BMS small-molecule inhibitors and PD-L1. A detailed mechanism of the interaction process between small-molecule inhibitors and PD-L1 was proposed and validated by molecular dynamics simulations.
The BMS small-molecule inhibitors tended to interact with one PD-L1 monomer first and further formed dimer with the other monomer for an advantage of stability. The results of binding free energy and water occupancy calculation revealed the key stability factors for ligand-induced PD-L1 dimers including the hydrophobic contribution of Ile54, Tyr56, Met115, Ala121, and Tyr123 on both monomers and the water bridges involved in ALys124. The unbinding pathway prediction also indicated that the tunnel formed by PD-L1 dimers tended to be stable upon the getting away of BMS-inhibitors. The R-group QSAR model suggested that the substituents at R2, R4, R6, and R7 had a significant effect on the inhibition activities of BMS inhibitors. The structural modification with these substituent positions tended to be an effective way to improve the inhibition activities of BMS inhibitors. Taken together, this study would provide a comprehensive view of the inhibition mechanism for BMS small-molecule inhibitors and guide the further development of more potential small-molecule inhibitors targeting PD-L1.
Data Availability Statement
The raw data supporting the conclusions of this manuscript will be made available by the authors, without undue reservation, to any qualified researcher.
Author Contributions
DS, HL, and XY designed the research. DS and XY were responsible for the writing and revising of the manuscript. XA, QB, SZ, and ZB were responsible the main data analysis in the manuscript.
Funding
This work was supported by the National Natural Science Foundation of China (Grant Nos. 21475054 and 2175060).
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fchem.2019.00764/full#supplementary-material
Figure S1. The interactions diagrams between the monomer conformation A of PD-L1 (the initial crystal structure and three respective dynamics structures) and BMS-8 in the monomer system of replica 2.
Figure S2. The interaction diagrams between the monomer conformation B of PD-L1 (the initial crystal structure and three respective dynamics structures) and BMS-8 in the monomer system of replica 2.
Figure S3. The interaction diagrams between the monomer conformation A of PD-L1 (the initial crystal structure and three respective dynamics structures) and BMS-1166 in the monomer system of replica 2.
Figure S4. The interaction diagrams between the monomer conformation B of PD-L1 (the initial crystal structure and three respective dynamics structures) and BMS-1166 in the monomer system of replica 2.
Figure S5. The linear correlation between experimental pIC50 and the absolute values of the docking scores.
Figure S6. The distance and residue contribution analysis of the binding poses of eight representative small-molecule inhibitors. (A) The respective and average distance between the small-molecule inhibitor and the residues on PD-L1 dimer. (B) The respective energy contribution of residues on PD-L1 dimer when interacting with the small-molecule inhibitor.
Figure S7. The binding pose analysis of eight representative small-molecule inhibitors. (A) The surrounding residues of the substituent groups at R1 to R3 for eight representative small-molecule inhibitors. (B–I) The surrounding residues of the substituent groups at R4 to R7 for small-molecule inhibitor with NO. of 4, 101, 102, 103, 104, 108, 119, 110, respectively.
Table S1. The detailed structural and activity information for 110 BMS small-molecule inhibitors.
References
Abdel-Magid, A. F. (2015). Inhibitors of the PD-1/PD-L1 pathway can mobilize the immune system: an innovative potential therapy for cancer and chronic infections. ACS Med. Chem. Lett. 6, 489–490. doi: 10.1021/acsmedchemlett.5b00148
Adhikari, N., Halder, A. K., Saha, A., Das Saha, K., and Jha, T. (2015). Structural findings of phenylindoles as cytotoxic antimitotic agents in human breast cancer cell lines through multiple validated QSAR studies. Toxicol. In Vitro 29, 1392–1404. doi: 10.1016/j.tiv.2015.05.017
Barducci, A., Bonomi, M., and Parrinello, M. (2011). Metadynamics. Wiley Interdiscip. Rev. Comput. Mol. 1, 826–843. doi: 10.1002/wcms.31
Bellissent-Funel, M. C., Hassanali, A., Havenith, M., Henchman, R., Pohl, P., Sterpone, F., et al. (2016). Water determines the structure and dynamics of proteins. Chem. Rev. 116, 7673–7697. doi: 10.1021/acs.chemrev.5b00664
Bernardi, R. C., Melo, M. C., and Schulten, K. (2015). Enhanced sampling techniques in molecular dynamics simulations of biological systems. Biochim. Biophys. Acta. Gen. Sub. 1850, 872–877. doi: 10.1016/j.bbagen.2014.10.019
Bonomi, M., Branduardi, D., Bussi, G., Camilloni, C., Provasi, D., Raiteri, P., et al. (2009). PLUMED: a portable plugin for free-energy calculations with molecular dynamics. Comput. Phys. Commun. 180, 1961–1972. doi: 10.1016/j.cpc.2009.05.011
Casasnovas, R., Limongelli, V., Tiwary, P., Carloni, P., and Parrinello, M. (2017). Unbinding kinetics of a p38 MAP kinase type II inhibitor from metadynamics simulations. J. Am. Chem. Soc. 139, 4780–4788. doi: 10.1021/jacs.6b12950
Case, D., Babin, V., Berryman, J., Betz, R., Cai, Q., Cerutti, D., et al. (2014). AMBER 14. San Francisco, CA: University of California.
Cavalli, A., Spitaleri, A., Saladino, G., and Gervasio, F. L. (2015). Investigating drug-target association and dissociation mechanisms using metadynamics-based algorithms. Acc. Chem. Res. 48, 277–285. doi: 10.1021/ar500356n
Chang, H. N., Liu, B. Y., Qi, Y. K., Zhou, Y., Chen, Y. P., Pan, K. M., et al. (2015). Blocking of the PD-1/PD-L1 interaction by a D-peptide antagonist for cancer immunotherapy. Angew. Chem. Int. Edit. 54, 11760–11764. doi: 10.1002/anie.201506225
Chen, F., Sun, H., Wang, J., Zhu, F., Liu, H., Wang, Z., et al. (2018). Assessing the performance of MM/PBSA and MM/GBSA methods. 8. Predicting binding free energies and poses of protein-RNA complexes. RNA 24, 1183–1194. doi: 10.1261/rna.065896.118
Duan, J., Dixon, S. L., Lowrie, J. F., and Sherman, W. (2010). Analysis and comparison of 2D fingerprints: insights into database screening performance using eight fingerprint methods. J. Mol. Graph. Model. 29, 157–170. doi: 10.1016/j.jmgm.2010.05.008
Gatalica, Z., Snyder, C., Maney, T., Ghazalpour, A., Holterman, D. A., Xiao, N., et al. (2014). Programmed cell death 1 (PD-1) and its ligand (PD-L1) in common cancers and their correlation with molecular cancer type. Cancer Epidem. Biomar. 23, 2965–2970. doi: 10.1158/1055-9965.EPI-14-0654
Gauto, D. F., Petruk, A. A., Modenutti, C. P., Blanco, J. I., Di Lella, S., and Marti, M. A. (2013). Solvent structure improves docking prediction in lectin-carbohydrate complexes. Glycobiology 23, 241–258. doi: 10.1093/glycob/cws147
Genheden, S., and Ryde, U. (2015). The MM/PBSA and MM/GBSA methods to estimate ligand-binding affinities. Expert Opin. Drug Discov. 10, 449–461. doi: 10.1517/17460441.2015.1032936
Guzik, K., Zak, K. M., Grudnik, P., Magiera, K., Musielak, B., Törner, R., et al. (2017). 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. 60, 5857–5867. doi: 10.1021/acs.jmedchem.7b00293
Hirons, L., Holliday, J. D., Jelfs, S. P., Willett, P., and Gedeck, P. (2005). Use of the R-group descriptor for alignment-free QSAR. QSAR Comb. Sci. 24, 611–619. doi: 10.1002/qsar.200510102
Holliday, J., Jelfs, S., Willett, P., and Gedeck, P. (2003). Calculation of intersubstituent similarity using R-group descriptors. J. Chem. Inf. Comput. Sci. 43, 406–411. doi: 10.1021/ci025589v
Kalé, L., Skeel, R., Bhandarkar, M., Brunner, R., Gursoy, A., Krawetz, N., et al. (1999). NAMD2: greater scalability for parallel molecular dynamics. J. Comput. Phys. 151, 283–312. doi: 10.1006/jcph.1999.6201
Kolarevic, A., Ilic, B. S., Anastassova, N., Mavrova, A. T., Yancheva, D., Kocic, G., et al. (2018). Benzimidazoles as novel deoxyribonuclease I inhibitors. J. Cell. Biochem. 119, 8937–8948. doi: 10.1002/jcb.27147
Lee, H. T., Lee, J. Y., Lim, H., Lee, S. H., Moon, Y. J., Pyo, H. J., et al. (2017). Molecular mechanism of PD-1/PD-L1 blockade via anti-PD-L1 antibodies atezolizumab and durvalumab. Sci. Rep. 7:5532. doi: 10.1038/s41598-017-06002-8
Lee, J. Y., Lee, H. T., Shin, W., Chae, J., Choi, J., Kim, S. H., et al. (2016). Structural basis of checkpoint blockade by monoclonal antibodies in cancer immunotherapy. Nat. Commun. 7:13354. doi: 10.1038/ncomms13354
Liu, A., Dong, L., Wei, X. L., Yang, X. H., Xiao, J. H., and Liu, Z. Q. (2016). Development of amino- and dimethylcarbamate-substituted resorcinol as programmed cell death-1 (PD-1) inhibitor. Eur. J. Pharm. Sci. 88, 50–58. doi: 10.1016/j.ejps.2016.03.023
Liu, K., Tan, S., Chai, Y., Chen, D., Song, H., Zhang, C. W., et al. (2016). Structural basis of anti-PD-L1 monoclonal antibody avelumab for tumor therapy. Cell Res. 27, 151–153. doi: 10.1038/cr.2016.102
Magiera-Mularz, K., Skalniak, L., Zak, K. M., Musielak, B., Rudzinska-Szostak, E., Berlicki, Ł., et al. (2017). Bioactive macrocyclic inhibitors of the PD-1/PD-L1 immune checkpoint. Angew. Chem. Int. Ed. 129, 13920–13923. doi: 10.1002/ange.201707707
Masetti, M., Cavalli, A., Recanatini, M., and Gervasio, F. L. (2009). Exploring complex protein–ligand recognition mechanisms with coarse metadynamics. J. Phys. Chem. B 113, 4807–4816. doi: 10.1021/jp803936q
Pascolutti, R., Sun, X., Kao, J., Maute, R. L., Ring, A. M., Bowman, G. R., et al. (2016). Structure and dynamics of PD-L1 and an ultra-high-affinity PD-1 receptor mutant. Structure 24, 1719–1728. doi: 10.1016/j.str.2016.06.026
Patel, S. P., and Kurzrock, R. (2015). PD-L1 expression as a predictive biomarker in cancer immunotherapy. Mol. Cancer Ther. 14, 847–856. doi: 10.1158/1535-7163.MCT-14-0983
Perry, E., Mills, J. J., Zhao, B., Wang, F., Sun, Q., Christov, P. P., et al. (2019). Fragment-based screening of programmed death ligand 1 (PD-L1). Bioorg. Med. Chem. Lett. 29, 786–790. doi: 10.1016/j.bmcl.2019.01.028
Pettersen, E. F., Goddard, T. D., Huang, C. C., Couch, G. S., Greenblatt, D. M., Meng, E. C., et al. (2004). UCSF chimera—a visualization system for exploratory research and analysis. J. Comput. Chem. 25, 1605–1612. doi: 10.1002/jcc.20084
Qin, M., Cao, Q., Zheng, S., Tian, Y., Zhang, H., Xie, J., et al. (2019). Discovery of [1,2,4]Triazolo[4,3-a]pyridines as potent inhibitors targeting the programmed cell death-1/programmed cell death-ligand 1 interaction. J. Med. Chem. 62, 4703–4715. doi: 10.1021/acs.jmedchem.9b00312
Romero, J. M., Trujillo, M., Estrin, D. A., Rabinovich, G. A., and Di Lella, S. (2016). Impact of human galectin-1 binding to saccharide ligands on dimer dissociation kinetics and structure. Glycobiology 26, 1317–1327. doi: 10.1093/glycob/cww052
Sastry, M., Lowrie, J. F., Dixon, S. L., and Sherman, W. (2010). Large-scale systematic analysis of 2D fingerprint methods and parameters to improve virtual screening enrichments. J. Chem. Inf. Model. 50, 771–784. doi: 10.1021/ci100062n
Sharma, P., and Allison, J. P. (2015a). The future of immune checkpoint therapy. Science 348, 56–61. doi: 10.1126/science.aaa8172
Sharma, P., and Allison, J. P. (2015b). Immune checkpoint targeting in cancer therapy: toward combination strategies with curative potential. Cell 161, 205–214. doi: 10.1016/j.cell.2015.03.030
Sharpe, A. H., Butte, M. J., and Oyama, S. (2011). Modulators of Immunoinhibitory Receptor pd-1, and Methods of Use Thereof. U.S. Patent Application.
Skalniak, L., Zak, K. M., Guzik, K., Magiera, K., Musielak, B., Pachota, M., et al. (2017). Small-molecule inhibitors of PD-1/PD-L1 immune checkpoint alleviate the PD-L1-induced exhaustion of T-cells. Oncotarget 8, 72167–72181. doi: 10.18632/oncotarget.20050
Sun, H., Duan, L., Chen, F., Liu, H., Wang, Z., Pan, P., et al. (2018). Assessing the performance of MM/PBSA and MM/GBSA methods. 7. Entropy effects on the performance of end-point binding free energy calculation approaches. Phys. Chem. Chem. Phys. 20, 14450-14460. doi: 10.1039/c7cp07623a
Sun, H., Li, Y., Shen, M., Li, D., Kang, Y., and Hou, T. (2017). Characterizing drug-target residence time with metadynamics: how to achieve dissociation rate efficiently without losing accuracy against time-consuming approaches. J. Chem. Inf. Model. 57, 1895–1906. doi: 10.1021/acs.jcim.7b00075
Tan, S., Chen, D., Liu, K., He, M., Song, H., Shi, Y., et al. (2016). Crystal clear: visualizing the intervention mechanism of the PD-1/PD-L1 interaction by two cancer therapeutic monoclonal antibodies. Protein Cell 7, 866–877. doi: 10.1007/s13238-016-0337-7
Tan, S., Liu, K., Chai, Y., Zhang, C. W., Gao, S., Gao, G. F., et al. (2017). Distinct PD-L1 binding characteristics of therapeutic monoclonal antibody durvalumab. Protein Cell 9, 135–139. doi: 10.1007/s13238-017-0412-8
Ts Mavrova, A., Dimov, S., Yancheva, D., Kolarevic, A., Ilic, B. S., Kocic, G., et al. (2018). Synthesis and DNase I inhibitory properties of some 5,6,7,8-tetrahydrobenzo[4,5]thieno[2,3-d]pyrimidines. Bioorg. Chem. 80, 693–705. doi: 10.1016/j.bioorg.2018.07.009
Xue, W., Jin, X., Ning, L., Wang, M., Liu, H., and Yao, X. (2013). Exploring the molecular mechanism of cross-resistance to HIV-1 integrase strand transfer inhibitors by molecular dynamics simulation and residue interaction network analysis. J. Chem. Inf. Model. 53, 210–222. doi: 10.1021/ci300541c
Zak, K. M., Grudnik, P., Guzik, K., Zieba, B. J., Musielak, B., Dömling, A., et al. (2016). Structural basis for small molecule targeting of the programmed death ligand 1 (PD-L1). Oncotarget 7, 30323–30335. doi: 10.18632/oncotarget.8730
Zak, K. M., Kitel, R., Przetocka, S., Golik, P., Guzik, K., Musielak, B., et al. (2015). Structure of the complex of human programmed death 1, PD-1, and its ligand PD-L1. Structure 23, 2341–2348. doi: 10.1016/j.str.2015.09.010
Zarganes-Tzitzikas, T., Konstantinidou, M., Gao, Y., Krzemien, D., Zak, K., Dubin, G., et al. (2016). Inhibitors of programmed cell death 1 (PD-1): a patent review (2010-2015). Expert Opin. Ther. Pat. 26, 973–977. doi: 10.1080/13543776.2016.1206527
Zhang, F., Qi, X., Wang, X., Wei, D., Wu, J., Feng, L., et al. (2017a). Structural basis of the therapeutic anti-PD-L1 antibody atezolizumab. Oncotarget 8, 90215–90224. doi: 10.18632/oncotarget.21652
Keywords: PD-L1, small-molecule inhibitors, molecular dynamics simulation, metadynamics simulation, R-group QSAR
Citation: Shi D, An X, Bai Q, Bing Z, Zhou S, Liu H and Yao X (2019) Computational Insight Into the Small Molecule Intervening PD-L1 Dimerization and the Potential Structure-Activity Relationship. Front. Chem. 7:764. doi: 10.3389/fchem.2019.00764
Received: 02 May 2019; Accepted: 24 October 2019;
Published: 12 November 2019.
Edited by:
Simone Brogi, University of Pisa, ItalyReviewed by:
Eduardo De Faria Franca, Federal University of Uberlandia, BrazilChanin Nantasenamat, Mahidol University, Thailand
Copyright © 2019 Shi, An, Bai, Bing, Zhou, Liu and Yao. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Huanxiang Liu, hxliu@lzu.edu.cn; Xiaojun Yao, xjyao@lzu.edu.cn