Investigation of protein-protein interactions and hot spot region between PD-1 and PD-L1 by fragment molecular orbital method

Inhibitors to interfere protein-protein interactions (PPI) between programmed cell death 1 (PD-1) and programmed death ligand-1 (PD-L1) block evasion of cancers from immune surveillance. Analyzing hot spot residues in PPI is important for small-molecule drug development. In order to find out hot spots on PPI interface in PD-1/PD-L1 complex, we analyzed PPI in PD-1/PD-L1 with a new analysis method, 3-dimensional scattered pair interactions energies (3D-SPIEs), which assorts significant interactions with fragment molecular orbital (FMO) method. By additionally analyzing PPI in PD-1/antibody and PD-L1/antibody complexes, and small-ligand interactions in PD-L1/peptide and PD-L1/small-molecule complexes, we narrowed down the hot spot region with 3D-SPIEs-based interaction map, which integrates PPI and small-ligand interactions. Based on the map, there are two hot spot regions in PPI of PD-1/PD-L1 and the first hot spot region is important for inhibitors. In particular, LY56, LE58, and LN66 in the first hot spot of PD-L1 are important for PD-L1-antibodies and small-inhibitors in common, while LM115 is important for small-inhibitors. Therefore, the 3D-SPIEs-based map would provide valuable information for designing new small-molecule inhibitors to inhibit PPI of PD-1/PD-L1 and the FMO/3D-SPIEs method provides an effectual tool to understand PPI and integrate PPI and small-ligand interactions at a quantum mechanical level.

approaches based on the molecular mechanics (MM) methods 14,15 . Nonetheless, a direct application of the present MM to PPI has limits, because MM methods depend on predetermined parameters optimized to reproduce average values and classical mechanics could not describe the molecular system precisely. The use of quantum mechanics (QM) comes to the fore, because QM could overcome the limitations of MM. Although QM could not be easily applied to the large biological systems, such as PPI, due to huge computational cost 15 , it provides the analysis of the components of molecular interaction in post Hartree-Fock level with energy decomposition analysis (EDA) method, developed by Morokuma et al. 16 . The EDA in QM can describe the important dispersion interactions and non-classical interactions, including CH-π, cation-π, and so on in PPI, which were not easy for MM to describe.
As a practical application of EDA, Kitaura et al. developed fragment molecular orbital (FMO) method for large molecules and molecular systems in 1999 17 , and pair interaction energy decomposition analysis (PIEDA) in 2007 18 . It has advantages of drastically reducing the computing time without compromising the accuracy of the results compared with traditional QM method, and of providing PIEDA as a practical application of EDA. Recent progress in method development and applications was well summarized in review articles of Fedorov et al. 19,20 It has been successfully applied to investigate the protein-ligand interactions based on PIEDA qualitatively and quantitatively 21,22 . In case of PPI as a large molecular system, there have been a few cases to apply FMO method, where a partner protein was defined as a fragment for facilitating PPI analysis and it required huge computational cost 23 .
In this work, we applied FMO method to investigate PPI between PD-1 and PD-L1. In order to find hot spot residues in PD-1/PD-L1 complex, we performed FMO calculation of wild type complex. For analyzing PPI with computational cost being low by using a-fragment-per-a-residue fragmentation scheme, we devised a 3-dimensional scattered pair interaction energies (3D-SPIEs) analysis tool. The 3D-SPIEs method sorts pair interaction energy (PIE) out by using distance between fragments and energy value information, so it provides a 3-dimensional scatter plot and important interactions in PPI. To narrow down hot spot region in PD-1/PD-L1, we performed FMO calculations of PD-1/antibody complexes, PD-L1/antibody complexes, PD-L1/peptide complexes, and PD-L1/small-molecule complexes. For integrating PPI in the antibody complexes and the small-ligand interactions in small-ligand complexes, we made 3D-SPIEs-based interaction map in PD-1/PD-L1. In order to validate PPI predictability of the FMO results, we compared them with the reported site-directed mutagenesis results. Finally, we summarized FMO/3D-SPIEs results as an interaction map, by integrating the information from PPI and small-ligand interactions, and found the hot spot region in PD-1/PD-L1 at a quantum-mechanical level.

