Combining Well-Tempered Metadynamics Simulation and SPR Assays to Characterize the Binding Mechanism of the Universal T-Lymphocyte Tetanus Toxin Epitope TT830-843

Abstract

Peptide TT830-843 from the tetanus toxin is a universal T-cell epitope. It helps in vaccination and induces T-cell activation. However, the fine molecular interaction between this antigen and the major histocompatibility complex (MHC) remains unknown. Molecular analysis of its interaction with murine MHC (H-2) was proposed to explore its immune response efficiency. Molecular dynamics simulations are important mechanisms for understanding the basis of protein-ligand interactions, and metadynamics is a useful technique for enhancing sampling in molecular dynamics. SPR (surface plasmon resonance) assays were used to validate whether the metadynamics results are in accordance with the experimental results. The peptide TT830-843 unbinding process was simulated, and the free energy surface reconstruction revealed a detailed conformational landscape. The simulation described the exiting path as a stepwise mechanism between progressive detachment states. We pointed out how the terminus regions act as anchors for binding and how the detachment mechanism includes the opening of α-helices to permit the peptide’s central region dissociation. The results indicated the peptide/H-2 receptor encounter occurs within a distance lesser than 27.5 Å, and the encounter can evolve to form a stable complex. SPR assays confirmed the complex peptide/H-2 as a thermodynamically stable system, exhibiting enough free energy to interact with TCR on the antigen-presenting cell surface. Therefore, combining in silico and in vitro assays provided significant evidence to support the peptide/H-2 complex formation.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

Clostridium tetani is a pathogenic gram-positive bacterium causing tetanus. This bacterium can be isolated from soil and animal intestinal tracts and, as such, can contaminate many surfaces and substances [1]. Clostridium tetani is a restricted anaerobic bacillus, a spore producer that allows it to survive under aerobic conditions and produces an exotoxin, a tetanus toxin, that attaches itself to the system nervousness, provoking the symptomatology of the disease [2].

The tetanus toxin is a peptidase tentoxilysin (Clan MA, M27 family) [3], which is produced by the anaerobic spore-forming bacteria C. tetani as a single-chain polypeptide that is 1315 amino acids (aa) in length [4]. The protein is endogenously processed to yield the 52 kDa light chain (457 aa) and the 98 kDa heavy chain (858 aa). A disulfide bridge forms a dimeric protein of 150 kDa, showing a toxic property due to its ability to bind to specific membrane receptors on presynaptic motor nerve cells. Therefore, this interaction causes interneuron discharge inhibition affecting the motor and autonomic nervous system [5–8].

The heavy chain fragment of the tetanus toxin is marked by its antigenic and immunogenic properties [9], which is related to effective, specific, and safe antitoxin vaccines [10]. This fragment stands out by its ability to induce the production of antibodies when used as a fusion partner for foreign antigens [11–13]. Tetanus toxin fragment adjuvant potency is related to the presence of promiscuous T-lymphocyte epitopes in the protein [14]. The tetanus toxoid protein contains several T-lymphocyte epitopes (p2, p21, p23, p30, and p32), which have been extensively studied. They were demonstrated to be universally immunogenic T-cell epitopes in both mice and humans. Among these peptides stands the p2 epitope TT830-843 (QYIKANSKFIGITE), located at the protein C-terminus region, which can bind major histocompatibility complex (MHC) proteins [14–16].

As a linear T-lymphocyte epitope, the TT830-843 sequence is related to T-lymphocyte activation and is used to improve vaccine potential. For example, it has been reported that a TCR (T-cell receptor) set can recognize the TT830-843 epitope presented by DRB1$\ast$1302 (human MHC class II DRB1$\ast$1302 allele). It indicates that TCRs with distinct CDR3s (complementary-determining region 3) in the conserved V beta-chain may not differ in how they recognize the ligand and may provide new insights into understanding the TCR/peptide/MHC supramolecular complex formation [17].

The TT830-843 sequence has been used as an adjuvant universal epitope of the immune system in different vaccine platforms due to its significant properties as a T-lymphocyte helper epitope. This peptide fragment enhances antigenicity and the overall immune response. In this context, the use of the TT830-843 epitope as a peptide conjugated to antigenic epitopes of interleukin-13 has been proposed to increase the efficacy of vaccine protective humoral immunity in murine asthma models [18]. Another example is type I diabetes in mice whose expression of the insulin B:9–23 sequence is required to develop anti-insulin autoimmunity [19, 20]. The sensibilization approach using the TT830-843 epitope as an adjuvant to modify insulin B:9–23, intranasally induced, provided significant suppression of diabetes [21]. Rodrigues-da-Silva et al. [22] demonstrated the fusion of the TT830-843 sequence in a synthetic construction against Plasmodium vivax-enhanced specific T- and B-cell responses to a vaccine candidate.

In addition, a study of the immunogenicity of hepatitis B virus epitope-based polypeptides to trigger a specific HLA I-restricted (human leukocyte antigen class I) CD8⁺ T-lymphocyte response was proposed by constructing mimetic peptides based on the introduction of the TT830-843 epitope to strengthen the T-lymphocyte response. This therapeutic construction improved the induction of CD8⁺ CTL-mediated (cytotoxic T-lymphocyte) cytotoxicity in HLA-A2⁺ human peripheral blood lymphocytes [23]. In Figure 1, we present a schematic model of the action of the tetanus toxin-derived peptide TT830-843.

[figure omitted; refer to PDF]

Even with all these data showing the enhanced immune response provided by the TT830-843 fragment, little is known regarding to the molecular mechanisms underlying this phenomenon. Bioinformatics tools can explain the molecular foundations of immunity and validate potential epitopes for vaccine candidates [24]. However, current in silico computational methods for the validation and investigation of T-cell epitope interactions are far from satisfactory and vary in degree of accuracy. For example, most of these current methods confirm experimentally predicted binding to MHC molecules of most peptides predicted, but with only ~10% of those shown to be immunogenic [25]. Therefore, we explore the feasibility of combining an in silico computational approach and an in vitro physicochemical assay.

In this work, the association of an in silico and a physicochemical approach was proposed to access the details of the molecular interaction between the TT830-843 peptide and the murine MHC (H-2) receptor. An in silico computational approach, well-tempered metadynamics (WTMetaD), was applied to enhance the sampling of the free energy configurational space. WTMetaD introduces a bias potential that acts on a select number of degrees of freedom, called collective variables (CVs) [26]. Additionally, the peptide binding affinity to the H-2 protein and the kinetic parameters of this association were accessed by physiochemical assays. Surface plasmon resonance (SPR) was used for this purpose, monitoring the changes in the refractive index at the surface of a carboxyl sensor chip (COOH5) [27].

2. Materials and Methods

2.1. Molecular Dynamics (MD) and Well-Tempered Metadynamics (WTMetaD)

MD simulation is an important method for understanding the physical basis of biological macromolecule structure and function. MD became an especially useful computational technique to simulate biological processes inside the cells [28]. With the understanding of the system’s internal motions and their implications, questions concerning particles’ conformational changes as a function of simulation time can be explored. Therefore, MD simulation results can be used to address questions about specific properties of biological entities that would be more difficult to address in real systems [28].

