Understanding the structural and energetic basis of PD-1 and monoclonal antibodies bound to PD-L1: A molecular modeling perspective

BBA - General Subjects 1862 (2018) 576–588

Contents lists available at ScienceDirect

BBA - General Subjects journal homepage: www.elsevier.com/locate/bbagen

Understanding the structural and energetic basis of PD-1 and monoclonal antibodies bound to PD-L1: A molecular modeling perspective ⁎

Danfeng Shia, Shuangyan Zhoua,b, Xuewei Liua, Chenxi Zhaoa, Huanxiang Liub, , Xiaojun Yaoa,c,

T ⁎⁎

a

State Key Laboratory of Applied Organic Chemistry and Department of Chemistry, Lanzhou University, Lanzhou 730000, China School of Pharmacy, Lanzhou University, Lanzhou 730000, China c State Key Laboratory of Quality Research in Chinese Medicine, Macau Institute for Applied Research in Medicine and Health, Macau University of Science and Technology, Taipa, Macau, China b

A R T I C L E I N F O

A B S T R A C T

Keywords: PD-L1 PD-1 Monoclonal antibody Molecular dynamics simulation Virtual alanine scanning mutagenesis

Background: The inhibitors blocking the interaction between programmed cell death protein 1(PD-1) and programmed death-ligand 1(PD-L1) can activate the immune response of T cell and eliminate cancer cells. The crystallographic studies have provided structural insights of the interactive interfaces between PD-L1 and its protein ligands. However, the hotspot residues on PD-L1 as well as structural and energetic basis for different protein ligands still need to be further investigated. Methods: Molecular modeling methods including molecular dynamics simulation, per-residue free energy decomposition, virtual alanine scanning mutagenesis and residue-residue contact analysis were used to qualitatively and quantitatively analyze the interactions between PD-L1 and different protein ligands. Results: The results of virtual alanine scanning mutagenesis suggest that Y56, Q66, M115, D122, Y123, R125 are the hotspot residues on PD-L1. The residue-residue contact analysis further shows that PD-1 interacts with PD-L1 mainly by F and G strands while monoclonal antibodies like avelumab and BMS-936559 mainly interact with PD-L1 by CDR2 and CDR3 loops of the heavy chain. Conclusions: A structurally similar β-hairpin peptide with 13 or 14 residues was extracted from each protein ligand and these β-hairpin peptides were found tightly binding to the putative hotspot residues on PD-L1. General significance: This study recognizes the hotspot residues on PD-L1 and uncovers the common structural and energetic basis of different protein ligands binding to PD-L1. These results will be valuable for the design of small molecule or peptide inhibitors targeting on PD-L1.

1. Introduction Recently, the signaling pathway controlled by programmed cell death protein 1(PD-1) and its cognate ligand programmed death-ligand 1 (PD-L1) have attracted great attention in the anti-cancer drug development [1–5]. It's recognized that the coinhibitory signal upon PD-1 and PD-L1 binding is the basis for the immune evasion of tumors and the blockade of this protein-protein interaction could activate the antitumor immune responses of T cell and eliminate tumor cells among cancer patients [6–10]. Several monoclonal antibodies such as nivolumab, pembrolizumab, BMS-936559 and avelumab, which could competitively block PD-1/PD-L1 interaction, are either approved by FDA or undergoing clinical research [10,11]. The inhibitor development targeting these immune checkpoints is emerging as a promising approach to fight against malignant neoplasms in cancer

⁎

immunotherapy [12–14]. Compared to the rapid development of monoclonal antibodies, the research of the small molecule or peptide inhibitors is still at the very beginning [15–18]. The difficulty of inhibitor discovery targeting on PD-1 or PD-L1 is due to a very hydrophobic, large and flat proteinprotein interaction (PPI) between PD-1 and PD-L1 [19–22]. According to the complex structures of murine PD-1/human PD-L1 (PDB code:3BIK) [20] and human PD-1/human PD-L1 (PDB code:4ZQK) [19], PD-1 interacted with PD-L1 through their CC'FG strands in the Ig variable (IgV) domain. The recent studies suggest that PD-L1 is a predictive biomarker in cancer immunotherapy for the close relationship between the PD-L1 expression quantity and the treatment response of anti-PD-1 directed therapy [23,24]. The discovery of the first class of small molecule inhibitors targeting PD-L1 with validated complex crystal structures [25–27] also shows the possibility of small molecule

Corresponding author. Correspondence to: X. Yao, State Key Laboratory of Applied Organic Chemistry and Department of Chemistry, Lanzhou University, Lanzhou 730000, China. E-mail addresses: [email protected] (H. Liu), [email protected] (X. Yao).

⁎⁎

https://doi.org/10.1016/j.bbagen.2017.11.022 Received 10 September 2017; Received in revised form 13 November 2017; Accepted 29 November 2017 Available online 02 December 2017 0304-4165/ © 2017 Elsevier B.V. All rights reserved.

BBA - General Subjects 1862 (2018) 576–588

D. Shi et al.

Fig. 1. The structure of the human PD-L1/human PD-1 complex. (a, b) The residue sequences of PD-L1 and PD-1 and the residues within secondary structures of helix and β-sheet are colored in yellow and cyan shadow respectively. The 11 mutant residues on PD-1 are labeled out in red. (c) The three dimensional crystal structures of human PD-1 and human PD-L1 complex. The secondary structures of helix and β-sheet are colored in yellow and cyan respectively.

interaction interface. A further understanding of the hotspot residues on PD-L1 and the contribution of the different domains on PD-1 or monoclonal antibodies will be important for the design and discovery of novel inhibitors targeting on PD-L1. Molecular dynamics simulations combined with Molecular Mechanics/Generalized Born Surface Area (MM-GBSA) method [31] have been widely used to study the thermodynamics effect of residue mutations [32–34], protein isoform difference [35,36] as well as the interfacial interaction in protein-protein systems [37–39]. Previous studies mainly focused on the interaction mechanism study of PD-1 binding with PD-L1 or PD-L2 [40,41]. For instance, Viricel et al. constructed the complex models of PD-1/PD-L1 and PD-1/PD-L2 complex by protein-protein docking method and further refined them using molecular dynamics simulation and binding free energy analysis, revealing that the rearrangements of the loops between the C'D strands and the BC strands in PD-1 appeared to play an important function in readjusting the beta sheets of the PD-1V domain to interact with the two different ligands [41]. The dynamics characteristics of the IgV domain of human PD-L1 have also been systematically evaluated by Ahmed et al. using both classical and accelerated molecular dynamics simulation [42]. According to their studies, the specific flexible regions

inhibition of PD-1/PD-L1 interaction. The reported crystal structures of the complex of PD-L1 and its protein ligands provided an initial structural basis for the study of protein-protein interactions, and the experiments have also evaluated the affinity of PD-L1 with different protein ligands [19,28–30]. However, these results cannot provide detailed information about the recognition process between PD-L1 and its protein ligands. Molecular dynamics simulation combined with binding free energy calculation has the following two advantages for the study of interaction between PD-L1 and its protein ligands: firstly, the crystal structures obtained in the experiment is often a single conformational state, which often ignored the potential conformational changes of the proteins. Molecular dynamics simulation can restore the conformational state of the protein complex, which reflects the dynamic characteristics of the protein better. Secondly, it can be seen from the crystal structures that the protein-protein interaction interface is often flat, wide and scattered, which make it difficult for peptide or small molecule design to choose the suitable binding site. The binding free energy calculation and contact analysis based on the molecular dynamics simulations can quantitatively evaluate the contribution of important residues or residue pairs and these quantitative information is helpful to distinguish the key residues or residue pairs on the protein-protein

