Computational discovery of small drug-like compounds as potential inhibitors of PD-1/PD-L1 interactions

Abstract The programmed cell death ligand protein 1 (PD-L1) is a strong immunosuppressive molecule that inactivates tumor-specific T cells by binding to the programmed cell death- 1 protein (PD-1). Cancer immunotherapy based on the monoclonal antibodies targeting the PD-1/PD-L1 pathway has demonstrated therapeutic responses without precedent over a wide range of cancers. However, the antibody-based immunotherapies have several limitations such as high production cost or the induction of severe immune-related adverse effects. Small-molecule inhibitors of the PD-1/PD-L1 pathway are a promising alternative or complementary therapeutic to antibodies. Currently, the field of developing anti-PD-1/PD-L1 small-molecule inhibitors is intensively explored. In the present study a pharmacophore model was generated based on previously developed compounds and their atomistic structures with the PD-L1 dimer. Structure-based affinity-based virtual screening of small-molecule inhibitors of the PD-1/PD-L1 pathway according to the pharmacophore model followed by a screening in terms of drug-likeness resulted in ten hit compounds of high affinity towards the PD-L1 dimer and the satisfaction to all of the drug-likeness rules. Molecular dynamics (MD) simulations showed that nine of ten compounds formed stable complexes with the PD-L1 dimer as evidenced by the analysis of MD trajectories. Molecular mechanics Poisson- Boltzmann surface area (MM-PBSA) calculation revealed very low binding energies (<-46 kcal/mol) for the interactions of these ligands with the PD-L1 dimer, suggesting that identified compounds may serve as good scaffolds for the design of novel agents of antitumor immunotherapy able to target the PD-1/PD-L1 interaction Communicated by Ramaswamy H. Sarma


Introduction
PD-1 (Programmed cell Death protein 1) and its ligand PD-L1 (Programmed cell Death Ligand protein 1) (PD-L1) are transmembrane immunosuppressive checkpoint proteins (Dermani et al., 2019). PD-1 is expressed on activated T cells, and PD-L1 is expressed on somatic tissue and antigen-presenting cells (Bellmunt et al., 2017;Patsoukis et al., 2020). The binding of PD-L1 to PD-1 causes a downregulation of T cell proliferative gene expression and disruption of T cell cytotoxic activity. The blockade of the interaction of PD-1 and PD-L1 can prevent immune evasion of tumor cells and significantly prolong the survival of cancer patients. Currently marketed drugs targeting PD-1/PD-L1 pathway are all monoclonal antibodies (mAbs) that have achieved great success in clinical trials (Chen et al., 2019, McBride et al, 2021Peters et al, 2017). Compared to PD-1 inhibitors, PD-L1 inhibitors can reduce the incidence of side effects resulting from immune disorders (Khunger et al., 2017;Kythreotou et al., 2018;Spagnuolo & Gridelli, 2018), because they are involved, in contrast to PD-1 inhibitors, in a limited number of signaling pathways. The FDA has approved three monoclonal IgG1 antibodies targeting PD-L1, Atezolizumab, Avelumab and Durvalumab (Deng et al., 2016;Kaufman et al, 2016;Stewart et al, 2015). Despite the remarkable therapeutic success of these antibodies, there are considerable drawbacks such as poor bioavailability, immunogenicity, and the high cost of large-scale production. Moreover, compared to peptides and small molecules, the immunogenicity of mAbs can result in severe immune-related adverse events (irAEs) in a number of cases (Naidoo et al., 2015). Small molecule inhibitors have the potential to overcome these limitations. Moreover, the oral bioavailability of small molecules is higher, and the synthesis of small molecules is easier . However, the development of small-molecule inhibitors (SMIs) of the PD-1/PD-L1 pathway is slow; only a few smallmolecule compounds with obvious inhibitory activities against the PD-1/PD-L1 interaction have been identified at the nonclinical study stage, and only one patented PD-1/PD-L1 inhibitor (CA-170) has undergone clinical trials (Perry et al., 2019;Zhan et al., 2016). Of note, compared to peptides, small molecules have advantages in terms of their oral and plasma stability. The study of small-molecule PD-L1 inhibitors has attracted significant interest. However, because of the complexity and plasticity of the PD-L1 surface, it is difficult to design active SMIs targeting PD-L1 (Guzik et al., 2019). Therefore, while many efforts have been made to develop such inhibitors, only a few of them have been reported and patented (Cheng et al. 2018;Kopalli et al., 2019;). First small-molecule inhibitors of PD-L1 have been discovered by scientists at Bristol-Mayers-Squibbs in 2016. More recently crystal structures of the complexes of PD-L1 with its small-molecule ligands have been provided (PDB IDs: 5J89, 5J8O, 5N2D, 5N2F, 5NIU, 5NIX) Skalniak et al., 2017;Zak et al., 2016). These crystal structures show that the interaction of PD-1/PD-L1 was blocked by the small molecules, which induced PD-L1 dimerization, making PD-L1 non-competent for binding to PD-1 (Zak et al., 2016).
These high-resolution data lay the foundation not only for understanding the mechanisms of small-molecule PD-1/PD-L1 inhibitors action, but also for developing novel efficient SMIs by direct methods of computer-aided drug design. In this study, an integrated computational approach to in silico drug discovery was carried out to discover small drug-like compounds able to show structural and functional mimicry of known small-molecule inhibitors of the PD-1/PD-L1 interaction. This computer-based approach included: (i) generation of a pharmacophore model representing 3 D arrangements of chemical functionalities that make ligands active towards the binding site of PD-L1 dimers; (ii) shape/ pharmacophore-based identification of the drug candidates by a virtual screening allowing one to search for small molecules based on their structural and chemical similarity to other active small molecules (Sunseri & Koes, 2016); (iii) elimination of compounds not satisfying drug-like properties; (iv) molecular docking of drug-like compounds to the PD-L1 dimer binding site; (v) calculation of affinities and ranking by affinities; (vi) molecular dynamics (MD) simulations of the identified compounds bound to the PD-L1 dimer binding site; (vii) affinity calculations of the optimized models; (viii) selection of the most promising compounds for biochemical assays. As a result, several compounds with expected high affinity binding to the PD-L1 dimer and good ADMET properties were identified. These compounds are suggested to form promising scaffolds for the development of novel and potent drugs against cancer.