However, with ongoing advances in computational simulations, standard MD methods often fail to adequately explore the configurational space to accurately evaluate proteins’ thermodynamics and kinetic properties [29]. Standard MD simulation requires a large amount of computational time to run and provide meaningful data for analysis. This situation occurs since high free energy barriers separate the relevant equilibrium configuration states, and the simulation tends to revisit the same energy minimum. For example, in rare events, like peptide ligands unbinding from MHC proteins, the system is trapped in configurational space local regions over the simulation time scale. This behavior is because there are significant high free energy barriers to be overcome, and the simulation cannot move from one stability state to another, i.e., from the binding state to the unbinding state. Therefore, in this case, standard MD simulations cannot reproduce biological processes in a feasible computational time [29]. Enhanced sampling techniques can be applied to address this issue and successfully perform biological rare phenomenon simulations.

WTMetaD belongs to a class of techniques that enhance the sampling of certain degrees of freedom, known as collective variables (CVs). WTMetaD facilitates the sampling of the configurational space by introducing a bias potential that acts on this selected number of degrees of freedom [26]. This method helps to accelerate computer simulations by adding a history-dependent potential to the system. One can say WTMetaD is a standard dynamic simulation in which it is imposed a harmonic restraint on a set of CVs. Therefore, WTMetaD can obtain accurate results in a relatively short simulation time [30].

WTMetaD is successfully applied to simulate rare events, predict binding affinity, and investigate intermolecular interactions occurred during the simulations. Additionally, enhanced sampling approaches permit to reconstruct the free energy surface (FES) associated with the protein interaction dynamics as a function of CVs [26, 30, 31]. WTMetaD operates a dimensional reduction of the degrees of freedom of a system. Energetic landscapes obtained from WTMetaD can be used to understand the intermolecular interactions occurring in the metastable states visited by the system. It also helps to explore many possible transition pathways between different free energy minima. Enhanced sampling technics have been applied to obtain meaningful data to study the dissociation mechanisms in protein complexes [31], becoming an important in silico method for biological research.

2.2. Initial Molecular Structure Preparation for Peptide/H-2D^b Complex

MHC class I starting structure was obtained from the Protein Data Bank (PDB) [32]. PDB entry 1jpf was used as a template to build the peptide/H-2D^b complex structure. PDB entry 1jpf has a resolution of 2.18 Å and represents the most extended peptide sequence (with 11 amino acids) whose structure has been determined for the H-2D^b haplotype [33]. Peptide TT830-843 was modeled into the H-2D^b receptor binding groove by using the psfgen package in VMD [34] and submitted the resulting coarse starting model to a proper refinement protocol.

The H-2D^b haplotype (receptor) structure in complex with the peptide TT830-843 (ligand) was submitted to the Rosetta FlexPepDock web server [35] for the refinement of the coarse starting model. Rosetta FlexPepDock is implemented in the Rosetta modeling suite framework [36] and performs a flexible peptide docking refinement protocol allowing full flexibility for the peptide and receptor sidechains. A standard refinement FlexPepDock protocol was applied (with default options selected) to model the peptide/H-2D^b complex. The FlexPepDock web server optimized the peptide conformation within the H-2D^b protein cleft, carrying out two hundred independent simulations: the first hundred in high-resolution mode and the remainder by applying the low-resolution preoptimization step, with a high-resolution posterior refinement. Resulting models were ranked based on Rosetta’s generic full-atom energy score [36], measuring the peptide/H-2D^b complex’s binding affinity. The best-ranked atomic coordinate pose was used (the lowest energy structure) (see Table S2 in Supplementary file 1) as the initial structure for WTMetaD simulation. The Rosetta FlexPepDock web server also calculates bb-RMSD (RMSD calculated for all peptide backbone atoms) from the starting conformation. Complex submitted to Rosetta FlexPepDock was obtained from a mutation in an existing PDB structure; consequently, the bb-RMSD calculated during the refinement does not correspond to a deviation from a crystallographic structure. Therefore, these conditions explain why we considered the Rosetta energy score to choose the initial structure for WTMetaD simulation. Top 10 resulting models obtained from the Rosetta FlexPepDock webserver are shown in Figure S2 of Supplementary file 1. The stability of the peptide initial structure obtained from the Rosetta FlexPepDock webserver is demonstrated in Figure S3 of Supplementary file 1.

2.3. Simulation Setup

All simulations were performed using NAMD 2.13 [37] with the CHARMMM27 force field [38]. Electrostatic interactions were evaluated using the Particle Mesh Ewald (PME) algorithm with a grid spacing of 1 Å. Nonbonded interactions were truncated using a cutoff of 12 Å and a switching function starting at a radius of 10 Å. Protein-peptide complex was immersed in an orthorhombic box using periodic boundary conditions. System was inserted within a 10 Å (for x- and y-axes) and 40 Å (along the z-axis) layer of water molecules, containing around 13,200 TIP3P water particles [39]. The complex structure was oriented to keep the peptide’s exiting direction aligned with the z-axis. To neutralize the system, Na⁺ counterions were added. The simulations ran in physiological pH, and the protonation states of the protein and the peptide residues were selected to conform to this physiological pH range. Therefore, the protonation states of all residues were assigned according to the pKa values of their side chains using the psf builder [34]. In particular, the protonation states for histidine residues were set as HSE (with their ϵ nitrogen protonated). Equations of motion were integrated using a velocity Verlet integration algorithm with a timestep of 2 fs in an NPT ensemble. The SHAKE algorithm was used to constrain covalent bonds. Energy minimization of the starting structure involved 15000 steps of the steepest descent method. Next, the system was equilibrated in a three-stage protocol: (i) heating up from 10 K to 310 K by increasing the temperature by 10 K for every 100 steps, and with the CA atom positions restrained using a harmonic potential with a force constant of 0.25 kcal/mol/Ų; (ii) 1 ns of water equilibration at 310 K, applying harmonic restraints to the CA atom positions (0.25 kcal/mol/Ų force constant) to allow water molecules to fully envelope the complexes; and (iii) peptide equilibration (1 ns of MD simulation) at 310 K, in which the peptide was free, keeping harmonic restraints only for CA atoms of the H-2D^b protein (0.25 kcal/mol/Ų force constant).

For the simulations in the NPT ensemble, the temperature was maintained at 310 K by Langevin dynamics, and the pressure was kept constant (1 atm) by the Langevin piston method. During WTMetaD simulation, the peptide and the two α-helical segments that flank the N-terminus peptide binding region (α1-helix residues 56-70 and α2-helix residues 155-175) were free to move, while other CA atom positions were restrained using a harmonic potential with the same force constant of 0.25 kcal/mol/Ų. WTMetaD simulations were performed for 16 ns, with Gaussians height set at 0.05 kcal/mol, a new hill added to the WTMetaD potential every 0.2 ps, and a biasing potential set at 3000 K. Collective variables (CVs) were defined as the distance between the centers of mass of the peptide and the H-2D^b protein, named CV_dist(CM-CM), and the distance between the centers of mass of the two N-terminus binding cleft α-helical segments (group of residues 56-70 and residues 155-175), named CV_{dist(α-helices)} (see Figure 2 for the visualization of these CVs). CV_dist(CM-CM) varied from 8 Å to 38 Å with a Gaussian width of 1.25331 Å, and CV_{dist(α-helices)} varied from 10 Å to 28 Å with a Gaussian width of 0.626657 Å. Simulation-saved potential mean force (PMF) maps every 1 ps.