Results
All X-ray structures in this work are summarized in Table 1. We used subscript notation 'L' before amino acid to indicate PD-L1, subscripts 'HC' and 'LC' to indicate heavy chain and light chain in antibody, and superscripts, ' A' , 'B' , 'C' , 'G' , 'H' , and 'L' , after water molecules to indicate chain arrangement.

3-dimensional scattered pair interaction energies (3D-SPIEs).
Analyzing protein-protein interactions (PPI) between a target protein and a partner protein can be correlated with analyzing pair interaction energy (PIE) in FMO calculations with a-fragment-per-a-residue fragmentation scheme. Because the calculated PIE often show unphysically large values primarily due to electrostatic interactions 24 , we used polarizable continuum model (PCM) for solvent screening effects and introduced distance information to sort out the meaningful PIEs from overestimated PIEs. As for distance information, we selected only PIEs within a specific distance (5.4 Å) between two fragments, which reflected the distance used for the approximate of electrostatic potential in FMO method 25 . The brief scheme of 3-dimensional scattered pair interaction energies (3D-SPIEs) method is shown in Fig. 1.
Two main hot spot regions in wild-type PD-1/PD-L1 complex. In order to find out amino acids on PPI interface that form hot spots in wild type PD-1/PD-L1 complex, we performed FMO calculation of the wild-type complex and analyzed it with 3D-SPIEs method. It detected 28 PPI between PD-1 and PD-L1, which are summarized in Supplementary Table S1, and 5 water-bridge interactions, which are summarized in Supplementary Table S2. The FMO results were in agreement with the previously reported 13 PPI and 3 water-bridge-interactions of Zak et al. 7 13 PPI: Y68: L D122, Q75: L D26, Q75: L R125, A132: L Q66, I134: L Y56, I134: L E58, I134: L Y123, E136: L R113, E136: L Y123, E136: L R125, N66: L D122, Q75: L I126, and K78: L T20, and 3 water-bridge-interactions: I134:HOH202: L Y56, I134:HOH202: L E58, and N66:HOH203: L A121. Zak et al. reported three hot spots on the PPI interface between PD-L1 and PD-1 7 the first hot spot was L Y56, L E58, L R113, L M115, and L Y123; the second was L M115, L A121, and L Y123; and the third was from L D122 to L R125. PPI in PD-1/PD-L1 complex analyzed with FMO method are delineated in Fig. 2. In the FMO results, there were only two hot spots on the PPI interface. In the second hot spot of Zak et al. 7 , only E136: L Y123 were strong, which can belong to the third hot spot. Therefore, we rearranged the reported three hot spot regions into two hot spot regions in PD-1/PD-L1 complex based on the FMO results. (see Fig. 2) The first hot spot region of PD-1 was comprised of residues from A132 to E136 located in the front sheet of PD-1, which corresponds to the first hot spot region of PD-L1, consisting of L Y56, L E58, L Q66, L R113, and L Y123 located in the opposite β strands to the front sheet of PD-1. The second hot spot region of PD-1 includes N66, Y68, M70, S73-K78, and E84, which corresponds to the second hot spot region of PD-L1, consisting of L A18-L T20, L D26, and L G120-L I126.

