Investigating the Role of the N-Terminal Loop of PD-1 in Binding Process Between PD-1 and Nivolumab via Molecular Dynamics Simulation

The blockade of immune checkpoints, such as programmed death receptor 1 (PD-1) and programmed death ligand 1 protein (PD-L1), is a promising therapeutic approach in cancer immunotherapy. Nivolumab, a humanized IgG4 antibody targeting PD-1, was approved by the US Food and Drug Administration for several cancers in 2014. Crystal structures of the nivolumab/PD-1 complex show that the epitope of PD-1 locates at the IgV domain (including the FG and BC loops) and the N-terminal loop. Although the N-terminal loop of PD-1 has been shown to play a dominant role in the complex interface of the static structure, its role in the dynamic binding process has not been illustrated clearly. Here, eight molecular systems were established for nivolumab/PD-1 complex, and long-time molecular dynamics simulations were performed for each. Results showed that the N-terminal loop of PD-1 prefers to bind with nivolumab to stabilize the interface between IgV and nivolumab. Furthermore, the binding of the N-terminal loop with nivolumab induces the rebinding between the IgV domain and nivolumab. Thus, we proposed a two-step binding model for the nivolumab/PD-1 binding, where the interface switches to a high-affinity state with the help of the N-terminal loop. This finding suggests that the N-terminal loop of PD-1 might be a potential target for anti-PD-1 antibody design, which could serve as an important gatekeeper for the anti-PD-1 antibody binding.


