Computational Insight Into the Small Molecule Intervening PD-L1 Dimerization and the Potential Structure-Activity Relationship

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(Lee et al., , 2017Liu 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  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-3biphenylyl)methanol scaffold were privileged for inducing the dimerization of PD-L1 and interacted with the hydrophobic tunnel formed by two PD-L1 monomers 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 . 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;. 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.

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  and 5NIX . 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 Frontiers in Chemistry | www.frontiersin.org 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 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.
Frontiers in Chemistry | www.frontiersin.org 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 (s i ) 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 (W 0 ) 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 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 IC 50 ) 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 smallmolecule 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.

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 A Tyr56 tended to be easily affected by the conformation of A Tyr56 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  a The occupancy of hydrogen bonds were analyzed through the last 20 ns trajectories and only hydrogen bonds with an occupancy more than 0.20 were extracted and shown. b,c The hydrogen bonds were determined by an acceptor-donor atom distance of <3.5 Å and acceptor H-donor angle of >120 • .   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.

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 G PB , G GB, and E exp 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 ( H PB , H GB ) 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.
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 nonpolar part and the polar part as shown in Figures 4A-D. It could be seen that BMS-8 and BMS-1166 mainly formed nonpolar 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 A Tyr56, A Met115, A Ala121, A Tyr123 and B Ile54, B Tyr56, B Gln66, B Met115, B Ala121 as shown in Figure 4E, while the key residues in BMS-1166 dimer system included A Ile54, A Tyr56, A Met115, A Ala121, A Asp122, A Tyr123, A Arg125 and B Ile54, B Tyr56, B Val76, B Met115, B Ala121, B Asp122 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 B Gln66 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 A Asp122. 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 A Ile54, A Arg125, B Val76, and B Asp122 could significantly enhance the inhibitory activities of BMS smallmolecule inhibitors.
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 A Asp122, A Tyr123, A Lys124, and B Gln66 were stable in both dimer systems.
Both ligands formed a strong water bridge with A Lys124 with an occupancy higher than 90%, which indicated that A Lys124 had a significant effect on the stabilities of the ligand conformations.

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.
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 A Tyr56 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.
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.

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 IC 50 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 pIC 50 and the experimental pIC 50 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 A Ile54, A Tyr56, B Asp122, B Tyr123.
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 B Ile54. 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 B Tyr56 and B Gln66, however, that of BMS-1166 mainly interacted with A Thr20 and A Asp122. 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 A Asp122, A Tyr123, A Lys124, A Arg125, B Tyr56, and B Gln63. 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.
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 B Asp122 (interacting with R1 to R3) and A Asp122, A Lys124, B Tyr56, B Gln66 (interacting with R4 to R7). The binding mode analysis of novel series of [1,2,4]triazolo[4,3a]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 B Tyr56 and π-anion interactions with A Asp122 (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 A Lys124. 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 smallmolecule 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.

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.   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.