Generation and validation of the pharmacophore model
In this study, the crystal structures of the PD-L1 dimer in complex with BMS-8 (PDB code 5J89), BMS-202 (PDB code 5J8O), BMS-1001 (PDB-code 5NIU), BMS-1166 (PDB code 6R3K) were used by the Pharmit server software (Sunseri & Koes, 2016) to generate the structure-based pharmacophore model for the PD-L1 dimer-small-molecule ligand high affinity interaction. The validation of the pharmacophore model was carried out using a directory of useful decoy À enhanced (DUD-E) dataset (Mysinger et al., 2012), which has been frequently applied as a benchmark to test proposed docking/scoring methodologies. For the validation, the set of 1957 decoys basing on 40 active compounds from US-patents with the established inhibitory activity in respect to PD-1/PD-L1 interaction was formed with the use of the DUD-E dataset (Mysinger et al., 2012). The chemical structures of the 40 active compounds are displayed in Figure S1. The screening of these ligands aided by ROC AUC (area under the curve, receiver operating characteristics curve) analysis was applied for estimation of the capability of the pharmacophore model to distinguish active compounds from decoy ligands. In addition to pharmacophoric points derived from ligand atoms, we used the atomic coordinates of the structural model of the complex of the PD-L1 dimer from the PDB-structure 5NIU with the US-patent WO2017 087777 compound #8, obtained using the Glide docking program , to define by the Pharmit server the 'excluded volumes' which the ligand is not allowed to occupy. This protein-ligand complex was chosen based on the results of computational experimentation (data not shown) as imposing the least conformational limitations on the structures of PD-L1 dimer-small molecule ligand complexes.