577

BBA - General Subjects 1862 (2018) 576–588

D. Shi et al.

Fig. 2. The residue interaction networks of interfacial residues in (a, b)PD-L1/PD-1(wild) and (c, d) PD-L1/PD-1(mutant) systems. The interfacial residues on PD-L1 and PD-1 are represented in gray and green sticks respectively in a zooming view. The residues of PD-L1 and PD-1 in the interaction networks are shown in gray nodes and green nodes respectively. The non-covalent interactions including hydrogen bond, van der Waals interaction, salt bridge and π-π stacking are depicted by red dash, gray solid, red solid and cyan solid respectively.

study [43,44] and the interaction network among the different domains of the protein [45–47]. Here, we investigated the protein-protein interfacial interaction by constructing a simplified network model using the initial crystal structures of complex systems. The topology file and node file for these complex were firstly generated in the web server RING [48] and then reconstructed into residue interaction network using Cytoscape 3.3.0 [49]. The crucial interfacial interacting residues between PD-L1 and its ligands were further extracted from the complex interaction network through the plugin RIN analyzer. Four interaction types, namely hydrogen bond, van der Waals, salt bridge and π-π stacking, were recognized and depicted using the relaxed distance thresholds in the web server.

such as C″D loop and BC loop as well as Met115 had significant effect on the conformational state of PD-L1 during the binding process with different ligands. In this study, molecular dynamics simulations were mainly performed to study the binding modes between PD-L1 and four different protein ligands namely wild-type PD-1, mutant PD-1 with 11 mutations (V64H, L65V, N66V, Y68H, M70E, N74G, K78T, C93A, L122V, A125V, A132I), avelumab and BMS-936559. The binding free energy decomposition and virtual alanine scanning mutagenesis provided the energy contribution of interfacial residues and helped to recognize the hotspot residues. We systematically compared the interfacial residues derived from crystal structures and MM-GBSA method and evaluated the importance of interfacial residue pairs in the four complex systems. It is worth mentioning that a structurally similar βhairpin peptide extracted from each ligand protein was found to closely interact with the putative hotspot residues. We expect that the results of this study would provide a good starting point for the development of cyclic peptides or peptide mimics inhibitors targeting on PD-L1.

2.2. Molecular dynamics simulation The four complex systems above were further used to perform molecular dynamics simulations. The initial structures of PD-1 and PDL1 were extracted by the residue sequence 33–146 and 18–132 from the crystal structures respectively. The lacking C′-D loop (sequence 85–92) in PD-L1/PD-1 (wild) system (Fig. 1) was constructed according to the corresponding domain on PD-L1/PD-1(mutant) system. The mutant PD1 with 11 residue mutations namely V64H, L65V, N66V, Y68H, M70E, N74G, K78T, C93A, L122V, A125V, A132I was tested to have a better affinity than the wild-type PD-1 experimentally, and in this study it was used to study the effect of mutations on the interfacial interaction between PD-L1 and PD-1. The initial structure of avelumab used in molecular dynamics simulation only contained Fv domain, while the initial structure of BMS-936559 contained both Fv and Fc domain. All the four complex systems were firstly prepared through structural inspection

2. Material and methods 2.1. Residue interaction network analysis on the crystal structures of PD-L1 in complex with different ligands The initial crystal structures of human PD-L1 in complex with wildtype PD-1, mutant PD-1 (V64H, L65V, N66V, Y68H, M70E, N74G, K78T, C93A, L122V, A125V, A132I), avelumab and BMS-936559 were derived from the Protein Data Bank (PDB codes:4ZQK [19], 5IUS [28], 5GRJ [29], 5GGT [30]). Recently, the residue interaction network analysis has been successfully applied in the protein-ligand interaction

578

BBA - General Subjects 1862 (2018) 576–588

D. Shi et al.

Fig. 3. The residue interaction networks of interfacial residues in (a, b) PD-L1/avelumab and (c, d)PD-L1/BMS-936559 systems. The interfacial residues on PD-L1, the heavy chain and the light chain of monoclonal antibody are represented in gray, blue and cyan sticks respectively in a zooming view. The residues of PD-L1, the heavy chain and the light chain of monoclonal antibody in the interaction networks are shown in gray, blue and cyan nodes respectively. The non-covalent interactions including hydrogen bond, van der Waals interaction, salt bridge and π-π stacking are depicted by red dash, gray solid, red solid and cyan solid respectively.

and optimization in Schrödinger software suite (Schrödinger, LLC: New York, NY, 2015). Then, the complex proteins were neutralized with Na+ ions and solvated in a rectangular box of TIP3P water box. The periodic boundary conditions were setup with all the solvents at least 10 Å away from the complex and the solvated systems were parameterized using the AMBER FF14SB force field [50]. The molecular dynamics simulations were performed in four steps. Firstly, energy minimization was performed to remove the local atomic collision in the systems. Both the descent steepest method and the conjugated gradient method with 5000 steps were adopted in the energy minimization. Then in the heat step, each system was gradually heated from 0 K to 300 K in the NVT ensemble with all the solute atoms constrained with a force constant of 2.0 kcal mol− 1 A− 2. After that, each system was equilibrated with the force constant decreasing from 2.0 to 0 kcal mol− 1 A− 2 in 1 ns. Finally, a production run of 200 ns was performed for each system in the NPT ensemble at 300 K with 1.0 atm pressure. The snapshots for all the trajectories were saved per 2 ps.

(2)

ΔH = ΔEgas + ΔEsol = ΔEpolar + ΔEnonpolar

(3)

ΔEgas = ΔEint + ΔEele + ΔEvdW

(4)

ΔHsol = ΔEele, sol + ΔEnonpl, sol

(5)

ΔEnonpl, sol =

γ ∗ΔSASA

(6)

ΔEpolar = ΔEele + ΔEele, sol

(7)

ΔEnonpl = ΔEvdW + ΔEnonpl, sol

(8)

500 snapshots were extracted from the last 20 ns trajectories and used for MM-GBSA calculation. The parameter settings during MMGBSA calculation were referred to the previous works published by our group [36–38]. The enthalpy (ΔH) was further divided into the polar energy contribution (ΔEpolar) and the non-polar energy contribution (ΔEnonpl) as shown in Eqs. (7) and (8). Then, the per-residue based decomposition was performed to identify the key residues in four complex systems. Finally, the entropy change (− TΔS) was calculated by 10 snapshots from the last 20 ns trajectory due to the expensive calculation cost.

2.3. MM-GBSA calculation For each complex system, the energy contributions were calculated for PD-L1/Ligand complex, PD-L1 and Ligand protein respectively, and the binding free energy was calculated according to the equation below:

ΔG = < GPD − L1/ Ligand − GPD − L1 − GLigand >

ΔG = ΔH − TΔS

2.4. Virtual alanine scanning mutagenesis

(1)

The virtual alanine scanning mutagenesis [51] was used to recognize the hotspot residues among the interfacial residues on PD-L1 by estimating the energy contribution of single residue in the PD-L1/PD1(wild) system. The basic process includes two steps, the generation of

where < > represents the average value for all the snapshots used for MM-GBSA calculation. Different energy terms can be estimated as follows:

579

BBA - General Subjects 1862 (2018) 576–588

D. Shi et al.

Fig. 4. The RMSDs and RMSFs of molecular dynamics trajectories for all the complex systems. (a) The RMSDs of all the non‑hydrogen atoms in PD-L1/ligand complex. (b) The RMSDs of the non-hydrogen atoms on PD-L1. (c) The RMSDs of the non-hydrogen atoms on different protein ligands. (d) The RMSFs of residues on PD-1(wild) and PD-1(mutant). (e) The RMSFs of residues on the heavy chain of avelumab and BMS-936559. (f) The RMSFs of residues on the light chain of avelumab and BMS-936559.

alanine during virtual alanine scanning mutagenesis. The contribution of the interfacial residues were evaluated based on the 500 snapshots of PD-L1/PD-1(wild) system and the enthalpy contribution difference between the wild type and the mutant was calculated as follows:

mutated snapshots and the calculation of the binding free energy difference between the wild-type and mutant complex. Because of the structural limitation of glycine and alanine, only the non-glycine and non-alanine residues among the interfacial residues were mutated to

ΔΔH = ΔHmutant − ΔHwild

Table 1 The calculated (MM-GBSA) and experimental Gibbs free energies, as well as the buried surface areas for PD-1(wild), PD-1(mutant), avelumab and BMS-936559 bound to PD-L1. Terms

PD-1(wild)

PD-1(mutant)

Avelumab

BMS-936559

ΔEvdw,gasa ΔEele,gasb ΔEnonpl,solc ΔEpolar,sold ΔEnonpolar,totale ΔEpolar,totalf ΔHg -TΔSh ΔGcali ΔGexpj Buried Area(Å2)

− 97.82 − 196.21 − 14.00 248.12 − 111.82 51.91 − 59.91 48.81 − 11.10 − 5.7 1970

−93.97 −291.73 −14.27 332.44 −108.19 40.71 −67.53 53.72 −13.81 −10.3 ND

− 81.71 − 243.08 − 11.22 276.94 − 92.93 33.86 − 59.07 42.25 − 16.82 − 14.1 1856

− 76.96 − 161.09 − 9.11 210.18 − 86.07 49.09 − 36.97 ND ND ND 1349

(9)

2.5. Residue-residue contact analysis The residue-residue contact analysis has been used to recognize the key residue pairs in protein folding [52] or uncover the protein motions upon substrate binding [53]. It's assumed that these residue-residue contacts exist as long as the distance between two sidechain atoms from different residues was within 3.5 Å. In this study, we performed residue-residue contact analysis to detect the key interactive residue pairs between PD-L1 and its four protein ligands. According to the frequency of occurrence, the stable contact between two residues can be recognized and the recognized residue pair was suggested to have an obvious contribution for the stability of the formation of interactions. All the calculations were performed in software VMD 1.91.

a

Van der Waals energy. Electrostatic energy. c Non-polar solvation free energy. d Polar solvation free energy. e Total non-polar binding free energy. f Total polar binding free energy. g Enthalpy change. h Entropy change. i Calculated Gibbs free energy. j Experimental Gibbs free energy derived from ref. 28, 29 and translated by the equation ΔG = − RT ln (1/Kd) at 298.15 K. All the energy terms are in Kcal/mol. ND means “not determined”. b

3. Results and discussion 3.1. The interfacial residue recognition from the crystal structures The simplified residue interaction networks of the interfacial residues between PD-L1 and its protein ligands are constructed by a series of nodes and edges. As shown in Fig. 2a, the interfacial residues between PD-L1 and PD-1(wild) are composed of PD-L1F19, PD-L1I54, PDL1Y56, PD-L1Q66, PD-L1R113, PD-L1M115, PD-L1S117, PD-L1A121, PD580

BBA - General Subjects 1862 (2018) 576–588

D. Shi et al.

Fig. 5. The per-residue energy contribution spectrums of PD-L1 in complex systems: (a) PD-L1/PD-1(wild), (b) PD-L1/PD-1(mutant), (c) PD-L1/avelumab, (d) PD-L1/BMS-936559. For each complex system, the molecular surface of PD-L1 is shown in white, while the recognized key residues with interaction energies less than − 2 Kcal/mol are colored in red.

L1D122, PD-L1Y123, PD-L1K124, PD-L1R125

and PD-1 V64, PD-1 N66, PDT76, PD-1D77, PD-1 K78, PD-1I126, PD-1 L128, PD1 K131, PD-1A132, PD-1I134, PD-1E136. According to the residue comparison between PD-L1/PD-1(wild) and PD-L1/PD-1(mutant) systems in Fig. 2a and b, the mutations on PD-1(mutant) including V64H, N66V, Y68H, M70E, K78T, A132I are involved in the interfacial interactions between PD-L1 and PD-1(mutant). Although with these mutations, most interfacial residues on PD-L1 are retained except PD-L1A18, PD-L1F19 and PD-L1T20 at the N-terminal. As for the PD-L1/monoclonal antibody complex systems (Fig. 3a, b), PD-L1 interacts with both heavy chains and light chains of the Fv domains on monoclonal antibodies and heavy chains seem to have bigger contribution than light chains, which can be obviously seen from the number of the involved residues. The interfacial residues of PD-L1 upon avelumab binding include PD-L1Y56, PDL1E58, PD-L1E60, PD-L1D61, PD-L1K62, PD-L1N63, PD-L1Q66, PD-L1V76, PDL1R113, PD-L1M115, PD-L1S117, PD-L1Y123. And the interfacial residues of PD-L1 upon BMS-936559 binding include PD-L1D49, PD-L1A51, PD-L1I54,

PD-L1Y56,

PD-L1H69, PD-L1M115, PD-L1S117, PDA comparison of the interfacial residues on PD-L1 for four complex systems (Table S1) shows that the residues PD-L1Y56, PD-L1Q66, PD-L1M115, PD-L1S117 are common for all the complex systems, which suggests that the interfacial sites on PD-L1 are partly conserved when binding with different protein ligands. Considering these interactions between PD-L1 and its protein ligands are intimately affected by the effect of solvents and the structural flexibility of protein ligand, the conformational changes of the complex systems are further analyzed by molecular dynamics simulations.

1Y68, PD-1Q75, PD-1

PD-L1Q66,

PD-L1V68,

L1G119, PD-L1G120, PD-L1A121.

3.2. Stability of the simulated systems 200 ns molecular dynamics simulations were performed to explore the equilibrated states for all the four complex systems. As shown in Fig. 4a, the root-mean-square deviations (RMSDs) of the non‑hydrogen atoms in complex show that all the systems reach a convergence. The

581

BBA - General Subjects 1862 (2018) 576–588

D. Shi et al.

