Inhibition of the checkpoint protein PD-1 by the therapeutic antibody pembrolizumab outlined by quantum chemistry

Much of the recent excitement in the cancer immunotherapy approach has been generated by the recognition that immune checkpoint proteins, like the receptor PD-1, can be blocked by antibody-based drugs with profound effects. Promising clinical data have already been released pointing to the efficiency of the drug pembrolizumab to block the PD-1 pathway, triggering the T-lymphocytes to destroy the cancer cells. Thus, a deep understanding of this drug/receptor complex is essential for the improvement of new drugs targeting the protein PD-1. In this context, by employing quantum chemistry methods based on the Density Functional Theory (DFT), we investigate in silico the binding energy features of the receptor PD-1 in complex with its drug inhibitor. Our computational results give a better understanding of the binding mechanisms, being also an efficient alternative towards the development of antibody-based drugs, pointing to new treatments for cancer therapy.


Results
The programmed cell death protein 1 (PD-1) is an important regulator for the immune tolerance and T cell exhaustion, being recently emerged as a key target in the treatment of several types of cancer. It is expressed after the T cell activation and binds to the ligants PD-L1 and PD-L2, suppressing immune response against autoantigens and playing an important role in the maintenance of peripheral immune tolerance 13 . However, the ligant PD-L1 is often overexpressed in different tumor including lymphoma, melanoma, non-small-cell lung cancer and other types, making the PD-1/PD-L1 signaling pathway crucial in dampening the immune surveillance of the tumor 29 .
In this context, the target of the PD-1/PD-L1 interaction with monoclonal antibodies has demonstrated to be an important strategy for the control and eradication of several types of cancers. The treatment with pembrolizumab was approved by US-FDA in 2014 for advanced melanoma, and recently it has also been approved to an increasing number of cancer types, such as Hodgkin lymphoma and non-small-cell lung cancer. It is believed that the weak frequency in humans and the induction of cell activation, characteristics of the members of IgG4 subclass, are some of the main basis for the success of this therapeutic compound 23 . As we can see from Fig. 2, the pembrolizumab recognition surface on PD-1 is filled by intermolecular direct and water-mediated hydrogen bonds (HBs and wHBs, respectively), non-conventional hydrogen bonds (nHBs), salt bridges (SB) and hydrophobic contacts. Therefore, mapping the relevant interactions among the pembrolizumab/PD-1 surface complex is highly important to the rationale drug development of new compounds. In this work, we employed the MFCC/DFT scheme to assess the relative energetic contribution of the individual as well as pair-interaction contribution of each residue from the pembrolizumab/PD-1 complex. From now on we will describe all the steps in this recognition process.
All other residues shown attractive interactions (negative energies), with the most intense one being observed for the heavy-chain amino-acids Y33 HC (ε 20   Pembrolizumab heavy-chain/PD-1 receptor interaction energy. Among the 408 residue-residue pairs analyzed here, 260 were pembrolizumab heavy-chain/PD-1 interactions. It is a reflection of the proximity between the HC/LC pembrolizumab with the receptor PD-1 represented in crystallographic structures, where pembrolizumab Fab heavy-chain fragment is closer to the PD-1 receptor than the light-chain one 26,30 . Besides a higher number of pairs, the sum of the energetic interaction between the pembrolizumab heavy-chain and the PD-1 receptor amino-acids shows also the higher value, accounting for -142.50 kcal mol −1 (−138.33 kcal mol −1 ) for the dielectric constant ε 20 (ε 40 ).
Although Y33 HC has been shown to be the most energetic pembrolizumab amino-acid residue, it only interacts with 15 residues from the PD-1 receptor. Figure 4a depicts the highest interaction energies calculated to Y33 HC /PD-1 residues. As one can see, the strongest energy of Y33 HC is mainly related to its binding with three residues: K78 (ε 20 Figure 5a displays how close the residues Y33 HC , Y101 HC , R102 HC and Y103 HC are to the PD-1 receptor. Analyzing it one can understand the reason why these amino-acids present some of the largest number of pair-interaction and binding energies. All of them are involved in a network of hydrogen bonds with the PD-1 residues from CC'FG β-strands and some of its loops, mainly those belonging to the C'D loop, which was described to intrude into the complementary determining region (CDR) of the pembrolizumab drug, a variable portion present in some igG molecules responsible for the recognition of specific antigens 30 . Figure 5b depicts some of these interactions, starting with the residue Y33 HC . Here one can see that it forms two direct hydrogen bonds with the residues K78 charged amine group and Q88 side-chain oxygen. Besides, the residue Y33 HC also makes a water-mediated hydrogen bond with the K78 carbonyl group of the main-chain, and is involved in non-polar contacts with P89. The two hydrogen bonds between the residues Y33 HC -K78 give the major contribution to their high energy.
The residues Y101 HC , R102 HC and Y103 HC interaction network is shown in Fig. 5c-e. Figure 5c depicts 3 directs, 2 water-mediated and 1 non-conventional hydrogen bonds. The strongest interaction of this residue is related to the pair Y101 HC -K78, where it creates a direct hydrogen bond with charged side-chain amine and a non-conventional hydrogen bond with amine of the main-chain. Three hydrogen bonds are formed with the residue T76, two with its oxygen atoms (carbonyl and hydroxyl) and a water-mediated with the hydroxyl group. It is also involved in one wHB with the residue Y68 side-chain hydroxyl and in some non-polar contacts with the residues D77 and N66. The residues R102 HC and Y103 HC are mainly involved in weak interactions. The latter form a σ-π interaction with the residue V64, a non-conventional hydrogen bond with the residue P83 and some non-polar contacts with the residues L128 and D85, while the former makes hydrogen bonds with A132 oxygen atom and N66 side-chain nitrogen, besides non-conventional hydrogen bonds with the side-chain of the residues I126 and K78.
The MFCC scheme yields not only the individual pembrolizumab/PD-1 interactions but also important information from the pair-interactions that could be otherwise missing, namely:  The residues M105 HC and D108 HC make hydrogen bonds with the residue K131. The main-chain oxygen carbonyl of M105 HC makes a direct hydrogen bonds with the charged amine of this lysine residue, while the residue D108 is engaged in a water-mediated hydrogen bond with the same molecular group, yet presenting its negatively charged carboxyl group (D108) in a closer distance to positively charged amine group (see Fig. 5f). The pair Y35 HC -S87 binding energy is governed by a hydrogen bond formed between the hydroxyl group of this tyrosine residue (Y35 HC ) and the main-chain carbonyl from S87. The pair R99 HC -D85 is described for some authors as the only salt bridge formed between the pembrolizumab drug and the PD-1 receptor 26,30 . Pembrolizumab light-chain /PD-1 receptor interaction energy. The pembrolizumab light-chain is a little more distant of the PD-1 receptor than the heavy-chain one, as it can be inferred from its contribution to the total binding energy (ε 20 : −44.67 kcal mol −1 ; ε 40 : −47.75 kcal mol −1 ), as well as from the number of pairs formed with the receptor (148). This is consistent with previous crystallographic results 25,26 . According to Fig. 6, the most intense binding energies for the LC amino-acids is associated to the pair E59 LC -K131 (ε 20 : −11.90 kcal mol −1 ; ε 40 : −9.70 kcal mol −1 ), being followed by the pairs S95 LC -S87 (ε 20 : −6.00 kcal mol −1 ; ε 40    7 depicts some interactions made by the pembrolizumab light-chain residues linked to the PD-1 receptor. As one can see, the residue Y36 LC is surrounded by some charged residues (E84, D85 and R86) from the PD-1 receptor (Fig. 7a). However, it only makes a single direct hydrogen bond with the main-chain carbonyl of the residue E84, all the other interactions being non-conventional hydrogen bonds (P83 and R86), and water-mediated hydrogen bonds. It would be inconsistent with the position occupied by this residue among the most energetic one of the pembrolizumab light-chain, if the number of pairs formed with Y36 LC (10) would not overcome it. The residue making the biggest number of pair-interaction is Y34 LC , namely 19. Notwithstanding, the residue Y34 LC shows a binding energy less than the Y36 LC one, due to an unfavorable interaction with the residue F63 (ε 20 : 5.27 kcal mol −1 ; ε 40 : 5.30 kcal mol −1 ).
Although only five pair-interactions were calculated with the residue S95 LC within a radius of 8.0 Å, it shows one of the highest interaction energies of the drug's light-chain. Similar to the residue Y36 LC , the majority of the surrounding residues make weak interactions, excluding S87 which forms a water-mediated hydrogen bond  through the hydroxyl group from both side-chains (see Fig. 7b). It is also important to notice the attractive binding energy found to the R96 LC -R86 pair, even with the positively charged guanidine group from both amino-acids, assuming a conformation where they are very close. These two arginine residues are in a T-shaped stacking interaction, which favors an attractive bind 30 . Finally, the E59 LC -K131 pair forms the second salt bridge of the complex drug/receptor, depicting the highest individual interaction energy in the drug's light-chain. This interaction is due to the negatively charged carboxyl group of the residue E59 LC , and the positively charged amine of the residue K131 side-chain. It gives us the idea of the dynamic process which govern the interaction between the pembrolizumab Fab fragment and the extracellular region of the PD-1 receptor, such as the formation of a new salt bridge between the residues D108 HC -K131 whose opposite charges from the side-chain are very close (3.6 Å).
For completeness, we display in Fig. 8 the electrostatic potential isosurface with projected electron densities for the pembrolizumab amino-acids bound to some of the most important residues at the binding pocket site.

Total binding energy.
To evaluate the binding interaction energies through fragment-based quantum mechanics method, it is important to take into account every significant attractive and repulsive amino-acid residue which can influence this mechanism. Therefore, instead of taking an arbitrary region of the binding site, we performed a search for an optimal binding pocket radius (r) in which a variation less than 10% of the sequential pocket radius could be observed after a radius increase. For this task, the binding pocket radius r is varied from 2.0 Å (−32.04 kcal mol −1 for ε = 20 and −27.29 kcal mol −1 for ε = 40) to 8.0 Å (−187.17 kcal mol −1 for ε = 20 and −186.08 kcal mol −1 for ε = 40) in order to determine the best value of r, found to be 6.0 Å corresponding to an energy of −179.94 (−179.60) kcal mol −1 for ε = 20 (40), from which the convergence was achieved (see Fig. 9). This result is in agreement with previous works, stating that after 6.0-7.0 Å the molecular interactions are weak 31,32 .

Discussion
Cancer immunotherapy research is trying to overcome the cancer's ability to resist the immune responses by stimulating the body's own mechanisms to remain effective against the disease. Early clinical results using blocking agents against the human cell surface receptor PD-1 and its ligants PDL-1 and PDL-2, point to unprecedented rates of long lasting anti-tumor activity in patients with metastatic cancers of different histologies 33 . From the pharmaceutical point of view, however, the clinical development of the PD-1 pathway blockers requires an understanding of the signals that induce expression of its ligands within the tumor.
However, there are some patients showing resistance to the blockage of the receptor PD-1 34 . Over stimulation of immune responses, on the other hand, may lead to damage normal, healthy tissue. Therefore, to be more efficient, the design and production of such a pharmaceutical drug require a precise knowledge not only of its biochemical structure, but also its binding mechanism with the receptor PD-1, leading to a better understanding of how this complex influences the cancer tumor and its anti-angiogenesis therapies as well as its microenvironment. Computational tools based on quantum chemistry may be a promising step in this route, further increased by combination with other anti-cancer therapies.
Within this context, we presented here an in silico quantum biochemistry calculation of the electronic structure of the complex pembrolizumab Fab drug/PD-1 receptor in order to map its recognition surface, searching for the binding interactions that stabilize it. The quantum MFCC approach used here is a route to investigate accurately large biological systems with low computational cost, and has been applied previously to describe molecular interactions at the quantum level related to the collagen stability 35 , central nervous system disorders 36 , and breast cancer 37 , to cite just a few. There is no other quantum mechanical study relating the binding of the pembrolizumab drug at the PD-1 receptor, albeit its high pharmaceutical relevance.
After a convergence study on the size of the binding pocket sphere, we have considered all ligand-residue interactions within a pocket radius of 8.0 Å from the ligand. In addition, we have established the residues that play the most important roles on the binding affinity, namely: Y33 HC > R102 HC > F103 HC > Y101 HC > E59 LC > M105 HC > N52 HC > Y36 LC > S95 LC > D104 HC . We also highlighted the energetic and structural differences between the pembrolizumab's heavy-chain and light-chain besides unveiling the relevance to take into account the water molecules and weak interactions in the design of this therapeutic compounds. The strong attractive character observed in the amino-acid residues Y101 HC , R102 HC , F103 HC , D104 HC , M105 HC and D108 HC could give the base for the development of an improved drug targeting the PD-1/PD-L1 pathway by activating an immune resistance mechanism in response to endogenous anti-tumor activity, avoiding an immunosuppressive tumor microenvironment. It is also remarkable the relevance of the residues K78, D85, P89, I126 and L128 of the receptor PD-1. These amino-acids have been shown making pairs with almost all pembrolizumab residues, and mutational studies have displayed their relevance to the PD-1/ligand interaction 38 .
Some of these residues binding energies spectra are consistent with previous experimental studies 26,39 . Our results show also that in the Pembrolizumab/PD-1 complex, the residue D85 is associated with the strongest interaction energy regime, forming pairs with the residues R99 HC and Y36 LC through a salt bridge and hydrophobic contacts, respectively, while the residue R86 has no interaction whatsoever with any of the residues from the ligand, in agreement with the experimental results 30 . As shown by Fessas et al. 40 , one of the most relevant structural aspect of the Pembrolizumab/PD-1 complex is the interaction of the former with residues from the C'D loop of PD-1 (amino-acids P83-S93), emphasized in our previous section (Results), excluding the amino-acids D100 HC , Y101 HC , Y53 LC , Y57 LC and E59 LC , commonly showing one energetically relevant pair. These achievements confirm the relevance of this loop for anchoring Pembrolizumab and, likely, its ligand-based compounds.
Despite the absence of mutation experimental results, outcomes based on the crystallographic structures of the PD-1 protein mutants bound to the receptors PD-L1 and PD-L2 could be an useful background to understand the energetic behavior depicted in this work. Recently, Lázár-Molnár performed a mutagens study taking into account the wild-type human PD-1 receptor and a mutant form of the protein A132L 41 . Mutation of the amino-acids K78A in both proteins (wild-type and mutant) was responsible for the complete loss of the PD-L1 and PD-L2 interaction bind energy. Notwithstanding, K78, a residue from the C' strand, is seen forming strong interaction pairs with a number of amino-acids from the Pembrolizumab, including the most energetic heavy-chain residues Y33 HC , D100 HC , Y101 HC , R102 HC and F103 HC . Also, although the amino-acid I126 makes non-conventional hydrogen bonds only with the heavy-chain residue R102 HC , it is responsible for one of its strongest interaction energies. Meanwhile, the residues I134, L128 and E136, although not presenting a large number of energetically strong pairs, are related to R102 HC and Y103 HC , which are the second and third most energetic Pembrolizumab heavy-chain residues, thus defining, along with the residue I126, a region of stabilization of the ligand. Taken together, it becomes clear the relevance of these residues for the inhibition mechanism of the Pembrolizumab.
To close this section, we would like to propose two pharmacophores that may be useful for the rational design of ligands targeting the receptor PD-1, looking for the development of antibody-based drugs. As one can note, the Pembrolizumab/PD-1 binding surface is mainly governed by a set of hydrophobic contacts strengthened by strong hydrogen bonds. Five residues from the PD-1 receptor can be observed to form two hot-spots required for the inhibition of the PD-1/PD-ligands. The first one is in the C'D loop and C' sheet, while the second one is in the FG loop and F sheet, respectively. The former is composed by the residues K78, D85 and P89, with their most distinctive feature being the steric complementarity between the residue K78 and the hydrophobic/hydrophilic heavy-chain amino-acids Y33 HC , Y101 HC , R102 HC and F103 HC . They are coupled to flexibilize the C'D loop, permitting the residue D85 to interact with the light-chain amino-acids S36 LC and S95 LC . This pocket is ruled by hydrogen bonds and is, likely, one of the first docking points on the PD-1 surface due to its more external position, when compared to other motifs. The latter is located just nearby and has a most hydrophobic character, being composed by the residues I126 and L128. These residues are seen to accommodate the heavy-chain amino-acids R102 HC and F103 HC by anchoring them through hydrophobic interactions and non-conventional h-bonds. These proposed pharmacophores should provide a solid point for the design of new antibodies or small molecules targeting the receptor PD-1.
In summary, we have provided a quantum mechanics description of the Pembrolizumab/PD-1 complex at a molecular and energetic level by decomposing both proteins into amino-acid fragments and calculating the residue-residue interaction pair through an MFCC-based method within the density functional theory framework. The use of quantum mechanics in our calculations provided a good characterization of important molecular features by including an accurate electronic description and enhancing the differences of interaction intensity among the amino-acid pairs. Since Pembrolizumab has shown great improvement in the treatment of several cancer types, our results could favor the development of new compounds targeting the PD-1/PD-ligands, in special those based on the hot-spots, by avoiding structural features which are in minor relationship with this binding scenario. Notwithstanding, the quantum chemistry computational methods used in this work emerged as a simple and efficient alternative to unveil the drug's amino-acids residues that play the most important role on the binding affinity of the receptor-ligand complex. No doubt, taking into account the small cost/benefit of the operation, in silico approach, not only in clinical oncology but also in pharmaceutical medicine in general, is an important initial step toward the development of efficient cancer pharmaceutical drugs.

Methods
The calculations performed in this study have taken full advantage of the first X-ray crystal structure of human programed cell death protein 1 (PD-1) solved in complex with a Fab fragment of pembrolizumab (PDB ID: 5GGS) at 2.00 Å of resolution. The protonation state of the receptor PD-1 and its ligand was adjusted according to the results obtained from the PROPKA 3.1 package 42 . Afterwards, atoms not resolved by X-ray diffraction and therefore absent in the crystallographic files, as hydrogens and some amino-acid side-chains, were added to the structures and submitted to a classical geometry optimization fixing the other atoms. This optimization was performed using the classical force field CHARMm22 (Chemistry at Harvard Molecular Mechanics 22 forcefield), which is especially parametrized for organic molecules 43 . The calculations were carried out with convergence tolerances set to 10 −5 kcal mol −1 (total energy variation), 0.001 kcal mol −1 Å −1 (RMS gradient) and 10 −5 Å (maximum atomic displacement).
The interaction energy between the amino-acids from the pembrolizumab (at position R i ) and PD-1 (at position R j ) was calculated by using the formalism of MFCC fragmentation scheme proposed by Zhang and Zhang 44,45 , and modified by Rodrigues and co-workers to calculate protein-protein interactions 28 , within a DFT framework. For each amino acid residue of interest at position R i , we draw an imaginary sphere with radius equal to 8.0 Å, and evaluate the interaction energy EI(R i -R j ) with each residue at position R j , considering at least one atom inside the sphere. Hence, the scheme to compute the residue-residue interaction energies follows the MFCC approach, in which the interaction energy between the residue R i and the other residues R j was calculated according to: where the C i and C i* (C j and C j* ) caps are respectively the neighbour residues R i−1 and R i+1 (R j−1 and R j+1 ) attached to the reference residue R i (R j ). The term E(C i R i C i* − C j R j C j* ) corresponds to the total energy of the fragment comprises both residues as well as their capped residues. The second (third) term, E( , gives the total energy of the system formed by the capped residue R i (R j ) and the hydrogenated caps of R j (R i ). E(C i C i* − C j C j* ) is the total energy of the system formed by the caps only. In order to achieve the structural stability of the complex promoted by interactions with extended hydration network, all water molecules forming hydrogen bonds with a particular residue were included for completeness in the fragments. The analysis of the binding scenario were based in the criteria adopted by Bezerra et al. 35 .
Energetic calculations for each residue were performed using the Gaussian G09 code 46 , within the density functional theory (DFT) formalism, with the generalized gradient approximation (GGA) functional B97D (for a review see ref. 47 ). It was proposed by Grimme 48 and include dispersion terms in order to improve the description of the non-covalent interactions, being recently used together with the D3 dispersion correction under vacuum and COSMO continuum solvation model to calculate the binding energy of some protein ligand complexes within a MFCC scheme 49 . We selected the 6-311 + G(d,p) basis set for the calculations of the expanded Kohn-Sham orbitals for all electrons.
The representation of molecular environment is a necessary step for the theoretical study of molecular properties. Recently, models of implicit solvation have been used in the MFCC scheme by Ourique et al. 50 and Dantas et al. 51 to calculate binding affinity in protein-ligand complexes, by assuming different values of the dielectric constant (ε) to take into account the effect of the protein and solvent environment in the evaluation of the electrostatic energies. Here, we have employed the MFCC approach together with the CPCM (Conductor-like Polarizable Continuum Model) continuum solvation model 52 with dielectric constants ε considered to be 20 (ε 20 and 40 (ε 40 ), respectively) to increase the similarity with the protein environment and estimate energy effects, such as the electrostatic polarization promoted by the solvent. Data availability. All other relevant data are available from the corresponding authors upon reasonable request.