Pharmacophore-based virtual screening (PBVS)
Pharmacophore-based virtual screening (PBVS) of the library of small-molecule compounds in terms of high binding affinity to the PD-L1 dimer and according to the generated pharmacophore model was applied to identify new chemical agents able to stabilize the PD-L1 dimer with high efficiency. PBVS was performed in the PubChem molecular library containing over 90 million chemical structures and about 450 millions of conformations using the PHARMIT web server, resulting in an identification of a set of compounds satisfying the pharmacophore model. Then, the PD-L1 dimer -ligand complexes were subjected to a fast free energy minimization using the minimization program Smina of AutoDock Vina (Forli et al., 2016) implemented in Pharmit, and the complexes with energies lower than À8 kcal/mol were saved for further screening.
With SwissADME, Bioavailability Radar tool taking into account six physicochemical properties, such as lipophilicity, size, polarity, solubility, flexibility and saturation, were used at first as a rapid approach to exclude compounds incompatible with an acceptable pharmacokinetics profile. Next, additional pharmacokinetics parameters, such as skin permeation (Kp), passive human gastrointestinal absorption (HIA), bloodbrain barrier (BBB) permeation,binding to the permeability glycoprotein (P-gp) and cytochromes P450 (CYP) were predicted as restrictors. Additionally, five different rule-based filters, with diverse range of properties inside of which the molecule is defined as drug-like., were used to exclude the compounds incompatible with oral bioavailability. These were the Lipinski filter (Lipinski et al., 2001), Ghose filter (Ghose et al., 1999), Veber filter (Veber et al., 2002), Egan filter (Egan et al., 2000), Muegge filter (Muegge et al., 2001). Multiple estimations allow consensus views or selection.
Besides, a list of 105 fragments identified in (Brenk et al., 2008) was used to exclude the molecules containing the fragments to be putatively toxic, chemically reactive, metabolically unstable or to bear properties responsible for poor pharmacokinetics 2.4. Molecular docking based screening Subsequent screening of the compounds selected at the pharmacophore-based stage was carried out by molecular docking. The rationale behind performing molecular docking is to make a systematic prediction of the ideal pose or conformation of a ligand in a protein's binding site, which could be taken further for molecular dynamics simulation studies. The Protein Preparation Wizard module (Madhavi Sastry et al., 2013) of Schr€ odinger software package, version 2017-2 (Schr€ odinger Release 2017-2) was used to prepare the active PD-L1 dimer state (PDB codes: 5NIU, 6NM7). In so doing, the optimization of ionization state, for predicting pKa values in proteins (pH 7.0 ± 2.0) and orientations of side chain functional groups (e.g., hydroxy group in Ser, Thr and Tyr; amino group in Asn and Gln), as well as the addition of missing hydrogens, removal of water molecules and assignment of the right bond order was performed. A restrained minimization of the complex was then carried out (cutoff_ root mean square deviation (RMSD) 0.3 Å), so as to avoid any steric clashes. The compounds saved after the pharmacophore filtering were prepared using the LigPrep module of Schr€ odinger software package, version 2017-2 (Madhavi Sastry et al., 2013; Schr€ odinger Release 2017-2).
In the SIFD protocol, docking simulations were performed as three steps: initial docking, receptor refinement, and redocking. Ring conformational sampling with a 2.5 kcal/mol energy barrier, as well as a non-planar conformation penalty on amide bonds was applied to the IFD protocol. The scaling for both receptor and ligand was set at 0.5 with a maximum of 10 allowable poses per ligand. Residues within 5 Å of the docked ligand were further refined using Prime Refinement algorithm implemented in Maestro v11.2. Prime energy was used to rank the refined protein-ligand complexes. The receptor structures within 30 kcal/mol of the minimum energy structure were submitted for a final round of Glide docking and scoring. Each ligand was re-docked into every single refined low-energy receptor structure in the subsequent second docking step using the default Glide XP settings. For HTVS, only one pose was considered, whereas for the II-IV steps, 10 poses were generated for each molecule (with all parameters at their default values and by employing the OPLS3 force field) (Jorgensen et al., 1996) [54]. The induced-fit docking by the RosettaLigand program that employs the Monte Carlo minimization protocol (Davis & Baker, 2009;Lemmon & Meiler, 2012). was applied to take into account the induced fit processes of the receptor backbone and side chains.
The RosettaLigand induced fit docking (RLIFD) comprised three stages, which progressed from low-resolution conformational sampling and scoring to full atom optimization using all-atom energy function. In the first, low-resolution stage, the end state of the receptor/ligand complex that was obtained with Schr€ odinger Induced Fit Docking was used as a starting configuration of each complex, where ligand was perturbed 500 times and rotated 1,000 times in a random direction. The second, high-resolution stage employed the Monte Carlo minimization protocol in which the ligand position and orientation were randomly perturbed by a small deviation (0.05 Å and 3 ); receptor side chains were repacked using a rotamer library. The ligand position, orientation, and torsions and protein sidechain torsions were simultaneously optimized using quasi-Newton minimization; and the end result was accepted or rejected based on the Metropolis criterion. Scoring used the full-atom Rosetta energy function with softened van der Waals repulsion (Davis & Baker, 2009;Lemmon & Meiler, 2012). The sidechain rotamers were searched simultaneously during 'full repack' cycles and one at a time in the 'rotamer trials' cycles. The full repack made 1000 random rotamer substitutions at random positions and accepted or rejected each based on the Metropolis criterion. Rotamer trials chose the single best rotamer at a random position in the context of the current state of the rest of the system, with the positions visited once each in random order. 12 cycles of rotamer trials and a full repack after every three cycles were performed. The third and final stage was a more stringent, gradient-based minimization of the ligand position, orientation, and torsions and receptor torsions for both side chains and the backbone. Binding site for the Glide docking phase was that of the complexes PD-L1 dimer-BMS-1001 (5NIU) and PD-L1 dimer-BMS-36 (6NM7) (Zak et al., 2016). For the I-III steps, the docking poses were subjected to a pose filtering by Glide docking scores (DG Glide ) computed using Maestro (Schr€ odinger, LLC:, 2017). Molecules with a DG Glide value lower than À9, À11, À13.5 kcal/mol were saved after HTVS, SP and XP docking, respectively for their use in subsequent steps. For the hit selection after the Schr€ odinger IFD and the RosettaLigand IFD, the compounds were discarded if the MM/GBSA DGbind was higher than À30 kcal/mol or RosettaLigand Interface score (Isc) (calculated for the lowest total energy (E tot )) of each complex was higher than À12 REU ('Rosetta Energy Units'). Two previously reported PD-L1 inhibitors BMS-8 and BMS-1166 were also submitted to all docking stages, as reference compounds and the choice of cutoff values was made in accordance with those for the PD-L1 dimer -BMS-8 reference complex which possesses the worst IC 50 value (146 nM) among BMS-inhibitors of PD-L1 currently chartacterized experimentally (Guzik et al., 2019).