PD
Hot spot region which five PD-L1-antibodies occupied in PD-L1. In order to narrow down hot spot regions in PD-1/PD-L1 complex, we performed FMO calculations of five antibody complexes in complex with PD-L1 and analyzed them with 3D-SPIEs method. There were the reported five antibody complexes in complex with PD-L1: BMS-936559, Avelumab, KN035, Atezolizumab, and Durvalumb.  Supplementary Table S9, and 17 water-bridge interactions, which are summarized in Supplementary Table S10.
The FMO results were in agreement with the published site-directed mutagenesis studies 31,32 . The FMO results showed 21 PPI related with the mutagenesis data, where the mutations on L Y56, L E58, L D61, L N63, L Q66, L R113, L M115, L Y123, and L R125 in PD-L1 made the lower binding affinity to KN035. In particular, it was reported that I54A in PD-L1 decreased binding affinity between PD-L1 and KN035. The FMO results showed 4 intra-protein interactions with L I54 in PD-L1: L I54: L L50, L I54: L H69, L I54: L G70, and L I54: L Y118. Because L I54 in β-sheet 4 may stabilize the structure of PD-L1 by interacting with L L50 in loop, L H69 next to β-sheet 5, L G70 next to β-sheet 6, and L Y118 next to β-sheet 3, I54A mutation may make β-sheet flexible and the stability unstable.
PD-L1/Atezolizumab (PDB ID: 5XXY). In PD-L1/Atezolizumab complex, the FMO results detected 40 PPI, which are summarized in Supplementary Table S11. The FMO results were in agreement with the previously The common hot spot region in PD-L1 against five PD-L1-antibodies and PD-1. In order to investigate the common hot spot regions among five PD-L1-antibodies and PD-1, we illustrated the FMO results in Fig. 3. The five PD-L1-antibodies have PPI with L Y56 and L E58 in common, while four of them had PPI with L D61, L Q66, L R113, L Y123, and L R125 in PD-L1. In the wild-type complex, L D61 had no PPI between PD-1 and PD-L1, but had intra-protein interactions with L R113 (−14.636 kcal/mol), which may contribute to the stability in the first hot spot in PD-L1. Some antibodies had PPI with L E60, L H69, and L D122. The FMO results showed that most PPI between PD-L1-antibodies and PD-L1 focused on the PPI on the first hot spot between PD-1 and PD-L1.

Hot spot region which small ligand inhibitors occupied in PD-L1.
In order to narrow down hot spot regions in PD-1/PD-L1 complex, we performed FMO calculations of two peptide complexes in complex with PD-L1 and four small-molecule complexes in complex with two PD-L1 proteins, and analyzed them with 3D-SPIEs method in case of peptide inhibitors and traditional bar graph method in case of small-molecule inhibitors.
PD-L1/peptide-57. In PD-L1/peptide-57 complex, the FMO results detected 9 strong interactions between peptide-57 and PD-L1 with 9 residues and 6 water-bridge interactions: 9 residues: L Y56, L E58, L N63, L Q66, L E71, L E72, L D73, L R113, and L M115 are summarized in Supplementary Table S14: and 6 water-bridge interactions are summarized in Supplementary Table S15. The FMO results were in agreement in with the published mutagenesis data, where peptides which lack W8 and W10 residues had relatively low binding affinity 34  Four small molecule inhibitors in complex with two PD-L1 proteins. There were the reported four nonpeptidic small-molecule inhibitors in complex with two PD-L1 proteins by Bristol-Myers-Squibb (BMS-8, BMS-37, BMS-200, and BMS-202). Because it was reported that these compounds made dimerization of two PD-L1 proteins, we used superscripts ' A' and 'B' to annotate protein molecules depending on the chain arrangement in the crystal structures.
The FMO results of BMS compounds and two PD-L1 proteins are illustrated in Fig. 4, and their interactions are organized in Supplementary Table S18-S21. In BMS-8/PD-L1 complex, the FMO results detected 10 strong interactions with PD-L1 with 9 residues and one water molecule: 9 residues: L Y56 A , L E58 A , L Q66 A , L M115 A , L D122 A , L Y56 B , L M115 B , L S117 B , and L D122 B : and one water molecule: HOH206. In BMS-37/PD-L1 complex, the FMO results detected 10 strong interactions with PD-L1 with 9 residues and one water molecule: 9 residues: L T20 A , L Y56 A , L M115 A , L D122 A , L Y56 B , L E58 B , L Q66 B , L M115 B , and L D122 B : and one water molecule: HOH227. In BMS-200/PD-L1 complex, the FMO results detected 16 strong interactions with PD-L1 with 12 residues and four water molecules: 12 residues: L T20 A , L Y56 A , L M115 A , L S117 A , L D122 A , L Y56 B , L E58 B , L Q66 B , L M115 B , L A121 B , L D122 B , and L Y123 B : and four water molecules: HOH229, HOH250, HOH306, and HOH316. In BMS-202/PD-L1 complex, the FMO results detected 13 strong interactions with PD-L1 with 10 residues and three water molecule: 10 residues: L T20 A , L Y56 A , L M115 A , L S117 A , L D122 A , L Y56 B , L E58 B , L Q66 B , L D73 B , L M115 B , and L D122 B : and three water molecules: HOH301, HOH316, and HOH325.