[figure omitted; refer to PDF]

Initially, the peptide was exposed to continuous rearrangements seeking favorable interactions, implying a transition from basins A to B. Consequently, in basin B, the complex adopts a new energy conformation that is ~1.15 kcal/mol lower in energy than that in basin A.

The deepest minimum, basin B, corresponds to the peptide in its binding site, in which some persistent H-bond interactions with H-2D^b occurred (see Table 1). In this basin, N-terminus peptide residue P1 formed an H-bond network with residues GLU63, LYS66, TYR159, and TRP167 in the cleft. Additionally, the peptide formed relevant H-bond interactions between residue P2 and residues GLU9, GLU63, and LYS66. On the other hand, the C-terminus peptide residue P14 interacted with the α2-helix through a significant H-bond interaction network, mainly with residue LYS146. Some other relevant interactions between the peptide’s central part and H-2D^b receptor were observed, namely, interactions with residues P4 and P6 (see Figure 4, box B).

The third energy basin (basin C) is ~7.31 kcal/mol higher in energy than basin B. B-to-C transition occurred after opening of the α-helices, and it might account for the dissociation of the peptide’s central part from H-2D^b. While the simulation was filling basin C, the interactions established between the N-terminus peptide residue P1 and residues GLU63 and TRP167 were conserved. Moreover, the C-terminus peptide residue P14 still maintained interaction with residue LYS146 of the α2-helix. Therefore, at this stage of the undocking pathway, the peptide was moved through the cleft of H-2D^b. The peptide stabilized at a position in which P1 maintained the H-bond with residue GLU63 of the α1-helix and residue TRP167 of the α2-helix and in which P14 maintained H-bond interaction with residue LYS146 of the α2-helix (see Figure 4, box C). WTMetaD overcame this energy barrier and pushed out the system from basin B to basin C. During the B-to-C transition, distance between α-helices reached values ~18 Å (CV_{dist(α-helices)}) to facilitate the dissociation of residues P2, P4, and P6 from the cleft.

Basin D was reached just after the peptide N-terminus region dissociation. In this metastable state, only C-terminus peptide residue P14 interacted with the receptor, maintaining H-bond interactions with the α2-helix of the H-2D^b protein, mainly with residue LYS146. After this stage of the undocking pathway, the ending of these H-bond interactions permitted the completed disruption of the peptide/H-2D^b complex (see Figure 4, box D).

At the last minimum (basin E), the distance between α-helices decreased to values close to 16 Å (CV_{dist(α-helices)}), and the distance between centers of mass of H-2D^b and the peptide achieved ~22.5 Å in CV_dist(CM-CM) (see Figure 4, box E). It indicated that the peptide reached the fully hydrated state. This state corresponds to the complete disruption of all interactions between the peptide and the H-2D^b receptor. At the end of the undocking pathway, the free energy of dissociation calculated as the difference between the energies in basins B (docking position) and E (undocking position) had a mean value of 8.23 kcal/mol. The gif animated image in Supplementary file 4 shows the projection of the FES plot of the configurations in CV space sampled during the WTMetaD simulation.

Time evolution of CV_dist(CM-CM) mostly describes the peptide exiting from the H-2D^b receptor. Therefore, the FES was projected onto the CV_dist(CM-CM) coordinate to measure each basin’s free energy value on this surface. In Figure 5, the projection was plotted by integrating the Boltzmann factor over the CV_{dist(α-helices)} (Equation (1)). This projection considered the energy values for each pair of CVs [44]. Free energy values for basins A, B, C, D, and E are 1.15, 0, 7.31, 5.72, and 5.10 kcal/mol, respectively.$$\begin{array}{}\text{(1)}& -\beta {\omega }_{\gamma }Ƶ=\ln \frac{\int e^{-\beta \omega Ƶ,Ƴ}dƳ}{\int e^{-\beta \omega Ƶ,Ƴ}dƵdƳ}.\end{array}$$

[figure omitted; refer to PDF]

In Equation (1), $Ƴ=\text{C}{\text{V}}_{\text{dist}\alpha -\text{helices}}$, $Ƶ=\text{C}{\text{V}}_{\text{dist}\text{CM}-\text{CM}}$, and $\beta =1/k_{b}T$, where $k_b$ is the Boltzmann constant and $T=310\text{K}$, $\omega$ is the function that gives the free energy for each pair of ($Ƶ$, $Ƴ$), and ${\omega }_{Ƴ}$ is the same function but fixing the value of $Ƴ$.

Binding contacts were monitored by analyzing the intermolecular interactions between the peptide and the cleft of H-2D^b. The “half-lives” of the interaction energy between residues from the peptide and the H-2D^b receptor were monitored. The first half-life was defined as a factor of the simulation time required for the interaction energy to reduce 50% of its initial value. In the same way, the second half-life is the factor required to reduce 75% of the initial value. Initial interaction energy contributions per residue are different. Therefore, half-life factors were weighted and rescaled according to Equation (2) to have values between 0 and 1.$$\begin{array}{}\text{(2)}& x_{\text{new}}=\frac{x-x_{\min }}{x_{\max }-x_{\min }}.\end{array}$$

In Equation (2), $x_{\text{new}}$ is the new weighted initial energy calculated for each peptide residue, $x$ is the original value of each residue’s initial energy, and $x_{\min }$ and $x_{\max }$, respectively, are the minimum and maximum values of the set of initial energies. In Figure 6, the half-life factors indicate the residues of the peptide that play an important role in the unbinding trajectory.

[figure omitted; refer to PDF]

One more aspect to consider about the interaction between the peptide and the H-2D^b receptor is the analysis of the long-range electrostatic energy as a function of the distance between the centers of mass of the peptide and H-2D^b, i.e., the CV_dist(CM-CM). Our analysis suggested a minimum distance between the ligand and receptor in which the long-range electrostatic attraction can promote an association of them. Our analysis detected three regions of interaction for nonbonded energy. The first region (labelled R₁ in Figure 7), with the lowest nonbonded energy values, corresponds to a distance between the peptide and H-2D^b of up to 18.5 Å. In this region, the nonbonded energy can promote the association of the peptide and the receptor. The second region (labelled R₂ in Figure 7) shows that the interaction intensity was reduced to values higher than the median of nonbonded energy obtained for the first region (-347.72 kcal/mol). In the second region, contacts between the H-2D^b protein and the peptide occur to form a precomplex. In this region, molecular diffusion plays a decisive role for binding, driving the complex formation process [45]. Finally, the third region (labelled R₃ in Figure 7) is where no influence of nonbonded energy was observed, with the intermolecular interaction energy approaching zero (less than 10% of the initial intermolecular interaction energy, i.e., -47.63 kcal/mol). Therefore, beyond the limit of 27.5 Å, no influence of the H-2D^b receptor over the peptide TT830-843 was observed.

[figure omitted; refer to PDF]