INTRODUCTION
Cancer, the leading cause of death worldwide, constitutes a considerable burden to society (Bray et al., 2018;Miller et al., 2019). Cancer immunotherapy, that is, harnessing the immune system to battle tumors, has attained remarkable achievement in cancer treatment. Cancer immunotherapy comprises various treatment approaches, including antitumor monoclonal antibodies, cancer vaccines, and antibodies that block immune inhibitory pathways (Mellman et al., 2011;Couzin-Frankel, 2013). Among these treatments, blockade of immune checkpoints is the most promising approach to activate therapeutic antitumor immunity (Ribas and Wolchok, 2018). The 2018 Nobel Prize in Physiology or Medicine was awarded to James P. Allison and Tasuku Honjo for their pioneering discoveries that led to the development of immune checkpoint inhibitors, which block the inhibitory action of T cell molecules, including programmed death receptor-1 (PD-1) and cytotoxic T lymphocyte-associated protein 4 (CTLA-4) (Pardoll, 2012;Guo, 2018;Zaidi and Jaffee, 2018).
PD-1 is an immune checkpoint receptor of the CD28 family expressed in tumor and immune cells (Keir et al., 2008). It is a 288 amino acid type I transmembrane receptor, and its ectodomain consists of four domains, namely, signal peptide, N-terminal loop, extracellular immunoglobulin variable (IgV) domain, and stalk region. The blockade of the interaction between PD-1 and its ligand programmed death ligand 1 protein (PD-L1) was observed to restore the attenuated immune response and lead to increased antitumor and antiviral activities (Hirano et al., 2005;Zitvogel and Kroemer, 2012;Sharma and Allison, 2015). Several PD-1/PD-L1 pathway inhibitors, including pembrolizumab, nivolumab targeting PD-1 and atezolizumab, durvalumab, and avelumab targeting PD-L1, have been approved by the US Food and Drug Administration (FDA) for the treatment of multiple cancers to date (Callahan et al., 2016;Ivashko and Kolesar, 2016;Leventakos and Mansfield, 2016;Bellmunt et al., 2017;Kim, 2017;Muller et al., 2017;Rittmeyer et al., 2017;Sidaway, 2017;Syed, 2017). Humanized IgG4 antibody nivolumab received the most attention and was approved for the treatment of melanoma, metastatic non-small-cell lung cancer, renal cell carcinoma, and Hodgkin lymphoma in 2014.
Two crystal structures of the nivolumab/PD-1complex have been reported in 2015 and 2016, respectively, providing interface information at the atomic level (Figures 1A,E; Lee et al., 2016;Tan et al., 2017). Three loops of PD-1 provide a flexible platform for nivolumab binding, including the N-terminal loop, the FG loop, and the BC loop. The FG and BC loops locate on the IgV domain. Previous studies showed that the FG loop (i.e., PD−1 PRO130 and PD−1 LYS131) is a binding site for PD-L1 as well as a "hot spot" for several immune checkpoint blockade monoclonal antibodies, such as GY-5 and GY-14 (Zak et al., 2015;Chen et al., 2019a). They clearly suggest a steric clash blockade mechanism of nivolumab. Unexpectedly, the N-terminal loop of PD-1 is far from the interface of PD-1/PD-L1 complex but contributes the majority of hydrogen bonds (H-bonds) to the binding of nivolumab and PD-1. Experiments further proved that the truncation of the N-terminal loop of PD-1 would abolish the nivolumab binding (Tan et al., 2017). Thus, the N-terminal loop of PD-1 plays a dominant role in the complexation of PD-1 and nivolumab. In spite of this, how the N-terminal loop regulates the dynamic binding process has not been answered clearly. Molecular recognition is a dynamic process, and the binding of a ligand to its receptor is regarded not as a single, frozen structure but rather a macromolecule in constant motion (Moroni et al., 2015). Considering the two-site mode (the N-terminal loop and the IgV domain) at the interface of the nivolumab/PD-1 complex, we assumed a two-step binding model, where the N-terminal loop will help switch the binding to a stronger state whether it comes across nivolumab firstly or not.
Molecular dynamics (MD) simulation is well suited for studying the dynamics of proteins (Liu et al., 2018). We have used MD simulations to elucidate conformational selection and induced fit mechanisms in the binding of PD-1 and PD-L1 via MD simulations, and found that the CC' loop of PD-1 is flexible and switches from an open form to a close one to stabilize the PD-1/PD-L1 complex . MD simulation is also proven to be a useful tool to detect the "hot-spots" in the complex interface, such as the 6B4/GPIba complex, the 10B12/GPVI complex, the PD-1/PD-L1 complex, and the PD-1/pembrolizumab complex (Fang et al., 2012;Liu et al., 2016;. Therefore, MD simulations were adopted here to investigate the role of the N-terminal loop of PD-1 in the dynamic binding process between PD-1 and nivolumab. Two crystal structures of the nivolumab/PD-1 complex were used to build eight molecular systems with different binding states to mimic the scenarios with or without the N-terminal loop, and the N-terminal loop binds firstly or not (Figure 1). The results show that the N-terminal loop of PD-1 prefers to bind with nivolumab to stabilize the complex interface between the IgV domain (i.e., FG loop and BC loop) and nivolumab. The binding of the N-terminal loop with nivolumab also induces the rebinding between the IgV domain and nivolumab. These findings suggest a two-step binding model, in which the interface of nivolumab/PD-1complex switches to a stronger binding state with the help of the N-terminal loop of PD-1.

System Setup
Eight nivolumab/PD-1 complex structures were set up for MD simulations (Figure 1). First, two crystal structures of the nivolumab/PD-1 complex with accession codes of 5GGR and 5WT9 were downloaded from the PDB, and were designated as Complex I and Complex II respectively (Figures 1A,E). Both structures include a long N-terminal loop at the complex interface, but in different lengths ( PD−1 SER27-PD−1 ASN33 for complex I, PD−1 LEU25-PD−1 ASN33 for complex II). Actually, Complex II has an intact N-terminal loop and Complex I only lacks two residues, because the residues before PD−1 LEU25 belong to the signal peptide, which will be post-translationally removed and cannot be secreted. Second, the N-terminal loops of Complex I and II were cut off, and the remaining structures were designated as Complex I-N-truncated and Complex II-Ntruncated (Figures 1B,F). Third, the N-terminal loop of Complex I and Complex II was rotated backward at the interface with 90 • to dissociate from nivolumab, with the IgV domain of PD-1 and nivolumab fixed. These two structures were used to mimic the scenario where the IgV domain of PD-1 binds to nivolumab at the first step and designated as Complex I-N-rotated and Complex II-N-rotated (Figures 1C,G). Finally, the IgV domain of Complex I ( PD−1 PRO34-PD−1 ARG143) and Complex II ( PD−1 PRO34-PD−1 LEU142) was rotated backward at the interface with 90 • to dissociate from nivolumab, with the N-terminal loop and nivolumab fixed. These two structures were used to mimic the scenario where the N-terminal loop of PD-1 binds to nivolumab at the first step and designated as Complex I-IgV-rotated and Complex II-IgV-rotated (Figures 1D,H).
The N-and C-terminal residues of each complex were acetylated and amidated, respectively, to mimic the continuation of the protein chains. The missing residues of PD-1 in each complex structure were modeled by the SWISS-MODEL server, with the free PD-1 structures with PDB accession codes 3RRQ and 2M2D as templates (Bordoli et al., 2008;Biasini et al., 2014). The protonation state of each protein residue at neutral pH was determined with the software PROPKA (Bas et al., 2008). Each complex was first solvated with TIP3P water molecules in a rectangular box with walls at least 15Å away from any protein atom. Then, Na + and Cl − ions were added to neutralize the systems at a 150 mM salt concentration. The final system each consists of ∼35,500 water molecules, ∼100 Na + and ∼100 Cl − ions (Supplementary Table S1).