Molecular dynamics (MD) simulations
Whereas molecular docking allows the structure-based-knowledge-based guesses about most likely binding modes, the dynamical nature and stability of the noncovalent intermolecular interactions are not considered. To assess the consistency, stability and suitability of the hits extracted from virtual screening, molecular dynamics (MD) simulations are frequently used. In this work, MD simulations were carried out using the GROMACS À2019.4 software (Pronk et al., 2013) with the CHARMM36 all-atom force field (Huang & MacKerell, 2013). The protein was solvated using the TIP3P water model (Jorgensen et al., 1983) in a dodecahedron box of 130 nm 3 volume and counterions were added to keep systems neutral. Periodic boundary conditions were applied and Lennard-Jones interactions were truncated at 12 Å with a force switch smoothing function from 10 to 12 Å. The integration time step was 2.0 fs. The nonbonded interaction lists were generated with a cutoff distance of 16 Å. The V-rescale thermostat (Bussi et al., 2007) was used to maintain the temperature at 300 K and the Parrinello-Rahman barostat for keeping the pressure at 1 bar (Parrinello & Rahman, 1981). Electrostatic interactions were calculated explicitly at a distance smaller than 1.0 nm, long-range electrostatic interactions were treated by particle-mesh Ewald summation at every step (Darden et al., 1993). After a 50 000-step steepest descent minimization with maximum tolerated force of 2.39 kcal/mol/Å, the entire systems were then progressively heated up to 300 K on a time scale of 100 ps (under NVT ensemble), followed by a 100-ps isothermal-isobaric ensemble (NPT) density equilibration. After a subsequent 1-ns NPT simulation as equilibration, the production simulations were run for 300 ns in the NPT ensemble with the Verlet leap-frog algorithm coupled with the particle mesh-Ewald method for long-range electrostatics and Verlet cut-off with the distance of 1 nm for short-range interactions (Essmann et al., 1995). GPU-acceleration was used in all MD simulations. MD simulations in duplicate were undertaken for each selected complex starting with the same RosettaLigand final structures and using in each case two different initial velocity assignments (seeds). The values of root mean square deviations (RMSD) for protein backbone in respect to starting structures and for ligands in respect to the protein binding site together with the number of hydrogen bonds formed between protein and ligands were monitored along the whole MD simulations to understand the relative stabilities of the complexes. With ligand RMSDs, the residues participating in key interactions with ligands in structural models obtained by docking calculations were grouped in sets of binding site residues for their use as reference structural groups in the ligand RMSD calculations. These residues were determined with the use of the protein-ligand interaction profiler PLIP (Salentin et al., 2015). To explore the variations of protein flexibility, the values of root mean square fluctuation (RMSF) reflecting the time-averaged fluctuations of amino acid residues along the protein were calculated.

Analysis of interaction modes and protein-ligand binding affinity calculation
Five different approaches, three knowledge-based and two physics-based, were used to calculate the binding affinity of ligands to the PD-L1 dimer. These were (i) Glide-based protein-ligand binding affinity evaluation  (v) the affinity evaluation by the physics-based MM/PBSA method (Kollman et al., 2000;Kumari et al., 2014). When using the MM/PBSA method, 100 snapshots in the range between 150 and 200 ns at 0.5 ns-intervals were considered, with the solute dielectric constant (e solute ) set to 2. General characterization of protein-ligand interactions was carried out using protein-ligand interaction profiler PLIP (Salentin et al., 2015) applied to the representative (lowest energy) structures of the complexes. The time evolution of the noncovalent interactions of the receptor À ligand complexes were analyzed over 200 ns of MD simulation using PyContact software (Scheurer et al., 2018) In this study, two compounds of BMS-series, BMS-8 and BMS-1166, were taken as the reference ones, as possessing binding affinities in the nanomolar range with BMS-1166 being the best inhibitor in the series in this range (IC 50 ¼ 1.4 nM) and BMS-8 being the worst one (IC 50 ¼146 nM)