Binding simulation between the TT830-843 epitope with the H-2D^b protein assessed by in silico assay was confirmed by in vitro assay (see Figure 8). Real-time interactions were measured by surface plasmon resonance methodology, showing the parameters for complex formation. Kinetics of the interaction between DimerX and the peptide on the sensor chip after immobilization of the complexes was evaluated by two consecutive applications of the SensiQ Pioneer biosensor, an initial injection of the DimerX protein and a final injection of the peptide. DimerX was bound to the COOH5 chip by the Fc region and exhibited an average binding rate of 4548 RU/s (Figure 8(a)). Therefore, the α1 and α2 domains of the H-2D^b protein were free to interact with the peptide in solution, as previously shown to identify H-2 epitopes and detect specific CD8⁺ T-lymphocytes in the immune response of Leishmania (Leishmania) amazonensis cysteine proteinase B [46, 47].

[figures omitted; refer to PDF]

Interactions between the peptide and H-2D^b were performed at pH 7.2 at the temperature of 37°C. They showed $6.7\pm 0.5\times {10}^3$ RUmax values in three round assays (Figure 8(b)). Receptor-ligand binding kinetics were also accessed to determine the association ($k_a=9\pm 1\times {10}^5{\text{M}}^{-1}{\text{s}}^{-1}$) and dissociation ($k_d=0.31\pm 0.02{\text{s}}^{-1}$) rate constants, being possible to determine the equilibrium constant as $K_D=344\pm 60\text{nM}$. Results of SPR assays showed the reaction between immobilized H-2D^b and the epitope follows pseudo-first-order kinetics since the ligand concentration is constant, and the analyte concentration is in excess [48–50].

4. Discussion

Peptide TT830-843 is a universal epitope used as a control and adjuvant in acquired immunity studies. This peptide’s immunological properties should be a subject of fine molecular interaction studies since it may reveal the nuances of the protective efficiency of the immune response. This work showed the relevant structural and kinetic aspects of this linear epitope necessary to interact with H-2D^b proteins.

Exploring the peptide unbinding pathway from the H-2D^b receptor and reconstructing the FES assigned to this process, relevant aspects of the mechanisms involved in the peptide exiting were exhibited. Reconstructed FES reveals two distinct energy minima within the bound region, corresponding to different binding modes, labeled as A and B in Figure 4. During the transition from basins A to B, the peptide was exposed to continuous rearrangements seeking favorable interactions. These conformational rearrangements resulted in a protein-ligand complex with a tighter binding conformation. It allowed the complex to accommodate the peptide binding and led the peptide to adopt a low-energy conformation alternative state (~1.15 kcal/mol lower in energy). At this stage, the peptide’s central part bulged out until finding a favorable pose. Peptide length can explain this movement. This behavior was already reported for other haplotypes, assigning the peptide length to unpredictably find the exact binding mode for a peptide 14 aa in length [51]. A PDB entry in which the peptide is 11 aa long was used as a scaffold structure to build the complex; therefore, in the first step of WTMetaD, the system sought for a more favorable peptide initial position within the cleft.

In basin B (corresponding to the docked position), H-bond interactions formed between N-terminus peptide residue P1 and residues GLU63, LYS66, TYR159, and TRP167 stabilized this end of the bound peptide. Those H-bond interactions were previously described in the literature [31, 52, 53]. In particular, the influence of residue GLU63 on peptide binding is a key point to understand the effects of inducing cytokines to promote a T-cell response. GLU63 carboxylic moiety and the peptide P1 residue formed a stable H-bond interaction, indicating GLU63 plays an important role in binding affinity. Experimentally, GLU63 was also reported as a relevant residue for peptide binding, inducing the expression of specific cytokines that influence the balance between types of T-cell responses [54]. A second anchor occurred at residue P2, forming H-bond interactions with residues GLU9, GLU63, LYS66, and TYR159. Atoms of the P2 residue formed a persistent H-bond interaction with the polymorphic glutamic acid GLU9 of the beta-sheet floor in the H-2D^b receptor. Overall, WTMetaD revealed a persistent H-bond interaction network between peptide residues P1 and P2 and glutamic acid residues in the receptor’s groove. It suggests a functional role for glutamic acid residues in the cleft in TT830-843 peptide binding, notably GLU63 and GLU9.

The most persistent H-bond interaction occurred between the C-terminus peptide residue P14 and the α2-helix residue LYS146. Undocking pathway was completed only after the disruption of this interaction. Therefore, C-terminus residue P14 was the most important anchor for peptide binding, being a critical component along the undocking pathway.

A transition mechanism for the peptide dissociation was identified. Our analysis of the unbinding trajectory showed that the peptide rested in some stable state. The reconstructed FES consists of five basins, within which there are some transition events connected by transition states. During the transition from basins A to B, the peptide’s central part bulged flexibly out of the groove to compensate for length effects. Basin B corresponds to the peptide bounding state. In this state, the peptide stabilized, forming a persistent network of H-bond interactions. Basin C corresponds to an intermediary state visited by the peptide exiting from the cleft. Basin C was visited just after the opening of the α-helices, which permitted the dissociation of the peptide’s central part but maintained some important H-bond interactions with H-2D^b. Therefore, the peptide’s N- and C-termini anchored the binding to the receptor, with its central part bulged out of the groove. In basin D, only the C-terminus residue P14 anchored the peptide, maintaining H-bond interactions only with the H-2D^bα2-helix. In basin E, the distance between the centers of mass of the peptide and H-2D^b indicated that the peptide achieved a fully hydrated state.

WTMetaD simulation provided an approach to reconstruct the underlying FES for the peptide unbinding process, which allowed deriving relevant details about this energetic profile (see Figure 5). In FES, the global energy minimum (basin B) corresponds to the docked state at ~11 Å for CV_dist(CM-CM). Transition from basins A to B had to overcome a barrier of ~3 kcal/mol, seeking a tighter peptide binding conformation. After sampling the global minimum, the simulation had to overcome three energy barriers, totalizing ~12 kcal/mol. The simulation overcame these barriers, escaping from basin B. It permitted the opening of the α-helices to allow the dissociation of the peptide’s central part (transition from basins B to C). After sampling basin C, the simulation reached basin D, the last step of the dissociation process. Finally, the peptide had to pass the last barrier to escape from basin D and fall into the unbound state E. Therefore, the basins were separated by free energy barriers, which were larger than the thermal energy at 310 K ($k_{B}T=0.62\text{kcal}/\text{mol}$). It implies that each barrier crossing would take a time scale inaccessible running a standard MD simulation. The mean of the free energy of dissociation ($\Delta G_{\text{diss}}=8.23\text{kcal}/\text{mol}$) indicated an external driving force should be applied to dissociate the peptide from the receptor. Therefore, the complex peptide/H-2D^b is thermodynamically stable.

In silico mutations were introduced using the psfgen plugin of VMD [34] to evaluate the effects of the substitution of P14 for ALA and GLY residues in WTMetaD simulation. The alanine mutation resulted in a loss of the H-bond association rate with LYS146, but no substantial affinity changes relative to the whole dissociation energy between the peptide and H-2D^b were observed ($\Delta G_{\text{diss}}=0.58\text{kcal}/\text{mol}$). This effect of alanine substitution retains the beta carbon and introduces a small mutant side chain for the P14 residue. This mutation increased the number of H-bond interactions formed between backbone residues of mutant P14 alanine and the cleft of H-2D^b (see Table S3 in Supplementary file 1). The glycine mutation removes the beta carbon and can cause flexibility and conformational changes in the peptide. We observed an absence of highly persistent H-bound interactions between the mutant P14 glycine and the H-2D^b F pocket, mainly due to the loss of persistence with LYS146 (see Table S4 in Supplementary file 1). It suggests a negative impact on the peptide/H-2D^b complex stability and can explain the lower free energy of dissociation estimated for this mutant peptide ($\Delta G_{\text{diss}}=5.16\text{kcal}/\text{mol}$).