1(wild) system, consistent with the results from the reported literatures [28]. As for monoclonal antibodies complex systems, the van der Waals interaction (ΔEvdw,gas) and electrostatic interaction (ΔEele,gas) in PD-L1/ BMS-936559 were both weaker than PD-L1/avelumab probably due to the less buried surface. The discussion above suggests the MM-GBSA method could still be a reliable tool to evaluate the interfacial interaction between PD-L1 and its ligands. 3.4. The key residues recognition on PD-L1 through per-residue energy decomposition and virtual alanine scanning mutagenesis The energy contribution of the residues on PD-L1 for all the 4 systems were firstly evaluated by per-residue free energy decomposition. The residues with the energy contribution value lower than − 2 Kcal/ mol are labeled out as key residues. As shown in Fig. 5a, the key residues on PD-L1 in PD-L1/PD-1(wild) system are mainly located at the N-terminal (PD-L1A18, PD-L1F19), C strand (PD-L1I54, PD-L1Y56), F strand (PD-L1M115) and G strand (PD-L1A121, PD-L1Y123, PD-L1R125). We compared the key residues in the other three complex systems to PD-L1/PD1(wild) system. In PD-L1/PD-1(mutant) system (Fig. 5b), the energy contributions of PD-L1A18, PD-L1F19, PD-L1I54, PD-L1A121 are obviously weakened while the energy contributions of PD-L1Q66, PD-L1D122, PDL1R125 are greatly enhanced. Specially, PD-L1R125 tends to have the strongest energy contribution in PDL1/PD-1(mutant) system. As for PDL1/avelumab system (Fig. 5c), the key residues are mainly located in C strand (PD-L1I54, PD-L1Y56), C-C′ loop (PD-L1D61, PD-L1N63), C′-D loop (PD-L1V76) and F strand (PD-L1M115), indicating that avelumab does not interact closely with the G strand which is important for PD-1 binding. As shown in Fig. 5d, the key residues are mainly located at B-C loop (PDL1A51), C strand (PD-L1I54), C′ strand (PD-L1H69), F strand (PD-L1M115), G strand (PD-L1G120, PD-L1A121) suggesting that BMS-936559 interacts with PD-L1 at a similar site to PD-1(wild) and PD-1(mutant). According to the distribution of key residues in all the four complex systems (Fig. 5), it can be seen that C strand and F strand are the main interactive domains of PD-L1 in all the 4 systems. We further compared the recognized residues of PD-L1 derived from crystal structures and MM-GBSA energies in all the 4 systems. As shown in Table S1, the recognized key residues are almost completely included by the interfacial residues derived from crystal structures except for PDL1A18 in PD-L1/PD-1(wild) system and PD-L1I54 in PD-L1/avelumab system. Among the common interfacial residues PD-L1Y56, PD-L1Q66, PDL1M115, PD-L1S117, only PD-L1M115 was recognized as key residue in all the four systems; PD-L1Y56 appeared as key residue in three systems (not in PD-L1/BMS-936559 system); PD-L1Q66 was recognized as key residue in PD-L1/PD-1(mutant) system; PD-L1S117 didn't appear as key residue in all the four systems. The importance of PD-L1M115 and PD-L1Y56 was also elucidated by the crystal structures of BMS inhibitors as the side chain of PD-L1Y56 and PD-L1M115 had obvious rearrangements upon BMS-202 binding comparing to PD-L1/PD-1(wild) complex [54]. Recently, the work by Ahmed et al. also proposed that the side chain of Met115 turned to control the open or closure states of the cryptic pocket formed by residues Ile54, Tyr56, Met115 and Ser117 by monitoring its side-chain conformations at the apo states of PD-L1 [42]. The obvious energy contributions of M115 in four different complex systems further validated this mechanism upon the binding of different protein ligands. The polar and non-polar energy contributions of all the recognized key residues in four systems are shown in Fig. 6a and b. It can be seen that the non-polar energy contributions are the main driving force for all the four systems and involved in almost all the recognized interfacial residues while the polar energy contributions are only dominant in PDL1/PD-1(wild) and PD-L1/PD-1(mutant) systems and involved in PDL1F19, PD-L1Q66, PD-L1A121, PD-L1D122 and PD-L1R125. As shown in Table 2, the hydrogen bonds are compared between PD-L1/PD-1(wild) and PD-L1/PD-1(mutant) systems and it can be seen that more hydrogen bonds are formed in PD-L1/PD-1(mutant) than PD-L1/PD-

Fig. 6. The polar energy and non-polar energy contributions of the key residues on PD-L1 in all the four complex systems.

RMSDs of the IgV domain of PD-L1 in Fig. 4b show few fluctuations during the simulation implying that the PD-L1 have strong conformational stability upon the binding of ligands in all the systems. The conformational stability of the IgV domain of human PD-L1 consisted with the studies by Ahmed et al. that limited structural differences between the PD-L1 conformations occurred in both the existing crystal structures and those populated during our MD simulations following a conformational selection mechanism [42]. However, the RMSDs for different protein ligands in Fig. 4c show that there are obvious flexibility difference among different protein ligands. So we further analyzed the root mean square fluctuations (RMSFs) of residues in different protein ligands. The comparison between PD-1(wild) and PD-1(mutant) in Fig. 4d shows that the main difference occur on the C′-D loop of PD-1 (residue 85 to 95). The root-mean-square fluctuations in Fig. 4e and f also show that the Fc domain seems to have a larger flexibility than Fv especially on the light chains. Comparing to avelumab, the initial structure of BMS-936559 contain an extra domain of Fc. The influence of Fc probably account for the reason why BMS-936559 has a relatively larger RMSD than BMS-936559 in Fig. 4c. According to the analysis above, all the 4 systems reach final equilibrium along the 200 ns trajectory. 3.3. The calculated binding free energy between PD-L1 and its protein ligands The binding free energies between PD-1 and its protein ligands in four complex systems were calculated by MM-GBSA method and listed in Table 1. It's found that the total Gibbs free energies have a good consistency between the calculated values (ΔGcal) and the experimental values (ΔGexp). According to the values of different energy terms, enthalpy contribution (ΔH) are the main driving force for the formation of protein-protein interaction. The enthalpy contribution (ΔH) were then split into non-polar and polar sections. For all the 4 systems, the nonpolar interaction contribution (ΔEnonpl,total) was favorable for ligand binding while the polar interaction contribution (ΔEpolar,total) was adverse mainly owing to the effect of solvent (ΔEpolar,sol). The non-polar and polar interaction contribution in the gas phase were further evaluated by the Van der Waals interaction (ΔEvdw,gas) and electrostatic interaction (ΔEele,gas). It can be seen that the electrostatic interaction (ΔEele,gas) was stronger in PD-L1/PD-1(mutant) system than PD-L1/PD582

BBA - General Subjects 1862 (2018) 576–588

D. Shi et al.