Generation and validation of the pharmacophore model
The Pharmit server software (Sunseri & Koes, 2016) was used to generate the structure-based pharmacophore model for high-affinity protein-ligand complexes between the PD-L1 dimer and small-molecule ligands. The best obtained pharmacophore model consists of three aromatic rings, one hydrogen bond donor, one hydrogen bond acceptor, and two hydrophobic points (Figure 1). Key parameters of the obtained pharmacophore model are given in Table 1. Hydrophobic regions H2 and H4 correspond to a hydrophobic channel which is formed upon PD-L1 dimerization. Aromatic rings R1 and R3 are localized in this channel and feature the biphenyl scaffold which is typical of all known PD-1/PD-L1 inhibitors. The channel dimension imposes strong restrictions on the conformational flexibility of the ligand portion localized in the channel with pharmacophore tolerance radii 1.5 and 1 Å for aromatic rings and hydrophobic regions, respectively. To allow for a greater structural variability of peripheral cyclic groups, the tolerance radii for R6 and D7 were increased to 2 and 2.5 Å, respectively. The validation of the obtained pharmacophore model was performed using a directory of useful decoy À enhanced (DUD-E) dataset. In the process of screening of test set of 1997 ligands including 40 active compounds and 1957 DUD-E generated decoys the obtained pharmacophore model could identify 35 active compounds from 40 active compounds and 44 active compounds from 1957 decoys. Thus, the sensitivity, specificity and the proportion of active compounds were 87.5%, 97.2%, and 79.5%, respectively. The ROC AUC has value of 0.92 (see Figure S1 in Supplementary Information file). Together, this suggests a high prognostic power of the obtained pharmacophore model.
The generation of pharmacophore models is a proven useful tool for drug design. Recently, an alternative pharmacophore model has been derived for the search of inhibitors of PD-1/PD-L1 interaction (Mej ıas & Guirola, 2019). Unlike our model, which is based solely on small-molecule ligands inducing PD-L1 dimerization, the pharmacophore model of Mejias and Guirola, which includes seven hydrophobic features and one polar feature, was learned, in addition to small-molecule ligands, on proteins and peptides that disrupt directly the PD-1/PD-L1 interaction. Despite a clear benefit of such approach in the search of direct inhibitors of the PD-1/PD-L1 interaction, the use of such multidirectional restrictors can narrow the search base in respect to the compounds that use solely the strategy of the PD-L1 dimerization induction

Virtual screening
The shape/pharmacophore-based virtual screening of the PubChem molecular library with the use of the Pharmit webserver followed by rapid energy minimization with the smina Pharmit program resulted in 174, 000 small-molecule compounds satisfying the developed pharmacophore model and forming high-affinity (score values < À8 kcal/mol) complexes with the PD-L1 dimer. These compounds were selected and saved for additional structure refinement and analysis from the standpoint of drug-like properties.

QikProp screening
The use of the QikProp (Schr€ odinger, LLC, New York, NY, 2020) filtering resulted in about 70, 000 compounds as druglike ones for further consideration.

Molecular docking based identification of promising hits
HTVS and Standard Protocol (SP) Glide Docking algorithms implemented in Maestro Schr€ odinger were used to screen over 70, 000 compounds obtained at the previous step. In so doing the compounds with the Glide score for PD-L1 dimersmall-molecule ligand interaction higher than À9 kcal/mol were rejected. At the next step of standard precision docking, the compounds with the Glide score higher than À11 kcal/mol were rejected. At the extra-precision docking step, about 100 top-score structures with the Glide score less than À14 kcal/mol were selected and saved for following consideration.

Screening by SwissADME
The four-stage SwissADME filtering was applied to reject the compounds inacceptable in respect to drug-likeness from the set of $100 compounds that were saved after docking. After the application of SwissADME screening, 10 compounds were selected (shown in Figure 2); list of compounds with the values of the most important SwissADME descriptors are shown in Tables S1a and S1b). The results of SwissADME Bioavailability Radar for six physicochemical properties, such as lipophilicity, size, polarity, solubility, flexibility and saturation, are shown in Figure 3. The Boiled-Egg diagram of intuitive evaluation of passive gastrointestinal absorption (HIA) and brain penetration (BBB) in function of the position of the molecules in the WLOGPversus-TPSA referential (Daina & Zoete, 2016) is shown in Figure 4. According to the diagram, all ten compounds that have been selected during the precedent rejection stage (compounds IN1-IN10) are well-absorbed, not accessing the brain and P-gp positive.  Table 2). Interestingly, that despite the fact that the application of RosettaLigand enhances conformational sampling, no lowering of MM/GBSA-based DGbind was observed in a number of cases. Furthermore, a significant elevation in DGbind took place for the case of compound IN7. The complexes of compounds IN1, IN2, IN3, IN4, IN5, IN6,  IN7, IN8, IN9, IN10, BMS-8 and BMS-1166 with the PD-L1 dimer were subjected to 200 ns MD simulations. MD simulations for each of the 12 complexes were performed by two replicas corresponding to two different velocity seeds. The trajectories were analyzed for stability and compared for evaluation of reproducibility. In addition, the complete details of conformations of protein-ligand complexes were observed to reconfirm the results of docking. The stability of all the complexes was studied by root mean square deviation (RMSD), root-mean-square fluctuation (RMSF), and the time dependence of the number of hydrogen bonds formed between the PD-L1 dimer and ligands. In addition, the lifetime analysis of the MD trajectories was used to evaluate the stability of the complexes and convergence of the replicas.