Interaction kinetics are important to analyze (un)binding processes [45], and the half-life factors can help to examine the role of some peptide residues. Kinetic analysis introduces the time component of the observation and reveals other aspects of the interaction between the ligand and its receptor. For instance, in drug discovery research, the kinetic of the drug-receptor binding process can be as important or even more important than affinity in determining drug efficacy, particularly when the pharmacological duration effect is a significant component of in vivo efficacy [45, 55–57]. A kinetic analysis was retrieved from the values of half-life factors (described in Figure 6) and the persistence of H-bond interactions in Table 1. Such analysis examined the role of some residues to maintain the peptide binding. It is related to the time during which the peptide remains bound. In Figure 6, the half-life factors are the highest for N- and C-termini, indicating that these regions play an important role in the unbinding trajectory. Specifically, the C-terminus residue P14 has the highest half-life factor and has the longest persistence of H-bond interactions (Table 1). Therefore, it suggests an important role of this residue to drive peptide binding, mainly when interacting with residue LYS146 in the cleft. Standard in vitro or in vivo assays cannot identify the peptide residues’ atomic rearrangements within the cleft and their interaction kinetics. The approach applied in this work allowed to simulate the peptide dissociation on computationally tractable time scales. Therefore, it can be a useful method to observe the interaction kinetics and, consequently, shed light on the molecular mechanism of the peptide exiting. We hypothesize that residues with longer residence time on the receptor are kinetically decisive for peptide binding.

Protein-protein binding sites can reveal an amino acid residue network that communicates structurally and energetically with one another. This site-to-site communication contributes to protein binding, and interactions over long distances provide long-range communication between such binding sites [58]. For binding to occur, the approximation of the two biomolecules should happen at a suitable distance and orientation, for which the molecular diffusion acts on both the ligand and the receptor to form a precomplex. A precomplex is defined as a transitional state necessary for the binding of two molecules. Long-range electrostatic energy can promote this encounter, and it can evolve to a stable state to form a complex [59, 60].

In previous work, our results suggested a behavior pattern for peptides that bind H-2 receptors [31]. Those results were obtained assuming a preset condition for the parameters used to perform enhanced sampling simulations. Computing the nonbonded energy between the peptide and H-2D^b, three interaction regions were detected. Our results indicated that an encounter between the peptide and the H-2D^b receptor can happen at a distance ~27.5 Å to form a precomplex. Within this distance, molecular diffusion can favor the peptide binding, and the process can evolve to form a stable complex. We suggest that this feature can be part of the peptide loading mechanism that occurs in endosomes for cross-presentation. In this case, peptide loading is considered a diffusional step with an entropic barrier due to the penalty associated with peptide rotational and translational entropy loss. Within a certain distance, the peptide can overcome entropic barriers and form a diffusional encounter complex that evolves to bind the receptor.

In SPR assays, although the kinetic findings are in accordance with predicting good protein-protein binding affinity [61], further comments are needed. SPR assay conditions do not reflect the cellular physiochemical conditions in which these interactions occur. However, the measured $k_a$ value indicates that the peptide can form a complex with the H-2D^b protein in antigen-presenting cells during the immune response. This fast and spontaneous binding can be one reason for preserving this fragment from complete enzymatic hydrolysis by the proteasome system. Low $k_d$ values found indicate that the peptide maintains a more stable complex with the H-2D^b protein to persist long enough to interact with T-lymphocytes.

Additionally, as the recombinant H-2D^b protein is supplemented with recombinant beta2-microglobulin (β2M), SPR data can be related to the supramolecular complex stability for T-lymphocyte recognition by this epitope [62]. An important experimental consideration for these assays was the Gibbs free energy of the formed complex: $\Delta \text{G}=-8.3\text{kcal}/\text{mol}$. $\Delta \text{G}$ negative values indicated spontaneous binding for supramolecular complex formation under pH 7.2 and 37°C. Therefore, the peptide/H-2D^b complex presented enough free energy to interact with a TCR. Estimated energy suggests the complex should be stable on the antigen-presenting cell surface, whose pH is close to neutral, as described for epitopes from the C-terminus extension of cysteine proteinase B from L. (L.) amazonensis [46].

5. Conclusion

Using an in silico method, some key residues for peptide binding to H-2D^b receptors were identified. Glutamic acids in the cleft, notably GLU9 and GLU63, play an essential role in peptide binding. It also pointed out that the C-terminus residue P14 is the most important anchor for the peptide interaction with the H-2D^b receptor.

During the simulation, the exiting from the cleft is characterized by a stepwise mechanism between progressively detached states until the full dissociation of the peptide. An important functional feature for peptide detachment was observed: the opening of α-helices to permit the peptide’s central part dissociates from the cleft. Additionally, some relevant interactions that occurred during the dissociation process were described: in particular, the H-bond interactions occurred between the receptor and C- and N-terminus peptide residues. These unbinding pathway functional features were detected running a relatively short simulation. This in silico assay was crucial to describe important functional features, which has been impossible to investigate using experimental techniques. Therefore, we assumed WTMetaD as a valuable method to confirm interactions among multiple residues of a protein complex. This method can be used to investigate binding contacts that have a functional implication.

The approach applied in this work helped to understand the peptide dissociation mechanistic basis and indicate whether a complex is a thermodynamically stable system. To date, we do not know any experimental results of calculating the free energy of dissociation for complexes formed between MHC class I receptors and 14 residue-long peptides. Therefore, although previous work has used this approach to calculate and compare the free energy of dissociation for 8-10 amino acid peptides in complexes with MHC class I [31], further studies should be conducted to evaluate the effect of longer peptides in determining the free energy of dissociation in complexes with 14 residue-long peptides.

FES projection onto the CV_dist(CM-CM) revealed the energetic profile of the peptide exiting, providing an understanding of how the unbinding process overcame the surface’s energy barriers.

We suggested a minimal distance between the peptide and the H-2D^b receptor to favor the peptide binding. We can understand this feature as a diffusional step that can be part of the peptide loading mechanism, in which the peptide must overcome an entropic barrier to load into the receptor. Future work can also consider this feature to figure out the peptide presentation to the immune system and how they compete to bind to MHC class I molecules.

In silico assays were confirmed by additional surface biosensing experiments that use the DimerX (H-2D^b haplotype). Data gathered in these assays reinforce that the peptide can form a stable complex with H-2D^b. This complex exhibits enough free energy to interact with a TCR on the antigen-presenting cell surface. Therefore, the DimerX used in SPR assays in combination with WTMetaD simulation indicates a promising approach to study real-time interactions between ligands and receptors. In our specific case, the combination of in silico and in vitro assays provided significant evidence supporting the formation of a stable complex between the peptide TT830-843 and H-2D^b receptor.

Acknowledgments