interactions. These results are consistent with the structural and dynamics studies by Roberta et al. showing that these mutations act in a synergistic manner to enhance the polar interactions between PD-L1 and PD-1(mutant). The recognized key residues on PD-L1 in all the four systems including PD-L1F19, PD-L1I54, PD-L1Y56, PD-L1D61, PD-L1N63, PD-L1Q66, PDL1H69, PD-L1V76, PD-L1M115, PD-L1D122, PD-L1Y123, PD-L1R125(except PDL1A18, PD-L1A51, PD-L1G120, PD-L1A121) were further validated by virtual alanine scanning mutagenesis. Although the residue PD-L1A121 were ignored during the calculation, the obvious energy contributions of PD-L1A121 in PD-L1/PD-1(wild) and PD-L1/BMS-936559 systems still suggested the importance of this residue. Otherwise, the residue mutation results by Wang et al. [55] and crystal structure analysis by Lin et al. [20] also confirmed the importance of residues 121 to 124 in murine PD-L1. The residue energy contribution in Fig. 5 showed that K124 didn't show obvious energy contribution for all the complex systems and was omitted for virtual alanine scanning mutagenesis. As shown in Fig. 7, six residues on PD-L1 namely PD-L1Y56, PD-L1Q66, PDL1M115, PD-L1D122, PD-L1Y123, PD-L1R125 turn out to be the hotspot residues with an energy contribution higher than 2 Kcal/mol. Among the recognized hotspot residues, PD-L1Y56, PD-L1M115, PD-L1D122, PDL1Y123 and PD-L1R125 have been identified as pharmacophores in the human PD-L1/human PD-1(wild) complex crystal structures [19]. Comparing to the reported combination of anchors(PD-L1Y56, PD-L1R113, PD-L1A121, PD-L1D122 and PD-L1Y123) recognized by Li et al. [26], all the other residues except for PD-L1R113 have also been validated by the virtual alanine scanning mutagenesis. Although PD-L1R113 appears as interfacial residue in the crystal structures of PD-L1/PD-1(wild), PD-L1/ PD-1(wild) and PD-1/avelumab systems, the results of per-residue energy decomposition show that PD-L1R113 doesn't have obvious energy contributions for all the four systems. As for PD-L1Q66, it appears as interfacial residues in all the four complex systems (Table S1), shows obvious energy contributions in the PD-L1/PD-1(mutant) system. Recently, the structural basis study of KN035, a novel PD-L1 nanobody, shows that PD-L1Q66 forms three hydrogen bonds upon KN035 binding and the replacement of PD-L1Q66 with alanine decreases the KN035 binding by 162-fold [56]. These results suggest that PD-L1Q66 have significant potential to be hotspot residues. So the identified hotspot residues can be used as potential binding site for the further small molecule or peptide inhibitor design.

Table 2 The hydrogen bonds at the interfaces of PD-L1/PD-1(wild), PD-L1/PD-1(mutant) systems. Acceptor PD-L1/PD-1(wild) PD-1D85@OD1 PD-1P130@O PD-1S87@HG PD-1S87@HB3 PD-1E136@OE2 PD-L1D122@OD1 PD-L1A121@O PD-L1/PD1(mutant) PD-1P130@O PD-1E70@OE1 PD-1E136@OE2 PD-1E136@OE1 PD-1 K131@HA PD-1E70@OE1 PD-1E136@OE1 PD-1E136@OE2 PD-L1D122@OD1 PD-L1Q66@OE1 PD-L1D122@OD2 PD-L1D122@CG PD-L1D122@OD1 PD-L1Q66@HE21 PD-L1A121@O

Donor

Distancea (Å)

Anglea (°)

Occupancyb

PD-L1F19@N

2.87 3.02 3.29 3.30 2.94 2.66 2.87

156.96 151.63 140.82 156.85 150.09 166.11 156.19

99.90% 91.71% 68.01% 56.24% 50.23% 99.95% 99.29%

2.92 2.84 3.01 3.06 3.14 2.96 2.87 2.94 2.86 3.03 2.80 3.29 3.09 3.25 3.05

156.72 152.70 146.11 144.91 130.93 145.05 153.04 147.56 160.33 163.73 152.63 157.87 141.27 142.48 138.99

98.71% 87.49% 81.51% 76.32% 71.17% 70.26% 69.42% 54.23% 99.88% 97.05% 88.16% 78.38% 72.43% 61.88% 55.35%

PD-L1Q66@NE2 PD-L1F19@CB PD-L1F19@CB PD-L1R125@NH1 PD-1Y68@OH PD-1N66@ND2

PD-L1Q66@NE2 PD-L1R25@NH1 PD-L1R125@NE PD-L1R125@NE PD-L1Q66@NE2 PD-L1R125@NH2 PD-L1R125@NH1 PD-L1R125@NH1 PD-1H68@NE2 PD-1I132@N

T78@OG1 T78@OG1 PD-1 T8@OG1 PD-1 K131@CA PD-1H64@NE2 PD-1 PD-1

a The hydrogen bonds are determined by an acceptor-donor atom distance of less than 3.5 Å and acceptor H-donor angle of greater than 120°. b To evaluate the stability and the strength of the hydrogen bond, only hydrogen bonds with occupancy more than 50% of the simulation time were listed.

1(wild) system. The hydrogen bonds between PD-L1F19 and PD-1D85, PD1S87 only exist in PD-L1/PD-1(wild) system owing to the conformational flexibility of C′-D loop in wild-type PD-1. In addition, the hydrogen bonding between PD-L1A121 and PD-1N66 in PD-L1/PD-1(wild) system is regenerated between PD-L1A121 and PD-1H64 owing to the residue mutations of V64H and N66V in mutant PD-1. Furthermore, the hydrogen bonds associated with PD-L1Q66, PD-L1D122 and PD-L1R125 in PD-L1/PD-1(wild) system are significantly enhanced by forming multiple hydrogen bonds with the mutant residues including PD-1Y68H, PD1M70E, PD-1K78T, PD-1A132I in PD-L1/PD-1(mutant) system, suggesting that these mutations are the main cause of the enhancement of polar

Fig. 7. The virtual alanine scanning mutagenesis of the key residues on PD-L1. The residues with the energy contribution higher than 2 Kcal/mol are picked out as hotspot residues. The electrostatic surface of hotspot residues are labeled out on PD-L1.

583

BBA - General Subjects 1862 (2018) 576–588

D. Shi et al.

Fig. 8. The energy cumulative curves for (a) PD-1(wild) and PD-1(mutant), (b) the heavy chains of monoclonal antibodies and (c) the light chains of monoclonal antibodies. The C-C′ domain, C′-D loop and F-G domain refer to the residue sequences 62–84, 85–95, 120–140 respectively. The CDR domains of monoclonal antibodies are highlighted by grey boxes as follows: HCDR1 (residue sequence: 25–34), HCDR2 (residue sequence: 51–61), HCDR3 (residue sequence: 98–108), LCDR1 (residue sequence: 25–34), LCDR2 (residue sequence: 49–63), LCDR3 (residue sequence: 90–100). The corresponding energy values (Kcal/mol) are labeled out at the position 40,80 and 118 in the heavy chains and at the position 40,80 and 109 in the light chains.

the main contacts with the F-G domain on both wild-type PD-1 and mutant PD-1. The interacting residue pairs contain PD-L1I54-PD-1L128, PD-L1Y56-PD-1A132, PD-L1Y56-PD-1I134, PD-L1Q66-PD-1A132, PDL1M115-PD-1I126, PD-L1M115-PD-1I134, PD-L1Y123-PD-1I126, PDL1Y123-PD-1I134, PD-L1Y123-PD-1E136 in both systems. According to the contact analysis in Table S2, except for the significant polar energy contribution of PD-L1Y123-PD-1E136, the other residue pairs mainly provide non-polar energy contributions. Otherwise, the contact of the residue pairs involved in residue mutations including PD-L1D122-PD1V66, PD-L1Y123-PD-1V66, PD-L1R125-PD-1E70 in PD-L1/PD-1(wild) system, are apparently enhanced in PD-L1/PD-1(mutant) system implying that N66 V, M70E in PD-1 (mutant) significantly affect the contact with the residues PD-L1D122, PD-L1Y123 and PD-L1R125. As for the monoclonal antibodies, the heavy chain and light chain of Fv domain both contain three complementary determining regions (CDRs), namely HCDR1, HCDR2, HCDR3, LCDR1, LCDR2, LCDR3. As shown in Fig. 8b and c, the interfacial interaction energy contributions