Hot spot region in PD-1/PD-L1 integrated with 3D-SPIEs-based interaction map. To integrate
information between PPI and small-molecule-protein interactions, we made 3D-SPIEs-based interaction map. (see Figs 3 and 4). In the wild-type complex, we found that there were two hot spot regions. In the first hot spot, (2019) 9:16727 | https://doi.org/10.1038/s41598-019-53216-z www.nature.com/scientificreports www.nature.com/scientificreports/ most of PD-L1-antibodies, peptides, and small-molecule inhibitors have common significant interactions with L Y56, L E58, and L Q66. In the second hot spot, most of PD-L1-antibodies had interactions with L D122 and L R125, while most of peptides and small-molecule inhibitors had interactions with L D122. According to the map, the first hot spot might be the most important hot spot on PPI interface in the PD-1/PD-L1 complex. Notably, L M115, located between the first and the second hot spot regions on PPI interface, had significant interactions with all of peptides and small-molecule inhibitors, while there was no PPI in wild-type complex, while two of the five PD-L1-antibodies had PPI with L M115. In intra-protein interactions in PD-L1 of the wild-type, L M115 interacted with L W57 (−17.157 kcal/mol), L Y123 (−7.154 kcal/mol), and L K124 (−16.566 kcal/mol). It indicated that L M115 may be also important for maintenance of the first hot spot region in PD-L1 and play a role of pivot, connecting the first and the second hot spots, for small ligands.