MD Simulations
VMD 1.9.3 program was used for visualization, modeling, data analysis, and conformation presentation (Humphrey et al., 1996). NAMD 2.11 program with CHARMM36 all-atom force field was used for simulations (Phillips et al., 2005;Best et al., 2012). The cMAP correction for protein backbone, a time step of 2 fs, and a periodic boundary condition were applied in the simulations. The particle mesh Ewald method and a smooth cutoff of 12 Å were employed to calculate the full electrostatic and van der Waals forces. First, each system was energy-minimized for 5,000 steps with all protein atoms fixed and for another 5,000 steps with all atoms free. Next, each system was heated gradually from 0 to 310 K in 1 ns. Then, equilibrium simulation of 100 ns was performed thrice (named Equ1, Equ2, and Equ3) for each system. A 310 K heat bath was manipulated using the Langevin thermostat, and a 1 atm pressure was controlled using the Langevin piston method during equilibriums.

Data Analysis
Transient complex formation usually relies on H-bonds (Kar et al., 2012), and they are the dominant linkers at the interface of nivolumab/PD-1 complex (Lee et al., 2016;Tan et al., 2017). Therefore, H-bonds across the interface in simulations were detected using VMD software with in-house scripts. An H-bond was defined if the donor-acceptor distance and bonding angle were smaller than 3.5 Å and 30 • , respectively. The survival ratio of an H-bond was defined as the percentage of bond survival time.
The root mean square deviation (RMSD) of heavy atoms was used to illustrate the stability of the structures as well as the conformational changes of the N-terminal loop and IgV domain of PD-1. When analyzing the RMSD for the N-terminal loop and IgV domain of PD-1, the structures of nivolumab were aligned. Buried solvent-accessible surface area (SASA) of each complex, representing the interface area, was calculated using VMD software with handwritten scripts. The interaction energy, mainly including van der Waals and electrostatic energy, was calculated using NAMD Energy plugin (version 1.4) provided in VMD 1.9.3.