Molecular dynamics simulations of the complexes of the PD-L1 dimers with identified hits
The RMSD values of the protein Ca atoms in respect to starting structures and the ligands relative to protein binding site residues are illustrated in Figures 5 and 6, respectively.  Figure 6A) together with the results of the lifetime analysis of replica trajectories given in Figure S5 suggests that two replicas in this case correspond to the location, most of the time, of the ligand in two different binding sites. Even greater divergence between two replicas took place in the case of IN3. As may be seen from Figure 6, while IN3 behavior is well stabilized upon first replica simulation and keep their binding pose close to that of the starting structure, the situation was drastically changing with the second replica, when the ligand ran away from the protein. The consideration of additional replicas (data not Figure 3. The results of the application of the SwissADME Bioavailability Radar to compounds I-X. The SwissADME Bioavailability Radar enables a first glance at the drug-likeness/bioavailability of a molecule. The pink area represents the optimal range for each properties (lipophilicity: XLOGP3 between À 0.7 and þ 5.0, size: MW between 150 and 500 g/mol, polarity: TPSA between 20 and 130 Å 2 , solubility: log S not higher than 6, saturation: fraction of carbons in the sp3 hybridization not less than 0.25, and flexibility: no more than 9 rotatable bonds. I stands for PubChem-10412069, II stands for PubChem-11418589, III stands for PubChem-18753825, IV stands for PubChem-18754056, V stands for PubChem-54675254, VI stands for PubChem-118588236, VII stands for PubChem-128364226, VIII stands for PubChem-129208611, IX stands for PubChem-136060551, X stands for PubChem-136060618. shown) was not capable to reverse the trend of alternation of the stable behavior with the ligand dissociation.
In order to determine the ligand-binding effect on the stability of protein structure, RMSF analysis for protein backbone atoms was carried out for all 12 complexes above over the 100-200 ns MD simulation range (depicted in Figure S3). in all 12 ligand-protein complexes were showing stability by maintaining the RMSF values less than 0.25 nm, suggesting that these residues are rather rigid in the presence of ligands. Greater fluctuations were typical for the PD-L1 C''D loop (residues 74-85, more than 2.7 nm) and BC loop (residues 44-52; more than 3.5 nm) in both PD-L1 monomers for both replicas, suggesting that these loops are most flexible parts of PD-L1.
The analysis of time dependence of the numbers of Hbonds formed between the ligands and protein (shown in Figure S4) shows that the predicted H-bonding pattern of all ten compounds is well correlated with MD analysis. Except for the first replica of compound IN9 and second replicas of compounds IN6 and IN3, at least one hydrogen bond was formed through entire MD simulations for all considered complexes of the PD-L1 dimer ( Figure S4), with the most stable H-bonding taking place upon the binding of compound IN8.

Calculation of binding affinities of identified compounds towards the PD-L1 dimer
The lowest energy structures for the complexes of compounds IN1, IN2, IN3, IN4, IN5, IN6, IN7, IN8, IN9, IN10, BMS-8 and BMS-1166 with the PD-L1 dimer that were obtained with the use of the MM-PBSA approach and 100 MD-snapshots in the range between 150 and 200 ns were taken as the representative structures of the trimeric PD-L1 dimer-ligand complexes (shown in Figure 7). The obtained binding affinities are given in Table 3. The individual components of interactions are presented in Table S25 of Supplementary Information File. In addition to the physics-based MM-PBSA analysis, we carried out knowledge-based analysis using protein-ligand affinity evaluations by Glide  and Prodigy-Lig (Vangone et al., 2019). Corresponding values for the complexes of the PD-L1 dimer with BMS-8 and BMS-1166 were taken as a positive control. The estimated affinities are given in Table 3. The results of the analysis of the binding modes of the predicted compounds, namely hydrophobic interactions, hydrogen bonds, salt bridges, p-stacking, and p-cation interactions using the protein-ligand interaction profiler PLIP (Salentin et al., 2015) show that common to all 12 compounds are hydrophobic interactions with Ile54A/B, Tyr56A/B, Met115A/B, Ala121A/B and Tyr123A/B of PD-L1. In many cases, the lead compounds formed hydrogen bonds with Gln66 (compounds IN1, IN5, IN6, IN8, IN9), Ala121 (compounds IN1, IN5, IN6, IN8, IN9, IN10) Asp122 (compounds INI,  IN6, IN7, IN8). One case of p-stacking with Tyr56A took place. Notably, the same p-stacking with Tyr56A was detected for the interaction of BMS-8 with the PD-L1 dimer. The detailed results of this analysis are given in Tables S2-S24  The binding energies for all protein-ligand complexes in our study were calculated for the last 50 ns of MD trajectories using 100 snapshots at 0.5 ns intervals. All 10 complexes showed negative binding energies (Table 3) indicating that all the complexes were energetically stable. However, there was not ideal correlation between binding energies obtained using different physics-based and knowledge-based metrics. So, IN5 showed the lowest binding energy in terms of physics-based MM-PBSA (-99.4 kcal/mol) but was only ninth in terms of knowledge-based Prodigy metrics (-18.8 kcal/mol). In addition, IN5 exhibited lower binding energies towards the PD-L1 dimer in terms of the MM-PBSA evaluation than those of reference Figure 4. The diagram of "intuitive evaluation" of passive gastrointestinal absorption (HIA) and brain penetration (BBB) in function of the position of the molecules in the WLOGP-versus-TPSA referential for compounds I-X with the use of Boiled-Egg diagram. The white region is for high probability of passive absorption by the gastrointestinal tract, and the yellow region (yolk) is for high probability of brain penetration. The points are coloured in blue if predicted as actively effluxed by P-gp (PGPþ) and in red if predicted as non-substrate of Pgp (PGPÀ). All studied compounds (1-10) are predicted as well-absorbed but not accessing the brain (in the white) and PGPþ (blue dot). compounds BMS-8 and BMS-1166, but had worse affinity according to Prodigy scoring. The analysis of individual components of the IN5 interaction (Table S25) shows a significant contribution of electrostatic interactions in the low energy of the complex. Of note, in doing so, the intermolecular hydrogen bonds in the lowest energy structure (representative structure) are not formed (see Table S8), suggesting significant contribution of other electrostatic components.
The comparison of the Prodigy DG with individual components of DG MM-PBSA (Table 3) shows a good correlation with the van der Waals component of the latter (R-coefficient is 0.71), suggesting a satisfactory performance of Prodigy metrics in ranking by affinity if solvation is slight. Because none of the metrics is the ultimate in accuracy, whereas all ten compounds exhibit good stabilization of the PD-L1 dimer, there is good reason to consider all these leads in further optimization and biochemical assays.
All ten selected inhibitors are consistent with the pharmacophore model. All of them possess the H-bond acceptor, which closes the hydrophobic channel from the solution (see   pharmacophore model by R6 and D7, respectively (shown in Figure 1 and Table 1).
The detection of non-covalent PD-L1 dimer-ligand contacts with the use of the protein-ligand interaction profiler (PLIP) service (Salentin et al., 2015) shows that all lead compounds possessing a significant binding affinity to the PD-L1 dimer bind at the same binding site where most of them form noncovalent contacts with Tyr56, Ala121, and Tyr 123 of both PD-L1 monomers, as well as with Met115 of the monomer A (Tables S2-S25 in the Supplementary Information File). The 3 D views for the representative (MMPBSA-based lowest energy structures) and 2 D schematic diagrams for the non-covalent PD-L1 dimer-ligand contacts are given in Figure 7.
In this study, two replicas of MD simulations were considered. While for the most complexes of examined ligands with the PD-L1 dimer, DRMSD 1,2 between the lowest energy structures of two replicas were less than 2.5 Å (see Table S26 of Supplementary Information), suggesting that the systems in both replicas are in close proximity to the same energy minima, in the case with IN3 we observed both the strong binding of the ligand with one replica and the detachment of the ligand from the protein in other case not allowing unambiguous conclusion about prediction reliability and longer simulations are necessary to clarify the situation.
Previously, two hot regions on PD-L1 have been detected as key ones in the interaction of the PD-L1 dimer with BMSligands: the first one (referred to as P1) including Tyr56, Glu58, Arg113, Met115, Tyr123 and the second one (P2) containing Met115, Ala121, Tyr123 (Paciotti et al., 2020;Zak et al., 2016). As it follows, from the analysis of the individual components of the interactions of the identified compounds with the PD-L1 dimer, given in Table S27, electrostatic interactions in most cases are significantly compensated by desolvation energy, resulting in a crucial role of hydrophobic interactions in a high protein-ligand affinity. In most cases, Ile54, Tyr56, Met115, Ala121, Tyr123 are involved in the hydrophobic contacts of the identified ligands with the PD-L1 dimer, similarly to hydrophobic interactions of known BMS-ligands (Zak et al., 2016, Paciotti et al., 2020.
3. 9. The used protocol combining SIFD, RLIFD with MD predicts the PD-L1 dimer-BMS-8 complex structure with high precision The results that were obtained in this study for the PD-L1 dimer-BMS-8 complex were compared with the X-ray structure for this complex (Zak et al., 2016; PDB ID: 5J8O), The structural superposition of the X-ray structure on the representative structures (the lowest energy structures) of the PD-L1 dimer-BMS-8 obtained by MM/PBSA analysis for replicas 1 and 2 is shown in Figure 8. The RMSDs of the replicas 1 and 2 PD-L1 dimer-BMS-8 structures and 5J8O were 0.97 and 1.1 Å, respectively, thus suggesting a high precision of the prediction.

Conclusion
Involvement of PD-L1 in an immune evasion of tumor cells makes it an important drug target. In a number of studies it has been shown that some ligands can induce dimerization of PD-L1 monomers resulting in the prevention of PD-L1 interaction with PD-1 and thus preventing PD-L1-mediated immunosuppression. Therefore, identification of compounds which can effciently and stably bind to PD-L1 was undertaken through computer aided drug design (CADD) study. In this paper, we presented a combined study in which the generation of pharmacophore model, virtual screening, molecular docking and screening of compounds in terms of their drug-likeness by QikProp and SwissADME services were used in order to identify the promising candidates as potent inhibitors of the PD-1/PD-L1 interaction. In this process, a large collection of over 90 000, 000 compounds was evaluated against the pharmacophore filtering which helped to identify 174, 000 compounds with potential features to induce PD-L1 dimerization. These molecules were then submitted to drug-likeness screening, followed by molecular docking in an increasing degree of precision, shortlisting ten compounds. Under static conditions, they exhibited a significant degree of binding affinity and important molecular recognitions with the PD-L1 dimer. To evaluate the binding of the identified hits to the PD-L1 dimer under dynamical conditions, MD simulations in duplicate were undertaken for each of the ten complexes. System stability and binding energy analyses helped to identify among these ten compounds nine ones as promising candidates for inhibiting the interaction of PD-L1 with PD-1. Binding free energy calculations by MM/PBSA showed low binding energy values less than À35.4 kcal/mol. Taking into account the results of the physicochemical characterizations and ADME/Tox predictions, one might expect that these nine compounds could be considered as promising scaffolds for the development of novel and potent inhibitors of PD-1/PD-L1 interaction. These molecules can be taken up further for in vitro evaluation of their PD-L1 inhibiting effect. Previously, the compounds IN2, IN3,  IN4, IN5, IN6, IN8 and IN9 have been patented, but none of them has been claimed as an antitumor agent or inhibitor of the PD-1/PD-L1 interaction. It is clear that, despite promising in silico profile, the analyzed compounds are only intermediate points for the development of new highly potent drug candidates. In this connection, these compounds should further go through a lead optimization, iterative process of altering the molecule structure to identify their chemical modifications with improved immunotherapeutical potency and improved ADMET parameters. For this purpose, the modern QSAR methods commonly used as a lead optimization approach in drug discovery may be applied (Golbraikh et al., 2017;Kuseva et al., 2019;Schultz et al., 2018). Thereafter the inhibitory activity of these lead compounds should be further tested in vitro and animal studies. In addition to QSAR, the integrative approaches combining in silico studies with experimental assays may be also beneficial (Amir et al., 2021). The generation of pharmacophore models is a proven useful tool for drug design. Recently, an alternative pharmacophore model has been derived for the search of inhibitors of PD-1/PD-L1 interaction (Mejhas & Guirola, 2019). Unlike our model, which is based solely on small-molecule ligands inducing PD-L1 dimerization, the pharmacophore model of Mejhas and Guirola, which includes seven hydrophobic features and one polar feature, was learned, in addition to small-molecule ligands, on proteins and peptides that disrupt directly the PD-1/PD-L1 interaction. Despite a clear benefit of such approach in the search of direct inhibitors of the PD-1/PD-L1 interaction, the use of such multidirectional restrictors can narrow the search base in respect to the compounds that use solely the strategy of the PD-L1 dimerization induction.
In this study, a combination of induced fit docking using two IFD docking approaches, Schr€ odinger IFD and RosettaLigand IFD, with molecular dynamic simulations was used. In the control study of the PD-L1 dimer-BMS-8 complex, we were able to perform the prediction of complex structure with high accuracy, thus suggesting a considerable promise of taking into account protein flexibility in protein-ligand interaction by such computationally low-cost approach.

Author contributions
V. V conceived of the idea and V. V. and V.U. planned the work for this publication; V. U., A.D., and V. V. carried out the simulations; V. V. and V. U. analyzed the data and wrote the paper. All authors read and approved the final version.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Disclosure statement
No potential conflict of interest was reported by the authors.