3.5. The residue-residue contacts between PD-L1 and its protein ligands In order to figure out the structural and energy basis for protein ligands, we applied the residue-residue contact analysis to recognize the key residue pairs between PD-L1 and its protein ligands. From the analysis above, it can be seen that the interfacial residues of PD-L1 turn out to be partly conservative when binding with different protein ligands. The conformational superposition and cumulative energy analysis of PD-1(wild) and PD-1(mutant) (Fig. 8a) show that PD-L1 interact with PD-1(wild) mainly on C′-D loop and F-G domain, while interacting with PD-1(mutant) mainly on C-C′ domain and F-G domain. The interactions with F-G domain tend to be conversed for PD-1(wild) and PD1(mutant) with energies of −18.03 and − 15.87 Kcal/mol respectively. However, the contacts of the C′-D loop and C-C′ domain are quite different between PD-L1/PD-1(wild) and PD-L1/PD-1(mutant) systems owing to the conformational change of these domains. The distribution of the residue-residue contact maps in Fig. 9 also show that PD-L1 has

584

BBA - General Subjects 1862 (2018) 576–588

D. Shi et al.

Fig. 9. The residue-residue contact maps between the residues on PD-L1 and the residues on (a) PD-1(wild) and (c) PD-1(mutant). The occurrence frequency of the contact between two residues are calculated throughout the last 20 ns trajectories. The key residue pairs in (b)PD-L1/PD-1(wild) and (d)PD-L1/PD-1(mutant) systems are shown in sticks with residues on PDL1 colored gray and residues on PD-1(wild) or PD-1(mutant) colored green.

of the heavy chain and the light chain are − 32.05 and −11.04 Kcal/ mol in avelumab, and − 16.94 and −3.15 Kcal/mol in BMS-936559 respectively. For both monoclonal antibodies, the interfacial interaction of the heavy chain are much stronger than the light chain and the HCDR2 and HCDR3 domains turn to dominate the interactions with PDL1. The further structural superposition of monoclonal antibodies also shows that the HCDR2 domain in both monoclonal antibodies seem to interact with the same domain on PD-L1. The residue-residue contact analysis in Fig. 10 shows that avelumab and BMS-936559 interact with PD-L1 mainly by their HCDR2 and HCDR3 domains, which is consistent with results of the interfacial interaction energy contributions (Fig. 8b). In PD-L1/avelumab system, the residues PD-L1D61, PD-L1N63 and PDL1V76 mainly interact with residues HL101 to HT106 on HCDR3, while the other four residues PD-L1I54, PD-L1Y56, PD-L1M115 and PD-L1A121 form close contact with residues HY52 to HF59 on HCDR2, which is consistent with the results from the crystal structure analysis [29]. PDL1D61 is the most interactive residues by forming residue-residue contacts with HT103, HV104, HT105, HT106, LY34, LY93, LR99 and the energy contributions of PD-L1D61-HT103, PD-L1D61-HV104, PDL1D61-HT105, PD-L1D61-HT106, PD-L1D61-LR99 are mainly polar energy contributions (Table S2). However, these polar interactions between residue pairs are largely compensated by the polar solvation effect, which explains the reason why PD-L1D61 has the most non-polar

contribution but no polar contribution for avelumab binding (Fig. 6b). As for PD-L1/BMS-936559 system, the residues PD-L1A51, PD-L1I54, PDL1H69 on PD-L1 form the main contacts with HCDR3 in BMS-936559, while PD-L1Y56, PD-L1M115, PD-L1G120, PD-L1A121, PD-L1Y123 forms the main contacts with HCDR2. PD-L1H69 on PD-L1 has the most non-polar contributions by forming contact with LW94, HP107, HS106 in BMS936559 (Table S2 and Fig. 10c). A comparison of the residue pairs between two monoclonal antibody systems reveals that the HCDR2 on both antibodies form contacts with PD-L1Y56, PD-L1M115 and PD-L1A121 on PD-L1. These results show a good agreement with the results from the crystal structures [29,30] and will provide important guidance for structural modification of monoclonal antibodies.

3.6. The potential drug design strategy targeting on PD-L1 Based on the analysis above, two β-hairpin peptides (124–136) with 13 residues were extracted from the F-G domain of PD-1(wild) and PD1(mutant), and two β-hairpin peptides (48–61) with 14 residues were extracted from the HCDR2 domain of avelumab and BMS-936559 as shown in Table 3. The residue energy contributions and the corresponding binding mode are shown in Fig. 11, suggesting that all the selected β-hairpin peptides have close contact with the putative hotspot domains (PD-L1Y56, PD-L1Q66, PD-L1M115, PD-L1D122, PD-L1Y123, PD-L1

585

BBA - General Subjects 1862 (2018) 576–588

D. Shi et al.

Fig. 10. The residue-residue contact maps of the residues on PD-L1 and the residues on (a) avelumab and (c) BMS-936559. The occurrence frequency of the contact between two residues are calculated throughout the last 20 ns trajectories. The key residue pairs in (b) PD-L1/avelumab and (d) PD-L1/BMS-936559 systems are shown in sticks with residues on PD-L1 colored gray, residues on the heavy chains colored green and residues on the light chains colored cyan.

interaction network analysis and binding free energy decomposition and then the hotspot residues were validated by virtual alanine scanning mutagenesis. According to our computational results, the putative hotspot residues on PD-L1 include PD-L1Y56, PD-L1Q66, PD-L1M115, PDL1Asp122, PD-L1Tyr123, PD-L1R125, which are mainly located on the CC'FG strands of PD-L1. The following cumulative contributions analysis and residue-residue contact analysis show that wild-type PD-1 and the mutant PD-1 mainly interact with PD-L1 by their F-G domains while monoclonal antibodies like avelumab and BMS-936559 mainly interact with PD-L1 by their HCDR2 and HCDR3 domains. A structurally similar β-hairpin domain with 13 or 14 residues was extracted from each ligand protein and turned out to bind with the hotspot residues in all the four complex systems. These peptide sequences can be used as a good starting point for rational design of cyclic peptides or peptide mimics targeting PD-L1.

Table 3 The residue sequences and energy contributions for the extracted β-hairpin peptides. Protein ligand

Domain

Number

Sequence

Contribution (Kcal/mol)

PD-1(wild) PD-1(mutant) Avelumab BMS-936559

F-G F-G HCDR2 HCDR2

124–136(1–13) 124–136(1–13) 48–61(1–14) 48–61(1–14)

GAISLAPKAQIKE GVISLAPKIQIKE VSSIYPSGGITFYA MGGIIPIFGKAHYA

− 18.04 − 18.02 − 11.38 − 9.43

R125) despite their different residue sequences or conformational states. And these β-hairpin peptides can be set as good starting point for peptide inhibitor design targeting on PD-L1. 4. Conclusions

Transparency document In this work, we analyzed the interfacial interactions between PD-L1 and four different protein ligands and recognized the hotspot residues on PD-L1 through molecular dynamics simulations. The key interfacial residues on PD-L1 surface are were firstly recognized by residue

The http://dx.doi.org/10.1016/j.bbagen.2017.11.022 associated with this article can be found, in online version.

586

BBA - General Subjects 1862 (2018) 576–588

D. Shi et al.

Fig. 11. (a) The per-residue energy contributions of the extracted β-hairpin peptides and the binding modes of the β-hairpin peptides extracted from (b) PD-1(wild), (c) PD-1(mutant), (d) avelumab and (e) BMS-936559. The residues on peptides are shown in cyan sticks while the hotspot residues on PD-L1 are shown in red surface.