We are grateful to the technical support of Fundação Oswaldo Cruz and Instituto Oswaldo Cruz platforms: Water Platform Grade Reagent Type I and II and Surface Resonance Plasmonic Platform (RPT03E), Conselho Nacional de Desenvolvimento Científico e Tecnológico-Brasil (CNPq), Coordenação de Aperfeiçoamento de Pessoal de Nível Superior-Brasil (CAPES), and FAPERJ. This study was financed in part by Conselho Nacional de Desenvolvimento Científico e Tecnológico-Brasil (CNPq: 301744/2019-0), Coordenação de Aperfeiçoamento de Pessoal de Nível Superior-Brasil (CAPES) (Finance Code 001), and Fundação Carlos Chagas Filho de Amparo à Pesquisa (FAPERJ) (processo: E-26/010.002021/2019).

References

[1] World Health Organization: Tetanus, 2020. https://www.who.int/immunization/diseases/tetanus/en/

[2] T. V. Murphy, B. A. Slade, K. R. Broder, K. Kretsinger, T. Tiwari, P. M. Joyce, J. K. Iskander, K. Brown, J. S. Moran, "Advisory Committee on Immunization Practices (ACIP) Centers for Disease Control and Prevention (CDC). Prevention of pertussis, tetanus, and diphtheria among pregnant and postpartum women and their infants. Recommendations of the Advisory Committee on Immunization Practices (ACIP)," Morbidity and Mortality Weekly Report (MMWR), vol. 57 no. RR-4, 2008.

[3] MEROPS, The peptidase database, 2020. https://www.ebi.ac.uk/merops/cgi-bin/pepsum?id=M27.001;type=P

[4] U. Eisel, W. Jarausch, K. Goretzki, A. Henschen, J. Engels, U. Weller, M. Hudel, E. Habermann, H. Niemann, "Tetanus toxin: primary structure, expression in E. coli, and homology with botulinum toxins," The EMBO Journal, vol. 5 no. 10, pp. 2495-2502, DOI: 10.1002/j.1460-2075.1986.tb04527.x, 1986.

[5] J. H. Kerr, J. L. Corbett, C. Prys-Roberts, A. C. Smith, J. M. Spalding, "Involvement of the sympathetic nervous system in tetanus: studies on 82 cases," Lancet, vol. 2 no. 7562, pp. 236-241, DOI: 10.1016/s0140-6736(68)92350-7, 1968.

[6] D. J. Kanarek, B. Kaufman, S. Zwi, "Severe sympathetic hyperactivity associated with tetanus arch," Internal Medicine, vol. 132 no. 4, pp. 602-604, 1973.

[7] V. M. Hollow, G. M. Clarke, "Autonomic manifestations of tetanus," Anaesthesia and Intensive Care, vol. 3 no. 2, pp. 142-147, DOI: 10.1177/0310057X7500300211, 1975.

[8] N. Buchanan, G. W. Cane, G. Wolfson, M. de Andrade, "Autonomic dysfunction in tetanus: the effects of a variety of therapeutic agents, with special reference to morphine," Intensive Care Medicine, vol. 5 no. 2, pp. 65-68, DOI: 10.1007/BF01686048, 1979.

[9] R. F. Coico, G. Sunshine, "Immunology: a short course," In: Annals of Biomedical Engineering, 2009.

[10] T. R. Poulsen, P. J. Meijer, A. Jensen, L. S. Nielsen, O. S. Andersen, "Kinetic, affinity, and diversity limits of human polyclonal antibody responses against tetanus toxoid," Journal of Immunology, vol. 179 no. 6, pp. 3841-3850, DOI: 10.4049/jimmunol.179.6.3841, 2007.

[11] M. B. Spellerberg, D. Zhu, A. Thompsett, C. A. King, T. J. Hamblin, F. K. Stevenson, "DNA vaccines against lymphoma: promotion of anti-idiotypic antibody responses induced by single chain Fv genes by fusion to tetanus toxin fragment C," The Journal of Immunology, vol. 159 no. 4, pp. 1885-1892, 1997.

[12] C. M. Khan, B. Villarreal-Ramos, R. J. Pierce, G. Riveau, R. Demarco de Hormaeche, H. McNeill, T. Ali, N. Fairweather, S. Chatfield, A. Capron, "Construction, expression, and immunogenicity of the Schistosoma mansoni P28 glutathione S-transferase as a genetic fusion to tetanus toxin fragment C in a live aro attenuated vaccine strain of Salmonella," Proceedings of the National Academy of Sciences, vol. 91 no. 23, pp. 11261-11265, DOI: 10.1073/pnas.91.23.11261, 1994.

[13] C. M. Khan, B. Villarreal-Ramos, R. J. Pierce, R. Demarco de Hormaeche, H. McNeill, T. Ali, S. Chatfield, A. Capron, G. Dougan, C. E. Hormaeche, "Construction, expression, and immunogenicity of multiple tandem copies of the Schistosoma mansoni peptide 115-131 of the P28 glutathione S-transferase expressed as C-terminal fusions to tetanus toxin fragment C in a live aro-attenuated vaccine strain of Salmonella," The Journal of Immunology, vol. 153 no. 12, pp. 5634-5642, 1994.

[14] P. Panina-Bordignon, A. Tan, A. Termijtelen, S. Demotz, G. Corradin, A. Lanzavecchia, "Universally immunogenic T cell epitopes: promiscuous binding to human MHC class II and promiscuous recognition by T cells," European Journal of Immunology, vol. 19 no. 12, pp. 2237-2242, DOI: 10.1002/eji.1830191209, 1989.

[15] P. T. Kaumaya, S. Kobs-Conrad, Y. H. Seo, H. Lee, A. M. Vanbuskirk, N. Feng, J. F. Sheridan, V. Stevens, "Peptide vaccines incorporating a 'promiscuous' T-cell epitope bypass certain haplotype restricted immune responses and provide broad spectrum immunogenicity," Journal of Molecular Recognition, vol. 6 no. 2, pp. 81-94, DOI: 10.1002/jmr.300060206, 1993.

[16] J. I. Jang, J. S. Kim, J. S. Eom, H. G. Kim, B. H. Kim, S. Lim, I. S. Bang, Y. K. Park, "Expression and delivery of tetanus toxin fragment C fused to the N-terminal domain of SipB enhances specific immune responses in mice," Microbiology and Immunology, vol. 56 no. 9, pp. 595-604, DOI: 10.1111/j.1348-0421.2012.00480.x, 2012.

[17] B. Boitel, U. Blank, D. Mège, G. Corradin, J. Sidney, A. Sette, O. Acuto, "Strong similarities in antigen fine specificity among DRB1 ∗ 1302-restricted tetanus toxin tt830-843-specific TCRs in spite of highly heterogeneous CDR3," Journal of Immunology, vol. 154 no. 7, pp. 3245-3355, 1995.

[18] F. Wu, Y. Huang, P. Zhang, C. Wang, Y. Tian, L. Lu, G. He, L. Yang, "Interleukin-13 peptide vaccine induces protective humoral immunity in murine asthma models," Oncotarget, vol. 9 no. 6, article 19950, pp. 6678-6690, DOI: 10.18632/oncotarget.19950, 2018.