Discussion
Given the successful prospects in clinical studies of monoclonal antibodies targeting the PD-1/PD-L1 axis, several agents to block the PD-1/PD-L1 axis have been developed. The understanding of PPI between PD-1 and PD-L1 is vital for successful structure-based drug design. Even though the crystal structures of PD-1/PD-L1 complex were solved, it has been still a challenge to find small-molecule inhibitors. Visual inspection or traditional molecular mechanic approaches do not always give a clear answer, especially in protein-protein interactions. In the present study, we made the 3D-SPIEs-based interaction map in PD-1/PD-L1, which integrated PPI and small-molecule-protein interactions at a quantum-mechanical level. In the map, it indicated that the first hot spot region is important for inhibiting PPI in PD-1/PD-L1, which all antibodies and small-ligands focused on. The integrated information can be used in the design of new compounds in hit-to-lead or lead optimization stages for targeting PD-1/PD-L1 in drug discovery.
Water molecules play a crucial role in the protein-protein interactions, though they can be quickly exchanged with bulk solvent 35 . Here, we analyzed water molecules which were resolved in the crystal structures. Using the www.nature.com/scientificreports www.nature.com/scientificreports/ FMO method to analyze the water-bridge interactions in protein-protein system can allow us to identify water molecules which are energetically favorable or unfavorable. And this information can be used to design several agents that can interact with certain water molecules or displace the interactions. However, identification and verification of water molecules depends on the resolution of X-ray crystal structures. The FMO method is suitable for analyzing the roles of already predicted, or identified, water molecules which have been determined by other computational methods 15 .
FMO/PIEDA method has been a practical tool for rational structure-based drug design, because it offers a precise and comprehensive information on the protein-ligand interactions 15,19 . The information from the FMO/ PIEDA method allows us to distinguish whether the interactions are strong, moderate, weak or not, and whether they are attractive or repulsive between the ligand and its neighboring residues. Furthermore, it can be used as an efficient ab initio tool for design, assessment and filtering of new compounds, which may reduce the effort and the cost of chemical synthesis in drug discovery.
Even though the FMO method was successfully applied to the protein-protein interaction of influenza hemagglutinin 19,24,36 , it has not been generally employed in other PPI cases. Although the analysis scheme, summing all PIEs over all residues of a partner protein, facilitated the analysis of PPI with a bar graph, it required the huge computational cost. Here, we introduced the 3D scattered pair interaction energies (3D-SPIEs) to sort and analyze PIEs between fragments, which provided two advantages. Firstly, it allowed the FMO method to be applied easily to protein-protein systems due to lower computational cost from a-fragment-per-a-residue scheme. It would provide a protocol to analyze PPI, which can be applied to protein-RNA and protein-DNA complexes in future. Secondly, most protein-protein interactions, detected by the FMO/3D-SPIEs, method were in agreement with the published reports and mutagenesis data. Not only providing the reported PPI, it had additional advantage of distinguishing inter-protein (PPI) and intra-protein (non-PPI) interactions. In some cases, the alanine scanning experiments could be inconclusive, because the mutations could affect the binding free energy by a mechanism unrelated with PPI 37 . Some mutations might destabilize the unbound state of protein and alter its conformation, so hot spots identified by alanine scanning experiments could be false positives 37 . It indicated that the application of FMO method to PPI can predict site-directed mutagenesis study results, and furthermore provide intra-protein interactions affecting stability. It would provide the starting points to modify proteins and nucleic-acids as a powerful tool in protein-engineering and nucleic-acid-engineering for the higher binding affinity.
In summary, the outcome of this study provided three points. Firstly, as a protocol, the FMO/3D-SPIEs method can be applied to analyze hot spot region of PPI. Secondly, in analyzing case-by-case PPI, the FMO/3D-SPIEs results had a qualitative correlation with site-directed mutagenesis results, and additionally provided intra-protein interactions for stability. Thirdly, the integrated information with PPI and small-ligand interactions in PD-1/PD-L1, the FMO/3D-SPIEs results provided valuable information for hot spot regions in PD-1/PD-L1, and designing new potential small-molecule targeting PD-1/PD-L1. Therefore, the FMO/3D-SPIEs method allowed the FMO method to be used to understand not only the nature of protein-ligand complexes, but also protein-protein complexes. The FMO method can be hereby used to analyze hot spots in PPI, protein-RNA, and protein-DNA, and provide new direction to hit-to-lead, lead optimization, protein engineering, and nucleic-acid engineering in future.  Two-body FMO method was applied to PD-1/PD-L1, PD-1/mAb, PD-L1/mAb, and PD-L1/small-ligand complexes to investigate key interactions at the second order Møller-Plesset perturbation (MP2) 43 and polarizable continuum model (PCM) 44 with 6-31G** basis set (FMO-MP2/6-31G**/PCM level). The implicit solvation model (PCM) was employed with explicitly expressed water molecules present in the X-ray crystal structures. In FMO calculations of PD-1/mAb and PD-L1/mAb, only Fv domains of mAb and crystallographic water molecules around the region were included. All inputs files were prepared with an in-house code in compliance with the hybrid orbital projection (HOP) 45 scheme fragmentation. Each residue and water molecule was defined as a fragment and two cysteine residues forming the disulfide bridge were defined as a fragment (a-fragment-per-a-residue fragmentation scheme).
Two body FMO calculation consists of four sequential steps on (i) fragmentation, (ii) fragment (monomer) self-consistent field (SCF) calculation, (iii) fragment pair (dimer) SCF calculation, and iv) total property, including energy and gradient, evaluation 15 .
(i) Fragmentation is a starting point for FMO method. In a-fragment-per-a-residue fragmentation scheme, each residue of protein, ligand and water molecule can be defined as a fragment except two cysteine residues forming disulfide bond. It should be noted that each fragment, representing amino acid residues, is different from the normal assignment for the amino acid residues 45,46 . In FMO calculations, all the residues in protein are divided at the sp3 bond of Cα position, not peptide bond, based on the hybrid orbital projection (HOP) 45 scheme, where the fragmentation is performed at an atom, not midway of bond 47 in Fig. 1.
If the main-chain carbonyl oxygen contributes to hydrogen bond, the interaction will be shifted to the next amino acid by HOP fragmentation, because the carbonyl group of the amino acid belongs to the next amino acid in the scheme. The scheme was introduced for reducing computational cost significantly and correcting errors from projection operator 19,[45][46][47] . (ii) All the MOs on each fragment (monomer) are optimized by self-consistent-field (SCF) 48 theory in the whole electrostatic field and all the electron densities are self-consistently solved through self-consistent-charge (SCC) iterations 17,24 . (iii) Schrödinger equations for fragment pair (dimer) are self-consistently solved in the same way as for the second step, independently 24 . Dimer fragments are made of two monomer fragments. The difference between the third step (dimer) and the second step (monomer) is that the Hamiltonian operator of dimer includes the electrostatic potential generated by surrounding N-2 monomers and that of monomer includes it generated by surrounding N-1 monomers, where N means the number of fragments. (iv) In total property evaluation step, all MO calculations for monomer and dimer are pieced together to generate whole picture of the system. FMO provides the pair interaction energies (PIEs) between fragments. In pair interaction energy decomposition analysis (PIEDA), PIE or ∆E ij int , between fragments i and j, indicates a sum of four energies including the electrostatic ∆ ( ) 3D Scattered pair interaction energies (3D-SPIEs). The main advantage of FMO is that it can account for the energy and structural fitness that stabilizes the complex by calculating the interaction energy (PIE) between the fragments, which make up the protein-protein or protein-ligand complexes. Therefore, it is necessary to sort and analyze the information obtained from the FMO calculations in an easy-to-understand way.
The FMO results consist of three-dimensional data, whose elements are pair interaction energies and names of two fragments. In traditional and powerful way, the FMO results can be analyzed with two-dimensional bar graphs 23,36 , because name of ligand fragment can be omitted in protein-ligand system. However, names of residue fragments in PPI cannot be omitted because they include crucial information. For facilitating analysis of PPI and representing the results with two-dimensional bar graphs, a partner protein can be defined as a fragment, which requires huge computational cost. In order to reduce computational cost and obtain interaction information of all residue pairs in PPI system, we used a-fragment-per-a-residue fragmentation scheme and devised 3-dimensional scattered pair interaction energies (3D-SPIEs) analysis method. (see Fig. 1) In 3D-SPIEs, we sorted out PIE with distance between fragments and energy value information to find important PPI information from all fragment interaction information. We calculated the distance between fragments with single-linkage clustering, where the distance between two fragments equals to the distance between the two atoms that are nearest apart from each other. To make 3D-SPIEs, we used 3D scatter plot graph implemented in plotly 49 . In plots from 3D-SPIEs, the interactions are marked as dots which have 5 components: x, y, z, energy values, and color. The x and y axes indicate name of fragments (e.g. the residues' name of PD-1 and PD-L1), z axis indicates the distance between two fragments, and the dots have color according to their PIEs values. And then, PIEDA information is shown when dots are hovered.