Acknowledgements

[10] G. Abril-Rodriguez, A. Ribas, SnapShot: immune checkpoint inhibitors, Cancer Cell 31 (2017) 848. [11] D.B. Page, M.A. Postow, M.K. Callahan, J.P. Allison, J.D. Wolchok, Immune modulation in cancer with antibodies, Annu. Rev. Med. 65 (2014) 185–202. [12] P. Sharma, J.P. Allison, The future of immune checkpoint therapy, Science 348 (2015) 56–61. [13] S.L. Topalian, D.M. Weiner, G.J. Fau, D.M. Pardoll, Cancer immunotherapy comes of age, J. Clin. Oncol. 29 (2011) 4828–4836. [14] I. Mellman, G. Coukos, G. Fau, G. Dranoff, Cancer immunotherapy comes of age, Nature 480 (2011) 480–489. [15] M.M. Zhan, X.Q. Hu, X.X. Liu, B.F. Ruan, J. Xu, C. Liao, From monoclonal antibodies to small molecules: the development of inhibitors targeting the PD-1/PD-L1 pathway, Drug Discov. Today 21 (2016) 1027–1036. [16] T. Zarganes-Tzitzikas, M. Konstantinidou, Y. Gao, D. Krzemien, K. Zak, G. Dubin, T.A. Holak, A. Domling, Inhibitors of programmed cell death 1 (PD-1): a patent review (2010–2015), Expert Opin. Ther. Pat. 26 (2016) 973–977. [17] H. Weinmann, Cancer immunotherapy: selected targets and small-molecule modulators, ChemMedChem 11 (2016) 450–466. [18] J.L. Adams, J. Smothers, R. Srinivasan, A. Hoos, Big opportunities for small molecules in immuno-oncology, Nat. Rev. Drug Discov. 14 (2015) 603–622. [19] K.M. Zak, R. Kitel, S. Przetocka, P. Golik, K. Guzik, B. Musielak, A. Domling, G. Dubin, T.A. Holak, Structure of the complex of human programmed death 1, PD1, and its ligand PD-L1, Structure 23 (2015) 2341–2348. [20] D.Y. Lin, Y. Tanaka, M. Iwasaki, A.G. Gittis, H.P. Su, B. Mikami, T. Okazaki, T. Honjo, N. Minato, D.N. Garboczi, The PD-1/PD-L1 complex resembles the antigen-binding Fv domains of antibodies and T cell receptors, Proc. Natl. Acad. Sci. U. S. A. 105 (2008) 3011–3016. [21] M.R. Arkin, Y. Tang, J.A. Wells, Small-molecule inhibitors of protein-protein interactions: progressing toward the reality, Chem. Biol. 21 (2014) 1102–1114. [22] M.R. Arkin, J.A. Wells, Small-molecule inhibitors of protein-protein interactions: progressing towards the dream, Nat. Rev. Drug Discov. 3 (2004) 301–317. [23] S.P. Patel, R. Kurzrock, PD-L1 expression as a predictive biomarker in cancer immunotherapy, Mol. Cancer Ther. 14 (2015) 847–856. [24] Z. Gatalica, C. Snyder, T. Maney, A. Ghazalpour, D.A. Holterman, N. Xiao, P. Overberg, I. Rose, G.D. Basu, S. Vranic, H.T. Lynch, D.D. Von Hoff, O. Hamid, Programmed cell death 1 (PD-1) and its ligand (PD-L1) in common cancers and their correlation with molecular cancer type, Cancer Epidem. Biomar. 23 (2014) 2965–2970. [25] K. Guzik, K.M. Zak, P. Grudnik, K. Magiera, B. Musielak, R. Törner, L. Skalniak,

This work was supported by the National Natural Science Foundation of China (Grant No. 21475054) and the Fundamental Research Funds for the Central Universities (Grant No. lzujbky-2014191). Appendix A. Supplementary data Supplementary data to this article can be found online at https:// doi.org/10.1016/j.bbagen.2017.11.022. References [1] V.A. Boussiotis, Molecular and biochemical aspects of the PD-1 checkpoint pathway, N. Engl. J. Med. 375 (2016) 1767–1778. [2] K.C. Ohaegbulam, A. Assal, E. Lazar-Molnar, Y. Yao, X. Zang, Human cancer immunotherapy with antibodies to the PD-1 and PD-L1 pathway, Trends Mol. Med. 21 (2015) 24–33. [3] C.G. Drake, E.J. Lipson, J.R. Brahmer, Breathing new life into immunotherapy: review of melanoma, lung and kidney cancer, Nat. Rev. Clin. Oncol. 11 (2014) 24–37. [4] M. Saresella, V. Rainone, N.M. Al-Daghri, M. Clerici, D. Trabattoni, The PD-1/PD-L1 pathway in human pathology, Curr. Mol. Med. 12 (2012) 259–267. [5] L. Chen, Co-inhibitory molecules of the B7-CD28 family in the control of T-cell immunity, Nat. Rev. Immunol. 4 (2004) 336–347. [6] P. Sharma, J.P. Allison, Immune checkpoint targeting in cancer therapy: toward combination strategies with curative potential, Cell 161 (2015) 205–214. [7] D.M. Pardoll, The blockade of immune checkpoints in cancer immunotherapy, Nat. Rev. Cancer 12 (2012) 252–264. [8] G.J. Freeman, Structures of PD-1 with its ligands: sideways and dancing cheek to cheek, Proc. Natl. Acad. Sci. U. S. A. 105 (2008) 10275–10276. [9] A. Hoos, Development of immuno-oncology drugs - from CTLA4 to PD1 to the next generations, Nat. Rev. Drug Discov. 15 (2016) 235–247.

587

BBA - General Subjects 1862 (2018) 576–588

D. Shi et al.

[26]

[27]

[28]

[29]

[30]

[31]

[32]

[33]

[34]

[35]

[36]

[37]

[38]

[39]

[40]