[19] E. Simone, D. Daniel, N. Schloot, P. Gottlieb, S. Babu, E. Kawasaki, D. Wegmann, G. S. Eisenbarth, "T cell receptor restriction of diabetogenic autoimmune NOD T cells," Proceedings of the National Academy of Sciences, vol. 94 no. 6, pp. 2518-2521, DOI: 10.1073/pnas.94.6.2518, 1997.

[20] M. Nakayama, T. Castoe, T. Sosinowski, X. L. He, K. Johnson, K. Haskins, D. A. A. Vignali, L. Gapin, D. Pollock, G. S. Eisenbarth, "Germline TRAV5D-4 T-cell receptor sequence targets a primary insulin peptide of NOD mice," Diabetes, vol. 61 no. 4, pp. 857-865, DOI: 10.2337/db11-1113, 2012.

[21] M. Kobayashi, N. Abiru, T. Arakawa, K. Fukushima, H. Zhou, E. Kawasaki, H. Yamasaki, E. Liu, D. Miao, F. S. Wong, G. S. Eisenbarth, K. Eguchi, "Altered B:9–23 insulin, when administered Intranasally with cholera toxin adjuvant, suppresses the expression of insulin autoantibodies and prevents diabetes," The Journal of Immunology, vol. 179 no. 4, pp. 2082-2088, DOI: 10.4049/jimmunol.179.4.2082, 2007.

[22] R. N. Rodrigues-da-Silva, D. Correa-Moreira, I. F. Soares, P. M. de-Luca, P. R. R. Totino, F. N. Morgado, M. G. Oliveira Henriques, A. L. Peixoto Candea, B. Singh, M. R. Galinski, A. Moreno, J. Oliveira-Ferreira, J. C. Lima-Junior, "Immunogenicity of synthetic peptide constructs based on PvMSP9 E795-A808 , a linear B-cell epitope of the P . vivax Merozoite Surface Protein-9," Vaccine, vol. 37 no. 2, pp. 306-313, DOI: 10.1016/j.vaccine.2018.10.016, 2019.

[23] T. D. Shi, Y. Z. Wu, Z. C. Jia, L. Y. Zou, W. Zhou, "Therapeutic polypeptides based on HBV core 18-27 epitope can induce CD8+ CTL-mediated cytotoxicity in HLA-A2+ human PBMCs," World Journal of Gastroenterology, vol. 10 no. 13, pp. 1902-1906, DOI: 10.3748/wjg.v10.i13.1902, 2004.

[24] E. S. Bergmann-Leitner, S. Chaudhury, N. J. Steers, M. Sabato, V. Delvecchio, A. S. Wallqvist, C. F. Ockenhouse, E. Angov, "Computational and experimental validation of B and T-cell epitopes of the in vivo immune response to a novel malarial antigen," PLoS One, vol. 8, article e71610 no. 8,DOI: 10.1371/journal.pone.0071610, 2013.

[25] J. L. Sanchez-Trincado, M. Gomez-Perosanz, P. A. Reche, "Fundamentals and methods for T- and B-cell epitope prediction," Journal of Immunology Research, vol. 2017,DOI: 10.1155/2017/2680160, 2017.

[26] A. Barducci, G. Bussi, M. Parrinello, "Well-tempered metadynamics: a smoothly converging and tunable free-energy method," Physical Review Letters, vol. 100 no. 2, article 020603,DOI: 10.1103/PhysRevLett.100.020603, 2008.

[27] R. B. M. Schasfoort, "Chapter 1: introduction to surface plasmon resonance," Handbook of Surface Plasmon Resonance, vol. no. 2,DOI: 10.1039/9781788010283-00001, 2017.

[28] M. Karplus, J. A. McCammon, "Molecular dynamics simulations of biomolecules," Nature Structural Biology, vol. 9 no. 9, pp. 646-652, DOI: 10.1038/nsb0902-646, 2002.

[29] S. A. Adcock, J. A. McCammon, "Molecular dynamics: survey of methods for simulating the activity of proteins," Chemical reviews, vol. 106 no. 5, pp. 1589-1615, DOI: 10.1021/cr040426m, 2006.

[30] D. Branduardi, G. Bussi, M. Parrinello, "Metadynamics with adaptive Gaussians," Journal of Chemical Theory and Computation, vol. 8 no. 7, pp. 2247-2254, DOI: 10.1021/ct3002464, 2012.

[31] A. A. M. L. Brandt, P. R. Batista, F. Souza-Silva, C. R. Alves, E. R. Caffarena, "Exploring the unbinding ofLeishmania(L.)amazonensisCPB derived-epitopes from H2 MHC class I proteins," Proteins, vol. 84 no. 4, pp. 473-487, DOI: 10.1002/prot.24994, 2016.

[32] H. M. Berman, J. Westbrook, Z. Feng, G. Gilliland, T. N. Bhat, H. Weissig, I. N. Shindyalov, P. E. Bourne, "The Protein Data Bank," Nucleic Acids Research, vol. 28 no. 1, pp. 235-242, DOI: 10.1093/nar/28.1.235, 2000.

[33] C. Ciatto, A. C. Tissot, M. Tschopp, G. Capitani, F. Pecorari, A. Pluckthun, M. G. Grütter, "Zooming in on the hydrophobic ridge of H-2D b : implications for the conformational variability of bound peptides 1," Journal of molecular biology, vol. 312 no. 5, pp. 1059-1071, DOI: 10.1006/jmbi.2001.5016, 2001.

[34] W. Humphrey, A. Dalke, K. Schulten, "VMD: visual molecular dynamics," Journal of molecular graphics, vol. 14 no. 1, pp. 33-38, DOI: 10.1016/0263-7855(96)00018-5, 1996.

[35] B. Raveh, N. London, O. Schueler-Furman, "Sub-angstrom modeling of complexes between flexible peptides and globular proteins," Proteins, vol. 78 no. 9, pp. 2029-2040, DOI: 10.1002/prot.22716, 2010.

[36] R. Das, D. Baker, "Macromolecular modeling with Rosetta," Annual Review of Biochemistry, vol. 77 no. 1, pp. 363-382, DOI: 10.1146/annurev.biochem.77.062906.171838, 2008.

[37] J. C. Phillips, R. Braun, W. Wang, J. Gumbart, E. Tajkhorshid, E. Villa, C. Chipot, R. D. Skeel, L. Kalé, K. Schulten, "Scalable molecular dynamics with NAMD," Journal of computational chemistry, vol. 26 no. 16, pp. 1781-1802, DOI: 10.1002/jcc.20289, 2005.

[38] A. D. MacKerell, D. Bashford, M. Bellott, R. L. Dunbrack, J. D. Evanseck, M. J. Field, S. Fischer, J. Gao, H. Guo, S. Ha, D. Joseph-McCarthy, L. Kuchnir, K. Kuczera, F. T. Lau, C. Mattos, S. Michnick, T. Ngo, D. T. Nguyen, B. Prodhom, W. E. Reiher, B. Roux, M. Schlenkrich, J. C. Smith, R. Stote, J. Straub, M. Watanabe, J. Wiórkiewicz-Kuczera, D. Yin, M. Karplus, "All-atom empirical potential for molecular modeling and dynamics studies of proteins," The Journal of Physical Chemistry. B, vol. 102 no. 18, pp. 3586-3616, DOI: 10.1021/jp973084f, 1998.