Interface Analysis of the Nivolumab/PD-1 Complex
To describe the dynamic picture of the interface of the nivolumab/PD-1 complex, two crystal structures of the complex with accession codes 5GGR and 5WT9 were downloaded from the PDB database, and designated as Complex I and Complex II (Figures 1A,E), respectively. They were solvated with TIP3P water molecules in a rectangular box. Equilibrium simulation of 100 ns was performed thrice (named Equ1, Equ2 and Equ3) for each complex after an energy minimization of 10,000 steps. The RMSD of heavy atoms showed that these two complexes had reached a local minimum after 20 ns (Supplementary Figures S1, S2).
Variations of buried SASA, interaction energy, and number of H-bonds were recorded for each complex to evaluate the stability of the complex interface (Figures 2, 3 H-bonds with a maximum survival ratio of above 0.2 were sorted out to recognize the interaction residues because they are proposed as the dominant linkers across the interface (Lee et al., 2016;Tan et al., 2017; Table 1). In total 15 H-bonds were detected for Complex I, and 20 for Complex II. Interaction residues on PD-1 of these two complexes all located on three loops, including the N-terminal loop ( PD−1 ASP26, PD−1 PRO28, PD−1 ASP29, and PD−1 ARG30), the BC loop ( PD−1 THR59, PD−1 SER60, and PD−1 GLU61), and the FG loop ( PD−1 ALA129, PD−1 PRO130, and PD−1 LYS131) (Figures 2D-F, 3D-F). The FG loop formed five H-bonds in both the interfaces of Complex I and Complex II with the highest survival ratio of 0.68. The BC loop formed four and six H-bonds at the interface of Complex I and Complex II, respectively, with the highest survival ratio of 0.53. Although the numbers of H-bonds were similar to those of the FG loop, the interface formed by the BC loop was more vulnerable because it only appeared in two of three runs for both Complex I and Complex II. The average survival ratios of the H-bonds involved with the FG loop were higher than those with the BC loop, while the standard deviations with the FG loop were lower than those with the BC loop ( Table 1). As the FG loop is a binding region for PD-L1, these results clearly suggest a stable steric clash blockade mechanism of nivolumab.
In the crystal structures, the N-terminal loop of PD-1 is not a binding site for PD-L1 but forms a major interface with nivolumab. In our simulations, it involved in six H-bonds at the interface of Complex I and contributed nine H-bonds to the interface of Complex II, with the highest survival ratio of 0.82. The standard deviations of the survival ratios of the H-bonds formed by the N-terminal loop were a little high. This might be due to the high flexibility of the long N-terminal loop, which consists of eight residues ( PD−1 SER27 to PD−1 PRO34) in Complex I and 10 residues ( PD−1 LEU25 to PD−1 PRO34) in Complex II, and can only be constrained   by one side. The PD−1 ASP26 located at the N-terminal loop formed two H-bonds with nivolumab in Complex II, but it was missing in Complex I. Thus, the N-terminal loop in Complex II formed a larger interface with nivolumab. These MD results showed that interactions between the N-terminus of PD-1 and nivolumab are definite and stable on the nanosecond time scale we simulated.

Truncation of the N-Terminal Loop of PD-1 Impairs the Interface Between PD-1 and Nivolumab
The N-terminal loop of PD-1 greatly contributes to the interface of the nivolumab/PD-1 complex, and mutagenesis study revealed that the cut-off of the N-terminal loop abolished the binding between PD-1 and nivolumab. How will the truncation of the N-terminal loop impair the interfaces? To answer this question, we cut off the N-terminal loops in Complex I and II, and designated them as Complex I-N-truncated and Complex II-Ntruncated (Figures 1B,F). These two systems were simulated with the same scenario as before. The RMSD of heavy atoms showed that these two complexes reached a local minimum after 20 ns (Supplementary Figures S3, S4).
Variations of buried SASA, interaction energy, and number of H-bonds along the simulation time were analyzed, as shown in Figures 5A-C, 6A-C, and their distributions are shown in Figure 4. The buried SASA of Complex I-N-truncated greatly decreased to around 900 Å 2 , and its interaction energy and number of H-bonds dropped to -180 kcal/mol and 3, respectively. The buried SASA of Complex II-N-truncated greatly decreased to around 1,000 Å 2 , and its interaction energy and number of H-bonds dropped to -180 kcal/mol and 3, respectively. These results indicate that the truncation of the N-terminal loop seriously impairs the binding between PD-1 and nivolumab (Figures 5D-F, 6D-F). Moreover, the binding strength of Complex I-N-truncated is similar to that of Complex II-N-truncated, implying that Complex II is more stable than Complex I because it has a longer N-terminal loop ( PD−1 LEU25 to PD−1 PRO34) than Complex I ( PD−1 SER27 to PD−1 PRO34). This conclusion was further proved by the H-bonds with a survival ratio of above 0.2, where nine H-bonds  for Complex I-N-truncated, and ten for Complex II-N-truncated were found with the interaction residues contributed by the FG and BC loops of PD-1, similar to those in Complex I and Complex II ( Table 2).
Although the interface of the nivolumab/PD-1 was seriously impaired by cutting off the N-terminal loop of PD-1, the dissociation was not observed. A possible reason is that the energy barrier involved in the FG and BC loops of PD-1 was too high to overcome within our simulation time. However, analysis of accessibility of water molecules around the FG and BC loops revealed one major difference before and after deletion of the N-terminal loop. As shown in Figures 7A,B, the difference lies in extra water accessibility near the BC loop, where about nine and ten more water molecules entered within 4 Å of the BC loops of Complex I-N-truncated and Complex II-N-truncated. For the FG loop, the water accessibility showed no significant change. This result suggests the interface involved in the BC loop is protected from water attack by the N-terminal loop. Deleting

The N-Terminal Loop of PD-1 Prefers to Interact With Nivolumab to Stabilize the Complex Interface Further
Two binding regions were mapped on PD-1 for nivolumab, namely, the N-terminal loop and the IgV domain (including the FG loop and the BC loop), implying the possibility of a twostep binding process. Therefore, four additional complexes were built to verify this hypothesis. First, the N-terminal loop of PD-1 in Complex I and Complex II was rotated backward against the interface to dissociate from nivolumab to mimic the scenario where the IgV domain of PD-1 binds to nivolumab at the first step, designated as Complex I-N-rotated and Complex II-Nrotated, respectively (Figures 1C,G). Second, the IgV domain of PD-1 in Complex I and Complex II was rotated backward against the interface to dissociate from nivolumab to mimic the scenario where the N-terminal loop of PD-1 binds to nivolumab at the first step, designated as Complex I-IgV-rotated and Complex II-IgV-rotated, respectively (Figures 1D,H).
Similarly, Complex I-N-rotated and Complex II-N-rotated were simulated for 100 ns thrice after an energy minimization of 10,000 steps. The RMSD of heavy atoms showed that these two complexes reached a local minimum after 20 ns (Supplementary  Figures S5, S6). Buried SASA, interaction energy, and number of H-bonds of Complex I-N-rotated and Complex II-N-rotated are shown in Figures 8A-C, 9A-C, respectively, and their distributions are demonstrated in Figure 4. For Complex I-Nrotated, its buried SASA decreased to around 800 Å 2 during Equ1 and Equ2, and its interaction energy and number of H-bonds decreased to around -180 kcal/mol and 3, respectively, which was close to the binding strength of Complex I-N-truncated. However, its buried SASA increased to around 1,400 Å 2 during Equ3 with interaction energy and number of H-bonds fluctuating around −280 kcal/mol and 7, respectively, which was close to the binding strength of the initial Complex I.
Were these changes induced by the N-terminal loop or the IgV domain of PD-1? To answer this question, the number of intermolecular H-bonds formed by the N-terminal loop and the IgV domain of PD-1 was calculated (Figures 8D,E). H-bonds with a survival ratio of above 0.2 are listed in Table 3. The results clearly demonstrated that the increase of the complex interface was mainly caused by the N-terminal loop. The RMSD of the N-terminal loop in relative to its initial conformation in Complex I further confirmed that it returned back toward nivolumab in Equ3 ( Figure 8F). The N-terminal loop bound to nivolumab with six H-bonds (Bond No. 1-6 in Table 3) formed by PD−1 ARG30 and PD−1 PRO28 after 10 ns in Equ3, but it kept free in Equ1 and Equ2 until the end of simulations (Figures 8G-I).
The buried SASA of Complex II-N-rotated increased to around 1500 Å 2 in all three runs (Figures 9A-C). Its interaction energy decreased to around −350 kcal/mol in Equ1 and Equ3, and to around −250 kcal/mol in Equ2. The number of H-bonds of Complex II-N-rotated increased to 4 in Equ1 and Equ2, but to 6 in Equ3. Next, the number of H-bonds formed by the N-terminal loop as well as its RMSD in relative to its  initial conformation in Complex II were calculated, as shown in Figures 9D-F and H-bonds with a survival ratio of above 0.2 are listed in Table 3. The results reveal that the N-terminal loop rebuilt the complex interface and formed stable H-bonds with nivolumab in all three runs after 20 ns, especially in Equ3. PD−1 ARG30 and PD−1 ASP26 on the N-terminal loop formed two H-bonds with ASN31 H and LYS57 H of nivolumab (Bond No. 16 and 22 in Table 3) in Equ1, whereas PD−1 LEU25 formed one bond with TYR53 H in Equ2 (Bond No. 17 in Table 3). The interface between the N-terminal loop of PD-1 and the nivolumab in Equ3 was most stable, with five H-bonds formed by PD−1 ARG30 with ASN31 H and ASP100 H of nivolumab (Bond No. 15,(18)(19)(20)(21) Table 3 and Figures 9G-I).
Overall, on the nanosecond time scale, the N-terminal loop of PD-1 prefers to interact with nivolumab to stabilize the complex interface further. Interfaces of Complex I-N-rotated and Complex II-N-rotated tend to be stronger with the help of the N-terminal loop of PD-1. The binding strength indexes showed bimodal distributions, especially for the interaction energy of Complex II-N-rotated (red lines in Figure 4).

Binding of the N-Terminal Loop With Nivolumab Could Induce the Rebinding of the IgV Domain With Nivolumab
For the IgV-rotated complexes, the RMSD of heavy atoms showed that Complex I-IgV-rotated fluctuated more violently than Complex II-IgV-rotated in three runs (Supplementary  Figures S7, S8). Buried SASA, interaction energy, and number of H-bonds of Complex I-IgV-rotated and Complex II-IgV-rotated are shown in Figures 10A-C, 11A-C, respectively, and their distributions are demonstrated in Figure 4. For all three runs of Complex I-IgV-rotated, the buried SASA fluctuated around 700 Å 2 , the interaction energy fluctuated around −70 kcal/mol, and the number of H-bonds fluctuated around 3 (Figures 10A-C). The interface of complex I-IgV-rotated was mainly contributed by the N-terminal loop of the PD-1 within our simulation time (Figures 10D,E), with seven H-bonds formed by PD−1 PRO28, PD−1 ASP29, and PD−1 ARG30 with ASN31 H , GLY33 H , TYR53 H , and ASN99 H of nivolumab   Table 4 and Figures 10G-I). The RMSD of the IgV domain in relative to its *The name of the residues with H or L indicating that the residues are on the heavy or the light chain of nivolumab. # Equ1, Equ2, and Equ3 donate three runs. Max represents the maximum value of three survival ratios of a bond. Ave represents the average value of three survival ratios of a bond. Std represents the standard deviation of three survival ratios of a bond.
initial conformation in Complex I showed that the IgV domain of PD-1 remained free in all three runs, which caused great fluctuations ( Figure 10F). By contrast, for all three runs of Complex II-IgV-rotated, the buried SASA increased to nearly 1,600 Å 2 , the interaction energy decreased to nearly −300 kcal/mol, and the number of H-bonds increased to around 6. The number of H-bonds formed by the N-terminal loop and the IgV domain of PD-1, as well as the RMSD of the IgV domain in relative to its initial conformation in Complex II are shown in Figures 11D-F. It can be seen that the N-terminal loop of PD-1 maintained a stable interface with nivolumab on the nanosecond time scale, and changes were mainly caused by the IgV domain, which got close to and formed stable interface with nivolumab in all three runs. The FG loop ( PD−1 LYS131) and BC loop ( PD−1 THR59) of PD-1 interacted firmly with nivolumab (ASN31 H , ASN99 H , and ASP101 H ) by forming four bonds in Equ1 26 in Table 4 and Figure 11G). PD−1 LYS131 of FG loop and PD−1 GLU61 of BC loop formed five H-bonds with ASP50 L and THR28 H of nivolumab in Equ2   Table 4 and Figure 11H). The FG loop (i.e., PD−1 LYS131 and PD−1 PRO130) formed two H-bonds with ASP101 H and THR56 L of nivolumab in Equ3 (Bond No. 19-20, Figure 11I). Thus, in our simulations, binding of the N-terminal loop with nivolumab could further induce the interaction of the IgV domain with nivolumab, which would switch the interface of nivolumab/PD-1 complex to a stronger binding state (purple lines in Figures 4D-F).