[41] C. Viricel, M. Ahmed, K. Barakat, Human PD-1 binds differently to its human ligands: a comprehensive modeling study, J. Mol. Graph. Model. 57 (2015) 131–142. [42] M. Ahmed, K. Barakat, The too many faces of PD-L1: a comprehensive conformational analysis study, Biochemistry 56 (2017) 5428–5439. [43] W. Xue, Y. Ban, H. Liu, X. Yao, Computational study on the drug resistance mechanism against HCV NS3/4A protease inhibitors vaniprevir and MK-5172 by the combination use of molecular dynamics simulation, residue interaction network, and substrate envelope analysis, J. Chem. Inf. Model. 54 (2014) 621–633. [44] W. Xue, X. Jin, L. Ning, M. Wang, H. Liu, X. Yao, Exploring the molecular mechanism of cross-resistance to HIV-1 integrase strand transfer inhibitors by molecular dynamics simulation and residue interaction network analysis, J. Chem. Inf. Model. 53 (2013) 210–222. [45] G. Hultqvist, S.R. Haq, A.S. Punekar, C.N. Chi, A. Engstrom, A. Bach, K. Stromgaard, M. Selmer, S. Gianni, P. Jemth, Energetic pathway sampling in a protein interaction domain, Structure 21 (2013) 1193–1202. [46] L. Xie, J. Li, L. Xie, P.E. Bourne, Drug discovery using chemical systems biology: identification of the protein-ligand binding network to explain the side effects of CETP inhibitors, PLoS Comput. Biol. 5 (2009) e1000387. [47] Z. Hu, D. Bowen, W.M. Southerland, A. del Sol, Y. Pan, R. Nussinov, B. Ma, Ligand binding and circular permutation modify residue interaction network in DHFR, PLoS Comput. Biol. 3 (2007) e117. [48] D. Piovesan, G. Minervini, S.C. Tosatto, The RING 2.0 web server for high quality residue interaction networks, Nucleic Acids Res. 44 (2016) W367–374. [49] P. Shannon, A. Markiel, O. Ozier, N.S. Baliga, J.T. Wang, D. Ramage, N. Amin, B. Schwikowski, T. Ideker, Cytoscape: a software environment for integrated models of biomolecular interaction networks, Genome Res. 13 (2003) 2498–2504. [50] D. Case, V. Babin, J. Berryman, R. Betz, Q. Cai, D. Cerutti, T. Cheatham Iii, T. Darden, R. Duke, H. Gohlke, A. Goetz, S. Gusarov, N. Homeyer, P. Janowski, J. Kaus, I. Kolossváry, A. Kovalenko, T. Lee, S. LeGrand, T. Luchko, R. Luo, B. Madej, K. Merz, F. Paesani, D. Roe, A. Roitberg, C.R. Sagui Salomon-Ferrer, G. Seabra, C. Simmerling, W. Smith, J. Swails, R. Walker, J. Wang, R. Wolf, X. Wu, P. Kollman, AMBER14, University of California, San Francisco, 2014. [51] T. Kortemme, D.E. Kim, D. Baker, Computational alanine scanning of protein-protein interfaces, Sci. Signal. 2004 (2004) pl2. [52] L. Ning, J. Guo, Q. Bai, N. Jin, H. Liu, X. Yao, Structural diversity and initial oligomerization of PrP106-126 studied by replica-exchange and conventional m3lecular dynamics simulations, PLoS One 9 (2014) e87266. [53] U. Doshi, M.J. Holliday, E.Z. Eisenmesser, D. Hamelberg, Dynamical network of residue–residue contacts reveals coupled allosteric effects in recognition, catalysis, and mutation, Proc. Natl. Acad. Sci. U. S. A. 113 (2016) 4735–4740. [54] K.M. Zak, P. Grudnik, K. Guzik, B.J. Zieba, B. Musielak, A. Dömling, G. Dubin, T.A. Holak, Structural Basis for Small Molecule Targeting of the Programmed Death Ligand 1 (PD-L1), Oncotarget, Vol. 7 2016, pp. 30323–30335. [55] S. Wang, J. Bajorath, D.B. Flies, H. Dong, T. Honjo, L. Chen, Molecular modeling and functional mapping of B7-H1 and B7-DC uncouple costimulatory function from PD-1 interaction, J. Exp. Med. 197 (2003) 1083–1091. [56] F. Zhang, H. Wei, X. Wang, Y. Bai, P. Wang, J. Wu, X. Jiang, Y. Wang, H. Cai, T. Xu, A. Zhou, Structural basis of a novel PD-L1 nanobody for immune checkpoint blockade, Cell Discovery 3 (2017) 17004.

A. Dömling, G. Dubin, T.A. Holak, Small-molecule inhibitors of the programmed cell Death-1/programmed death-ligand 1 (PD-1/PD-L1) interaction via transientlyinduced protein states and dimerization of PD-L1, J. Med. Chem. 60 (2017) 5857–5867. Q. Li, L. Quan, J. Lyu, Z. He, X. Wang, J. Meng, Z. Zhao, L. Zhu, X. Liu, H. Li, Discovery of peptide inhibitors targeting human programmed death 1 (PD-1) receptor, Oncotarget 40 (2016) 64967–64976. A.F. Abdel-Magid, Inhibitors of the PD-1/PD-L1 pathway can mobilize the immune system: an innovative potential therapy for cancer and chronic infections, ACS Med. Chem. Lett. 6 (2015) 489–490. R. Pascolutti, X. Sun, J. Kao, R.L. Maute, A.M. Ring, G.R. Bowman, A.C. Kruse, Structure and dynamics of PD-L1 and an ultra-high-affinity PD-1 receptor mutant, Structure 24 (2016) 1719–1728. K. Liu, S. Tan, Y. Chai, D. Chen, H. Song, C.W. Zhang, Y. Shi, J. Liu, W. Tan, J. Lyu, S. Gao, J. Yan, J. Qi, G.F. Gao, Structural basis of anti-PD-L1 monoclonal antibody avelumab for tumor therapy, Cell Res. 27 (2017) 151–153. J.Y. Lee, H.T. Lee, W. Shin, J. Chae, J. Choi, S.H. Kim, H. Lim, T. Won Heo, K.Y. Park, Y.J. Lee, S.E. Ryu, J.Y. Son, J.U. Lee, Y.S. Heo, Structural basis of checkpoint blockade by monoclonal antibodies in cancer immunotherapy, Nat. Commun. 7 (2016) 13354. I. Massova, P.A. Kollman, Combined molecular mechanical and continuum solvent approach (MM-PBSA/GBSA) to predict ligand binding, Perspect. Drug Discov. 18 (2000) 113–135. H.-Y. Sun, F.-Q. Ji, A molecular dynamics investigation on the crizotinib resistance mechanism of C1156Y mutation in ALK, Biochem. Biophys. Res. Commun. 423 (2012) 319–324. H. Liu, X. Yao, C. Wang, J. Han, In silico identification of the potential drug resistance sites over 2009 influenza A (H1N1) virus neuraminidase, Mol. Pharm. 7 (2010) 894–904. J. Zhang, T. Hou, W. Wang, J.S. Liu, Detecting and understanding combinatorial mutation patterns responsible for HIV drug resistance, Proc. Natl. Acad. Sci. U. S. A. 107 (2010) 1321–1326. J. Zhu, P. Pan, Y. Li, M. Wang, D. Li, B. Cao, X. Mao, T. Hou, Theoretical studies on beta and delta isoform-specific binding mechanisms of phosphoinositide 3-kinase inhibitors, Mol. BioSyst. 10 (2014) 454–466. Y. Yang, Y. Shen, S. Li, N. Jin, H. Liu, X. Yao, Molecular dynamics and free energy studies on aurora kinase A and its mutant bound with MLN8054: insight into molecular mechanism of subtype selectivity, Mol. BioSyst. 8 (2012) 3049–3060. Y. Yang, H. Liu, X. Yao, Understanding the molecular basis of MK2-p38alpha signaling complex assembly: insights into protein-protein interaction by molecular dynamics and free energy studies, Mol. BioSyst. 8 (2012) 2106–2118. J. Guo, X. Wang, H. Sun, H. Liu, X. Yao, The molecular basis of IGF-II/IGF2R recognition: a combined molecular dynamics simulation, free-energy calculation and computational alanine scanning study, J. Mol. Model. 18 (2012) 1421–1430. H. Liu, X. Yao, Molecular basis of the interaction for an essential subunit PA-PB1 in influenza virus RNA polymerase: insights from molecular dynamics simulation and free energy calculation, Mol. Pharm. 7 (2010) 75–85. W.P. Liu, B. Huang, Y.S. Kuang, G.J. Liu, Molecular dynamics simulations elucidate conformational selection and induced fit mechanisms in the binding of PD-1 and PD-L1, Mol. BioSyst. 13 (2017) 892–900.

588