[39] W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey, M. L. Klein, "Comparison of simple potential functions for simulating liquid water," The Journal of Chemical Physics, vol. 79 no. 2, pp. 926-935, DOI: 10.1063/1.445869, 1983.

[40] R. Yoshida, "MHC class I recognition by monocyte-/macrophage-specific receptors," Advances in immunology, vol. 124, pp. 207-247, DOI: 10.1016/B978-0-12-800147-9.00007-8, 2014.

[41] E. Seong, T. L. Saunders, C. L. Stewart, M. Burmeister, "To knockout in 129 or in C57BL/6: that is the question," Trends in genetics, vol. 20 no. 2, pp. 59-62, DOI: 10.1016/j.tig.2003.12.006, 2004.

[42] A. K. Bronowska, "Thermodynamics of Ligand-Protein Interactions: Implications for Molecular Design," Thermodynamics-Interaction Studies-Solids, Liquids and Gases,DOI: 10.5772/19447, 2011.

[43] M. Zacharias, S. Springer, "Conformational flexibility of the MHC class I α 1 - α 2 domain in peptide bound and free states: a molecular dynamics simulation study," Biophysical Journal, vol. 87 no. 4, pp. 2203-2214, DOI: 10.1529/biophysj.104.044743, 2004.

[44] J. P. M. Jämbeck, A. P. Lyubartsev, "Exploring the free energy landscape of solutes embedded in lipid bilayers," The Journal of Physical Chemistry Letters, vol. 4 no. 11, pp. 1781-1787, DOI: 10.1021/jz4007993, 2013.

[45] M. Bernetti, A. Cavalli, L. Mollica, "Protein-ligand (un)binding kinetics as a new paradigm for drug discovery at the crossroad between experiments and modelling," Medchemcomm, vol. 8 no. 3, pp. 534-550, DOI: 10.1039/c6md00581k, 2017.

[46] F. Souza-Silva, B. A. Pereira, L. C. Finkelstein, V. Zucolotto, E. R. Caffarena, C. R. Alves, "Dynamic identification of H2 epitopes from Leishmania (Leishmania) amazonensis cysteine proteinase B with potential immune activity during murine infection," Journal of Molecular Recognition, vol. 27 no. 2, pp. 98-105, DOI: 10.1002/jmr.2334, 2014.

[47] L. S. Da Silveira-Júnior, F. Souza-Silva, B. A. Pereira, L. Cysne-Finkelstein, G. B. Cavalcanti Júnior, C. R. Alves, "Exploring the association of surface plasmon resonance with recombinant MHC:Ig hybrid protein as a tool for detecting T lymphocytes in mice infected with Leishmania ( Leishmania ) amazonensis," BioMed Research International, vol. 2017,DOI: 10.1155/2017/9089748, 2017.

[48] P. Bjorquist, S. Bostrom, "Determination of the kinetic constants of tissue factor/factor VII/factor VIIA and antithrombin/heparin using surface plasmon resonance," Thrombosis Research, vol. 85 no. 3, pp. 225-236, DOI: 10.1016/S0049-3848(97)00007-8, 1997.

[49] R. W. Glaser, "Antigen-antibody binding and mass transport by convection and diffusion to a surface: a two-dimensional computer model of binding and dissociation kinetics," Analytical Biochemistry, vol. 213 no. 1, pp. 152-161, DOI: 10.1006/abio.1993.1399, 1993.

[50] D. J. Oshannessy, M. Brighamburke, K. K. Soneson, P. Hensley, I. Brooks, "Determination of rate and equilibrium binding constants for macromolecular interactions using surface plasmon resonance: use of nonlinear least squares analysis methods," Analytical biochemistry, vol. 212 no. 2, pp. 457-468, DOI: 10.1006/abio.1993.1355, 1993.

[51] M. Probst-Kepper, H. J. Hecht, H. Herrmann, V. Janke, F. Ocklenburg, J. Klempnauer, B. J. van den Eynde, S. Weiss, "Conformational restraints and flexibility of 14-meric peptides in complex with HLA-B ∗ 3501," Journal of Immunology, vol. 173 no. 9, pp. 5610-5616, DOI: 10.4049/jimmunol.173.9.5610, 2004.

[52] D. Hudrisier, H. Mazarguil, F. Laval, M. B. Oldstone, J. E. Gairin, "Binding of Viral Antigens to Major Histocompatibility Complex Class I H-2D b Molecules Is Controlled by Dominant Negative Elements at Peptide Non-anchor Residues:," The Journal of biological chemistry, vol. 271 no. 30, pp. 17829-17836, DOI: 10.1074/jbc.271.30.17829, 1996.

[53] A. C. Young, S. G. Nathenson, J. C. Sacchettini, "Structural studies of class I major histocompatibility complex proteins: insights into antigen presentation," The FASEB journal, vol. 9 no. 1, pp. 26-36, DOI: 10.1096/fasebj.9.1.7821756, 1995.

[54] B. A. S. Pereira, F. S. Silva, K. M. Rebello, M. Marin-Villa, Y. M. Traub-Cseko, T. C. B. Andrade, A. L. Bertho, E. R. Caffarena, C. R. Alves, "In silico predicted epitopes from the COOH-terminal extension of cysteine proteinase B inducing distinct immune responses during Leishmania (Leishmania) amazonensis experimental murine infection," BMC Immunology, vol. 12 no. 1,DOI: 10.1186/1471-2172-12-44, 2011.

[55] R. A. Copeland, D. L. Pompliano, T. D. Meek, "Drug-target residence time and its implications for lead optimization," Nature Reviews. Drug Discovery, vol. 5 no. 9, pp. 730-739, DOI: 10.1038/nrd2082, 2006.

[56] A. C. Pan, D. W. Borhani, R. O. Dror, D. E. Shaw, "Molecular determinants of drug-receptor binding kinetics," Drug Discovery Today, vol. 18 no. 13-14, pp. 667-673, DOI: 10.1016/j.drudis.2013.02.007, 2013.

[57] D. C. Swinney, "The role of binding kinetics in therapeutically useful drug action," Current Opinion in Drug Discovery & Development, vol. 12 no. 1, pp. 31-39, 2009.

[58] E. J. Sundberg, "The molecular architecture of protein–protein binding sites," Computational Protein-Protein Interactions, pp. 43-55, DOI: 10.1201/9781420070071.ch3, 2009.

[59] R. R. Gabdoulline, R. C. Wade, "Biomolecular diffusional association," Current Opinion in Structural Biology, vol. 12 no. 2, pp. 204-213, DOI: 10.1016/s0959-440x(02)00311-1, 2002.

[60] J. Janin, "Kinetics and thermodynamics of protein-protein interactions," Protein-Protein Recognition, 2000.

[61] V. G. Veresov, "Predicting protein-protein binding affinity by in silico docking," JSM Chemistry, vol. 4 no. 1, 2016.

[62] L. Ø. Pedersen, A. Stryhn, T. L. Holtet, M. Etzerodt, J. Gerwien, M. H. Nissen, H. C. Thøgersen, S. Buus, "The interaction of beta 2-microglobulin (beta 2m) with mouse class I major histocompatibility antigens and its ability to support peptide binding. A comparison of human and mouse beta 2m," European journal of immunology, vol. 25 no. 6, pp. 1609-1616, DOI: 10.1002/eji.1830250621, 1995.