DISCUSSION
The N-terminal loop of PD-1 has not attracted much attention in recent years because it is not a binding region for PD-L1. However, it was proven critical for the binding between PD-1 and nivolumab, which is a humanized IgG4 antibody approved by the FDA for several cancers. Crystal structures of nivolumab/PD-1 complex showed that the interaction residues of PD-1 locate on the N-terminal loop and IgV domain. The N-terminal loop of PD-1 greatly contributes to the complex formation, and we believe that it should also play an important role in the dynamic molecular recognition process. As dynamics of terminal loops are hard to predict based on crystal structure alone, eight molecular systems of nivolumab/PD-1 complex with different binding states were established, and long-time MD simulations of three replicas were performed for each of them, with the total simulated time of 2.4 µs. Our results proposed a two-step binding mode, in which the N-terminal loop of PD-1 could switch the complex interface into a stronger binding state. When the IgV domain binds to nivolumab first, the N-terminal loop of PD-1 prefers to interact with nivolumab to stabilize the complex interface. When the N-terminal loop is occupied by nivolumab, it could further induce the binding between the IgV domain (i.e., the FG and BC loops) of PD-1 and nivolumab.
The present results provided a detailed picture on the dynamic properties of the N-terminal loop of PD-1 in molecular interactions. Although this is the first time to systematically study the function of the N-terminal loop of PD-1, the N-terminal loops of other proteins have been revealed to have similar regulatory functions. For example, the N-terminus of model protein thaumatin serves as a major binding for cisplatin fragments (Russo Krauss et al., 2016). The N-terminal loop residues of beta-amyloid plays a key role in its interactions with integrin receptor and cell surface (Venkatasubramaniam et al., 2014). The N-terminal loop region of A1 domain in von Willebrand factor could stabilize A1A2A3 complex and modulate platelet activation under shear stress (Ju et al., 2013). Surface-exposed loops, generally as the most flexible parts of a protein structure, are not mere connectors but also have the potential to interact with solvent, ligands, and other biomolecules (Papaleo et al., 2016).
Since loop regions are too flexible to be resolved by crystallography, our simulations pave the way for investigating the binding mechanism between PD-1 and nivolumab. With the proposed two-step binding mode, nivolumab might be at least twice as likely to bind PD-1 as other antibodies with only one binding site. Furthermore, due to the high flexibility and mobility of the N-terminal loop, it can greatly facilitate the scanning efficiency and thus increase the probability of PD-1-nivolumab binding. This is of great importance for molecular interactions in the crowded intracellular environment.
Besides the binding process, our work also revealed the role of the N-terminal loop in augmenting the PD-1-nivolumab residence time, which is defined as the reciprocal of the dissociation rate constant. An abundance of experimental data suggests that high-affinity drug interactions with macromolecular targets generally rely on multistep binding and dissociation described by the two-step, induced-fit model (Copeland, 2016). Here, we show that the dissociation trajectory for the PD-1/nivolumab complex probably involves a retrograde induced-fit mechanism, that the N-terminal loop of PD-1 is able to rebind the dissociated nivolumab and IgV domain before nivolumab is completely released from PD-1 (Figure 11). Thus, the retrograde induced-fit mechanism creates multiple kinetic and structural barriers to nivolumab dissociation. This might partially explain why nivolumab has a nearly 10-fold higher affinity than that of pembrolizumab (Kd = 3.06 vs. 29 pM) (Fessas et al., 2017), another humanized anti-PD-1 monoclonal antibody approved by FDA, although the epitope of nivolumab is distinctly smaller than that of pembrolizumab (buried surface areas = 1,487 vs. 2,126 Å 2 ). This model suggests that the N-terminal loop of PD-1 may be viewed as a "kinetic gatekeeper" that guides the docking of nivolumab onto the IgV domain and prevents nivolumab from dissociating.
Apart from the role of pulling the IgV domain back to nivolumab, the N-terminal loop may also increase the nivolumab residence time by shielding the IgV-nivolumab H-bonds from water. In previously reported studies, the relationship between long residence time and accessibility of water has been established (Schmidtke et al., 2011;Magarkar et al., 2019). In this study, binding of the N-terminal loop with nivolumab can create an environment of lower dielectric constant around the BC loop, shielded from bulk solvent (Figure 7). Such a solvent-shielded environment might result in the strengthening of non-covalent forces between the BC loop and nivolumab, such as hydrogen bonds, hydrophobic forces and van der Waals forces.
The function of the N-terminal loop of PD-1 for nivolumab requires its structure to be as complete as possible. The N-terminal loop of PD-1 in Complex I-N-rotated only bound with nivolumab in one run (Equ3 for complex I-N-rotated, Figure 8I), and it did not induce the rebinding of the IgV domain with nivolumab in the simulations (Figure 10). By contrast, the N-terminal loop of PD-1 interacted with nivolumab in all three runs of Complex II-N-rotated (Figure 9) and successfully induced the binding of the IgV domain with nivolumab in the simulations of Complex II-IgV-rotated (Figure 11). The reason is that the N-terminal loop of PD-1 contained 10 residues ( PD−1 LEU25 to PD−1 PRO34) in Complex II, but it only eight residues ( PD−1 SER27 to PD−1 PRO34) in Complex I. The missing residues in Complex I, such as PD−1 LEU25 and PD−1 ASP26, could also form stable H-bonds with nivolumab    Table 4). Therefore, the N-terminal loop of PD-1 in Complex II-N-rotated is more likely to be captured by nivolumab. Moreover, the interface between the N-terminal loop of PD-1 and nivolumab in Complex II-IgV-rotated (around 800 Å 2 and -140 kcal/mol) was stronger than that in Complex I-IgVrotated (around 700 Å 2 and −70 kcal/mol), which is more beneficial in pulling the IgV domain back to nivolumab.
Despite the use of massive computational resources and highly precise models (full atomic representation and detailed force field), plain all-atom MD simulation is still insufficient for adequately exploring biomolecular structural dynamics. Multiple evidences indicate that long simulations cannot address how to catch the possible transition paths, which are still rare during µslong MD due to inherent stochasticity and high-energy barriers. Nevertheless, we may predict some approximate behavior from simulations that suffer from sampling inefficiencies, in certain conditions e.g., upon introducing mutations or relaxation after removing ligands (Orellana, 2019). Here, we employed different starting configurations and multiple short simulations to enhances the configuration space sampling to better probe the conformational changes (Knapp et al., 2018;Chen et al., 2019b,c; TABLE 4 | Summary of survival ratio and involved residues of H-bonds in the interface of complex I-IgV-rotated and complex II-IgV-rotated. *The name of the residues with H or L indicating that the residues are on the heavy or the light chain of nivolumab. # Equ1, Equ2, and Equ3 donate three runs. Max represents the maximum value of three survival ratios of a bond. Ave represents the average value of three survival ratios of a bond. Std represents the standard deviation of three survival ratios of a bond. Chen et al., 2020). Fortunately, the conformational transitions of the nivolumab/PD-1 complex we concerned have been observed in a simulation time of 100 ns. For example, the N-terminal loop of PD-1 rotated back toward the interface and interacted with nivolumab in all three runs of Complex I-N-rotated (Figure 9). Moreover, the IgV domain of PD-1 also rotated back toward the interface and interacted with nivolumab with the help of the N-terminal loop of PD-1 in the simulations of Complex II-IgV-rotated (Figure 11). This work provides useful dynamics information on the role of the N-terminal loop in the molecular recognition process between PD-1 and nivolumab. Previous research usually focused on the "hot-spot" of PD-1, such as the FG loop and the C'D loop. Our results suggest that the N-terminal loop of PD-1, which acts as an important gatekeeper for the binding of nivolumab and PD-1, should also be considered in the anti-PD-1 blockade antibody design. We are hopeful that the results presented in this study will ultimately provide a theoretical framework to understand the structural landscape of N-terminal loop of PD-1. In general, this opens a new opportunity for medicinal biologists or chemists to optimize affinity for antibodies, if such gatekeepers can be identified.

DATA AVAILABILITY STATEMENT
All datasets presented in this study are included in the article/Supplementary Material.